From 227a2ac9080fe4a8a8c8f4546acf980838864acb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Chapoton?= Date: Mon, 17 Jun 2024 16:20:46 +0200 Subject: [PATCH 1/3] cleanup for binary recurrence sequences --- .../combinat/binary_recurrence_sequences.py | 378 +++++++++--------- 1 file changed, 181 insertions(+), 197 deletions(-) diff --git a/src/sage/combinat/binary_recurrence_sequences.py b/src/sage/combinat/binary_recurrence_sequences.py index f9c9cdbf2a0..ebe6570cc8b 100644 --- a/src/sage/combinat/binary_recurrence_sequences.py +++ b/src/sage/combinat/binary_recurrence_sequences.py @@ -123,7 +123,6 @@ def __init__(self, b, c, u0=0, u1=1): sage: R = BinaryRecurrenceSequence(1,1) sage: loads(R.dumps()) == R True - """ self.b = b self.c = c @@ -133,7 +132,7 @@ def __init__(self, b, c, u0=0, u1=1): self._PGoodness = {} # dictionary to cache primes that are "good" by some prime power self._ell = 1 # variable that keeps track of the last prime power to be used as a goodness - def __repr__(self): + def __repr__(self) -> str: """ Give string representation of the class. @@ -143,11 +142,10 @@ def __repr__(self): sage: R Binary recurrence sequence defined by: u_n = 3 * u_{n-1} + 3 * u_{n-2}; With initial conditions: u_0 = 2, and u_1 = 1 - """ return 'Binary recurrence sequence defined by: u_n = ' + str(self.b) + ' * u_{n-1} + ' + str(self.c) + ' * u_{n-2};\nWith initial conditions: u_0 = ' + str(self.u0) + ', and u_1 = ' + str(self.u1) - def __eq__(self, other): + def __eq__(self, other) -> bool: """ Compare two binary recurrence sequences. @@ -196,7 +194,7 @@ def __call__(self, n, modulus=0): v = vector(R, [self.u0, self.u1]) return list(F**n * v)[0] - def is_degenerate(self): + def is_degenerate(self) -> bool: """ Decide whether the binary recurrence sequence is degenerate. @@ -243,13 +241,12 @@ def is_degenerate(self): if D.is_square(): A = sqrt(D) else: - K = QuadraticField(D, 'x') - A = K.gen() + A = QuadraticField(D, 'x').gen() - aa = (self.u1 - self.u0*(self.b + A)/2)/(A) #called `a` in Docstring - bb = (self.u1 - self.u0*(self.b - A)/2)/(A) #called `b` in Docstring + aa = (self.u1 - self.u0 * (self.b + A)/2)/(A) # called `a` in Docstring + bb = (self.u1 - self.u0 * (self.b - A)/2)/(A) # called `b` in Docstring - #(b+A)/2 is called alpha in Docstring, (b-A)/2 is called beta in Docstring + # (b+A)/2 is called alpha in Docstring, (b-A)/2 is called beta in Docstring if self.b != A: if ((self.b+A)/(self.b-A))**6 == 1: @@ -261,7 +258,7 @@ def is_degenerate(self): return True - def is_geometric(self): + def is_geometric(self) -> bool: """ Decide whether the binary recurrence sequence is geometric - ie a geometric sequence. @@ -279,14 +276,14 @@ def is_geometric(self): sage: S.is_geometric() True """ - #If [u_0, u_1]^T is an eigenvector for the incrementation matrix F = [[0,1],[c,b]], then the sequence - #is geometric, ie we can write u_n = a*r^n for some a and r. + # If [u_0, u_1]^T is an eigenvector for the incrementation matrix F = [[0,1],[c,b]], then the sequence + # is geometric, ie we can write u_n = a*r^n for some a and r. - #We decide if u0, u1, u2 = b*u1+c*u0 are in geometric progression by whether u1^2 = (b*u1+c*u0)*u0 + # We decide if u0, u1, u2 = b*u1+c*u0 are in geometric progression by whether u1^2 = (b*u1+c*u0)*u0 return (self.u1)**2 == (self.b*self.u1 + self.c*self.u0)*self.u0 - def is_quasigeometric(self): + def is_quasigeometric(self) -> bool: """ Decide whether the binary recurrence sequence is degenerate and similar to a geometric sequence, i.e. the union of multiple geometric sequences, or geometric after term ``u0``. @@ -323,14 +320,13 @@ def is_quasigeometric(self): if D.is_square(): A = sqrt(D) else: - K = QuadraticField(D, 'x') - A = K.gen() + A = QuadraticField(D, 'x').gen() if ((self.b+A)/(self.b-A))**6 == 1: return True return False - def is_arithmetic(self): + def is_arithmetic(self) -> bool: """ Decide whether the sequence is degenerate and an arithmetic sequence. @@ -347,7 +343,7 @@ def is_arithmetic(self): sage: S.is_arithmetic() True """ - return (self(1) - self(0) == self(2) - self(1) == self(3) - self(2)) + return self(1) - self(0) == self(2) - self(1) == self(3) - self(2) def period(self, m): """ @@ -398,7 +394,7 @@ def period(self, m): .. NOTE:: The answer is cached. """ - #If we have already computed the period mod m, then we return the stored value. + # If we have already computed the period mod m, then we return the stored value. if m in self._period_dict: return self._period_dict[m] @@ -410,17 +406,17 @@ def period(self, m): Fac = list(m.factor()) Periods = {} - #To compute the period mod m, we compute the least integer n such that A^n*w == w. This necessarily - #divides the order of A as a matrix in GL_2(Z/mZ). + # To compute the period mod m, we compute the least integer n such that A^n*w == w. This necessarily + # divides the order of A as a matrix in GL_2(Z/mZ). - #We compute the period modulo all distinct prime powers dividing m, and combine via the lcm. - #To compute the period mod p^e, we first compute the order mod p. Then the period mod p^e - #must divide p^{4e-4}*period(p), as the subgroup of matrices mod p^e, which reduce to - #the identity mod p is of order (p^{e-1})^4. So we compute the period mod p^e by successively - #multiplying the period mod p by powers of p. + # We compute the period modulo all distinct prime powers dividing m, and combine via the lcm. + # To compute the period mod p^e, we first compute the order mod p. Then the period mod p^e + # must divide p^{4e-4}*period(p), as the subgroup of matrices mod p^e, which reduce to + # the identity mod p is of order (p^{e-1})^4. So we compute the period mod p^e by successively + # multiplying the period mod p by powers of p. for p, e in Fac: - #first compute the period mod p + # first compute the period mod p if p in self._period_dict: perp = self._period_dict[p] else: @@ -429,29 +425,29 @@ def period(self, m): FF = F**(p-1) p1fac = list((p-1).factor()) - #The order of any matrix in GL_2(F_p) either divides p(p-1) or (p-1)(p+1). - #The order divides p-1 if it is diagonalizable. In any case, det(F^(p-1))=1, - #so if tr(F^(p-1)) = 2, then it must be triangular of the form [[1,a],[0,1]]. - #The order of the subgroup of matrices of this form is p, so the order must divide - #p(p-1) -- in fact it must be a multiple of p. If this is not the case, then the - #order divides (p-1)(p+1). As the period divides the order of the matrix in GL_2(F_p), - #these conditions hold for the period as well. + # The order of any matrix in GL_2(F_p) either divides p(p-1) or (p-1)(p+1). + # The order divides p-1 if it is diagonalizable. In any case, det(F^(p-1))=1, + # so if tr(F^(p-1)) = 2, then it must be triangular of the form [[1,a],[0,1]]. + # The order of the subgroup of matrices of this form is p, so the order must divide + # p(p-1) -- in fact it must be a multiple of p. If this is not the case, then the + # order divides (p-1)(p+1). As the period divides the order of the matrix in GL_2(F_p), + # these conditions hold for the period as well. - #check if the order divides (p-1) + # check if the order divides (p-1) if FF*v == v: M = p-1 Mfac = p1fac - #check if the trace is 2, then the order is a multiple of p dividing p*(p-1) + # check if the trace is 2, then the order is a multiple of p dividing p*(p-1) elif FF.trace() == 2: M = p-1 Mfac = p1fac - F = F**p #replace F by F^p as now we only need to determine the factor dividing (p-1) + F = F**p # replace F by F^p as now we only need to determine the factor dividing (p-1) - #otherwise it will divide (p+1)(p-1) + # otherwise it will divide (p+1)(p-1) else: M = (p+1)*(p-1) - p2fac = list((p+1).factor()) #factor the (p+1) and (p-1) terms separately and then combine for speed + p2fac = list((p+1).factor()) # factor the (p+1) and (p-1) terms separately and then combine for speed Mfac_dic = {} for i0, i1 in list(p1fac + p2fac): if i0 not in Mfac_dic: @@ -460,19 +456,15 @@ def period(self, m): Mfac_dic[i0] += i1 Mfac = list(Mfac_dic.items()) - #Now use a fast order algorithm to compute the period. We know that the period divides - #M = i_1*i_2*...*i_l where the i_j denote not necessarily distinct prime factors. As - #F^M*v == v, for each i_j, if F^(M/i_j)*v == v, then the period divides (M/i_j). After - #all factors have been iterated over, the result is the period mod p. + # Now use a fast order algorithm to compute the period. We know that the period divides + # M = i_1*i_2*...*i_l where the i_j denote not necessarily distinct prime factors. As + # F^M*v == v, for each i_j, if F^(M/i_j)*v == v, then the period divides (M/i_j). After + # all factors have been iterated over, the result is the period mod p. Mfac = list(Mfac) - C = [] - - #expand the list of prime factors so every factor is with multiplicity 1 - for i0, i1 in Mfac: - for _ in range(i1): - C.append(i0) + # expand the list of prime factors so every factor is with multiplicity 1 + C = [i0 for i0, i1 in Mfac for _ in range(i1)] Mfac = C n = M @@ -482,7 +474,7 @@ def period(self, m): n = b perp = n - #Now compute the period mod p^e by stepping up by multiples of p + # Now compute the period mod p^e by stepping up by multiples of p F = A.change_ring(Integers(p**e)) v = w.change_ring(Integers(p**e)) FF = F**perp @@ -498,12 +490,10 @@ def period(self, m): break Periods[p] = perpe - #take the lcm of the periods mod all distinct primes dividing m - period = 1 - for p in Periods: - period = lcm(Periods[p], period) + # take the lcm of the periods mod all distinct primes dividing m + period = lcm(Periods.values()) - self._period_dict[m] = period #cache the period mod m + self._period_dict[m] = period # cache the period mod m return period def pthpowers(self, p, Bound): @@ -569,17 +559,17 @@ def pthpowers(self, p, Bound): This function is primarily optimized in the range where ``Bound`` is much larger than ``p``. """ - #Thanks to Jesse Silliman for helpful conversations! + # Thanks to Jesse Silliman for helpful conversations! - #Reset the dictionary of good primes, as this depends on p + # Reset the dictionary of good primes, as this depends on p self._PGoodness = {} - #Starting lower bound on good primes + # Starting lower bound on good primes self._ell = 1 - #If the sequence is geometric, then the `n`th term is `a*r^n`. Thus the - #property of being a ``p`` th power is periodic mod ``p``. So there are either - #no ``p`` th powers if there are none in the first ``p`` terms, or many if there - #is at least one in the first ``p`` terms. + # If the sequence is geometric, then the `n`th term is `a*r^n`. Thus the + # property of being a ``p`` th power is periodic mod ``p``. So there are either + # no ``p`` th powers if there are none in the first ``p`` terms, or many if there + # is at least one in the first ``p`` terms. if self.is_geometric() or self.is_quasigeometric(): no_powers = True @@ -594,50 +584,50 @@ def pthpowers(self, p, Bound): else: raise ValueError("the degenerate binary recurrence sequence is geometric or quasigeometric and has many pth powers") - #If the sequence is degenerate without being geometric or quasigeometric, there - #may be many ``p`` th powers or no ``p`` th powers. + # If the sequence is degenerate without being geometric or quasigeometric, there + # may be many ``p`` th powers or no ``p`` th powers. elif (self.b**2+4*self.c) == 0: - #This is the case if the matrix F is not diagonalizable, ie b^2 +4c = 0, and alpha/beta = 1. + # This is the case if the matrix F is not diagonalizable, ie b^2 +4c = 0, and alpha/beta = 1. alpha = self.b/2 - #In this case, u_n = u_0*alpha^n + (u_1 - u_0*alpha)*n*alpha^(n-1) = alpha^(n-1)*(u_0 +n*(u_1 - u_0*alpha)), - #that is, it is a geometric term (alpha^(n-1)) times an arithmetic term (u_0 + n*(u_1-u_0*alpha)). + # In this case, u_n = u_0*alpha^n + (u_1 - u_0*alpha)*n*alpha^(n-1) = alpha^(n-1)*(u_0 +n*(u_1 - u_0*alpha)), + # that is, it is a geometric term (alpha^(n-1)) times an arithmetic term (u_0 + n*(u_1-u_0*alpha)). - #Look at classes n = k mod p, for k = 1,...,p. + # Look at classes n = k mod p, for k = 1,...,p. for k in range(1, p + 1): - #The linear equation alpha^(k-1)*u_0 + (k+pm)*(alpha^(k-1)*u1 - u0*alpha^k) - #must thus be a pth power. This is a linear equation in m, namely, A + B*m, where + # The linear equation alpha^(k-1)*u_0 + (k+pm)*(alpha^(k-1)*u1 - u0*alpha^k) + # must thus be a pth power. This is a linear equation in m, namely, A + B*m, where A = (alpha**(k-1)*self.u0 + k*(alpha**(k-1)*self.u1 - self.u0*alpha**k)) B = p*(alpha**(k-1)*self.u1 - self.u0*alpha**k) - #This linear equation represents a pth power iff A is a pth power mod B. + # This linear equation represents a pth power iff A is a pth power mod B. if _is_p_power_mod(A, p, B): raise ValueError("the degenerate binary recurrence sequence has many pth powers") return [] - #We find ``p`` th powers using an elementary sieve. Term `u_n` is a ``p`` th - #power if and only if it is a ``p`` th power modulo every prime `\\ell`. This condition - #gives nontrivial information if ``p`` divides the order of the multiplicative group of - #`\\Bold(F)_{\\ell}`, i.e. if `\\ell` is ` 1 \mod{p}`, as then only `1/p` terms are ``p`` th - #powers modulo `\\ell``. + # We find ``p`` th powers using an elementary sieve. Term `u_n` is a ``p`` th + # power if and only if it is a ``p`` th power modulo every prime `\\ell`. This condition + # gives nontrivial information if ``p`` divides the order of the multiplicative group of + # `\\Bold(F)_{\\ell}`, i.e. if `\\ell` is ` 1 \mod{p}`, as then only `1/p` terms are ``p`` th + # powers modulo `\\ell``. - #Thus, given such an `\\ell`, we get a set of necessary congruences for the index modulo the - #the period of the sequence mod `\\ell`. Then we intersect these congruences for many primes - #to get a tight list modulo a growing modulus. In order to keep this step manageable, we - #only use primes `\\ell` that have particularly smooth periods. + # Thus, given such an `\\ell`, we get a set of necessary congruences for the index modulo the + # the period of the sequence mod `\\ell`. Then we intersect these congruences for many primes + # to get a tight list modulo a growing modulus. In order to keep this step manageable, we + # only use primes `\\ell` that have particularly smooth periods. - #Some congruences in the list will remain as the modulus grows. If a congruence remains through - #7 rounds of increasing the modulus, then we check if this corresponds to a perfect power (if - #it does, we add it to our list of indices corresponding to ``p`` th powers). The rest of the congruences - #are transient and grow with the modulus. Once the smallest of these is greater than the bound, - #the list of known indices corresponding to ``p`` th powers is complete. + # Some congruences in the list will remain as the modulus grows. If a congruence remains through + # 7 rounds of increasing the modulus, then we check if this corresponds to a perfect power (if + # it does, we add it to our list of indices corresponding to ``p`` th powers). The rest of the congruences + # are transient and grow with the modulus. Once the smallest of these is greater than the bound, + # the list of known indices corresponding to ``p`` th powers is complete. else: @@ -656,50 +646,50 @@ def pthpowers(self, p, Bound): for n in range(Bound): # n is the index of the a0 - #Check whether a0 is a perfect power mod ell + # Check whether a0 is a perfect power mod ell if _is_p_power_mod(a0, p, ell): - #if a0 is a perfect power mod ell, check if nth term is ppower + # if a0 is a perfect power mod ell, check if nth term is ppower if _is_p_power(self(n), p): powers.append(n) - a0, a1 = a1, bf*a1 + cf*a0 #step up the variables + a0, a1 = a1, bf*a1 + cf*a0 # step up the variables else: - powers = [] #documents the indices of the sequence that provably correspond to pth powers - cong = [0] #list of necessary congruences on the index for it to correspond to pth powers - Possible_count = {} #keeps track of the number of rounds a congruence lasts in cong + powers = [] # documents the indices of the sequence that provably correspond to pth powers + cong = [0] # list of necessary congruences on the index for it to correspond to pth powers + Possible_count = {} # keeps track of the number of rounds a congruence lasts in cong - #These parameters are involved in how we choose primes to increase the modulus - qqold = 1 #we believe that we know complete information coming from primes good by qqold - M1 = 1 #we have congruences modulo M1, this may not be the tightest list - M2 = p #we want to move to have congruences mod M2 - qq = 1 #the largest prime power divisor of M1 is qq + # These parameters are involved in how we choose primes to increase the modulus + qqold = 1 # we believe that we know complete information coming from primes good by qqold + M1 = 1 # we have congruences modulo M1, this may not be the tightest list + M2 = p # we want to move to have congruences mod M2 + qq = 1 # the largest prime power divisor of M1 is qq - #This loop ups the modulus. + # This loop ups the modulus. while True: - #Try to get good data mod M2 + # Try to get good data mod M2 - #patience of how long we should search for a "good prime" - patience = 0.01 * _estimated_time(lcm(M2, p*next_prime_power(qq)), M1, len(cong), p) + # patience of how long we should search for a "good prime" + patience = 0.01 * _estimated_time(lcm(M2, p * next_prime_power(qq)), + M1, len(cong), p) tries = 0 - #This loop uses primes to get a small set of congruences mod M2. + # This loop uses primes to get a small set of congruences mod M2. while True: - #only proceed if took less than patience time to find the next good prime + # only proceed if took less than patience time to find the next good prime ell = _next_good_prime(p, self, qq, patience, qqold) if ell: - #gather congruence data for the sequence mod ell, which will be mod period(ell) = modu + # gather congruence data for the sequence mod ell, which will be mod period(ell) = modu cong1, modu = _find_cong1(p, self, ell) - CongNew = [] # makes a new list from cong that is now mod M = lcm(M1, modu) instead of M1 + # makes a new list from cong that is now mod M = lcm(M1, modu) instead of M1 M = lcm(M1, modu) - for k in range(M // M1): - for i in cong: - CongNew.append(k * M1 + i) + CongNew = [k * M1 + i for k in range(M // M1) + for i in cong] cong = set(CongNew) M1 = M @@ -728,15 +718,15 @@ def pthpowers(self, p, Bound): cong = list(cong) break - #Document how long each element of cong has been there + # Document how long each element of cong has been there for i in cong: if i in Possible_count: Possible_count[i] += 1 else: Possible_count[i] = 1 - #Check how long each element has persisted, if it is for at least 7 cycles, - #then we check to see if it is actually a perfect power + # Check how long each element has persisted, if it is for at least 7 cycles, + # then we check to see if it is actually a perfect power for i in Possible_count: if Possible_count[i] == 7: n = Integer(i) @@ -744,7 +734,7 @@ def pthpowers(self, p, Bound): if _is_p_power(self(n), p): powers.append(n) - #check for a contradiction + # check for a contradiction if len(cong) > len(powers): if cong[len(powers)] > Bound: break @@ -869,60 +859,56 @@ def _next_good_prime(p, R, qq, patience, qqold): False """ + # We are looking for pth powers in R. + # Our primes must be good by qq, but not qqold. + # We only allow patience number of iterations to find a good prime. - #We are looking for pth powers in R. - #Our primes must be good by qq, but not qqold. - #We only allow patience number of iterations to find a good prime. + # The variable _ell for R keeps track of the last "good" prime returned + # that was not found from the dictionary _PGoodness - #The variable _ell for R keeps track of the last "good" prime returned - #that was not found from the dictionary _PGoodness + # First, we check to see if we have already computed the goodness of a prime that fits + # our requirement of being good by qq but not by qqold. This is stored in the _PGoodness + # dictionary. - #First, we check to see if we have already computed the goodness of a prime that fits - #our requirement of being good by qq but not by qqold. This is stored in the _PGoodness - #dictionary. + # Then if we have, we return the smallest such prime and delete it from the list. If not, we + # search through patience number of primes R._ell to find one good by qq but not qqold. If it is + # not good by either qqold or qq, then we add this prime to R._PGoodness under its goodness. - #Then if we have, we return the smallest such prime and delete it from the list. If not, we - #search through patience number of primes R._ell to find one good by qq but not qqold. If it is - #not good by either qqold or qq, then we add this prime to R._PGoodness under its goodness. + # Possible_Primes keeps track of possible primes satisfying our goodness requirements we might return + # check to see if anything in R._PGoodness fits our goodness requirements + Possible_Primes = [item[0] for j, item in R._PGoodness.items() + if qqold < j <= qq and item] - #Possible_Primes keeps track of possible primes satisfying our goodness requirements we might return - Possible_Primes = [] - - #check to see if anything in R._PGoodness fits our goodness requirements - for j in R._PGoodness: - if (qqold < j <= qq) and len(R._PGoodness[j]): - Possible_Primes.append(R._PGoodness[j][0]) - - #If we found good primes, we take the smallest + # If we found good primes, we take the smallest if Possible_Primes: q = min(Possible_Primes) n = _goodness(q, R, p) - del R._PGoodness[n][0] #if we are going to use it, then we delete it from R._PGoodness + del R._PGoodness[n][0] # if we are going to use it, then we delete it from R._PGoodness return q - #If nothing is already stored in R._PGoodness, we start (from where we left off at R._ell) checking - #for good primes. We only tolerate patience number of tries before giving up. + # If nothing is already stored in R._PGoodness, we start (from where we left off at R._ell) checking + # for good primes. We only tolerate patience number of tries before giving up. else: i = 0 while i < patience: i += 1 R._ell = next_prime(R._ell) - #we require that R._ell is 1 mod p, so that p divides the order of the multiplicative - #group mod R._ell, so that not all elements of GF(R._ell) are pth powers. + # we require that R._ell is 1 mod p, so that p divides the order of the multiplicative + # group mod R._ell, so that not all elements of GF(R._ell) are pth powers. if R._ell % p == 1: - #requiring that b^2 + 4c is a square in GF(R._ell) ensures that the period mod R._ell - #divides R._ell - 1 - if legendre_symbol(R.b**2+4*R.c, R._ell) == 1: + # requiring that b^2 + 4c is a square in GF(R._ell) ensures that the period mod R._ell + # divides R._ell - 1 + if legendre_symbol(R.b**2 + 4 * R.c, R._ell) == 1: N = _goodness(R._ell, R, p) - #proceed only if R._ell satisfies the goodness requirements + # proceed only if R._ell satisfies the goodness requirements if qqold < N <= qq: return R._ell - #if we do not use the prime, we store it in R._PGoodness + # if we do not use the prime, we store it in R._PGoodness else: if N in R._PGoodness: R._PGoodness[N].append(R._ell) @@ -958,19 +944,17 @@ def _is_p_power_mod(a, p, N): False sage: sage.combinat.binary_recurrence_sequences._is_p_power_mod(2**3,3,29) True - """ - - #By the chinese remainder theorem, we can answer this question by examining whether - #a is a pth power mod q^e, for all distinct prime powers q^e dividing N. + # By the chinese remainder theorem, we can answer this question by examining whether + # a is a pth power mod q^e, for all distinct prime powers q^e dividing N. for q, e in N.factor(): - #If a = q^v*x, with + # If a = q^v*x, with v = a.valuation(q) - #then if v>=e, a is congruent to 0 mod q^e and is thus a pth power trivially. + # then if v>=e, a is congruent to 0 mod q^e and is thus a pth power trivially. if v >= e: continue @@ -979,58 +963,58 @@ def _is_p_power_mod(a, p, N): if v % p: return False - #in this cse it is a pth power if x is a pth power mod q^(e-v), so let x = aa, - #and (e-v) = ee: + # in this cse it is a pth power if x is a pth power mod q^(e-v), so let x = aa, + # and (e-v) = ee: - aa = a/q**v + aa = a / q**v ee = e - v - #The above steps are equivalent to the statement that we may assume a and qq are - #relatively prime, if we replace a with aa and e with ee. Now we must determine when - #aa is a pth power mod q^ee for (aa,q)=1. + # The above steps are equivalent to the statement that we may assume a and qq are + # relatively prime, if we replace a with aa and e with ee. Now we must determine when + # aa is a pth power mod q^ee for (aa,q)=1. - #If q != p, then by Hensel's lemma, we may lift a pth power mod q, to a pth power - #mod q^2, etc. + # If q != p, then by Hensel's lemma, we may lift a pth power mod q, to a pth power + # mod q^2, etc. if q != p: - #aa is necessarily a pth power mod q if p does not divide the order of the multiplicative - #group mod q, ie if q is not 1 mod p. + # aa is necessarily a pth power mod q if p does not divide the order of the multiplicative + # group mod q, ie if q is not 1 mod p. if q % p == 1: - #otherwise aa if a pth power mod q iff aa^(q-1)/p == 1 + # otherwise aa if a pth power mod q iff aa^(q-1)/p == 1 - if GF(q)(aa)**((q-1)/p) != 1: + if GF(q)(aa)**((q - 1) / p) != 1: return False - #If q = p and ee = 1, then everything is a pth power p by Fermat's little theorem. + # If q = p and ee = 1, then everything is a pth power p by Fermat's little theorem. elif ee > 1: - #We use the strong statement of Hensel's lemma, which implies that if p is odd - #and aa is a pth power mod p^2, then aa is a pth power mod any higher power of p + # We use the strong statement of Hensel's lemma, which implies that if p is odd + # and aa is a pth power mod p^2, then aa is a pth power mod any higher power of p if p % 2: - #ZZ/(p^2)ZZ^\times is abstractly isomorphic to ZZ/(p)ZZ cross ZZ/(p-1)ZZ. then - #aa is a pth power mod p^2 if (aa)^(p*(p-1)/p) == 1, ie if aa^(p-1) == 1. + # ZZ/(p^2)ZZ^\times is abstractly isomorphic to ZZ/(p)ZZ cross ZZ/(p-1)ZZ. then + # aa is a pth power mod p^2 if (aa)^(p*(p-1)/p) == 1, ie if aa^(p-1) == 1. - if Integers(p**2)(aa)**(p-1) != 1: + if Integers(p**2)(aa)**(p - 1) != 1: return False - #Otherwise, p=2. By the strong statement of Hensel's lemma, if aa is a pth power - #mod p^3, then it is a pth power mod higher powers of p. So we need only check if it - #is a pth power mod p^2 and p^3. + # Otherwise, p=2. By the strong statement of Hensel's lemma, if aa is a pth power + # mod p^3, then it is a pth power mod higher powers of p. So we need only check if it + # is a pth power mod p^2 and p^3. elif ee == 2: - #all odd squares a 1 mod 4 + # all odd squares a 1 mod 4 if aa % 4 != 1: return False - #all odd squares are 1 mod 8 + # all odd squares are 1 mod 8 elif aa % 8 != 1: return False @@ -1061,21 +1045,19 @@ def _estimated_time(M2, M1, length, p): sage: from sage.combinat.binary_recurrence_sequences import _estimated_time sage: _estimated_time(2**4*3**2*5*7*11*13*17, 2**4*3**2*5*7*11*13, 20, 7) # needs sage.symbolic 106.211159309421 - """ + # The heuristic run time of the CRT step to go from modulus M1 to M2 - #The heuristic run time of the CRT step to go from modulus M1 to M2 - - #length is the current length of cong + # length is the current length of cong - Q = p * log(M2) #Size of our primes. - NPrimes = log(M2/M1) / log(Q) #The number of primes + Q = p * log(M2) # Size of our primes. + NPrimes = log(M2 / M1) / log(Q) # The number of primes - return (length * (Q/p)**NPrimes).n() + return (length * (Q / p)**NPrimes).n() -#Find the list of necessary congruences for the index n of binary recurrence -#sequence R using the fact that the reduction mod ell must be a pth power +# Find the list of necessary congruences for the index n of binary recurrence +# sequence R using the fact that the reduction mod ell must be a pth power def _find_cong1(p, R, ell): """ Find the list of permissible indices `n` for which `u_n = y^p` mod ``ell``. @@ -1103,33 +1085,33 @@ def _find_cong1(p, R, ell): u1 = F(R.u1) bf, cf = F(R.b), F(R.c) a0 = u0 - a1 = u1 #a0 and a1 are variables for terms in sequence + a1 = u1 # a0 and a1 are variables for terms in sequence - #The set of pth powers mod ell + # The set of pth powers mod ell PPowers = set(i**p for i in F) - #The period of R mod ell + # The period of R mod ell modu = R.period(ell) - #cong1 keeps track of congruences mod modu for the sequence mod ell + # cong1 keeps track of congruences mod modu for the sequence mod ell cong1 = [] - for n in range(modu): # n is the index of the a0 + for n in range(modu): # n is the index of the a0 - #Check whether a0 is a perfect power mod ell + # Check whether a0 is a perfect power mod ell if a0 in PPowers: - #if a0 is a perfect power mod ell, add the index - #to the list of necessary congruences + # if a0 is a perfect power mod ell, add the index + # to the list of necessary congruences cong1.append(n) - a0, a1 = a1, bf*a1 + cf*a0 #step up the variables + a0, a1 = a1, bf * a1 + cf * a0 # step up the variables cong1.sort() return cong1, modu -def _is_p_power(a, p): +def _is_p_power(a, p) -> bool: """ Determine whether ``a`` is a perfect ``p`` th power. @@ -1141,16 +1123,18 @@ def _is_p_power(a, p): OUTPUT: - - True if ``a`` is a ``p`` th power; else False. + boolean, whether ``a`` is a ``p`` th power EXAMPLES:: - sage: sage.combinat.binary_recurrence_sequences._is_p_power(2**7, 7) # needs sage.symbolic + sage: from sage.combinat.binary_recurrence_sequences import _is_p_power + sage: _is_p_power(2**7, 7) True - sage: sage.combinat.binary_recurrence_sequences._is_p_power(2**7*3**2, 7) # needs sage.symbolic + sage: _is_p_power(2**7*3**2, 7) False """ - return int(a**(1/p))**p == a - # slower tentative ? - # _, test = Integer(a).nth_root(p, truncate_mode=True) - # return test + try: + Integer(a).nth_root(p) + except ValueError: + return False + return True From e050fd5fa9005b810296535ac728abc9cfcd013c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Chapoton?= Date: Mon, 17 Jun 2024 17:23:11 +0200 Subject: [PATCH 2/3] Update binary_recurrence_sequences.py fix typos --- src/sage/combinat/binary_recurrence_sequences.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sage/combinat/binary_recurrence_sequences.py b/src/sage/combinat/binary_recurrence_sequences.py index ebe6570cc8b..f473ce162f2 100644 --- a/src/sage/combinat/binary_recurrence_sequences.py +++ b/src/sage/combinat/binary_recurrence_sequences.py @@ -214,7 +214,7 @@ def is_degenerate(self) -> bool: - `F` is nondiagonalizable -- this corresponds to `\\alpha = \\beta`. This sequence will be the point-wise product of an arithmetic and geometric sequence. - - `F^k` is scaler, for some `k>1` -- this corresponds to `\\alpha/\\beta` a `k` th root of unity. This sequence is a union of several geometric sequences, and so we again call it ``quasigeometric``. + - `F^k` is scalar, for some `k>1` -- this corresponds to `\\alpha/\\beta` a `k` th root of unity. This sequence is a union of several geometric sequences, and so we again call it ``quasigeometric``. EXAMPLES:: @@ -289,7 +289,7 @@ def is_quasigeometric(self) -> bool: i.e. the union of multiple geometric sequences, or geometric after term ``u0``. If `\\alpha/\\beta` is a `k` th root of unity, where `k>1`, then necessarily `k = 2, 3, 4, 6`. - Then `F = [[0,1],[c,b]` is diagonalizable, and `F^k = [[\\alpha^k, 0], [0,\\beta^k]]` is scaler + Then `F = [[0,1],[c,b]` is diagonalizable, and `F^k = [[\\alpha^k, 0], [0,\\beta^k]]` is a scalar matrix. Thus for all values of `j` mod `k`, the `j` mod `k` terms of `u_n` form a geometric series. From 94ed5e30cf7ffa1c18a5e98eaa06f2e02793e1f4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Chapoton?= Date: Tue, 18 Jun 2024 08:55:17 +0200 Subject: [PATCH 3/3] Apply suggestions from code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Matthias Köppe --- src/sage/combinat/binary_recurrence_sequences.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sage/combinat/binary_recurrence_sequences.py b/src/sage/combinat/binary_recurrence_sequences.py index f473ce162f2..76765b38aab 100644 --- a/src/sage/combinat/binary_recurrence_sequences.py +++ b/src/sage/combinat/binary_recurrence_sequences.py @@ -289,7 +289,7 @@ def is_quasigeometric(self) -> bool: i.e. the union of multiple geometric sequences, or geometric after term ``u0``. If `\\alpha/\\beta` is a `k` th root of unity, where `k>1`, then necessarily `k = 2, 3, 4, 6`. - Then `F = [[0,1],[c,b]` is diagonalizable, and `F^k = [[\\alpha^k, 0], [0,\\beta^k]]` is a scalar + Then `F = [[0,1],[c,b]` is diagonalizable, and `F^k = [[\\alpha^k, 0], [0,\\beta^k]]` is a diagonal matrix. Thus for all values of `j` mod `k`, the `j` mod `k` terms of `u_n` form a geometric series. @@ -900,7 +900,7 @@ def _next_good_prime(p, R, qq, patience, qqold): # requiring that b^2 + 4c is a square in GF(R._ell) ensures that the period mod R._ell # divides R._ell - 1 - if legendre_symbol(R.b**2 + 4 * R.c, R._ell) == 1: + if legendre_symbol(R.b**2 + 4*R.c, R._ell) == 1: N = _goodness(R._ell, R, p)