attachment:NumerationClass.sage of days72

   1 #################################################################
   2 # NumerationClass.sage
   3 # Author: Attila Kovács
   4 # http://compalg.inf.elte.hu/~attila
   5 # Last modification: June 5, 2016
   6 #################################################################
   7 
   8 ##### Exception Handling #####
   9 class NumberSystemException(Exception):
  10     def __init__(self, value):
  11         self.value = value
  12     def __str__(self):
  13         return repr(self.value)
  14         
  15 class NumberSystemFullResidueSystemException(NumberSystemException):
  16     def __init__(self, value):
  17         self.value = value
  18         
  19 class NumberSystemExpansivityException(NumberSystemException):
  20     def __init__(self, value):
  21         self.value = value
  22         
  23 class NumberSystemUnitConditionException(NumberSystemException):
  24     def __init__(self, value):
  25         self.value = value
  26         
  27 class NumberSystemRegularityException(NumberSystemException):
  28     def __init__(self, value):
  29         self.value = value
  30         
  31 ##### The Number System Class #####
  32 class NumberSystem:
  33     ### Checking the Expansivity
  34     def checkExpansivity(self):
  35         for i in [abs(p) for p in self.base.eigenvalues()]:
  36             if i <= 1:
  37                 return false
  38         return true
  39     
  40     ### Checking the Unit Condition
  41     def checkUnitCondition(self):
  42         cp = self.base.charpoly()
  43         if abs(cp.subs(x=1)) == 1:
  44             return false
  45         return true
  46     
  47     ### Checking the Full Residue System Property and building the hashes of the digits
  48     def checkCRSPropertyAndBuildDigitsHashes(self):
  49         if len(self.digits) <> self.absDeterminant:
  50             raise NumberSystemFullResidueSystemException("The digit set must be a full residue system, it should have |det(M)| elements...")
  51         digitsList = []
  52         self.digitsHash = []
  53         for v in self.digits:
  54             res = 0
  55     	    i = self.dimension-1
  56     	    while i >= 0 and self.smithdiagonal[i] > 1:
  57                 s = 0
  58                 for j in range(self.dimension):
  59             	    s = s + (self.smU[i][j]*v[j] % self.smithdiagonal[i])
  60                 res = res * self.smithdiagonal[i] + (s % self.smithdiagonal[i])
  61                 i = i-1
  62             if res in digitsList:
  63                 raise NumberSystemFullResidueSystemException("The digit set must be a full residue system, there are congruent elements...")
  64             else:
  65         	    digitsList.append(res)
  66         for i in range(len(digitsList)):
  67             self.digitsHash.append(self.digits[digitsList.index(i)])
  68         return true
  69     
  70     ### Computes the symmetric modulus of a number
  71     def getSymm_mod(self, num, mod):
  72         p=Mod(num,mod)
  73         r = p.lift()
  74         if r*2 > p.modulus():
  75             r -= p.modulus()
  76         return r
  77     
  78     ### Computes the congruent class of a given vector v using the Adjoint method
  79     def getAdjCongruentClass(self, v):
  80         v1 = []
  81         for i in range(self.dimension):
  82             s = 0
  83             for j in range(self.dimension):
  84                 s = (s + self.adjointM[i,j]*v[j])
  85             s = self.getSymm_mod(s, self.base.det())
  86             v1.append(s)
  87         return v1
  88     
  89     ### Generates the j-canonical type digit set
  90     def generateCanonicalDigits(self, j):
  91         return [[a if b==j-1 else 0 for b in range(self.dimension)] for a in range(self.base.det().abs())]
  92     
  93     ### Generates the j-symmetric type digit set
  94     def generateSymmetricDigits(self, j):
  95         t=self.base.det().abs()
  96         return [[a-floor(t/2) if b==j-1 else 0 for b in range(self.dimension)] for a in range(t)]
  97     
  98     ### Generates the adjoint type digit set
  99     def generateAdjointDigits(self):
 100         ### first generating a complete residue system to cr_set
 101         insmU = self.smU.inverse()
 102         cr_set = []
 103         v = [0]*self.dimension
 104         j = 0
 105         finished = 0
 106         while finished == 0 and j < self.absDeterminant:
 107             cr_set.append((insmU * vector(v)).list())
 108             i = self.dimension-1
 109             while i >= 0 and v[i] == self.smithdiagonal[i] - 1:
 110                 v[i] = 0
 111                 i = i-1
 112             if i >= 0:
 113                 v[i] = v[i] + 1
 114                 j = j+1
 115             else:
 116                 finished = 1
 117         ### second producing the adjoint type complete residue system
 118         Bs = []
 119         for i in cr_set:
 120             if i == [0]*self.dimension:
 121                 Bs.append(i)
 122             else:
 123                 Bs.append([x/self.base.det() for x in self.base*vector(self.getAdjCongruentClass(i))])
 124         return(Bs)
 125     
 126     ### Computes the congruent element for a given point v
 127     def getCongruentElement(self,v):
 128         res = 0
 129         i = self.dimension-1
 130         while i>=0 and self.smithdiagonal[i] > 1:
 131             s = 0
 132             for j in range(self.dimension):
 133                 s = s + self.smU[i][j]*v[j]
 134             res = res * self.smithdiagonal[i] + (s % self.smithdiagonal[i])
 135             i = i-1
 136         return self.digitsHash[res]
 137     
 138     ### Computes the Phi function for a given point v
 139     def getPhi_fun(self,v):
 140         return((self.inverseM * (vector(v)-vector(self.getCongruentElement(v)))).list())
 141     
 142     ### Computes the Phi function for a given point v and gives back the congruent element as well
 143     def getPhi_fun_withDigit(self,v):
 144         w = self.getCongruentElement(v)
 145         return((self.inverseM * (vector(v)-vector(w))).list(),w)
 146     
 147     ### Computes the orbit from the starting point v
 148     def getOrbitFrom(self,v):
 149         R = [v]
 150         v1 = v
 151         v2 = self.getPhi_fun(v1)
 152         while not v2 in R:
 153             R.append(v2);
 154             v1 = v2
 155             v2 = self.getPhi_fun(v1);
 156         R.append(v2)
 157         return(R)
 158     
 159     ### Computes the set of points in the fraction set for plotting
 160     ### iterNum - the number of iterations, default is 7
 161     ### rgbcolor - the color, default is black
 162     ### flag - if it is 1 then the set of points of fractions are computed, if it is -1 then the opposite set, default is -1
 163     def getFractionsSetPlot(self, iterNum=7,rgbcolor=(0, 0, 0),flag=-1):
 164         K = [[0,0]]
 165         for i in range(1,iterNum):
 166             oldK = K[:]
 167             for d in self.digits:
 168                 for k in oldK:
 169                     newPoint = (self.inverseM^i*vector(d) + vector(k)).list()
 170                     if newPoint not in K:
 171                         K.append(newPoint)
 172         K = [(flag * vector(k)).list() for k in K]
 173         return points(K,rgbcolor=rgbcolor)
 174     
 175     ### Computes the covering box of the fraction set
 176     ### The output is the coordinates of the lower and the upper corners of the box
 177     ### info - output to the screen, which can be True or False
 178     def getCoveringBox(self, info=True, eps=0.01):
 179         X = matrix.identity(self.dimension)
 180         v1 = [0]*self.dimension; v2 = [0]*self.dimension; v3 = [0]*self.dimension; v4 = [0]*self.dimension
 181         while X.norm(Infinity) >= eps:
 182             X = X * self.inverseM
 183             l = [X * vector(i) for i in self.digits]
 184             for i in range(self.dimension):
 185                 y = 0; z = 0
 186                 for j in l:
 187                     y=min(j[i],y)
 188                     z=max(j[i],z)
 189                 v1[i] = y
 190                 v3[i] = z
 191             v2 = [v1[i]+v2[i] for i in range(len(v1))]
 192             v4 = [v3[i]+v4[i] for i in range(len(v3))]
 193         tempMultiplier = 1/(1-X.norm(Infinity))
 194         v3 = [-ceil(x * tempMultiplier) for x in v2]
 195         v1 = [-floor(x * tempMultiplier) for x in v4]
 196         if info == True:
 197             print "Vertices of the covering Box of the periodic points: "
 198         return (v1,v3)
 199     
 200     ### The decision method
 201     def decideMethodA(self):
 202         ##### Checking the points in the Box
 203         if (self.lowBox == [0]*self.dimension and self.upBox == [0]*self.dimension):
 204             self.lowBox,self.upBox = self.getCoveringBox(info=False)
 205         v = self.lowBox[:]
 206         while true:
 207             last = self.getOrbitFrom(v)[-1]
 208             if last != [0]*self.dimension:
 209                 print "Not GNS, one witness is: ", last
 210                 return
 211             i = 0
 212             while i < self.dimension and v[i] == self.upBox[i]:
 213                 v[i] = self.lowBox[i]
 214                 i = i + 1
 215             if i < self.dimension:
 216                 v[i] = v[i] + 1
 217             else:
 218                 print("Numeration System")
 219                 return
 220             
 221     ### The classification method of type A
 222     ### Gives back all the cycles
 223     def classifyMethodA(self):
 224         ##### Generating the points in the box
 225         if (self.lowBox == [0]*self.dimension and self.upBox == [0]*self.dimension):
 226             self.lowBox,self.upBox = self.getCoveringBox(info=False)
 227         v = self.lowBox[:]
 228         end = 0
 229         G = [];
 230         while end == 0:
 231             G.append(v[:])
 232             i = 0
 233             while i < self.dimension and v[i] == self.upBox[i]:
 234                 v[i] = self.lowBox[i]
 235                 i = i + 1
 236             if i < self.dimension:
 237                 v[i] = v[i]+1
 238             else:
 239                 end = 1
 240         ##### Determining the cycles
 241         Li = []
 242         while len(G) > 0:
 243             l = []
 244             a = G[0];
 245             while a in G:
 246                 G.remove(a)
 247                 l.append(a)
 248                 v = self.getPhi_fun(a)
 249                 a = v
 250             if a in l:
 251                 l.append(a)
 252                 Li.append(l[l.index(a):len(l)])
 253         print "The classification is: "
 254         return Li
 255     
 256     ### Constructor
 257     ### m - the operator
 258     ### digits - the list of digits
 259     ### jDim - for canonical or symmetric digits set construction
 260     ### safeInit - for safe initialization, default is false
 261     def __init__(self, m, digits, jDim=1, safeInit=false):
 262         self.base = m
 263         self.absDeterminant = abs(self.base.det())
 264         
 265         if self.absDeterminant == 0:
 266             raise NumberSystemRegularityException("The operator must be regular")
 267         if self.checkUnitCondition() == False:
 268             raise NumberSystemUnitConditionException("abs(det(M-I)) must be greater than one")
 269         if self.checkExpansivity() == False:
 270             raise NumberSystemExpansivityException("The operator must be expansive")
 271             
 272         self.dimension = self.base.nrows()
 273         self.inverseM = self.base.inverse()
 274         self.adjointM = self.base.adjoint()
 275         self.sm,self.smU,self.smV = self.base.smith_form()
 276         self.smithdiagonal = [self.sm[i][i] for i in range(self.sm.nrows())]
 277         self.lowBox = [0]*self.dimension; self.upBox = [0]*self.dimension
 278         
 279         if (digits=='symmetric' and jDim > 0 and jDim <= self.dimension):
 280             self.digits=self.generateSymmetricDigits(jDim)
 281         elif (digits=='canonical' and jDim > 0 and jDim <= self.dimension):
 282             self.digits=self.generateCanonicalDigits(jDim)
 283         elif digits=="adjoint":
 284             self.digits=self.generateAdjointDigits()
 285         else:
 286             self.digits = digits
 287             
 288         self.checkCRSPropertyAndBuildDigitsHashes()
 289         
 290