Ticket #11380: trac_11380_cont_frac_quadratic_2.patch

File trac_11380_cont_frac_quadratic_2.patch, 31.1 KB (added by mmasdeu, 2 years ago)

Second version

  • sage/rings/number_field/number_field_element_quadratic.pyx 

    # HG changeset patch
# User Marc Masdeu <masdeu@math.columbia.edu>
# Date 1306446759 14400
# Node ID 0c92a2edd152e4089378686528b1139a0c9f6fb2
# Parent  b97e6e8ba4a1155e42d5b8c355e065c76ca1add8
trac 11380: Fixed minor issues to make it (much) faster.

diff -r b97e6e8ba4a1 -r 0c92a2edd152 sage/rings/number_field/number_field_element_quadratic.pyx
    a b  
    14161416            P=QuadraticContinuedFraction(F,Nbound,Tbound) 
    14171417            P.solve() 
    14181418            F.__continued_fraction_cache=P 
    1419         print 'Done with precomputations...' 
     1419            print 'Done with precomputations...' 
    14201420        v=[] 
    14211421        knm1=self 
    14221422        kn=1 

  • deleted file sage/rings/number_field/quadratic_continued_fraction.py 

    diff -r b97e6e8ba4a1 -r 0c92a2edd152 sage/rings/number_field/quadratic_continued_fraction.py
    + -  
    1  
    2 ########################################## 
    3 # 
    4 # Added by mmasdeu 
    5 # Authors: Xevi Guitart and Marc Masdeu 
    6 # 
    7 ########################################## 
    8  
    9 from sage.matrix.constructor import Matrix 
    10 from sage.rings.all import RealField 
    11 from sage.structure.sage_object import SageObject 
    12 from sage.misc.all import union 
    13 from sage.functions.all import (floor,ceil) 
    14 from sage.rings.all import IntegerRing 
    15 from sage.rings.all import RationalField 
    16 from collections import deque 
    17  
    18 class Regions(SageObject): 
    19     def __init__(self,parent): 
    20         self._regions=dict([]) 
    21         self._F=parent._F 
    22         self._eps=parent._eps 
    23         self._epsto=[self._F(1)] 
    24         self.fundom_rep=parent.fundom_rep 
    25         self.embed=parent.embed 
    26         self._parent=parent 
    27         self._N=1 # Means that we have regions up to _N-1=0, so none 
    28         self._area=self._parent._Rw 
    29  
    30     def get_regions(self,B): 
    31         regions=self._regions 
    32         if(B>self._N): 
    33             F=self._F 
    34             eps=self._eps 
    35             newelts=dict([(IntegerRing()(n),[I.gens_reduced()[0] for I in v if I.is_principal()]) for n,v in self._F.ideals_of_bdd_norm(B-1).iteritems() if n>=self._N]) 
    36             self._N=B 
    37             for nn in range(len(self._epsto),B): 
    38                 self._epsto.append(self._epsto[nn-1]*eps) 
    39             for nn,v in newelts.iteritems(): 
    40                 assert not regions.has_key(nn) 
    41                 regions[nn]=[] 
    42                 for q in v: 
    43                     invden=1/q 
    44                     invden_red=self.fundom_rep(invden) 
    45                     q0=invden_red-invden 
    46                     m_invden=-invden 
    47                     m_invden_red=self.fundom_rep(m_invden) 
    48                     m_q0=m_invden_red-m_invden 
    49                     a=invden_red[0] 
    50                     b=invden_red[1] 
    51                     am=m_invden_red[0] 
    52                     bm=m_invden_red[1] 
    53  
    54                     regions[nn].extend([[invden_red,[q0,q],Region(self.embed(invden_red),RealField()(1/nn)),[a,b],nn],[m_invden_red,[m_q0,-q],Region(self.embed(m_invden_red),RealField()(1/nn)),[am,bm],nn]]) 
    55                     for jj in range(1,nn): 
    56                         x=self._epsto[jj]*invden 
    57                         x_red=self.fundom_rep(x) 
    58                         x_red_minus=self.fundom_rep(-x_red) 
    59  
    60                         regs=regions[nn] 
    61                         if(x_red==regs[0][0] or x_red_minus==regs[0][0]): 
    62                             continue 
    63  
    64                         q0_minus=x_red_minus+x 
    65                         q0=x_red-x 
    66                         q1=1/x 
    67                         a=x_red[0] 
    68                         b=x_red[1] 
    69                         am=x_red_minus[0] 
    70                         bm=x_red_minus[1] 
    71                         regions[nn].extend([[x_red,[q0,q1],Region(self.embed(x_red),RealField()(1/nn)),[a,b],nn],[x_red_minus,[q0_minus,-q1],Region(self.embed(x_red_minus),RealField()(1/nn)),[am,bm],nn]]) 
    72  
    73         else: 
    74             # We have all the required regions 
    75             pass 
    76         return [reg for nn,lst in regions.iteritems() if nn<B for reg in lst] 
    77  
    78 class Hole(SageObject): 
    79     def __init__(self,parent,v,depth=0): 
    80         self._parent=parent 
    81         if isinstance(v,list): 
    82             #Input in as [xmin,ymin,xmax,ymax] 
    83             if not (v[0]<v[2] and v[1]<v[3]): 
    84                 print 'The input should be of the form [xmin,ymin,xmax,ymax], with xmin<xmax and ymin<ymax.' 
    85                 assert(0) 
    86             (self.xmin,self.ymin,self.xmax,self.ymax)=v 
    87         elif(isinstance(v,Hole)): 
    88             self.xmin=v.xmin 
    89             self.ymin=v.ymin 
    90             self.xmax=v.xmax 
    91             self.ymax=v.ymax 
    92         self.size=max([self.xmax-self.xmin,self.ymax-self.ymin]) 
    93         self._corners=None 
    94         self._depth=depth 
    95  
    96     def depth(self): 
    97         return self._depth 
    98  
    99     def corners(self): 
    100         if self._corners==None: 
    101             self._corners=[(self.xmin,self.ymin),(self.xmin,self.ymax),(self.xmax,self.ymax),(self.xmax,self.ymin)] 
    102         return self._corners 
    103  
    104     def plot(self,color='blue'): 
    105         return plot(polygon(self.corners())) 
    106  
    107     def overlaps_fundamental_domain(self): 
    108         return any([self._parent.in_fundamental_domain(x) for x in self.corners()]) 
    109  
    110     def _repr_(self): 
    111         my_str='Hole with corners=('+str(self.xmin)+','+str(self.ymin)+'),('+str(self.xmax)+','+str(self.ymax)+')' 
    112         return my_str 
    113  
    114 class Region(SageObject): 
    115     def __init__(self,center,radius): 
    116         assert(isinstance(center,tuple)) 
    117         self._center=center 
    118         radius=abs(radius) 
    119         if(radius>1): 
    120             print radius 
    121             print center 
    122             assert(0) 
    123         self._radius=radius 
    124     def radius(self): 
    125         return self._radius 
    126     def center(self): 
    127         return self._center 
    128  
    129     def plot_center(self,color='red'): 
    130         return plot(point(self._center,color=color)) 
    131  
    132     def plot(self,color='red'): 
    133         r=2*self._radius 
    134         a=self._center[0] 
    135         b=self._center[1] 
    136         Pts=[(a,b+r),(a+r,b),(a,b-r),(a-r,b)] 
    137         return plot(polygon(Pts),color=color) 
    138  
    139     def contains_point(self,P): 
    140         return abs((P[0]-self._center[0])*(P[1]-self._center[1]))<self._radius 
    141  
    142     def contains_hole(self,h): 
    143         return all([self.contains_point(P) for P in h.corners()]) 
    144  
    145 class QuadraticContinuedFraction(SageObject): 
    146     def __init__(self,F,Nbound=50,Tbound=5): 
    147         self._F=F 
    148         self._disc=self._F.discriminant() 
    149         self._w=self._F.gen() 
    150         self._eps=self._F.units()[0] 
    151         self._Phi=self._F.real_embeddings(prec=300) 
    152  
    153         # We assume here that Phi orders the embeddings to that 
    154         # Phi[1] is the largest 
    155         self._Rw=self._Phi[1](self._w) 
    156         self._Rwconj=self._Phi[0](self._w) 
    157  
    158         self._used_regions=[] 
    159         self._Nboundmin=2 
    160         self._Nboundmax=Nbound 
    161         self._Nboundinc=1 
    162         self._Nbound=Nbound 
    163         self._epsto=[self._F(1)] 
    164         self._Dmax=self.embed(self._w) 
    165         self._Tbound=Tbound 
    166         self._ranget=sorted(range(-self._Tbound,self._Tbound+2),key=abs) 
    167         self._M=Matrix(RealField(),2,2,[1,self._Rw,1,self._Rwconj]).inverse() 
    168         self.initialize_fundom() 
    169         self._master_regs=Regions(self) 
    170         self._maxdepth=0 
    171     '''Takes an element in the number field F and returns 
    172     y in the fundamental domain such that such that (x-y) belongs to R''' 
    173     def fundom_rep(self,x,y=None): 
    174         if(y==None): 
    175             x1=x[0] 
    176             x2=x[1] 
    177         else: # Assume these are coordinates a,b of a+bw 
    178             x1=x 
    179             x2=y 
    180         d=x1.floor()+self._w*(x2.floor()) 
    181         return x-d 
    182  
    183     # Returns a function that, given a t-value, will return a range of a-values. The arguments specify the region that we want to cover 
    184     def rangea_gen(self,xmin,ymin,xmax,ymax): 
    185         wprime0=self._Rw 
    186         wprime1=self._Rwconj 
    187         def rangea(t): 
    188             tp0=t*wprime0 
    189             tp1=t*wprime1 
    190             return sorted(list(union(range(floor(xmin-1-tp0-wprime0),ceil(xmax-tp0)+1),range(floor(ymin-1-tp1-wprime1),ceil(ymax-tp1)+1))),key=abs) 
    191         return rangea 
    192  
    193     def in_fundamental_domain(self,P): 
    194         A=self._M*Matrix(RealField(),2,1,[P[0],P[1]]) 
    195         return A[0,0]>=0 and A[0,0]<1 and A[1,0]>=0 and A[1,0]<1 
    196  
    197     '''Given x in F (or in R) embeds it into R^2 
    198     ''' 
    199     def embed(self,x): 
    200         return (self._Phi[1](x),self._Phi[0](x)) 
    201  
    202     def embed_coords(self,a,b): 
    203         return ((a+self._Rw*b,a+self._Rwconj*b)) 
    204  
    205     '''Build the fundamental domain 
    206     ''' 
    207     def initialize_fundom(self): 
    208         F=self._F 
    209         Pts=[self.embed(a+b*self._w) for a in [0,1] for b in [0,1]] 
    210         xmin=min([x[0] for x in Pts]) 
    211         xmax=max([x[0] for x in Pts]) 
    212         ymin=min([x[1] for x in Pts]) 
    213         ymax=max([x[1] for x in Pts]) 
    214  
    215         self._aplusbomega=dict([]) 
    216         rangea=self.rangea_gen(xmin-1,ymin-1,xmax+1,ymax+1) 
    217         for b in self._ranget: 
    218             bw=b*self._w 
    219             for a in rangea(b): 
    220                 self._aplusbomega[(a,b)]=a+bw 
    221  
    222         dx=RealField()(0.499) 
    223         dy=dx 
    224         Nx=ceil((xmax-xmin)/dx) 
    225         Ny=ceil((ymax-ymin)/dy) 
    226         self._holes=deque([]) 
    227         for ii in range(Nx): 
    228             for jj in range(Ny): 
    229                 self._holes.append(Hole(self,[xmin+ii*dx,ymin+jj*dy,xmin+(ii+1)*dx,ymin+(jj+1)*dy])) 
    230  
    231  
    232     def plot_holes(self): 
    233         myplot=plot([]) 
    234         for h in self._holes: 
    235             myplot+=h.plot() 
    236         return myplot 
    237  
    238     def remaining_holes(self): 
    239         return len(self._holes) 
    240  
    241     def test_regions(self,corners,at_least_one,embed_coords,regions,Region,aplusbomega,ranget,rangea): 
    242         ii=0 
    243         NumRegions=len(regions) 
    244         #print "Looking at ",NumRegions," regions..." 
    245         atws=[[aplusbomega[(a,b)],embed_coords(a,b),[a,b]] for b in ranget for a in rangea(b)] 
    246         for ii in range(NumRegions): 
    247             reg=regions[ii] 
    248             for inc in atws: 
    249                 passed=[False for ii in range(4)] 
    250                 a=inc[1][0] 
    251                 b=inc[1][1] 
    252                 for ii,P in enumerate(corners): 
    253                     if reg[2].contains_point((P[0]-a,P[1]-b)): 
    254                         at_least_one[ii]=True 
    255                         passed[ii]=True 
    256                     else: 
    257                         break 
    258                 if(all(passed)): 
    259                     rparts=reg[3] 
    260                     incparts=inc[2] 
    261                     return [reg[0]+inc[0],[reg[1][0]+inc[0],reg[1][1]],Region(embed_coords(rparts[0]+incparts[0],rparts[1]+incparts[1]),reg[2].radius()),reg[3],reg[4]] 
    262         return None 
    263  
    264     # Given a hole, returns: 
    265     #  -2 if at least one of the points is not covered by the regions 
    266     #  -1 if no region covers all the points 
    267     #  0 if the a region can cover all the points 
    268     def evaluate_hole(self,h,verifying=False,maxdepth=-1): 
    269         if not h.overlaps_fundamental_domain(): 
    270             return 0 
    271         regs=[[],[],[],[]] 
    272         data=[] 
    273         used=self._used_regions 
    274         corners=h.corners() 
    275         atws=[] 
    276         w=self._w 
    277         if(verifying): 
    278             at_least_one=[False,False,False,False] 
    279  
    280         for reg in used: 
    281             passed=True 
    282             for ii,P in enumerate(corners): 
    283                 if not reg[2].contains_point(P): 
    284                     passed=False 
    285                     if(not verifying): 
    286                         break 
    287                 elif verifying: 
    288                     at_least_one[ii]=True 
    289             if(passed): 
    290                 return 0 
    291  
    292         if(verifying): 
    293             if(all(at_least_one)): 
    294                 return -1 
    295             else: 
    296                 if(h.depth()>maxdepth): 
    297                     return -2 
    298                 else: 
    299                     return -1 
    300  
    301         rangea=self.rangea_gen(h.xmin,h.ymin,h.xmax,h.ymax) 
    302         at_least_one=[False,False,False,False] 
    303  
    304         l=100 
    305         B=ceil(4/(l*(h.size**2))) 
    306         B=max([self._Nboundmin,B]) 
    307         #print 'B=',B,'Size ~ 10^',RealField(prec=10)((log(h.size)/log(10))) 
    308         regions=self._master_regs.get_regions(min([B,self._Nbound])) 
    309         nreg=self.test_regions(corners,at_least_one,embed_coords=self.embed_coords,regions=regions,Region=Region,aplusbomega=self._aplusbomega, ranget=self._ranget,rangea=rangea) 
    310         if(nreg!=None): 
    311             used.append(nreg) 
    312             return 0 
    313  
    314         return -1 
    315  
    316     def evaluate_pointxy(self,x,y): 
    317         P=(x,y) 
    318         for reg in self._used_regions: 
    319             if(reg[2].contains_point(P)): 
    320                 return reg 
    321         return -1 
    322  
    323     def evaluate_number(self,x): 
    324         y=self.fundom_rep(x) 
    325         d=x-y 
    326         P=self.embed(self.fundom_rep(x)) 
    327         for reg in self._used_regions: 
    328             if(reg[2].contains_point(P)): 
    329                 return [reg[1][0]+d,reg[1][1]] 
    330         return -1 
    331  
    332     def verify(self): 
    333         maxdepth=self._maxdepth 
    334         self.initialize_fundom() 
    335         percent=RealField(prec=20)(0) 
    336         holes=self._holes 
    337         self._minholes=len(holes) 
    338         while(len(holes)>0): 
    339             if(len(holes)-self._minholes>1000): 
    340                 assert(0) 
    341             h=holes.popleft() 
    342             result=self.evaluate_hole(h,verifying=True,maxdepth=maxdepth) 
    343             if(result==-2): 
    344                 print "The solution is incorrect." 
    345                 print "The ",h," is not covered by any region." 
    346                 return False 
    347             if(result==-1): 
    348                 xh=(h.xmax+h.xmin)/2 
    349                 yh=(h.ymax+h.ymin)/2 
    350                 depth=h.depth()+1 
    351                 self._maxdepth=max([self._maxdepth,depth]) 
    352                 h1=Hole(self,[h.xmin,h.ymin,xh,yh],depth=depth) 
    353                 h2=Hole(self,[h.xmin,yh,xh,h.ymax],depth=depth) 
    354                 h3=Hole(self,[xh,h.ymin,h.xmax,yh],depth=depth) 
    355                 h4=Hole(self,[xh,yh,h.xmax,h.ymax],depth=depth) 
    356                 holes.extendleft([h1,h2,h3,h4]) 
    357             else: 
    358                 percent+=100*(4**(-h.depth())) 
    359                 #print 'Finished hole ',self._minholes,' up to ',percent,'%' 
    360  
    361             if(len(holes)<self._minholes): 
    362                 self._minholes=len(holes) 
    363                 percent=RealField(prec=20)(0) 
    364                 #print "Remaining ",self._minholes,' holes.' 
    365         return True 
    366  
    367     def solve(self,verify=False): 
    368         percent=RealField(prec=20)(0) 
    369         holes=self._holes 
    370         self._minholes=len(holes) 
    371         while(len(holes)>0): 
    372             if(len(holes)-self._minholes>1000): 
    373                 assert(0) 
    374             h=holes.popleft() 
    375             result=self.evaluate_hole(h) 
    376             if(result==-2): 
    377                 self._Nboundmax=max([self._Nboundmax,self._Nbound]) 
    378                 self._holes.append(h) 
    379             if(result==-1): 
    380                 xh=(h.xmax+h.xmin)/2 
    381                 yh=(h.ymax+h.ymin)/2 
    382                 depth=h.depth()+1 
    383                 self._maxdepth=max([self._maxdepth,depth]) 
    384                 h1=Hole(self,[h.xmin,h.ymin,xh,yh],depth=depth) 
    385                 h2=Hole(self,[h.xmin,yh,xh,h.ymax],depth=depth) 
    386                 h3=Hole(self,[xh,h.ymin,h.xmax,yh],depth=depth) 
    387                 h4=Hole(self,[xh,yh,h.xmax,h.ymax],depth=depth) 
    388  
    389                 holes.extendleft([h1,h2,h3,h4]) 
    390             else: 
    391                 percent+=100*(4**(-h.depth())) 
    392                 #print 'Finished hole ',self._minholes,' up to ',percent,'%' 
    393  
    394             if(len(holes)<self._minholes): 
    395                 self._minholes=len(holes) 
    396                 percent=RealField(prec=20)(0) 
    397                 #print "Remaining ",self._minholes,' holes.' 
    398         if(verify==True): 
    399             assert(self.verify()) 
    400         return [self._disc,self._maxdepth,1+max([abs(self._F(reg[1][1]).norm()) for reg in self._used_regions]),len(self._used_regions),self._used_regions] 
    401  

  • new file sage/rings/number_field/quadratic_continued_fraction.pyx 

    diff -r b97e6e8ba4a1 -r 0c92a2edd152 sage/rings/number_field/quadratic_continued_fraction.pyx
    - +  
     1 
     2########################################## 
     3# 
     4# Added by mmasdeu 
     5# Authors: Xevi Guitart and Marc Masdeu 
     6# 
     7########################################## 
     8 
     9from sage.matrix.constructor import Matrix 
     10from sage.rings.all import RealField 
     11from sage.structure.sage_object import SageObject 
     12from sage.misc.all import union 
     13from sage.functions.all import (floor,ceil) 
     14from sage.rings.all import IntegerRing 
     15from sage.rings.all import RationalField 
     16from collections import deque 
     17from collections import namedtuple 
     18 
     19 
     20 
     21class Regions(SageObject): 
     22    def __init__(self,parent): 
     23        self._regions=dict([]) 
     24        self._F=parent._F 
     25        self._eps=parent._eps 
     26        self._epsto=[self._F(1)] 
     27        self.fundom_rep=parent.fundom_rep 
     28        self.embed=parent.embed 
     29        self._parent=parent 
     30        self._N=1 # Means that we have regions up to _N-1=0, so none 
     31        self._area=self._parent._Rw 
     32 
     33    def get_regions(self,B): 
     34        regions=self._regions 
     35        if(B>self._N): 
     36            F=self._F 
     37            eps=self._eps 
     38            newelts=dict([(IntegerRing()(n),[I.gens_reduced()[0] for I in v if I.is_principal()]) for n,v in self._F.ideals_of_bdd_norm(B-1).iteritems() if n>=self._N]) 
     39            self._N=B 
     40            for nn in range(len(self._epsto),B): 
     41                self._epsto.append(self._epsto[nn-1]*eps) 
     42            for nn,v in newelts.iteritems(): 
     43                assert not regions.has_key(nn) 
     44                regions[nn]=[] 
     45                for q in v: 
     46                    invden=1/q 
     47                    invden_red=self.fundom_rep(invden) 
     48                    q0=invden_red-invden 
     49                    m_invden_red=self.fundom_rep(-invden) 
     50                    m_q0=m_invden_red+invden 
     51 
     52                    a,b=self._parent._change_basis(invden_red) 
     53                    am,bm=self._parent._change_basis(m_invden_red) 
     54 
     55                    regions[nn].extend([[invden_red,[q0,q],Region(self.embed(invden_red),RealField()(1/nn)),[a,b],nn],[m_invden_red,[m_q0,-q],Region(self.embed(m_invden_red),RealField()(1/nn)),[am,bm],nn]]) 
     56 
     57                    for jj in range(1,nn): 
     58                        x=self._epsto[jj]*invden 
     59                        x_red=self.fundom_rep(x) 
     60                        x_red_minus=self.fundom_rep(-x_red) 
     61                        regs=regions[nn] 
     62                        if(x_red==regs[0][0] or x_red_minus==regs[0][0]): 
     63                            continue 
     64 
     65                        q0_minus=x_red_minus+x 
     66                        q0=x_red-x 
     67                        q1=1/x 
     68                        a,b=self._parent._change_basis(x_red) 
     69                        am,bm=self._parent._change_basis(x_red_minus) 
     70                        regions[nn].extend([[x_red,[q0,q1],Region(self.embed(x_red),RealField()(1/nn)),[a,b],nn],[x_red_minus,[q0_minus,-q1],Region(self.embed(x_red_minus),RealField()(1/nn)),[am,bm],nn]]) 
     71 
     72        else: 
     73            # We have all the required regions 
     74            pass 
     75        return [reg for nn,lst in regions.iteritems() if nn<B for reg in lst] 
     76 
     77class Hole(SageObject): 
     78    def __init__(self,parent,v,depth=0): 
     79        self._parent=parent 
     80        if isinstance(v,list): 
     81            #Input in as [xmin,ymin,xmax,ymax] 
     82            if not (v[0]<v[2] and v[1]<v[3]): 
     83                print 'The input should be of the form [xmin,ymin,xmax,ymax], with xmin<xmax and ymin<ymax.' 
     84                assert(0) 
     85            (self.xmin,self.ymin,self.xmax,self.ymax)=v 
     86        elif(isinstance(v,Hole)): 
     87            self.xmin=v.xmin 
     88            self.ymin=v.ymin 
     89            self.xmax=v.xmax 
     90            self.ymax=v.ymax 
     91        self.size=max([self.xmax-self.xmin,self.ymax-self.ymin]) 
     92        self._corners=None 
     93        self._depth=depth 
     94 
     95    def depth(self): 
     96        return self._depth 
     97 
     98    def corners(self): 
     99        if self._corners==None: 
     100            Pt=self._parent.Pointxy 
     101            self._corners=[Pt(self.xmin,self.ymin),Pt(self.xmin,self.ymax),Pt(self.xmax,self.ymax),Pt(self.xmax,self.ymin)] 
     102        return self._corners 
     103 
     104    def plot(self,color='blue'): 
     105        return plot(polygon(self.corners())) 
     106 
     107    def overlaps_fundamental_domain(self): 
     108        return any([self._parent.in_fundamental_domain(x) for x in self.corners()]) 
     109 
     110    def _repr_(self): 
     111        my_str='Hole with corners=('+str(self.xmin)+','+str(self.ymin)+'),('+str(self.xmax)+','+str(self.ymax)+')' 
     112        return my_str 
     113 
     114class Region(SageObject): 
     115    def __init__(self,center,radius): 
     116        #assert(isinstance(center,Pointxy)) 
     117        self._center=center 
     118        radius=abs(radius) 
     119        if(radius>1): 
     120            print radius 
     121            print center 
     122            assert(0) 
     123        self._radius=radius 
     124    def radius(self): 
     125        return self._radius 
     126    def center(self): 
     127        return self._center 
     128 
     129    def plot_center(self,color='red'): 
     130        return plot(point(self._center,color=color)) 
     131 
     132    def plot(self,color='red'): 
     133        r=2*self._radius 
     134        a=self._center[0] 
     135        b=self._center[1] 
     136        Pts=[(a,b+r),(a+r,b),(a,b-r),(a-r,b)] 
     137        return plot(polygon(Pts),color=color) 
     138 
     139    def contains_point(self,P): 
     140        return abs((P.x-self._center.x)*(P.y-self._center.y))<self._radius 
     141 
     142    def contains_hole(self,h): 
     143        return all([self.contains_point(P) for P in h.corners()]) 
     144 
     145class QuadraticContinuedFraction(SageObject): 
     146    def __init__(self,F,Nbound=50,Tbound=5): 
     147        self._F=F 
     148        self._solved=False 
     149        self._disc=self._F.discriminant() 
     150        self._w=self._F.maximal_order().ring_generators()[0] 
     151        self._eps=self._F.units()[0] 
     152        self._Phi=self._F.real_embeddings(prec=300) 
     153        a=F.gen() 
     154        self.Pointxy=namedtuple('Pointxy','x y') 
     155        if(self._disc%2==0): 
     156            self._changebasismatrix=1 
     157        else: 
     158            self._changebasismatrix=Matrix(IntegerRing(),2,2,[1,-1,0,2]) 
     159 
     160        # We assume here that Phi orders the embeddings to that 
     161        # Phi[1] is the largest 
     162        self._Rw=self._Phi[1](self._w) 
     163        self._Rwconj=self._Phi[0](self._w) 
     164 
     165        self._used_regions=[] 
     166        self._Nboundmin=2 
     167        self._Nboundmax=Nbound 
     168        self._Nboundinc=1 
     169        self._Nbound=Nbound 
     170        self._epsto=[self._F(1)] 
     171        self._Dmax=self.embed(self._w) 
     172        self._Tbound=Tbound 
     173        self._ranget=sorted(range(-self._Tbound,self._Tbound+2),key=abs) 
     174        self._M=Matrix(RealField(),2,2,[1,self._Rw,1,self._Rwconj]).inverse() 
     175        self.initialize_fundom() 
     176        self._master_regs=Regions(self) 
     177        self._maxdepth=0 
     178    '''Takes an element in the number field F and returns 
     179    y in the fundamental domain such that such that (x-y) belongs to R''' 
     180    def fundom_rep(self,x): 
     181        # Now xp contains a,b such that x=a+b sqrt(D) 
     182        y00,y10=self._change_basis(x) 
     183        d=y00.floor()+self._w*y10.floor() 
     184        return x-d 
     185 
     186    def _change_basis(self,x): 
     187        v=x.parts() 
     188        tmp=self._changebasismatrix*Matrix(RationalField(),2,1,v) 
     189        return tmp[0,0],tmp[1,0] 
     190 
     191    # Returns a function that, given a t-value, will return a range of a-values. The arguments specify the region that we want to cover 
     192    def rangea_gen(self,xmin,ymin,xmax,ymax): 
     193        wprime0=self._Rw 
     194        wprime1=self._Rwconj 
     195        def rangea(t): 
     196            tp0=t*wprime0 
     197            tp1=t*wprime1 
     198            return sorted(list(union(range(floor(xmin-1-tp0-wprime0),ceil(xmax-tp0)+1),range(floor(ymin-1-tp1-wprime1),ceil(ymax-tp1)+1))),key=abs) 
     199        return rangea 
     200 
     201    def in_fundamental_domain(self,P): 
     202        #assert(isinstance(P,Pointxy)) 
     203        A=self._M*Matrix(RealField(),2,1,[P.x,P.y]) 
     204        return A[0,0]>=0 and A[0,0]<1 and A[1,0]>=0 and A[1,0]<1 
     205 
     206    '''Given x in F (or in R) embeds it into R^2 
     207    ''' 
     208    def embed(self,x): 
     209        return self.Pointxy(self._Phi[1](x),self._Phi[0](x)) 
     210 
     211    def embed_coords(self,a,b): 
     212        return self.Pointxy(a+b*self._Rw,a+b*self._Rwconj) 
     213        #return ((ZZ(a)+self._Rw*ZZ(b),ZZ(a)+self._Rwconj*ZZ(b))) 
     214 
     215    '''Build the fundamental domain 
     216    ''' 
     217    def initialize_fundom(self): 
     218        F=self._F 
     219        Pts=[self.embed_coords(a,b) for a in [0,1] for b in [0,1]] 
     220        xmin=min([P.x for P in Pts]) 
     221        xmax=max([P.x for P in Pts]) 
     222        ymin=min([P.y for P in Pts]) 
     223        ymax=max([P.y for P in Pts]) 
     224 
     225        self._aplusbomega=dict([]) 
     226        rangea=self.rangea_gen(xmin-1,ymin-1,xmax+1,ymax+1) 
     227        for b in self._ranget: 
     228            bw=b*self._w 
     229            for a in rangea(b): 
     230                self._aplusbomega[(a,b)]=a+bw 
     231 
     232        dx=RealField()(0.49) 
     233        dy=dx 
     234        Nx=ceil((xmax-xmin)/dx) 
     235        Ny=ceil((ymax-ymin)/dy) 
     236        self._holes=deque([]) 
     237        for ii in range(Nx): 
     238            for jj in range(Ny): 
     239                self._holes.append(Hole(self,[xmin+ii*dx,ymin+jj*dy,xmin+(ii+1)*dx,ymin+(jj+1)*dy])) 
     240 
     241 
     242    def plot_holes(self): 
     243        myplot=plot([]) 
     244        for h in self._holes: 
     245            myplot+=h.plot() 
     246        return myplot 
     247 
     248    def remaining_holes(self): 
     249        return len(self._holes) 
     250 
     251    def test_regions(self,corners,at_least_one,embed_coords,regions,Region,aplusbomega,ranget,rangea): 
     252        ii=0 
     253        NumRegions=len(regions) 
     254        #print "Looking at ",NumRegions," regions..." 
     255        atws=[[aplusbomega[(a,b)],embed_coords(a,b),[a,b]] for b in ranget for a in rangea(b)] 
     256        for ii in range(NumRegions): 
     257            reg=regions[ii] 
     258            for inc in atws: 
     259                passed=[False for ii in range(4)] 
     260                a=inc[1].x 
     261                b=inc[1].y 
     262                for ii,P in enumerate(corners): 
     263                    if reg[2].contains_point(self.Pointxy(P.x-a,P.y-b)): 
     264                        at_least_one[ii]=True 
     265                        passed[ii]=True 
     266                    else: 
     267                        break 
     268                if(all(passed)): 
     269                    rparts=reg[3] 
     270                    incparts=inc[2] 
     271                    return [reg[0]+inc[0],[reg[1][0]+inc[0],reg[1][1]],Region(embed_coords(rparts[0]+incparts[0],rparts[1]+incparts[1]),reg[2].radius()),reg[3],reg[4]] 
     272        return None 
     273 
     274    # Given a hole, returns: 
     275    #  -2 if at least one of the points is not covered by the regions 
     276    #  -1 if no region covers all the points 
     277    #  0 if the a region can cover all the points 
     278    def evaluate_hole(self,h,verifying=False,maxdepth=-1): 
     279        if not h.overlaps_fundamental_domain(): 
     280            return 0 
     281        regs=[[],[],[],[]] 
     282        data=[] 
     283        used=self._used_regions 
     284        corners=h.corners() 
     285        atws=[] 
     286        w=self._w 
     287        if(verifying): 
     288            at_least_one=[False,False,False,False] 
     289 
     290        for reg in used: 
     291            passed=True 
     292            for ii,P in enumerate(corners): 
     293                if not reg[2].contains_point(P): 
     294                    passed=False 
     295                    if(not verifying): 
     296                        break 
     297                elif verifying: 
     298                    at_least_one[ii]=True 
     299            if(passed): 
     300                return 0 
     301 
     302        if(verifying): 
     303            if(all(at_least_one)): 
     304                return -1 
     305            else: 
     306                if(h.depth()>maxdepth): 
     307                    return -2 
     308                else: 
     309                    return -1 
     310 
     311        rangea=self.rangea_gen(h.xmin,h.ymin,h.xmax,h.ymax) 
     312        at_least_one=[False,False,False,False] 
     313 
     314        l=100 
     315        B=ceil(4/(l*(h.size**2))) 
     316        B=max([self._Nboundmin,B]) 
     317        #print 'B=',B,'Size ~ 10^',RealField(prec=10)((log(h.size)/log(10))) 
     318        regions=self._master_regs.get_regions(min([B,self._Nbound])) 
     319        nreg=self.test_regions(corners,at_least_one,embed_coords=self.embed_coords,regions=regions,Region=Region,aplusbomega=self._aplusbomega, ranget=self._ranget,rangea=rangea) 
     320        if(nreg!=None): 
     321            used.append(nreg) 
     322            return 0 
     323 
     324        return -1 
     325 
     326    def evaluate_pointxy(self,x,y): 
     327        P=self.Pointxy(x,y) 
     328        for reg in self._used_regions: 
     329            if(reg[2].contains_point(P)): 
     330                return reg 
     331        return -1 
     332 
     333    def evaluate_number(self,x): 
     334        y=self.fundom_rep(x) 
     335        d=x-y 
     336        P=self.embed(y) 
     337        for reg in self._used_regions: 
     338            if(reg[2].contains_point(P)): 
     339                return [reg[1][0]+d,reg[1][1]] 
     340        return -1 
     341 
     342    def verify(self): 
     343        maxdepth=self._maxdepth 
     344        self.initialize_fundom() 
     345        percent=RealField(prec=20)(0) 
     346        holes=self._holes 
     347        self._minholes=len(holes) 
     348        while(len(holes)>0): 
     349            if(len(holes)-self._minholes>1000): 
     350                assert(0) 
     351            h=holes.popleft() 
     352            result=self.evaluate_hole(h,verifying=True,maxdepth=maxdepth) 
     353            if(result==-2): 
     354                print "The solution is incorrect." 
     355                print "The ",h," is not covered by any region." 
     356                return False 
     357            if(result==-1): 
     358                xh=(h.xmax+h.xmin)/2 
     359                yh=(h.ymax+h.ymin)/2 
     360                depth=h.depth()+1 
     361                self._maxdepth=max([self._maxdepth,depth]) 
     362                h1=Hole(self,[h.xmin,h.ymin,xh,yh],depth=depth) 
     363                h2=Hole(self,[h.xmin,yh,xh,h.ymax],depth=depth) 
     364                h3=Hole(self,[xh,h.ymin,h.xmax,yh],depth=depth) 
     365                h4=Hole(self,[xh,yh,h.xmax,h.ymax],depth=depth) 
     366                holes.extendleft([h1,h2,h3,h4]) 
     367            else: 
     368                percent+=100*(4**(-h.depth())) 
     369                #print 'Finished hole ',self._minholes,' up to ',percent,'%' 
     370 
     371            if(len(holes)<self._minholes): 
     372                self._minholes=len(holes) 
     373                percent=RealField(prec=20)(0) 
     374                #print "Remaining ",self._minholes,' holes.' 
     375        return True 
     376 
     377    def solve(self,verify=False): 
     378        percent=RealField(prec=20)(0) 
     379        holes=self._holes 
     380        self._minholes=len(holes) 
     381        while(len(holes)>0): 
     382            if(len(holes)-self._minholes>1000): 
     383                assert(0) 
     384            h=holes.popleft() 
     385            result=self.evaluate_hole(h) 
     386            if(result==-2): 
     387                self._Nboundmax=max([self._Nboundmax,self._Nbound]) 
     388                self._holes.append(h) 
     389            if(result==-1): 
     390                xh=(h.xmax+h.xmin)/2 
     391                yh=(h.ymax+h.ymin)/2 
     392                depth=h.depth()+1 
     393                self._maxdepth=max([self._maxdepth,depth]) 
     394                h1=Hole(self,[h.xmin,h.ymin,xh,yh],depth=depth) 
     395                h2=Hole(self,[h.xmin,yh,xh,h.ymax],depth=depth) 
     396                h3=Hole(self,[xh,h.ymin,h.xmax,yh],depth=depth) 
     397                h4=Hole(self,[xh,yh,h.xmax,h.ymax],depth=depth) 
     398 
     399                holes.extendleft([h1,h2,h3,h4]) 
     400            else: 
     401                percent+=100*(4**(-h.depth())) 
     402                #print 'Finished hole ',self._minholes,' up to ',percent,'%' 
     403 
     404            if(len(holes)<self._minholes): 
     405                self._minholes=len(holes) 
     406                percent=RealField(prec=20)(0) 
     407                #print "Remaining ",self._minholes,' holes.' 
     408        if(verify==True): 
     409            assert(self.verify()) 
     410        self._solved=True 
     411        return [self._disc,self._maxdepth,1+max([abs(self._F(reg[1][1]).norm()) for reg in self._used_regions]),len(self._used_regions),self._used_regions] 
     412