(* This package, which contains over 200 functions and options relevant to computational number theory, accompanies the book A Course in Computational Number Theory, by David Bressoud and Stan Wagon of Macalester College. Copyright Springer-Verlag, New York, 1999. Not to be reproduced without permission. For more information contact: wagon@macalester.edu *) (* reference for pi(x): M. Deleglise and J. Rivat; Computing pi(x): the Meissel, Lehmer, Lagarias, Miller, Odlyzko method Math. Comp. 65 (1996) 235-45. *) If[Names["Explorer`*`*"] != {}, Remove["Explorer`*`*"]]; If[Names[#] =!= {}, Remove[#]] & /@ { "Global`ContinuantQ", "CNT`ContinuantQ", "Global`Euclid", "Global`Repunit","Global`CFFactor", "Global`SqrtModPrime"}; Off[General::spell, General::spell1]; Unprotect[PrimePi, Pi, UpArrow, Solve, Normal, PowerMod]; SetAttributes[Pi, Listable]; PrimePi[10^22] = 201467286689315906290; PrimePi[10^21] = 21127269486018731928; PrimePi[10^20] = 2220819602560918840; PrimePi[10^19] = 234057667276344607; PrimePi[10^18] = 24739954287740860; PrimePi[10^17] = 2623557157654233; PrimePi[10^16] = 279238341033925; PrimePi[10^15] = 29844570422669; PrimePi[10^13] = 346065536839; PrimePi[10^14] = 3204941750802; PrimePi[10^12] = 37607912018; PrimePi[10^11] = 4118054813; PrimePi[10^10] = 455052511; PrimePi[10.^21] = 21127269486018731928; PrimePi[10.^20] = 2220819602560918840; PrimePi[10.^19] = 234057667276344607; PrimePi[10.^18] = 24739954287740860; PrimePi[10.^17] = 2623557157654233; PrimePi[10.^16] = 279238341033925; PrimePi[10.^15] = 29844570422669; PrimePi[10.^14] = 3204941750802; PrimePi[10.^13] = 346065536839; PrimePi[10.^12] = 37607912018; PrimePi[10.^11] = 4118054813; PrimePi[10.^10] = 455052511; UpArrow::usage = "a\[UpArrow]b returns a^a^a... b times."; (a_) \[UpArrow] (n_) := CNT`Tower[a, n]; UpArrow[a_,n_] := CNT`Tower[a,n]; Global`SuperStar[Log][a_, n_] := CNT`LogStar[a, n]; DeclarePackage["StrongPseudoprimeData`", {"StrongPseudoprimes", "LucasPseudoprimesAllTypes"}]; DeclarePackage["NumberTheory`FactorIntegerECM`", "FactorIntegerECM"]; BeginPackage["CNT`",{ "NumberTheory`NumberTheoryFunctions`", "NumberTheory`ContinuedFractions`", "Utilities`FilterOptions`", "Algebra`InequalitySolve`", "Algebra`PolynomialPowerMod`"}]; Unprotect["Explorer`Primes`NextPrime", "NumberTheory`NumberTheoryFunctions`NextPrime", "NumberTheory`NumberTheoryFunctions`PreviousPrime" ]; v4Q = ($VersionNumber >= 4); Unprotect[CarmichaelLambda, ContinuedFraction, ClassNumber,ContinuedFractionForm, ChineseRemainder, ChineseRemainderTheorem, NextPrime, PrimitiveRoot, QuadraticRepresentation]; Clear[NumberTheory`NumberTheoryFunctions`ClassNumber]; AbundantQ::usage = "AbundantQ[n] gives True if n is an abundant number."; AdditionChain::usage = "AdditionChain[n, method] gives the addition chain for n according to any of the three methods: BinaryLeftToRight, BinaryRightToLeft, and FactorMethod."; AdmissibleQ::usage = "AdmissibleQ[c] gives True if the list of integers in c is an admissible constellation, meaning: there is no divisibility blockage to there being infinitely many occurrences in the primes."; AKSBPPolynomial::usage = "AKSBPPolynomial[n] gives, for n \[GreaterEqual] 7, a polynomial that is 0 if n is prime and seems to be nonzero otherwise. This is the AKSBPP (Agrawal, Kayal, Saxena, Bhattacharjee, Pandey) Conjecture for a polynomial time test for primality."; (*AliquotSequence::usage = "AliquotSequence[n] returns the sequence {n, s(n), s(s(n)), ..., 1, 0}, where s is the sum of the divisors of n, excluding n itself.";*) AmicableQ::usage = "AmicableQ[n] gives True if n is a member of a pair of amicable integers."; AssumedPrimeBound::usage = "AssumedPrimeBound is an option to PrattCertificate that causes primes under this bound to be deemed prime without need of certification."; AverageList::usage = "AverageList[f, x] gives the list of cumulative averages of f up to f[x]."; Average::usage = "Average[f, x] gives the average value of f[1],f[2],...,f[x]."; SummatoryList::usage = "SummatoryList[f, denom, {x, a, b, m:1}, start:0] forms the cumulative sum of f as x goes from a to b in steps of 1, keeping every mth entry only, and then dividing the entries by denom, with the appropriate x-value. If one wants to start at a specified value other than 0, that is done by the start argument."; Base::usage = "Base is an option to Repunit that sets the base."; BinaryLeftToRight::usage = "BinaryLeftToRight is a possible second argument to AdditionChain and a possible setting to the Method option of PowerAlgorithm."; BirthDate::usage = "BirthDate[n] returns the birthdate such that \!\((md)\^2 + y\^2 = n\) where m, d, and y are month, day, and year, respectively."; BinaryRightToLeft::usage = "BinaryRightToLeft is a possible second argument to AdditionChain and a possible setting to the Method option of PowerAlgorithm."; BLSPrimeProof::usage = "BLSPrimeProof[n] returns a proof of primality of n, which must be prime, based on the Brillhart-Lehmer-Selfridge certification procedure."; Buchstab\[Omega]::usage = "Buchstab\[Omega][u] gives Buchstab's \[Omega] function in symbolic form, defined for 1 \[LessEqual] u \[LessEqual] 5."; Buchstab\[Omega]To10::usage = "Buchstab\[Omega]To10[u] gives the Buchstab omega function defined numerically on the interval [1, 11], to 20 digits of precision. The first usage takes a little time."; CarmichaelLambda::usage = "CarmichaelLambda[n] gives the Carmichael function \[Lambda](n), defined as the smallest integer m such that \!\(k\^m\)\[Congruent]1 (mod n) for all k relatively prime to n. CarmichaelLambda[p, e] is a little faster for \!\(p\^e\), p prime."; CarmichaelQ::usage = "CarmichaelQ[n] gives True if n is a Carmichael number."; Certified::usage = "Certified is an option to RandomPrime that causes the result to be a certified prime. This is generally faster than asking for a mere probable prime. The default of Automatic returns a certified prime if no congruence conditions are specified."; CertifiedPrime::usage = "CertifiedPrime[d] generates a prime number having d digits that is certifiably prime. CertifiedPrime[d, {a, m}] generates a prime that is congruent to a (mod m), though if m is too large, the prime might have more than d digits."; CertificateTreeForm::usage = "CertificateTreeForm[data] takes a Pratt certificate in the list-form returned by PrattCertificate and turns it into a gridbox in the form of a tree."; CFFactor::usage = "CFFactor[n] attempts to factor n using the continued fraction method. It works best on integers having between 15 and 30 digits."; ChineseRemainder::usage = "ChineseRemainder[list1, list2] gives the minimal nonnegative integer solution of Mod[r, list2] == list1. The solution is unique modulo the LCM of list2."; ClassNumber::usage = "ClassNumber[d] gives the class number of d whenever d is negative and congruent to 0 or 1 (mod 4)."; (*FundamentalDiscriminantQ::usage = "xxx";*) Columns::usage = "Columns is an option to FunctionTable that sets the number of columns."; CompositeColor::usage = "CompositeColor is an option to GaussianPrimePlot that sets the color of the composites."; Conductor::usage = "Conductor[{a, b, c,...}] gives the least r such that r and all greater integers are representable as a x + b y + c z +... with nonnegative x, y, z,... . "; ConstellationEstimate::usage = "ConstellationEstimate[A, prec:6] returns the function, of R, say, that gives the estimate for the number of times the pattern A occurs in the primes up to R (or the Gaussian primes in the first octant of radius R). R can be either numeric or symbolic, and in the latter case the result is the symbolic function of R."; ContinuantQ::usage = "ContinuantQ[q_List, d_:1] gives the starting term of the Euclidean algorithm remainder sequence ending in d that has q has the list of quotients."; ContinuedFraction::usage = "ContinuedFraction[x] gives the continued fraction representation for the rational number x. ContinuedFraction[x, n] gives the continued fraction representation for the real number x to order n. Note that the order returned may be less than n if the continued fraction terminates before n steps. If x is a quadratic irrational then ContinuedFraction[x] returns a periodic representation of the continued fraction for x."; CRTGrid::usage = "CRTGrid[m, n] returns a graphic illustrating the mapping from pairs of integers mod m and mod n, respectively, to integers modulo the product m n."; DeficientQ::usage = "DeficientQ[n] gives True if n is a deficient number."; DiskBackground::usage = "DiskBackground is an option to ReachableGaussianPrimes that sets the color of the background disk."; DiskColor::usage = "DiskColor is an option to GaussianPrimePlot that sets the color of the background disk."; DivisorSigma1::usage = "DivisorSigma1[k, n] gives the sum of kth powers of the divisors of n that are congruent to 1 (modulo 4)."; DivisorSigma3::usage = "DivisorSigma3[k, n] gives the sum of kth powers of the divisors of n that are congruent to 3 (modulo 4)."; EratosthenesTable::usage = "EratosthenesTable[n] gives a table showing the sieve of Eratosthenes working on the integers from 2 to n."; Eratosthenes::usage = "Eratosthenes[n] returns a list of all primes up to n using the sieve. It is faster to use Prime[Range[PrimePi[n]]]."; Euclid::usage = "Euclid[n] is the product of the first n primes, plus 1."; EuclideanAlgorithmLength::usage = "EuclideanAlgorithmLength[a, b] gives the number of steps the Euclidean algorithm takes on the pair {a, b}."; EuclidSequence::usage = "EuclidSequence[n] gives the first n terms of the sequence defined by Euclid's proof that there are infinitely many prime numbers, where we start with 2 and use the smallest prime divisor of \!\(p\_1\)\!\(p\_2\)\[CenterDot]\[CenterDot]\[CenterDot]\!\(p\_r\)\[ThinSpace]+\[ThinSpace]1 at each step."; EulerPhiLines::usage = "EulerPhiLines is an option to VisualPowerGrid that shows lines marking off the columns in steps of \[Phi](n)."; EulerPhiUpperBound::usage = "EulerPhiUpperBound[n] gives a bound M such that if \[Phi][x] = n then x \[LessEqual] M."; EulerPseudoprimeQ::usage = "EulerPseudoprimeQ[a, n] returns True if the odd integer n is an Euler pseudoprime to base a (which entails that n be composite)."; EulerPseudoprimeTest::usage = "EulerPseudoprimeTest[a, n] returns True if the odd integer n passes the Euler pseudoprime test using base a."; EvaluateIntegral::usage = "EvaluateIntegral is an option to ConstellationEstimate that causes the integral in the formula to be evaluated in terms of the logarithmic integral, li."; ExtendedGCDValues::usage = "ExtendedGCDValues is an option to FullGCD that causes the sequences that generate the extended gcd coefficients to be included in the output."; ExponentForm::usage = ToString[StringForm[ "ExponentForm[n, d:10] writes n in the form `` ± a.",Superscript["d","e"]], StandardForm]; FactorBaseSize::usage = "FactorBaseSize is an option to CFFactor that sets the size of the factor base."; FactorForm::usage = "FactorForm[l] expresses a list l of the form returned by FactorInteger in typeset form. If l is an integer or Gaussian integer, it factors l first."; FactorMethod::usage = "FactorMethod is a possible second argument to AdditionChain and a possible setting to the Method option of PowerAlgorithm."; Fermat::usage = "Fermat is a possible setting for the Type option to PseudoprimeSearch, and causes the search to look for b-pseudoprimes where b is set by the FermatBase option."; FermatBase::usage = "FermatBase is an option to PseudoprimeSearch that specifies the base used when the type is either Fermat or Strong."; FermatNumber::usage = "FermatNumber[n] returns the Fermat number \!\(2\^\(2\^n\)+1 \)."; Fibonacci::usage = "Fibonacci[n] gives the nth Fibonacci number. Fibonacci[n, x] gives the nth Fibonacci polynomial, using x as the variable. Fibonacci is also a possible setting to the Type option to PseudoprimeSearch."; FirstSPSPWitness::usage = "FirstSPSPWitness[n, d:2] gives the first b \[GreaterEqual] d such that the composite number n is not a b-strong pseudoprime."; FormatOutputNicely::usage = "FormatOutputNicely is an option to ConstellationEstimate that applies TraditionalForm then StandardForm to the output, which improves the formatting."; FreePart::usage = "FreePart[n, p] returns the pair {m, r} where \!\(n = m\ p\^r\) and m is not divisible by p."; FrobeniusSolve::usage = "FrobeniusSolve[{a1,a2,a3,..}, b] finds one solution to x a1 + y a2 + c a3 + ... = b in nonnegative integers. It returns {} if there are no solutions. The coefficients ai must be relatively prime."; FrobeniusNonRepresentables::usage = "FrobeniusNonRepresentables[{a1, a2, a3,...}] finds the set of all b for which the equation x a1 + y a2 + c a3 + ... = b has no solution in nonnegative integers. The coefficients ai must be relatively prime."; FromCFSquared::usage = "FromCFSquared[{a,b,c,...{d,e,f...}}] gives the square of the quadratic irrational whose CF has preperiod a,b,c,... and periodic part {d,e,f,...}."; FullContinuant::usage = "FullContinuant[q, d_:1] gives the Euclidean algorithm remainder sequence ending in d that yields q as the quotients."; FullContinuedFraction::usage = "FullContinuedFraction[c] gives a table showing the development of the continued fraction of a rational or quadratic irrational."; FullGCD::usage = "FullGCD[m, n] gives the full table of Euclidean algorithm data starting from m and n. If one of m, n is complex, then the algorithm for the Gaussian integers is used."; FullPrattCertificate::usage = "FullPrattCertificate is an option to PrattCertificate that causes the classical Pratt certificate to be returned. When this is False, a more concise form of the certificate is returned, where only primes from a \"half-factorization\" of n - 1 are used."; FunctionHeading::usage = "FunctionHeading is an option to FunctionTable that specifies the heading of the column of function values."; FunctionTable::usage = "FunctionTable[f, {n0, n1}] gives a grid box of the values of f as n goes from n0 to n1. FunctionTable[f, n] uses values from 1 to n."; GaussianIntegerQ::usage = "GaussianInteger[e] gives True if e is a Gaussian integer."; GaussianNorm::usage = "GaussianNorm[z] returns the norm of the Gaussian integer z, which is \!\(a\^2 + b\^2\) if z = a + b I."; GaussianMod::usage = "GaussianMod[a, z] gives a special residue class of a modulo a, where a and z are Gaussian integers. The important distinction with Mod is that if z is a Gaussian prime that is not a rational prime, then GaussianMod[a, z] is an integer between 0 and GaussianNorm[z] - 1, inclusive."; GaussianPrimePlot::usage = "GaussianPrimePlot[n] creates an image of the Gaussian primes having absolute value at most \!\(n\^2\)."; GaussianPrimeQ::usage ="GaussianPrimeQ[u] returns True if the Gaussian integer u is a Gaussian prime; False otherwise."; GridBoxFix::usage = "GridBoxFix[m] is a utility that typesets the items in the matrix m so that a GridBox can be used on them."; HardyLittlewood::usage = "HardyLittlewood[A, prec:6] returns the Hardy-Littlewood constant for the set A of integers or Gaussian integers."; HarmonicRearrangement::usage = "HarmonicRearrangement[target] rearranges terms of the alternating harmonic series so that partial sums converge to target. The MaxTerms option controls the number of terms of the rearranged series. The output contains any of three items depending on the options ShowSeries, ShowPartialSums, and ShowSignChanges. The default is to show all three."; HeadingTextColor::usage = "HeadingTextColor is an option to PythagoreanTriples that controls the color of the text headings."; HighlightColor::usage = "HighlightColor is an option to IteratedPrimeGaps that sets the color of the first column."; HighlightFirstLevel::usage = "HighlightFirstLevel is an option to PrattCertificate that causes the first level to be colored red when the tree form of the certificate is called for."; HighlightFontColor::usage = "HighlightFontColor is an option to FullContinuedFraction that sets the color of the highlighted entries."; IntegerToString::usage = "IntegerToString[n] turns an integer that represents a string according to the coding of StringToInteger into the string."; IteratedPrimeGaps::usage = "IteratedPrimeGaps[n, i] returns a grid box showing i iterations of the absolute difference operation on the first n primes."; KnuthTrabbPardo::usage = "KnuthTrabbPardo[1, a] (resp., KnuthTrabbPardo[2, a]) gives the probability that the largest (resp., second-largest) prime factor of an integer n is less than \!\(n\^\(1/a\)\)."; (*KroneckerSymbol::usage = "KroneckerSymbol[d, n] gives the Kronecker symbol.";*) Labels::usage = "Labels is an option to FullGCD that controls whether the numbers are accompanied by descriptive labels."; LambdaExp::usage = "LambdaExp[m] is the least n such that LambdaLog[n] is m."; LambdaLog::usage = "LambdaLog[m] gives the number of times CarmichaelLambda must be applied to reduce n to 1."; LargestExponent::usage = "LargestExponent[n] gives the highest exponent of any of the primes appearing in the factorization of n."; LastFactorialDigit::usage = "LastFactorialDigit[n] gives the rightmost nonzero digit, in base 10, of n!."; LegendrePhi::usage = "LegendrePhi[n,a] gives the numbers of integers in [1, n] that are not divisible by any of the first a primes."; Sieved\[CapitalPhi]::usage = "Sieved\[CapitalPhi][x, z] gives the number of integers in [1, x] all of whose prime divisors are greater than or equal to z."; LegendrePiFull::usage = "LegendrePiFull[n] produces a table showing the recursion that underlies Legendre's algorithm for \[Pi](n). n must be under \!\(10\^8\)."; LegendreSymbol::usage = "LegendreSymbol[a, n] gives the Legendre symbol when n is an odd prime."; li::usage = "li[x] gives LogIntegral[x]."; LinearCongruenceSolve::usage = "LinearCongruenceSolve[{a, b}, m] gives all solutions to the congruence a x \[Congruent] b (mod m)."; LinearDiophantineSolve::usage = "LinearDiophantineSolve[\!\(\"{\\!\\(a\\_1\\),\\!\\(a\\_2\\),. . .,\\!\\(a\\_m\\)}\"\), c] returns a matrix whose columns are a particular solution together with an integral basis for the set of homogeneous solutions to \!\(\"\\!\\(a\\_1\\)\\!\\(x\\_1\\)+...+\\!\\(a\\_n\\)\\!\\(x\\_n\\)=c\"\). A dot product with {1, \!\(n\_1\), \!\(n\_2\),...,\!\(TraditionalForm\`n\_\(m - 1\)\)} will get the parametric form of the solution set. The SolutionType option can be used to see a set of explicit solutions."; LinearRecurrence::usage = "LinearRecurrence[{c0,c1,c2,...,cj}, {a0,a1,a2,...,aj},k,m:Infinity] returns the kth term (modulo m) of the sequence a(n) defined by the recurrence a(n) = c0 a(n-j) + c1 a(n-j+1) + ... with initial conditions a(0)=a0, a(1)=a1,... . If k is a list then the corresponding list of terms is given."; LogStar::usage = "LogStar[b, n] gives the greatest integer m such that Tower[b, m] \[LessEqual] n. "; Lucas::usage = "Lucas[{P, Q}, m][n] gives the nth Lucas pair {U, V}, modulo m, for the parameter values P and Q. Omitting m sets it to Infinity, which means no modular reduction is performed. Lucas is also a possible setting to the type option in PseudoprimeSearch."; LucasCertificate::usage = "LucasCertificate[p] gives a certificate of p's primality using a Lucas test in a recursive fashion. It requires a factorization of p+1."; LucasParameter::usage = "LucasParameter[n] finds a P value that can be paired with Q=1 to form a strong Lucas test."; LucasParameters::usage = "LucasParameters is an option to PseudoprimeSearch that specifies the parameters for the Lucas sequence. It can be a pair or one of MethodA, MethodAStar, MethodAStartAtMinusSeven, or Automatic. Automatic is the same as MethodAStar."; LucasRank::usage = "LucasRank[{P, Q}, n] gives the rank of the odd integer n in the Lucas sequence determined by P and Q."; LucasU::usage = "LucasU[{P, Q}, m][n] gives the nth term of the Lucas U-sequence modulo m, for the parameter values P and Q. Omitting m sets it to Infinity, which means no modular reduction is performed."; LucasV::usage = "LucasV[{P, Q}, m][n] gives the nth term of the Lucas V-sequence modulo m, for the parameter values P and Q. Omitting m sets it to Infinity, which means no modular reduction is performed. LucasV is also a possible setting to the type option in PseudoprimeSearch."; Lucas\[CapitalLambda]::usage = "Lucas\[CapitalLambda][{P, Q}, n] gives, for n odd, the least universal rank mod n for the Lucas rank determined by P and Q."; Lucas\[CapitalPhi]::usage = "Lucas\[CapitalPhi][{P, Q}, n] is the analog of the Euler \[Phi] function for the Lucas sequence determined by P and Q."; LucasPseudoprimeQ::usage = "LucasPseudoprimeQ[{P, Q}, n] returns True if n is a Lucas pseudoprime for the U-sequence determined by P and Q."; LucasPseudoprimeTest::usage = "LucasPseudoprimeTest[{P, Q}, n] returns True if n passes the Lucas pseudoprime test for the U-sequence determined by P and Q."; LucasVPseudoprimeQ::usage = "LucasVPseudoprimeQ[{P, Q}, n] returns True if n is a Lucas pseudoprime for the V-sequence determined by P and Q."; LucasVPseudoprimeTest::usage = "LucasVPseudoprimeTest[{P, Q}, n] returns True if n passes the Lucas pseudoprime test for the V-sequence determined by P and Q."; \[Lambda]::usage = "\[Lambda][n] gives the least universal exponent modulo n."; MatrixPowerMod::usage = "MatrixPowerMod[a, n, m] gives the nth power of the matrix a modulo m."; MaxTerms::usage = "MaxTerms is an option to HarmonicRearrangement that gives a number of terms that should be examined."; Mersenne::usage = "Mersenne[k] gives the Mersenne number \!\(2\^k - 1\)."; MersennePrimeQ::usage = "MersennePrimeQ[n] gives True if n is a Mersenne prime. The determination is by the Lucas-Lehmer test."; MethodA::usage = "MethodA[n, start:5] gives the Lucas sequence parameters {1, Q} that can be used in a Lucas pseudoprime test on n. The method finds d (= 1 - 4 Q) first, trying 5, -7, 9, -11 ,..until a nonresidue mod n is found. A second argument, one of -7, 9, -11, ..., can be given so that the search is started later. MethodA is also a possible setting for the LucasParameters option to PseudoprimeSearch."; MethodAStar::usage = "MethodAStar[n] gives the Lucas sequence parameters {1, Q} that can be used in a Lucas pseudoprime test on n. The method finds d (= 1 - 4 Q) first, trying 5, -7, 9, -11 ,..until a nonresidue mod n is found. But if Q turns out to be 1, the pair (5, 5) is used instead. MethodAStar is also a possible setting for the LucasParameters option to PseudoprimeSearch."; MethodAPositive::usage = "MethodAPositive[n, start:5] gives the Lucas sequence parameters {1, Q} that can be used in a Lucas pseudoprime test on n. The method finds d (= 1 - 4 Q) first, trying 5, -7, 9, -11 ,..until a value is found for which JacobiSymbol[d, n] = +1. A second argument, one of -7, 9, -11, ..., can be given so that the search is started later."; MethodAStarPositive::usage = "MethodAStarPositive[n] gives the Lucas sequence parameters {1, Q} that can be used in a Lucas pseudoprime test on n. The method finds d (= 1 - 4 Q) first, trying 5, -7, 9, -11 ,..until a value is found for which JacobiSymbol[d,n] = +1. But if Q turns out to be 1, the pair (5,5) is used instead. A second parameter, one of -7, 9, -11, ..., can be given so that the search is started later."; MethodAStartAtMinusSeven::usage = "MethodAStartAtMinusSeven is a possible setting to the LucasParameters option of PseudoprimeSearch."; MillerRabinPrimeQ::usage = "MillerRabinPrimeQ[n, d] runs d strong pseudoprime tests on n using random choices for the base. If d is omitted it is assumed to be 1. A GCD against the first 100 primes is run first, for speed."; RootsMod::usage = "RootsMod[f, x, n] gives the roots of f, a polynomial in x, modulo n, using the Cantor-Zassenhaus algorithm."; NextPrime::usage = "NextPrime[n] gives the smallest prime p such that p > n. NextPrime[n, {i, m}] adds the condition that p \[Congruent] i (mod m)."; Nonresidue::usage = "Nonresidue[p] finds c such that c is not a quadratic residue for p, which must be an odd prime."; NumberSought::usage = "NumberSought is an option to PseudoprimeSearch that causes the search to halt when the specified number of examples has been found."; OddPart::usage = ToString[StringForm["OddPart[n] gives the pair {m, r} where m is odd and n = m ``.", Superscript[2,"r"]], StandardForm]; OrderMod::usage = "OrderMod[a, m] gives the least positive integer for which \!\(a\^k\) \[Congruent] 1 (mod m). Such exists iff a and m are relatively prime. OrderMod[a, p, e] assumes m = \!\(p\^e\), and is a little faster."; PairwiseCoprimeQ::usage = "PairwiseCoprimeQ[l] returns true if l is a list of pairwise coprime integers."; PartialFactor::usage = "."; TimeLimit::usage = " "; PascalTriangle::usage = "PascalTriangle[n] creates a grid box displaying the first n rows of Pascal's triangle."; PellSolutions::usage="PellSolutions[d, n] gives the nth positive integer solution to \!\(x\^2 - d\ y\^2 = \(±1\)\). If the second argument is {m, n} then the solutions from the mth to the nth are given."; PellSolve::usage = "PellSolve[d] gives the 2\[Times]2 matrix whose powers yield (in the first row) all the positive integer solutions to \!\(x\^2 - d\ y\^2 = \(±1\)\)."; PerrinPseudoprimeQ::usage = "PerrinPseudoprimeQ[n] gives True if n is composite but a[n]\[Congruent]0 (mod n) where a[n] is the Perrin sequence 3,0,2,3,2,5,5,7,10,12,17,22,... ."; PhiExp::usage = "PhiExp[n] gives the first integer m such that \[Phi] must be iterated n times to reduce m to 1."; PhiInverse::usage = "PhiInverse[n] gives the set of x for which \[Phi][x] = n."; PhiLog::usage = "PhiLog[n] gives the number of times EulerPhi must be applied to reduce n to 1."; PhiMultiplicity::usage = "PhiMultiplicity[n] gives the number of values x for which \[Phi][x] = n."; PocklingtonPrimeWitness::usage = "PocklingtonPrimeWitness[n, help] gives a Pocklington witness to the primality of n, where it is assumed that help contains primes that are sufficient to factor n-1 up to its square root. If help is omitted, then n-1 is factored and all the primes are used."; PowerAlgorithm::usage = "PowerAlgorithm[a, n, m, opts] shows the workings of various methods of fast exponentiation to form \!\(a\^n\) modulo m. The method is controlled by the Method option. If m is omitted, no modular reduction takes place. Possible values of the Method option are BinaryLeftToRight, BinaryRightToLeft, and FactorMethod."; PowerGrid::usage = "PowerGrid[n] gives a gridbox showing the mod-n values of all powers \!\(a\^i\)."; PrattCertificate::usage = "PrattCertificate[p] returns a Pratt certificate of p's primality, where p is assumed to have returned True under PrimeQ."; (*PreviousPrime::usage = "PreviousPrime[n] gives the largest prime p such that p < n.";*) PrimeColor::usage = "PrimeColor is an option to GaussianPrimePlot that sets the color of the primes."; PrimeFactorial::usage = "PrimeFactorial[n] gives the product of the first n primes."; PrimeNumberRace::usage = "PrimeNumberRace[a, b, s:1] presents a plot of the mod-4 prime race over the interval [a, b]. The number of plotted points is determined by the step size s."; PrimePiMod4::usage = "PrimePiMod4[i][x] gives the number of primes less than or equal to x that are congruent to i (mod 4)."; (*PrimePowerQ::usage = "PrimePowerQ[n] gives True if n is a power of a prime.";*) PrimeQTest::usage = "PrimeQTest is an option to PrattCertificate that asks the PrimeQ test on the input to be performed. Setting this to False saves time in situations where it is known that the input passes this test."; PrimeZeta::usage = "PrimeZeta[s, prec:6] gives the prime zeta function, the sum of \!\(1/p\^s\) over prime p. The GaussianIntegers option causes this to compute the Gaussian prime zeta function."; PrimitiveRoot::usage = "PrimitiveRoot[n] gives a cyclic generator of the multiplicative group (mod n), when n is a prime power or twice a prime power."; PrimitiveRootColor::usage = "PrimitiveRootColor is an option to VisualPowerGrid that sets the color of the values that are primitive roots."; PrimitiveRoots::usage = "PrimitiveRoots[n] gives the set of primitive roots for n."; PrintImmediately::usage = "PrintImmediately is an option to PseudoprimeSearch that causes any found pseudoprimes to be printed immediately."; PrintTheCounts::usage = "PrintTheCounts is an option to ReachableGaussianPrimes that causes a count of the primes to be printed, but only if the UseLines option is set to False."; PseudoprimeQ::usage = "PseudoprimeQ[b, n] gives True if n is a b-pseudoprime (which entails that n be composite)."; PseudoprimeSearch::usage = "PseudoprimeSearch[n0, n1] gives the pseudoprimes betweem n0 (which defaults to 1 if omitted) and n1. The default type is Fermat pseudoprimes to base 2, but the Type option can be used to specify a variety of other types."; PseudoprimeTable::usage = "PseudoprimeTable[k, n] generates a table of the b-pseudoprimes that are less than equal to n, for the first k many prime values of b."; PythagoreanTriples::usage = "PythagoreanTriples[r, s] returns a table of s primitive Pythagorean triples starting with the rth. The ordering is inherited from the generators. Assumptions may be given as an option that is a logical condition on a, b, and c."; \[Phi]::usage = "\[Phi][n] gives EulerPhi[n]."; \[Pi]::usage = "Pi (\[Pi]) is the constant pi, with numerical value approximately equal to 3.14159. \[Pi][x] gives the number of primes less than or equal to x (which is also available as PrimePi[x])."; Quadratic::usage = "Quadratic is a possible setting for the Type option to PseudoprimeSearch, and causes the search to look for quadratic pseudoprimes where the Lucas parameters are set by the LucasParameters option."; QuadraticCongruenceSolve::usage = "QuadraticCongruenceSolve[{a, b, c}, m] gives all solutions to the congruence \!\(a\ x\^2 + b\ x + c \[Congruent] 0\ \((mod\ m)\)\)."; QuadraticConjugate::usage = "QuadraticConjugate[d] returns the quadratic conjugate of a quadratic irrational."; QuadraticFibonacci::usage="QuadraticFibonacci is a possible setting to the Type option to PseudoprimeSearch."; QuadraticIrrationalQ::usage = "QuadraticIrrationalQ[d] gives True if d is a real quadratic irrational."; QuadraticMod::usage = "QuadraticMod[c, n] returns the residue of a quadratic irrational c modulo the integer n. It avoids doing modular arithmetic on the radical in c."; QuadraticNormalForm::usage = "QuadraticNormalForm[c] gives {P, d, Q} that define c as \!\(\((P + \@d)\)/Q\) where Q divides \!\(d - P\^2\)."; QuadraticPowerMod::usage = "QuadraticPowerMod[a,n,p] returns the residue of the nth power of a quadratic irrational c modulo the integer n. It avoids doing modular arithmetic on the radical in c."; QuadraticPseudoprimeTest::usage = "QuadraticPseudoprimeTest[{P, Q}, n] returns True if n passes the Lucas pseudoprime test for both the U- and V-sequences determined by P and Q."; QuadraticPseudoprimeQ::usage = "QuadraticPseudoprimeQ[{P, Q}, n] returns True if n is composite but passes the Lucas pseudoprime test for both the U- and V-sequences determined by P and Q."; QuadraticRepresentation::usage= "QuadraticRepresentation[n, f, g] returns the set of all pairs of nonnegative integers {a, b} such that n = f a^2 + g b^2, for positive f and g."; QuadraticResidueQ::usage = "QuadraticResidueQ[b, p] returns True if p is prime and b is a quadratic residue mod p."; QuarticResidueSymbol::usage = "QuarticResidueSymbol[a, p] returns the value of the quartic residue symbol computed in the Gaussian integers. p must be a Gaussian prime different from 1+I and a can be any Gaussian integer."; RandomBases::usage = "RandomBases is an option to MillerRabinPrimeQ that, when True causes random bases to be used. When False, the bases are taken to be the primes starting from 2."; RandomPrime::usage = "RandomPrime[d] gives a d-digit random prime. RandomPrime[n, {i, m}] adds the condition that p \[Congruent] i (mod m)."; ReachableGaussianPrimes::usage = "ReachableGaussianPrimes[d, rad] displays the Gaussian primes within distance r from the origin that one can reach from 1 + I using steps of size at most d and stepping on the Gaussian primes only. The unreachable primes are shown in a contrasting style."; ReachablesOnly::usage="ReachablesOnly is an option to ReachableGaussianPrimes that causes the unreachable primes to not be shown. If UseLines is set to False, then ReachableOnly automatically becomes True."; ReachableStyle::usage = "ReachableStyle is an option to ReachableGaussianPrimes that sets the style of the reachable points."; ReducedResidues::usage = "ReducedResidues[n] gives the set of integers in {1,2,...,n} that are relatively prime to n."; ReducedResiduesOnly::usage = "ReducedResiduesOnly is an option to VisualPowerGrid that causes powers of reduced residues only to be shown."; Repunit::usage = "Repunit[n, d] is the n-digit integer consisting of only the digit d. The default base is 10."; RightSide::usage = "RightSide is an option to PellSolutions that specifies whether the right side of the equation is 1, -1, or {-1, +1}."; RSAEncode::usage = "RSAEncode[s, e, n] turns the message-string s into a number and then encodes that number by the RSA method using exponent e and modulus n. If s is too large, then a list of integers is returned."; RSADecode::usage = "RSADecode[m, k, n] decodes the integer (or list) m using the decoding inverse k and the modulus n."; SFTF::usage = "SFTF[expr] converts the expression to traditional form, and then converts that to standard form. This can often lead to traditional looking output, butin StandardForm."; Show2s::usage = "Show2s is an option to PrattCertificate that, when True, suppresses the 2s that occur in all the recursive factorizations."; ShowSmithValues::usage = "ShowSmithValues is an option to FullGCD that, for inputs p x where p is a prime congruent to 1 (mod 4) and x is a mod-p square root of -1, highlights the first two remainders under \!\(\@p\)."; ShowData::usage = "ShowData is an option to CFFactor that causes a summary of the algorithm's performance to be printed."; ShowNumbers::usage = "ShowNumbers is an option to VisualPowerGrid that causes numbers to be shown in the array in place of colored squares."; ShowPartialSums::usage = "ShowPartialSums is an option to HarmonicRearrangement that, when True, causes the partial sums of the rearranged series at the sign changes to be shown."; ShowCertificate::usage = "ShowCertificate is an option to CertifiedPrime that, when true, causes a primality certificate to be returned. Its first entry is the certified prime."; ShowPrimitiveRoots::usage = "ShowPrimitiveRoots is an option to VisualPowerGrid that highlights the primitive roots in red."; ShowSeries::usage = "ShowSeries is an option to HarmonicRearrangement that, when True, causes the rearranged series to be shown."; ShowSignChanges::usage = "ShowSignChanges is an option to HarmonicRearrangement that, when True, causes the positions of the sign changes in the rearranged series to be shown."; ShowSquares::usage = "ShowSquares is an option to PythagoreanTriples that causes the squares of the triple to be shown provided the grid box output is requested."; ShowWitnesses::usage = "ShowWitnesses is an option to PrattCertify that causes a list of the witnesses used to be returned along with the certificate."; SolutionType::usage = "SolutionType is an option to LinearDiophantineSolve that can take on the values All, Positive, or NonNegative. In the latter two cases, if there are finitely many solutions then the finite solution set is returned."; SpspSequence::usage = "SpspSequence[b, n] gives the b-sequence of powers used to decide whether n is a b-strong pseudoprime."; SquareColors::usage = "SquareColors is an option to VisualPowerGrid that specifies a list of n colors to use for the n possible values that occur in the table. The Automatic setting sets +1 to red, -1 to yellow, 0 to white, and other values to varying shades of blue."; SqrtModAll::usage = "SqrtModAll[a, n] returns the set of all mod-n square roots of a, in increasing order."; SqrtModPrime::usage = "SqrtModPrime[a, p] returns the smallest square root of a modulo p, where p is prime, or {} if none exists."; SqrtModPrimePower::usage = "SqrtModPrimePower[a, p, r] returns the set of all mod-\!\(p\^r\) square roots of a, in increasing order."; SRKCorrect::usage = "SRKCorrect[s] corrects the string s, assuming it has been protected by the Sethi, Rajaraman, and Kenjale base-37 method."; SRKProtect::usage = "SRKProtect[s] adds two check digits to the string s using the Sethi, Rajaraman, and Kenjale method. There must be 11 or fewer characters in s and they must be upper-case letters or numerals or a hyphen."; StartAtMinusSeven::usage = "StartAtMinusSeven is a possible setting for the Method option to PseudoprimeSearch, and sets the method to MethodA with a start at -7 instead of the usual 5."; StringToInteger::usage = "StringToInteger[s] turns the string s into an integer by replacing each character by a two-digit number."; Strong::usage = "Strong is a possible setting for the Type option to PseudoprimeSearch, and causes the search to look for strong b-pseudoprimes where b is set by the FermatBase option."; StrongQuadratic::usage = "StrongQuadratic is a possible setting for the Type option to PseudoprimeSearch, and causes the search to look for strong quadratic pseudoprimes (i.e., Lucas U and V pseudoprimes)."; StrongQuadraticPseudoprimeQ::usage = "StrongQuadraticPseudoprimeQ[{P, Q}, n] is True if n is composite and passes the strong quadratic pseudoprime test for parameters P, Q."; StrongQuadraticPseudoprimeTest::usage = "StrongQuadraticPseudoprimeTest[{P, Q}, n] is True if n passes the strong quadratic pseudoprime test for parameters {P, Q}."; StrongPseudoprimeQ::usage = "StrongPseudoprimeQ[b, n] gives True if n is a b-strong pseudoprime (which entails that n be composite)."; SumOfFourSquares::usage = "SumOfFourSquares[n] gives four integers so that n is the sum of their squares."; SumOfSquaresRNoCache::usage ="SumOfSquaresRNoCache[k, n] gives the number of representations of n as a sum of k squares, avoiding caching."; SumOfSquaresRPrimitive::usage = "SumOfSquaresRPrimitive[k, n] gives the number of primitive representations of n as a sum of k squares."; SurvivorCountOnly::usage = "SurvivorCountOnly is an option to EratosthenesTable that causes the number of survivors at each stage to be returned."; \[Sigma]::usage = ToString[StringForm["\[Sigma][n] gives the sum of divisors of n; ``[n] gives the sum of the kth powers. This is standard in Mathematica as DivisorSigma[k, n]. ``[n] gives the reduced function, where n is excluded. ``[n] gives the reduced `` value.", Subscript[\[Sigma], "k"], SuperMinus[\[Sigma]], Subsuperscript[\[Sigma], "k", "-"], Subscript[\[Sigma], "k"]], StandardForm]; TargetMargin::usage = "TargetMargin is an option to CFFactor that sets the excess over the factor base size in the number of squares sought. A setting of 10 will often be adequate, but the default is 30."; Tower::usage = "Tower[a,n] gives the exponential tower \!\(a\^\(a\^\(a\^\[AscendingEllipsis]\)\)\) with n terms."; TowerMod::usage = "TowerMod[a, n, m] gives a\[UpArrow]n modulo m, where a\[UpArrow]n is a^a^a... n times, bracketed top down."; TraceInterval::usage = "TraceInterval is an option to PseudoprimeSearch that causes a tracing message to be printed after every t search steps, where t is the value of the option."; TraceIntervalFactor::usage = "TraceIntervalFactor is an option to CFFactor that controls the tracing prints. It should be between 0 and 1."; Tree::usage = "Tree is an option to PrattCertificate that requests that the certificate be returned in the form of a tree."; Type::usage = "Type is an option to PseudoprimeSearch that specifies the type of pseudoprime that is being sought. Possible settings are Fermat, Strong, Lucas, LucasV, Quadratic, StrongQuadratic, Fibonacci, QuadraticFibonacci."; Unfactor::usage = "Unfactor[{{p, r}, {q, s}},...] computes \!\(\(p\^r\) \(q\^s\) \[Ellipsis]\)."; UnitsColor::usage = "UnitsColor is an option to GaussianPrimePlot that sets the color of the units."; UnreachableStyle::usage="UnreachableStyle is an option to ReachableGaussianPrimes that sets the style of the unreachable points."; UseArrows::usage = "UseArrows is an option to VisualPowerGrid that causes arrows to point to the primitive roots."; UseGridBox::usage = "UseGridBox is an option to FullGCD and UVSequence that causes the results to be shown in a GridBox."; UseLines::usage = "UseLines is an option to ReachableGaussianPrimes that shows the edges connecting the reachable primes."; UVSequence::usage = "UVSequence[{P, Q}, n] gives the Lucas sequence of U and V values for parameters P and Q and indices m, 2m, 4m, ..., \!\(2\^k\)m, where \!\(2\^k\)m is n ± 1."; VisualPowerGrid::usage = "VisualPowerGrid[{amax, imax}, n] shows a colored chart of all powers \!\(a\^i\) mod n where a runs from 1 to amax and i from 1 to imax. The default setting of SquareColors causes occurrences of +1 to be shown in red, -1 in yellow, 0 in white, and other nonzero values in shades of blue. VisualPowerGrid[n] takes the bounds on a and i to be n-1."; VisualSieve::usage = "VisualSieve[] generates a 5-frame movie showing the sieve of Eratosthenes up to 121, using the primes 2, 3, 5, 7, 11 to strike out composites."; ZetaGaussian::usage = "ZetaGaussian[s] gives the version of the Riemann zeta function appropriate for the Gaussian integers: the sum of 1/N(z)^s where z ranges over \[DoubleStruckCapitalZ][i], ignoring associates."; \[DoubleStruckCapitalZ]::usage = "\!\(\(\[DoubleStruckCapitalZ]\_n\%*\)\) gives the set of reduced residues modulo n."; If[!v4Q, FromContinuedFraction::usage = "FromContinuedFraction[list] reconstructs a number from the list of its continued fraction terms."; IntegerExponent::usage = "IntegerExponent[n,b] gives the highest k such that b^k divides n."; Convergents::usage = "Convergents[q] gives the list of continued fraction convergents for the rational number q. Convergents[x,n] gives the first n convergents for the irrational x."]; (*PeriodicCFForm::usage = "PeriodicCFForm[a, b] represents the continued fraction whose initial part is the list a and whose periodic part is b. Its output is formatted as << a | b >>."; *) ChineseRemainder::badlen = "The two lists do not have the same length."; Conductor::cantdo = "The algorithm cannot handle this case."; ContinuedFraction::badqua = "ContinuedFraction works symbolically (i.e., on single arguments) only if the input is rational or a positive quadratic irrational."; CRTGrid::badgcd = "The arguments to CRTGrid must be relatively prime."; FirstSPSPWitness::noncom = "`` is prime, and hence there are no witnesses to its compositeness."; IteratedPrimeGaps::bigitr = "The number of iterations cannot be more than the number of primes."; KnuthTrabbPardo::toobig = "This implementation handles only arguments up to 12."; LinearDiophantineSolve::infsol = "There are infinitely many solutions so a finite list cannot be given."; LinearDiophantineSolve::toomny = "The Assumptions option can be used only when the SolutionType option is set to Positive or NonNegative."; PrattCertificate::noprim = "`` is not prime."; PrattCertificate::notree = "The tree format will not work when the ShowWitnesses option is turned on"; PrattCertificate::nowitn = "The request for witnesses works only if the AssumedPrimeBound is set to its default, 2."; PrattCertificate::noprim = "`` is not prime."; PrattCertificate::notree = "The tree format will not work when the ShowWitnesses option is turned on."; PrattCertificate::nowitn = "The request for witnesses works only if the AssumedPrimeBound is set to its default, 2."; PseudoprimeSearch::toobig = "The search domain is too large. It must not involve integers larger than \!\(2\^31\)."; FrobeniusNonRepresentables::gcd = "This function requires that the coefficients be relatively prime."; FrobeniusSolve::gcd = "This function requires that the coefficients be relatively prime."; PartialFactor::ecm = "Trying ECM method."; Begin["`Private`"]; Unprotect[AliquotSequence]; Scan[(SetAttributes[#, Listable]) &, {AbundantQ, AdditionChain, AliquotSequence,AmicableQ,CFFactor,MethodA,MethodAStar, PrimitiveRoots, PrattCertificate,PrimePowerQ, EulerPhiLowerBound, SumOfFourSquares, PellSolve, EulerPhiUpperBound, ContinuedFraction, OrderMod,QuadraticIrrationalQ, CarmichaelLambda, LambdaLog, LambdaExp, Nonresidue,PhiLog, PhiExp, Convergents, Mersenne,Euclid, CarmichaelQ, Repunit, ExponentForm, PhiInverse, PhiMultiplicity, PseudoprimeQ, StrongPseudoprimeQ, Tower, TowerMod, LegendreSymbol, SqrtModAll, SqrtModPrimePower, SqrtModPrime, SpspSequence,LucasCertificate, GaussianNorm, CarmichaelLambda, GaussianMod, GaussianPrimeQ, GaussianIntegerQ, QuarticResidueSymbol, ClassNumber, KroneckerSymbol}]; If[!v4Q, SetAttributes[IntegerExponent, Listable]]; If[!v4Q, ChineseRemainder[x__] = ChineseRemainderTheorem[x]]; \[Phi] = EulerPhi; \[Pi][x_?NumericQ] = PrimePi[x]; \[Sigma][n_] := DivisorSigma[1,n]; Subscript[\[Sigma], k_Integer][n_] := DivisorSigma[k,n]; Subscript[SuperMinus[\[Sigma]], k_Integer][n_Integer] := DivisorSigma[k,n] - n^k; SuperMinus[Subscript[\[Sigma],k_]][n_Integer] := DivisorSigma[k,n] - n^k; SuperMinus[\[Sigma]][n_Integer] := DivisorSigma[1,n]-n; CarmichaelQ[n_Integer] := MemberQ[{561,1105,1729,2465, 2821,6601,8911,10585,15841,29341,41041,46657,52633, 62745,63973,75361,101101,115921,126217,162401,172081, 188461,252601,278545,294409,314821,334153,340561, 399001,410041,449065,488881,512461,530881,552721, 656601,658801,670033,748657,825265,838201,852841, 997633},n] /; n <= 1000000; CarmichaelQ[n_Integer] := (OddQ[n] && (!MemberQ[fi = FactorInteger[n], {_, _?(#>1&)}]) && Length[fi] >= 3 && Mod[n-1, CarmichaelLambda[n]] == 0 ); Euclid[n_] := 1 + Times @@ Prime[Range[n]]; Mersenne[n_] := 2^n - 1; MersennePrimeQ[3] = True; MersennePrimeQ[1] = False; MersennePrimeQ[n_Integer] := IntegerQ[m = Log[2, n+1]] && Nest[Mod[#^2 - 2, n] & , 4, m - 2] == 0; Options[Repunit] = {Base -> 10}; Repunit[n_, d_Integer:1, opts___Rule] := (b = Base /. {opts} /. Options[Repunit]; d (b^n - 1)/(b-1)); PhiLog[n_] := Length[FixedPointList[EulerPhi, n]] - 2; PhiExp[6] = 41; PhiExp[7] = 83; PhiExp[8] = 137; PhiExp[9] = 257; PhiExp[10] = 641; PhiExp[11] = 1097; PhiExp[12] = 17*137; PhiExp[13] = 17*257; PhiExp[14] = 41*257; PhiExp[n_] := Module[{i = 2^(n-1)+1}, While[PhiLog[i] != n, i++]; i]; Format[ExponentForm[n_, b_:10]] := Module[{e = N[Log[b, Abs[n]], 20], h, ff, nn = Abs[n]}, h = Round[e]; e = If[h < e + 0.4, h, h - 1]; Sign[n]*If[(ff = IntegerExponent[nn, b]) != 0 && nn/b^ff < Sqrt[nn], nn/b^ff*ToString[b]^ff, SFTF[ToString[b]^e + (nn - b^e)]]]; ExponentForm /: N[HoldPattern[ExponentForm[n_, b_:10]], prec_:6] := N[n,prec]; ExponentForm /: Round[HoldPattern[ExponentForm[n_, b_:10]]] := n; ExponentForm /: ToExpression[HoldPattern[ExponentForm[n_, b_:10]]] := n; (*ExponentForm[n_, b_:10] := Module[{e = N[Log[b,Abs[n]],20],h,ff,nn=Abs[n]}, h= Round[e]; e=If[hFalse}; FactorForm[n_Integer, opts___] := FactorForm[FactorInteger[n, opts]]; FactorForm[n_?GaussianIntegerQ, opts___] := FactorForm[FactorInteger[n, opts]]; FactorFormfor[l_List] := SequenceForm @@ (DisplayForm[RowBox[{#1, " "}]] & ) /@ MapThread[ If[#2 == 1, If[Head[#1] === Complex, StringJoin["(", ToString[#1], ")"], ToString[#1]], If[Head[#1] === Complex, SuperscriptBox[ StringJoin["(", ToString[#1], ")"], #2], SuperscriptBox[ToString[#1], ToString[#2]]]] & , Transpose[l]] /. {"(I)" -> "I", "(-I)" -> "-I"}; FactorFormfor[x_?Negative] := SequenceForm["-",FactorForm[FactorInteger[-x]]]; Format[HoldPattern[FactorForm[x_]]] := FactorFormfor[x]; FactorForm /: N[HoldPattern[FactorForm[l_]],prec_:16] := ( N[Unfactor[l], prec]) ; FactorForm /: Round[HoldPattern[FactorForm[l_]]] := ( Unfactor[l]); FactorForm /: ToExpression[HoldPattern[FactorForm[l_]]] := ( Unfactor[l]); If[!v4Q, Clear[ChineseRemainderTheorem]; (* the general CRT code is by Dan Lichtblau, Mike McGeachie and Craig Ortner, Macalester students *) PairwiseCoprimeQ[m_List] := Apply[LCM,m] ==Apply[ Times , m]; ChineseRemainder[a_List, m_List] := If[ Max[Abs[m]] < 3000000 && LCM @@ m == Times @@ m, CRTStandard[a, m], Catch[First[CRTlist[a, m, 1, Length[a]]]]] /; Length[a] == Length[m]; CRTStandard[a_, m_] := Module[{mm = Times @@ m, m1}, m1 = mm/m; Mod[a . (m1 * PowerMod[m1, -1, m]), mm]]; CRTpair[{a1_, a2_}, {m1_, m2_}] := Block[ {mm, d = GCD[m1, m2]}, If[Mod[a1 - a2, d] == 0, mm = m1/d; {Mod[a1 + (a2 - a1)*mm*PowerMod[mm, -1, m2/d], mm*m2], mm*m2}, Throw[{}]]]; CRTlist[a_, m_, first_, last_] /; last - first <= 16 := Block[{j, val = a[[first]], prod = m[[first]]}, Do[{val, prod} = CRTpair[{val, a[[j]]}, {prod, m[[j]]}], {j, first + 1, last}]; {val, prod}]; CRTlist[a_, m_, first_, last_] := Block[{len = Quotient[last - first, 2], c1, c2, m1, m2}, {c1, m1} = CRTlist[a, m, first, first + len]; {c2, m2} = CRTlist[a, m, first + len + 1, last]; CRTpair[{c1, c2}, {m1, m2}]]]; If[!v4Q, Convergents[q_Rational] := Module[{st = Take[FullGCD[Numerator[q], Denominator[q], UseGridBox->False, ExtendedGCDValues -> True], -2]}, -(1/Divide @@ (Drop[#, If[st[[1,3]]*st[[2,3]] == 0, 3, 2]] &) /@ st)]; Convergents[x_,m_] := Take[Convergents[ FromContinuedFraction[ContinuedFraction[x,m+1]]], m]]; QuadraticConjugate[d_] := d /. {Sqrt[n_] :> -Sqrt[n], 1/Sqrt[n_] :> -(1/Sqrt[n])}; Normal[PeriodicCFForm[initial_List:{}, periodic_List]] := Module[{raw, persoln, x, quadConj}, persoln = Solve[x == ( Normal[ContinuedFractionForm[Append[periodic, x]]]),x]; persoln = persoln[[If[persoln[[1,1,2]] < 0, 2, 1],1,2]]; raw=Together[ Normal[ContinuedFractionForm[Append[initial, persoln]]]]; quadConj = QuadraticConjugate[Denominator[raw]]; RootReduce[(Numerator[raw]*quadConj)/ Expand[Denominator[raw]*quadConj]]]; Format[PeriodicCFForm[i_:{}, p_]] := StringForm[ " \[LeftAngleBracket]\[LeftAngleBracket] `` | `` \[RightAngleBracket]\[RightAngleBracket]",Sequence@@(StringReplace[ StringDrop[StringDrop[ToString[#],1],-1],","->""]&/@ {i,p})]; If[!v4Q, qi[_Integer|_Rational|Sqrt[_Integer]|Sqrt[_Rational]] = True; qi[x:(_Plus|_Times)]:= Apply[And,Map[qi,Apply[List,x]]]; qi[_] = False; qi[_Integer|_Rational|Sqrt[_Integer]|Power[_Integer, -1/2] |Sqrt[_Rational]|Power[_Rational, -1/2]] = True; QuadraticIrrationalQ[n_] := Block[ {$RecursionLimit = Infinity}, If[ !NumericQ[n] || NumberQ[n], False, qi[RootReduce[n]]] ]]; (* Extract the (normalised) components of the quadratic irrational with general form (p + x*Sqrt[d])/q. The result lifts the coefficient field to the integers, so that the recurrence relations for p and q avoid rational arithmetic. By Mark Sonofriou. *) splitqi[ rt:(_. *Sqrt[_]) ]:= {0, Denominator[rt], Numerator[rt]}; splitqi[ qrecip_.*(p_. + (rt:(_. *Sqrt[_]))) ]:= {#*p, #^2, #*rt}& @ Denominator[qrecip]; QNF[1/Sqrt[n_Integer]] := {0, n, n}; QNF[a_Integer/Sqrt[n_Integer]]:={0,n a^2,Sign[a]n}; QNF[Sqrt[n_Integer]] :={0, n, 1}; QNF[c_] := Module[{nf = splitqi[RootReduce[c]]}, {Sign[nf[[3]]]Abs[nf[[2]]]nf[[1]], nf[[3]]^2 nf[[2]]^2, Sign[nf[[3]]] nf[[2]] Abs[nf[[2]]]}]; QuadraticNormalForm[a_] := QNF[RootReduce[a]]; (* Following implements an iterative method of searching for the periodic part, as described in Rosen, Chapter 10 *) CFSpecial[0,1,d_] := (* this is mark s's code *) Module[{a, a0, aend, cfl, p, q, k}, cfl = {}; p = 0; q = 1; a0 = IntegerPart[Sqrt[d]]; a = a0; aend = 2*a; (* Marker for end of cycle *) k=0; While[k++;If[Mod[k,100]==0,cfl=Flatten[cfl]]; True, p = a*q - p; q = Quotient[d - p^2, q]; a = Quotient[p + a0, q]; cfl = {cfl, a}; If[ a == aend, Break[] ]; ]; {a0, Flatten[cfl]} ]; CFSpecial[PP_, QQ_, d_, o_:Slow] := Module[{kk,a,P,Q}, a := Floor[(P + Sqrt[d])/Q]; {P,Q} = current = {PP,QQ}; past = {{P,Q}}; pasta = {a}; kk=0; While[P = a Q - P; Q = (d - P^2)/Q; Position[past, {P,Q}] == {}, kk++; past = {past, {P,Q}}; pasta = {pasta, a}; If[Mod[kk,100] == 0,pasta = Flatten[pasta], past = Partition[Flatten[past],2]]]; past = Partition[Flatten[past],2]; pos = Position[past, {P,Q}][[1,1]]; pasta = Flatten[pasta]; Append[pasta[[Range[pos-1]]], Drop[pasta, pos-1]]]; ReducedQuadraticIrrationalQ[c_?QuadraticIrrationalQ] := c > 1 && -1 < QuadraticConjugate[c] < 0; Options[FullGCD] = {UseGridBox -> True, ExtendedGCDValues -> False, Labels -> False, GaussianIntegers->False, FontSize->10, FontColor->GrayLevel[0], ShowSmithValues -> False,HighlightFontColor->RGBColor[1,0,0], Background->GrayLevel[0.9]}; FullGCD[a_, b_, opts___] := Module[{gridQ, data, giQ,n,ff}, {corQ,giQ,gridQ, fs, bg, fc, exQ, labQ, hcol} = {ShowSmithValues, GaussianIntegers, UseGridBox, FontSize, Background, FontColor, ExtendedGCDValues, Labels, HighlightFontColor} /. {opts} /. Options[FullGCD]; gridQ = gridQ || labQ; data = If[giQ || !FreeQ[{a,b}, Complex], Map[ If[gridQ, MakeBoxes, Identity], FullGCD[a, b, {a}, {}, GI],{2}], (ff=FullGCD[a, b, {a}, {}])]; If[corQ && Mod[b^2, a] == a-1, {cor1, cor2} = Select[ff[[1]], # < Sqrt[a]&,2]]; n = Length[data[[1]]]; If[gridQ, PrependTo[data[[1]], "Remainders"]; data[[2]] = Append[Join[{"Quotients", ""}, data[[2]]], ""]; PrependTo[data[[3]], "s"]; PrependTo[data[[4]], "t"]; If[labQ === True || labQ === Short, data = { Join[{data[[1,1]], MakeBoxes @ StringForm["`` = a = ``", Subscript["r",0],a], MakeBoxes @ StringForm["`` = b = ``", Subscript["r",1],b]}, MapIndexed[MakeBoxes @ StringForm["`` `` ``", If[aa=(labQ === Short && #2[[1]] < n-4), "", Subscript["r",#2[[1]]+1]],If[aa,"","="],#1]& , Drop[data[[1]], 3]]], Join[data[[2,{1,2}]], MapIndexed[ MakeBoxes @ StringForm["`` `` ``", If[aa = (labQ === Short && 1 < #2[[1]] < n-3), "", Subscript["q",#2[[1]]]], If[aa,"","="],#1]&, Drop[Drop[data[[2]],2],-1]], {data[[2,-1]]} ], Join[data[[3,{1}]], MapIndexed[ MakeBoxes @ StringForm["`` `` ``", If[aa=(labQ === Short && 2 < #2[[1]] < n-2), "", Subscript["s",#2[[1]]-1]],If[aa,"","="],#1] & , Rest[data[[3]]] ]], Join[data[[4,{1}]], MapIndexed[ MakeBoxes @ StringForm["`` `` ``", If[aa=(labQ === Short && 2 < #2[[1]] < n-2), "", Subscript["t",#2[[1]]-1]],If[aa,"","="],#1] & , Rest[data[[4]]] ]]} ]; cols=If[exQ, data, Drop[#,0]& /@ Drop[data, -2]]; colorrule=(x_ /; (Abs[x] == cor1 || Abs[x]==cor2)) :> StyleBox[x, FontColor->hcol, FontWeight->"Bold"]; DisplayForm[ StyleBox[GridBox[Transpose[If[Length[cols]==4, {cols[[1]]/. colorrule, cols[[2]],cols[[3]], cols[[4]]/. colorrule}, cols[[{1,2}]]]], FilterOptions[GridBox, opts], GridFrame -> 2, System`GridFrameMargins -> {{2, 3}, {0.5, 1}}, RowLines -> Join[{2, 1}, Flatten[Table[{1, 1, 1, 1, 2}, {1+Length[data[[1]]]/5}]]], ColumnLines -> 1, ColumnAlignments -> If[labQ === True, {Right, Right, Left, Left},Right], ColumnSpacings -> 3], Background -> bg, FontSize -> fs, FontColor->fc]], If[exQ, data, {Drop[data[[1]],1], data[[2]]}]]]; FullGCD[a_, b_, r_List, q_List] := FullGCD[b, Mod[a, b], {r, b}, {q, Quotient[a, b]}] /; b != 0; FullGCD[a_, 0, r_List, q_List] := Module[{s = {1, 0}, t = {0, 1}}, Scan[(t = Flatten[{t, t[[-2]] - #1*t[[-1]]}]; s = Flatten[{s, s[[-2]] - #1*s[[-1]]}]) & , Flatten[q]]; Flatten /@ {Append[r,0], q, s, t}]; FullGCD[a_, b_, r_List, q_List, GI] := Module[{ qu = Round[a/b]}, FullGCD[b, a - b qu, {r, b}, {q, qu}, GI]] /; b !=0; FullGCD[a_, 0, r_List, q_List, GI] := Module[ {s = {1, 0}, t = {0, 1}}, Scan[(t = Flatten[{t, t[[-2]] - # t[[-1]]}]; s = Flatten[{s, s[[-2]] - # s[[-1]]}]) & , Flatten[q]]; Flatten /@ {Append[r, 0], q, s, t}]; Options[FullContinuedFraction] = {FontSize->12, UseGridBox -> True, Background->RGBColor[1,1,.6], HighlightFontColor->RGBColor[1,0,0]}; FullContinuedFraction[c_Rational?Positive, opts___Rule] := Module[{fs,bg, gQ}, {fs, gQ,bg} = {FontSize, UseGridBox, Background} /. {opts} /. Options[FullContinuedFraction]; con = Convergents[c]; cf = ContinuedFraction[c]; If[cf[[1]] == 0, PrependTo[con, 0]]; If[!gQ, Return[Transpose[{cf, con}]]]; DisplayForm[StyleBox[ GridBox[Prepend[Transpose[ {cf, MakeBoxes /@ If[con[[1]]==0,Rest[con], con]}], {"Partial Quotients","Convergents"}], GridFrame->2, RowLines->1, ColumnLines->1], Background->bg, FontFamily->"Times",FontSize->fs]]]; FullContinuedFraction[c_?(NumericQ[#] && Positive[#] && Head[#] =!= Rational &), n_Integer, opts___Rule] := Module[{fs}, {fs,bg} = {FontSize, Background} /. {opts} /. Options[FullContinuedFraction]; cf = ContinuedFraction[c, n+1]; con = Convergents[FromContinuedFraction[cf]]; nn = Min[n, Length[cf]]; DisplayForm[StyleBox[GridBox[ Prepend[Transpose[ {Take[cf,nn], MakeBoxes /@ Take[If[cf[[1]]==0, Prepend[con,0], con],nn]}], {"Partial Quotients","Convergents"}], GridFrame->2, RowLines->1, ColumnLines->1], Background->bg, FontFamily->"Times", FontSize->fs]]]; fixrule = {FractionBox[z_,"1"] :> z, RowBox[{"0", "+", z_}] :> z}; FullContinuedFraction[c_?QuadraticIrrationalQ, opts___] := Module[{s,ss, bg,gQ,fs,P,Q,d,fc,jj}, {fs,gQ,bg, fc} = {FontSize,UseGridBox, Background, HighlightFontColor} /. {opts} /. Options[FullContinuedFraction]; {P, d, Q} = QuadraticNormalForm[c]; ans = CFSpecial[P,Q,d, Slow]; ans = {Drop[ans, -1], Flatten[ans[[-1]]]}; jj = Join[ans[[1]],ans[[2]], If[Length[ans[[2]]] == 1,ans[[2]], ans[[2,{1,2}]]]]; If[jj[[2]]== 1, AppendTo[jj, 1]]; conv = Convergents[Normal[ContinuedFractionForm[jj]]]; (*Return[{temp, conv, pasta, past}]; *) conv = Take[If[pasta[[1]] == 0, PrependTo[conv,0], conv], Length[pasta]+1]; {PP,QQ} = Last[past]; PP = Last[pasta] QQ - PP; QQ = (d - PP^2)/QQ; AppendTo[past, {PP,QQ}]; AppendTo[pasta, Floor[(PP+Sqrt[d])/QQ]]; {s,ss} = ToString /@ {PP,QQ}; If[!gQ, Return[{past, pasta, conv}]]; Do[pasta = ReplacePart[pasta,StyleBox[ (pasta[[i]]), FontWeight->"Bold", FontSize -> fs, FontColor -> fc], i], {i, Position[past, {PP,QQ}][[1,1]], Length[pasta]-1}]; DisplayForm[StyleBox[GridBox[Prepend[ Append[Transpose[ {(MakeBoxes /@ past) /. RowBox[{"{",RowBox[{s,",",ss}],"}"}] :> StyleBox[MakeBoxes[#]& @ {s,ss}, FontWeight -> "Bold",FontColor->fc], MakeBoxes /@ ((Together@ Expand[((#[[1]] +Sqrt[d])/ #[[2]])]; DisplayForm[ FractionBox[RowBox[{ToString[#[[1]]],"+",SqrtBox[d]}],ToString[ #[[2]]]] //. (fixrule)]) & /@ past), pasta, MakeBoxes /@ conv}], MakeBoxes /@ {"", "Periodic Form:", PeriodicCFForm @@ ans,""}], {"{X,Y}", "\[Alpha]", "Partial Quotients", "Convergents"}], GridFrame->2, GridFrameMargins->{{.5,.5},{1,.5}}, RowSpacings->Join[{2},Table[1, {Length[pasta]-1}],{3}], RowLines->Join[{2},Table[1, {Length[pasta]-1}],{2}], ColumnLines->1], FontFamily->"Times",FontSize->fs, Background->bg]]]; Options[PellSolutions] = {RightSide -> {-1, 1}}; PellSolutions[d_, {m_, n_}, opts___Rule] := (rhs = RightSide /. {opts} /. Options[PellSolutions]; mat = PellSolve[d]; rhs = Flatten[Sort[{rhs}]]; rough = First /@ NestList[#1 . mat & , mat, 2*n]; ss = Select[rough, MemberQ[rhs, #1^2 . {1, -d}] & ]; If[ss == {}, {}, Take[ss, {m, n}]]); PellSolutions[d_, n_Integer:1, opts___Rule] := (ans=PellSolutions[d, {n, n}, opts]; If[ans=={}, ans, ans[[1]]]); (* the old method, basic binary. MatrixPowerMod[mat_, n_, m_Integer] := Mod[Fold[ Mod[If[#2 == 1,#1.#1.mat,#1.#1],m] &, mat, Rest[IntegerDigits[n,2]]],m] MatrixPowerMod[mat_, n_, Infinity] := Fold[ If[#2 == 1,#1.#1.mat,#1.#1] &, mat, Rest[IntegerDigits[n,2]]]; *) (* new method; base 32, optimized for 100 digit numbers *) initializec[k_] := Module[{e}, Clear[c]; c[0] = {0, 0}; Do[c[i] = (e = IntegerExponent[i,2]; {e, i/2^e}), {i, 2^k-1}]]; initializec[5]; MatrixPowerMod[mat_, n_, m_]:= MatrixPowerMod[5][mat,n,m]; MatrixPowerMod[k_][mat_,n_,m_Integer]:=Module[{sq,odds,b}, sq[mmm_] :=Mod[mmm.mmm,m]; b = sq[mat]; odds[0] = IdentityMatrix[Length[mat]]; odds[1] = Mod[mat,m]; Do[odds[i] = Mod[odds[i-2].b,m], {i, 3, 2^k-1,2}]; Fold[Nest[sq,Mod[Nest[sq, #1,k-c[#2][[1]]].odds[c[#2][[2]]],m], c[#2][[1]]]&, Nest[sq,odds[c[#[[1]]][[2]]],c[#[[1]]][[1]]], Rest[#]]& [IntegerDigits[n,2^k]]]; MatrixPowerMod[k_][mat_,n_,Infinity]:=Module[{sq,odds,b}, sq[mmm_] := mmm.mmm; b = sq[mat]; odds[0] = IdentityMatrix[Length[mat]]; odds[1] = mat; Do[odds[i]=odds[i-2].b, {i, 3, 2^k-1,2}]; Fold[ Nest[sq,Nest[sq, #1,k-c[#2][[1]]]. odds[c[#2][[2]]],c[#2][[1]]]&, Nest[sq,odds[c[#[[1]]][[2]]],c[#[[1]]][[1]]], Rest[#]]& [IntegerDigits[n,2^k]]]; Clear[LinearRecurrence]; LinearRecurrence[coeffs_, init_, n_Integer, m_Integer] := Mod[First[MatrixPowerMod[ Append[RotateRight /@ Drop[IdentityMatrix[Length[coeffs]], -1], coeffs], n, m]] . init, m] /; n >= Length[coeffs]; LinearRecurrence[coeffs_, init_, n_Integer, m_Integer] := Mod[ init[[n+1]],m] /; n < Length[coeffs]; LinearRecurrence[coeffs_, init_, n_Integer, Infinity] := First[MatrixPower[ Append[RotateRight /@ Drop[IdentityMatrix[Length[coeffs]], -1], coeffs], n]] . init /; n >= Length[coeffs]; LinearRecurrence[coeffs_, init_, n_Integer, Infinity] := init[[n+1]] /; n < Length[coeffs]; LinearRecurrence[coeffs_, init_, n_Integer] := LinearRecurrence[coeffs, init, n, Infinity]; LinearRecurrence[coeffs_, init_, n_List, m_:Infinity] := LinearRecurrence[coeffs, init, #, m]& /@ n; coun[s_] := First@Rest@Position[Characters[s], _?(# != " "&)] weight[d_] := First[Rest[Flatten[(d /. (x_String :> coun[x]))[[1,1,1]]]]] ; EuclideanAlgorithmLength[a_,b_] := Module[{aa = a, bb = b, i = 0}, While[bb != 0, {aa,bb,i} = {bb, Mod[aa,bb],i+1}]; i]; SumOfFourSquares[n_] :=Module[{m=Mod[n,4], r, s, a, b}, Sort[ Which[ n==0,{0,0,0,0}, m==0, 2 SumOfFourSquares[n/4] , m==2, While[!PrimeQ[n-a^2-b^2], a= Random[Integer,{0, Sqrt[n]}]; b= Random[Integer,{0, Sqrt[n-a^2]}]]; Flatten[{a,b, Abs[SumOfSquaresRepresentations[2, n-a^2-b^2][[1]]]}], True,{r,s,a,b}= Sort[SumOfFourSquares[2n],Mod[#1,2] -1]; StrongPseudoprimeQ[b_,n_] := !PrimeQ[n] && spspTest[b,n] /; OddQ[n]&&n > 2; StrongPseudoprimeQ[b_,n_?EvenQ] := False; If[True || (* to fix student reported bug *)!v4Q, spspTest[b_, n_?OddQ] := Module[{bSeq, m, r}, {m,r} = OddPart[n-1]; bSeq = NestList[PowerMod[#,2,n] &,PowerMod[b, m, n],r]; ((Union[bSeq]=={1}) || MemberQ[Drop[bSeq,-1], n-1])], spspTest[b_, n_?OddQ] := Module[ {r = IntegerExponent[n - 1, 2], start}, start = PowerMod[b, (n - 1)/2^r, n]; start == 1 || start == n - 1 || n - 1 == NestWhile[ PowerMod[#, 2, n]&, start, # != n-1 &, 1, r-1]]]; OrderMod[a_, 1] := 1; OrderMod[a_, m_] := Select[ Divisors[EulerPhi[m]], PowerMod[a, #, m] == 1 &, 1] /. {x_} :> x; OrderMod[a_, p_, e_] := Select[ Divisors[p^(e-1) (p-1)], PowerMod[a, #, p^e] == 1 &, 1] /. {x_} :> x; (* PrimePowerQ[n_] := Length[FactorInteger[n]] == 1;*) \[Lambda][n_] := CarmichaelLambda[n]; CarmichaelLambda[2, 1] = 1; CarmichaelLambda[2, 2] = 2; CarmichaelLambda[2, n_]:= 2^(n-2); CarmichaelLambda[m_, n_]:= m^(n-1)*(m-1); If[!v4Q, CarmichaelLambda[0] = 0; CarmichaelLambda[1] = 1; CarmichaelLambda[n_Integer]:= Apply[LCM, Apply[CarmichaelLambda, FactorInteger[Abs[n]], 1]]]; badQ[n_] := OddPart[n] == {1, _?(# > 3&)}; LambdaLog[m_] := Length[ FixedPointList[CarmichaelLambda, m]]-2 Options[LambdaExp] = {TraceInterval -> Infinity}; LambdaExp[n_, start_Integer:1, opts___Rule] := Module[ {i = Max[2^(n-1), start]}, ti = TraceInterval /. {opts} /. Options[LambdaExp]; While[If[Mod[i, ti] == 0, Print[i]]; !((PrimeQ[i] ||i == 3^(n-1))&& LambdaLog[i] == n), i++]; i]; (* LambdaExp[16] = 14348907 *) (* 3.1 note: maybe OK in new version. There were power of two bugs. *) Clear[PrimitiveRoot]; PrimitiveRoot[p_Integer] := Module[{res}, res /; NumberQ[ res = pr[p] ] ]; pr[2] = 1; pr[3] = 2; pr[4] = 3; pr[n_?badQ] := $Failed; pr[n_]:= Module[{n2, res}, n2 = Quotient[n,2]; res = pr[n2]; If[SameQ[res,$Failed], $Failed, If[OddQ[res], res, res + n2]] ] /; Mod[n, 4] == 2; pr[p_?PrimeQ] := Module[{g, flist}, g = 2; flist = (p-1)/Map[ First, FactorInteger[p-1] ]; While[JacobiSymbol[g, p] == 1 || Scan[If[PowerMod[g, #1, p] == 1, Return[True]]&, flist], g++]; g ]; pr[n_] := Module[{q, fi, g}, fi = FactorInteger[n]; If[ Length[fi] > 1, Return[$Failed]]; q = fi[[1,1]]; g = pr[q]; If[PowerMod[g, q - 1, q^2] == 1, g += q]; g] /; GCD[PowerMod[2, n, n] - 2, n] > 1 && !badQ[n]; Q[{}] = 1; Q[{q_}] := q; Q[q_]:= q[[1]] Q[Rest[q]]+Q[Drop[q,2]]; Q[q_, d_:1] := d Q[q]; FullContinuant[q_, d_:1] := (ContinuantQ[#, d] & ) /@ Reverse[FoldList[Join, {}, List /@ Reverse[q]]]; MakeBoxes[OrderedContinuant[args___], StandardForm] := MakeBoxes[Plus[args]]; Format[OrderedContinuant[args___], InputForm] := Plus[args]; Format[OrderedContinuant[args___], TraditionalForm] := StandardForm[OrderedContinuant[args]]; ContinuantQ[q_, d_:1] := Module[{ans = Expand[Q[q, d]]}, If[ListQ[List @@ ans], OrderedContinuant @@ Sort[ List @@ ans,( Length[#1] < Length[#2] || ( Length[#1] == Length[#2] && Order[List@@#1 ,List @@#2]==1) ) &], ans]]; AdditionChain[1 | 0, FactorMethod] := {1}; AdditionChain[n_, FactorMethod] := Module[{p = FactorInteger[n][[1,1]]}, Join[AdditionChain[p - 1, FactorMethod], p*AdditionChain[If[n == p, 0, n/p], FactorMethod]]]; AdditionChain[n_, BinaryRightToLeft] := Union[FoldList[Plus, 1, Rest[ Flatten[2^(Position[Reverse[IntegerDigits[n, 2]], 1] - 1)]]], 2^Range[Log[2, n]]]; AdditionChain[n_, BinaryLeftToRight] := Fold[Flatten[{#1, 2*#1[[-1]], If[#2 == 1, 2*#1[[-1]] + 1, {}]}] & , {1}, Rest[IntegerDigits[n, 2]]]; factorMethod[a_, 1 | 0, m_] := Module[{ans = mod[a,m]}, AppendTo[ansdata, ans]; ans]; factorMethod[a_, n_, m_] := Module[{p = FactorInteger[n][[1,1]], ans}, AppendTo[fdata, If[n == p, n, {p, n/p}]]; ans = factorMethod[a factorMethod[a, p-1, m], If[n == p, 0, n/p], m]; ans] Options[PowerAlgorithm] = {Method -> BinaryLeftToRight, FontSize -> 10, Background->GrayLevel[0.9]}; mod[x_, Infinity] := x; mod[x_, m_] := Mod[x, m]; (* the -1 below is a kludgy way to get Infinity as the default and still control the type of the argument *) PowerAlgorithm[a_, n_, mm_Integer:-1, opts___Rule] := Module[{ans = 1, sq = a, m = mm /. -1 -> Infinity, e = n, sz, met, par, data}, ansdata = fdata = {}; {met,sz,bg}={Method,FontSize, Background} /. {opts} /. Options[PowerAlgorithm]; fo[s_] := StyleBox[s, FontSize -> sz+2, FontWeight->"Bold"]; bo[s_] := StyleBox[s, FontWeight -> "Bold", FontSize -> sz+2]; Which[met === BinaryRightToLeft, data = {fo /@ {"Exponent", "Parity", (MakeBoxes[#1, TraditionalForm] &)[ If[m == Infinity, "Repeated Squares", StringForm["Squares Mod ``", m]]], "Answer"}}; While[e != 1, AppendTo[data, {e, Mod[e, 2], sq, ans}]; If[Mod[e, 2] == 1, ans = mod[ans * sq, m]]; e = Quotient[e, 2]; sq = mod[sq*sq, m]]; AppendTo[data, {1, 1, sq,ans}]; AppendTo[data, {"", "", "",mod[sq ans,m]}]; DisplayForm[StyleBox[GridBox[Join[{data[[1]]}, GridBoxFix[ Rest[data]]], GridFrame -> 1, RowLines -> {2, 1}, ColumnAlignments -> Center, ColumnLines -> 1], FontSize -> sz, Background -> bg]], met === BinaryLeftToRight, col1 = Flatten[{1,(digs = Rest@IntegerDigits[n,2]) /. 1 -> {1, ""}}]; col2 = {a}; Scan[(cur = Last[col2]; temp = mod[cur cur, m]; col2 = Flatten[{col2, If[# == 0, temp, {temp, mod[temp a, m]}]}])&, digs]; DisplayForm[ StyleBox[GridBox[Prepend[ GridBoxFix[ Transpose[{col1, col2}]], {bo @ MakeBoxes @ StringForm["Binary Digits of ``",n], bo @ "Answer"}], GridFrame -> 1, RowLines -> {2, 1}, ColumnAlignments -> Center, ColumnLines -> 1], FontSize -> sz, Background -> bg]], met === FactorMethod, factorMethod[a,n,m]; {lf,la} = Length /@ {fdata, ansdata}; data = If[lf >= la, {fdata, Join[ansdata, Table["", {lf - la}]]}, {Join[fdata,Table["", {la - lf}]], ansdata} ]; DisplayForm[ StyleBox[ GridBox[Prepend[ GridBoxFix[ Transpose[data]], bo /@ {"The Recursive Calls","Answer"}], GridFrame -> 1, RowLines -> {2, 1}, ColumnAlignments -> Center, ColumnLines -> 1], FontSize -> sz, Background -> bg]]]]; LargestExponent[n_] := Max[Last /@ FactorInteger[n]]; LastFactorialDigit[0 | 1] := 1 ; LastFactorialDigit[n_] := Module[{rProd = 1, m = n, q = 0}, While[m > 0, rProd = Mod[rProd * Mod[m, 5]!, 10]; m = Quotient[m, 5]; q = Mod[q + m, 4] ]; Mod[(2^q /. 1 -> 6) * rProd, 10]]; Options[CRTGrid] = {FontSize -> 10, Background->GrayLevel[0.9], FontColor->GrayLevel[0]}; CRTGrid[m_, n_, opts___] := Module[{fs, hf = HoldForm[CRTGrid[m, n]]}, If[GCD[m, n] != 1, Message[CRTGrid::badgcd]; Return[hf]]; {bg,fs,fc} = {Background,FontSize, FontColor} /. {opts} /. Options[CRTGrid]; Show[Graphics[{{bg, Rectangle[{-1, -m}, {n, 1}]}, Table[Text[ToString[ChineseRemainder[{i, j}, {n, m}]], {i + 0.5, -j - 0.5}, TextStyle -> {FontSize -> fs,FontColor -> fc}], {i, 0, n - 1}, {j, 0, m - 1}], Table[Line[{{i, 1}, {i, -m}}], {i, 0, n}], Table[Line[{{-1, j}, {n, j}}], {j, -m, 0}], {Line[{{0, 0}, {-1, 1}, {n, 1}}], Line[{{-1, 0}, {-1, -m}}]}, {Thickness[0.004], Line[{{0, -m}, {0, 0}}], Line[{{0, 0}, {n, 0}}], Line[{{-1, -m}, {-1, 1}, {n, 1}, {n, -m}, {-1, -m}}]}, Table[Text[i, {i + 0.5, 0.5}, TextStyle :> {FontSize -> fs + 2, FontColor -> fc,FontWeight -> "Bold"}], {i, 0, n - 1}], Table[Text[j, {-0.5, -j - 0.5}, TextStyle :> {FontSize -> fs + 2, FontColor -> fc,FontWeight -> "Bold"}], {j, 0, m - 1}], Text[Subscript["\[DoubleStruckCapitalZ]", m], {-0.62, 0.25}, TextStyle -> FontSize -> fs + 2], Text[Subscript["\[DoubleStruckCapitalZ]", n], {-0.28, 0.71}, TextStyle -> FontSize -> fs + 2]}], AspectRatio -> Automatic, TextStyle -> FontColor -> fc, PlotRange -> {{-1.1, n+.1}, {-m - 0.2, 1.1}}]]; Options[VisualPowerGrid] = {SquareColors -> Automatic, PrimitiveRootColor -> RGBColor[1,0,0], EulerPhiLines -> False, FontSize ->9, ShowPrimitiveRoots -> False, UseArrows -> True, ShowNumbers->False, ReducedResiduesOnly->False}; Box[{x_, y_}, {x1_, y1_}, shade_, fsty_:{}] := Module[{fs = If[!ListQ[fsty], {fsty}, fsty]}, {{shade, Rectangle[{x, y}, {x1, y1}]}, If[fs === {None}, {},Join[ fs,{Line[{{x,y},{x1,y},{x1,y1},{x,y1},{x,y}}]}]]}]; VisualPowerGrid[ bounds_List:{}, n_Integer, opts___Rule] := Module[{large, font, data, a, i, a1, i1, pos}, {a1, i1} = If[bounds == {}, {n-1, n-1}, bounds]; {euQ,cols, fs, prQ, numQ, prc, arrQ, redQ} = {EulerPhiLines,SquareColors, FontSize, ShowPrimitiveRoots, ShowNumbers, PrimitiveRootColor, UseArrows, ReducedResiduesOnly} /. {opts} /. Options[VisualPowerGrid]; prs = If[prQ, PrimitiveRoots[n],{}]; rr=ReducedResidues[n]; e = EulerPhi[n]; e1= 1+If[!redQ,a1,e Quotient[a1,n] + Length[Select[rr, # If[numQ, Table[GrayLevel[1], {n+1}], Join[{GrayLevel[1], Hue[0]}, Table[Hue[(i-1)/If[n!=4, (n-4), n-3]/2, .6, .6], {i,1,n-3}], {RGBColor[1,1,0]}]]; data = Table[{a, i, PowerMod[a,i,n]}, {i, i1}, {a, a1}]; fontSmall = fs /. Automatic -> If[n <= 10, 12, 9]; fontLarge = fs+2; Show[Graphics[{ Box[{1, -a1}, {i1 + 1, 0}, cols[[2]], Thickness[0.002]]; data1=If[redQ, DeleteCases[data, {a_,_,_} /; GCD[a,n]!=1, Infinity],data]; data1 /. { {a_, i_, 0} :> {cols[[1]], Rectangle[{i, -pos[a]}, {i + 1, -pos[a] + 1}]}, {a_, i_, 1} :> {cols[[2]], Rectangle[{i, -pos[a]}, {i + 1, -pos[a] + 1}]}, {a_, i_, n - 1} :> {cols[[n]], Rectangle[{i, -pos[a]}, {i+1, -pos[a]+1}]}, {a_,i_,(x_)?(#>1 &)} :> {RGBColor[0,0,x/(n-2)]; cols[[x+1]], Rectangle[{i,-pos[a]},{i+1,-pos[a]+1}]}}, If[numQ, data /. {a_, i_, nn_} :> Text[ToString@nn,{i+.5,-a+.5}],{}], Thickness[0.002], Table[Line[{{1, y},{i1+1,y}}], {y, -If[redQ,e1,a1], 0}], Table[Line[{{x,-If[redQ,e1,a1]},{x,0}}], {x, 1, i1+1}], Text["a", {1 - 2 i1/20, 0}, {1, 1}, TextStyle :> {FontSize->fontLarge}], If[euQ, {Thickness[0.006],GrayLevel[1],Table[ Line[{{x,-If[redQ,e1,a1]},{x,0}}],{x,e+1, i1, e}]} , {}], {Thickness[0.002], Table[Line[{{1, y},{i1+1,y}}], {y, -If[redQ,e1,a1], 0}]}, Table[ {pp = MemberQ[prs, Mod[a,n]]; col = If[pp, prc, GrayLevel[0]]; extra = " \[MediumSpace]"; If[arrQ && pp, extra = "\[Rule]\[MediumSpace]"]; If[redQ && GCD[a,n]!=1, {}, Text[StringJoin[extra,ToString[a]], {1-n/200, -pos[a]+0.5}, {1,0}, TextStyle -> {FontWeight->If[pp,"Bold","Plain"], FontColor -> col, FontSize->fontSmall}]]}, {a, a1}], Table[Text[ToString[i], {i+0.5, 1/45}, {0,-1}], {i, i1}], Text["i", {1-i1/20, 1/45}, {1, -1}, TextStyle :> { FontSize->fontLarge}]}], FilterOptions[Graphics, {opts}], PlotRange ->All, TextStyle->{FontFamily->"Courier", FontSize-> fontSmall}, AspectRatio -> Automatic]]; Options[PowerGrid] = {FontSize ->12}; PowerGrid[n_, opts___] := Module[{fs}, fs = FontSize /. {opts} /. Options[PowerGrid]; DisplayForm[StyleBox[ GridBox[Prepend[ MapIndexed[Prepend[#1, #2[[1]]]&, Table[ PowerMod[a,i,n] /. {n-1-> -1, 1 -> StyleBox["1", FontWeight->"Bold"]}, {a,1,n-1},{i,1, n-1}]]/. {z_/;(z==1 ||z==2),x__} :> { StringJoin["a = ",ToString[z]],x}, Prepend[ Range[n-1]," i ="]], RowLines->{2,0.6}, RowAlignments->"Center", RowSpacings->{2,1.8}, ColumnSpacings->{1.8,1}, ColumnAlignments->Center, ColumnLines->Prepend[Flatten@ Table[{0,0,0,0,1}, {Quotient[n,5]}],2], GridFrame->2, System`GridFrameMargins->{{1,0.6},{0.6,0.8}}], FontSize->fs, FontFamily->"Times", Background->GrayLevel[0.9]]]]; fix1[l_, m_] := (tt=Flatten@Transpose[{l, " "& /@ l}]; Join[sp=Table[" ", {(m- 2 Length@l)/2}],tt,sp]); Options[PascalTriangle] = {FontSize -> 10, Background->GrayLevel[0.9]}; bl[s_] := StyleBox[ToString[s], FontColor->GrayLevel[0]]; PascalTriangle[n_, opts___] := Module[{fs}, {fs,bg} = {FontSize,Background} /. {opts} /. Options[PascalTriangle]; DisplayForm@StyleBox[GridBox[{{StyleBox[ GridBox[MapIndexed[Prepend[#1, StyleBox[#2[[1]]-1, FontWeight ->"Bold", FontColor->GrayLevel[0]]]& , data = Table[fix1[#,2n+2]& [bl/@Binomial[i,Range[0,i]]],{i,0,n}]], GridFrame->0, System`GridFrameMargins -> {{0,1}, {1,1}}, RowLines->Join[{0},Flatten@Table[{0,0,0,0,1},{1+n/5}]], ColumnLines->{1,0}, ColumnAlignments->{Right, Center}, ColumnSpacings->{1,0}, System`ColumnWidths->{2, ((Ceiling@Log[10,1+Max@Cases[data, _Integer, Infinity]])/ 1.6)}], Background->bg,FontColor->GrayLevel[0.5], FontSize->fs, FontFamily -> "Times"]}}, GridFrame->1, System`GridFrameMargins->{{.2,.2},{1,1}}], FontColor->GrayLevel[0], Background->bg]]; pad[x_, n_] := Join[x, Table["", {n - Length[x]}]]; Options[PseudoprimeTable] = {FontSize -> 12}; PseudoprimeTable[bmax_, nmax_, opts___Rule] := Module[ {data}, fs = FontSize /. {opts} /. Options[PseudoprimeTable]; DisplayForm[StyleBox[GridBox[data = Table[{Prime[i], Select[Range[nmax], PseudoprimeQ[Prime[i], #1] & ]}, {i, bmax}]; max = Max[Length /@ Last /@ data]; data = data /. {z_, y_List} :> Join[{z}, pad[y, max]]; Prepend[GridBoxFix[data], pad[{StyleBox["b", FontWeight -> "Bold"], ""}, max + 1]], ColumnAlignments -> Left, GridFrame -> 2, RowLines -> Prepend[Flatten[Table[{0,0,0,0,1},{5}]], 2], ColumnLines -> {2, 1}, ColumnAlignments -> Right], FontSize -> fs, Background -> GrayLevel[0.9]]]]; Options[FunctionTable] = {Background -> GrayLevel[0.9], FontSize -> 12, Columns -> 5, FunctionHeading -> Automatic}; FunctionTable[f_, n_Integer, opts___Rule] := FunctionTable[f, {1, n}, opts]; FunctionTable[fcn_, {n0_, n1_}, opts___] := ( {fs, bg, cols, fh} = {FontSize, Background, Columns, FunctionHeading} /. {opts} /. Options[FunctionTable]; fh = (fh /. Automatic -> MakeBoxes[#] & )[TraditionalForm[HoldForm[fcn["n"]]]]; ht = 1 + Quotient[n1-n0, cols]; DisplayForm[StyleBox[GridBox[Prepend[Transpose[Flatten[Table[ {r = Range[n0 + ht*(i - 1), n0 + ht*i-1], (StringJoin[ToString[fcn[#]], " "] & ) /@ r}, {i, cols}], 1]], Flatten[Table[{"n", fh}, {i, cols}]]], GridFrame -> 2, RowLines -> Flatten[Table[{2, 1, 1, 1, 1}, {ht/5+1}]], ColumnAlignments -> Right, ColumnLines -> Flatten[Table[{1, 2}, {cols}]], (FilterOptions[GridBox, opts])], Background -> bg, FontSize -> fs, FontFamily -> Times]]); Average[fcn_, x_] := (sum=0.; Do[sum+=fcn[i],{i,x}]; sum/x); AverageList[fcn_, x_] := Rest[ FoldList[Plus, 0., fcn /@ Range[x]]]/Range[x]; SummatoryList[f_, pow_, {x_, a_, b_, step_:1}, start_:0] := Module[{s = start, data = {}}, Do[s += N[f]; If[Mod[x, step] == 0, data = {data, {x, s}}], {x, a, b}]; Cases[data, {xx_Integer, y_} :> {xx, N[y/pow /. x -> xx]}, Infinity]]; Unfactor[n_] := Times @@ Apply[Power, Round@n, {1}] ; EulerPhiUpperBound[n_] := (c = 1; p = 1; While[p <= n, c++; p *= Prime[c] - 1]; Floor[n/Apply[Times,1 - 1/Prime[Range[c-1]]]]); EulerPhiInverseDivisors[n_] := Module[{fi = FactorInteger[n], pr}, pr = First /@ FactorInteger[n]; ({#, 1 + If[!MemberQ[pr, #], 0, Cases[fi, {#, y_} :> y][[1]]]} & ) /@ Union[Select[Divisors[n] + 1, PrimeQ]]]; PhiCandidates[l_,_, fi_] := Table[ {{2, k}, Max[1,2^(k-1)]}, {k, 0, 1+fi[[1,-1]]}] /; Length[l] == 1; (* fi is the factorization of max; both are passed and used *) (* the use of logs below seems silly...one could keep a running list of the prime expoinents with each candidate and use that perhaps...but maybe not worth effort now *) (* Two PhiCandidatess are used, as the initial selection can be stricter *) PhiCandidates[l_,max_, fi_] := Module[{ r = PhiCandidates[Drop[l, -1],max, fi], p, e, pn,ss}, {p, e} = Last[l]; Join[r, Flatten[ Select[r, Mod[max,#[[-1]] * (p-1)] ==0 &] /. xx_List?(MatrixQ[Drop[#,-1], IntegerQ]&) :> (pn = xx[[-1]] * (p-1);xr = Drop[xx, -1]; Table[Join[xr, {{p, e1},pn p^(e1-1)}], {e1,(cc=Cases[fi, {p, _}]; If[cc == {}, 1,cc=cc[[1]]; ccc= cc[[1]]^cc[[2]]; 1+cc[[2]] - Log[p,GCD[xx[[-1]],ccc]]])}]), 1]]]; PhiCandidates2[l_,max_, fi_] := Module[{ r = PhiCandidates[Drop[l, -1],max, fi], p, e, pn,ss}, {p, e} = Last[l]; Join[Cases[r, {__,max}], Select[r, Mod[max,#[[-1]] * (p-1)] ==0 &] /. xx_List?(MatrixQ[Drop[#,-1], IntegerQ]&) :> (top = (cc=Cases[fi, {p, _}]; If[cc == {}, 1,cc=cc[[1]]; ccc= cc[[1]]^cc[[2]]; 1+cc[[2]] - Log[p,GCD[xx[[-1]],ccc]]]); Join[Drop[xx, -1], {{p, top},xx[[-1]] * (p-1) * p^(top-1)}])]] ; PhiInverse[1] = {1, 2}; PhiInverse[n_?OddQ] := {}; PhiInverse[n_] := (l = EulerPhiInverseDivisors[n]; Sort[Unfactor /@ (Cases[ DeleteCases[ DeleteCases[PhiCandidates2[ l,n, FactorInteger[n]], {2,0}, Infinity], {}], {x__,n} :> {x}])]); PhiMultiplicity[n_] := Length[PhiInverse[n]]; Solve[EulerPhi[x_Symbol] == n_Integer, x_] := {x->#}& /@ PhiInverse[n]; Solve[EulerPhi[x_Symbol] == n_Integer] := {x->#}& /@ PhiInverse[n]; AmicableQ[n_] := (s = SuperMinus[\[Sigma]][n]) != n && SuperMinus[\[Sigma]][s] == n; DeficientQ[n_] := SuperMinus[\[Sigma]][n] < n; AbundantQ[n_] := SuperMinus[\[Sigma]][n] > n; AliquotSequence[n_, max_:Infinity] := (ans = {n}; m = SuperMinus[\[Sigma]][n]; While[m < max && FreeQ[ans, m], ans = {ans, m}; m = SuperMinus[\[Sigma]][m]]; Flatten[ans]); PrimeFactorial[n_] := Times @@ Prime[Range[n]]; Options[LinearDiophantineSolve] = {SolutionType -> All, Assumptions -> None}; LinearDiophantineSolve[{a_Integer}, c_Integer, opts___Rule] := Module[{type, hf = HoldForm[LinearDiophantineSolve[{a}, c, opts]]}, {type,ass} = {SolutionType, Assumptions} /. {opts} /. Options[LinearDiophantineSolve]; If[ass =!= None && type === All, Message[LinearDiophantineSolve::toomny]; Return[hf]]; ass = ass /. None -> (True &); Which[Mod[c, a] != 0, {}, type === All, {{c/a}}, type === Positive && Positive[c/a], Select[{{c/a}}, ass], type === NonNegative && NonNegative[c/a], Select[{{c/a}}, ass], True, {}]]; (* the 2-vble case is handled by a standard formula *) LinearDiophantineSolve[{a_Integer, b_Integer}, c_, opts___Rule] := Module[ {e = ExtendedGCD[a,b], type, l, an, ass, hf = HoldForm[LinearDiophantineSolve[{a,b},c,opts]]}, {type, ass} = {SolutionType, Assumptions} /. {opts} /. Options[LinearDiophantineSolve]; If[ass =!= None && type === All, Message[LinearDiophantineSolve::toomny]; Return[hf]]; ass = ass /. None -> (True &); If[IntegerQ[c], If[Mod[c, e[[1]]] != 0, Return[{}]]]; Which[ type === NonNegative, Select[LDSNonneg[{a,b}, c], ass], type === Positive, l = LDSNonneg[{a,b}, c]; If[l === $Failed, l, Select[l, (ass[#] && !MemberQ[#,0])&]], type === All, an = Transpose[{c * e[[2]], {b, -a}}/e[[1]]]; an . {{1,0}, {If[Head[an[[1,1]]]===Integer && Head[an[[1,2]]]===Integer, -Round[an[[1,1]]/an[[1,2]]],0],1}}]]; LDSNonneg[vec_, c_] := Module[{ans}, ans = LinearDiophantineSolve[vec, c]; vbles = Table[Unique[x], {Length[ans] - 1}]; parsoln = ans . Prepend[vbles, 1]; is = InequalitySolve[And @@ Thread[parsoln >= 0], vbles]; pars = Check[InequalityToSet[is, vbles], Message[LinearDiophantineSolve::infsol];$Failed]; If[pars === $Failed, pars, Sort[(parsoln /. Thread[vbles -> #] &) /@ pars]]]; (* the general case is handled by replacing the last two variables, a y + b z by a single vble w and using recursion. The matrix produced recrusively is then modified. Note that symbolic computation is used here as symbols are used in the r.h.s. of the inductive case *) LinearDiophantineSolve[vec_List, c_Integer, opts___Rule] := Module[ {recur, type, l, cc, new, v,ass, hf = HoldForm[LinearDiophantineSolve[vec,c,opts]]}, {type,ass} = {SolutionType,Assumptions} /. {opts} /. Options[LinearDiophantineSolve]; If[ass =!= None && type === All, Message[LinearDiophantineSolve::toomny]; Return[hf]]; ass = ass /. None -> (True &); Which[ type === NonNegative, Select[LDSNonneg[vec, c], ass], type === Positive, l = LDSNonneg[vec, c]; If[l === $Failed, l, Select[l, (ass[#] && !MemberQ[#,0])&]], type === All, If[Mod[c, GCD @@ vec] != 0, Return[{}]]; v = Table[Unique[n], {Length[vec]-2}]; recur = LinearDiophantineSolve[Append[Drop[vec,-2], GCD @@ Take[vec,-2]],c]; new = Expand[LinearDiophantineSolve[Take[vec, -2] / GCD@@Take[vec,-2], recur[[-1]] . Prepend[v,1]]]; Join[Append[#,0] & /@ Drop[recur,-1], Table[Prepend[Append[ cc = Coefficient[new[[i,1]],#]& /@ v, new[[i,2]]], new[[i,1]] - cc.v], {i,2}]]]] /; Length[vec] > 2 Clear[InequalityToSetv4]; InequalityToSet = If[$VersionNumber < 4, InequalityToSetv3, InequalityToSetv4]; InequalityToSetv4[HoldPattern[Or[a__]], vbles_] := Module[{temp}, temp = (InequalityToSetv4[#1, Intersection[vbles, Cases[#, _Symbol, Infinity]]] &) /@ {a}; If[MemberQ[temp, $Failed], $Failed, Union @@ temp]]; InequalityToSetv4[HoldPattern[And[a_]], vbles_] := Module[{ineq, vble, first, ans}, ineq = Select[List @@ a, Length[Intersection[Cases[#1, _Symbol], vbles]] == 1 &, 1][[1]]; vble = Intersection[Cases[ineq, _Symbol], vbles][[1]]; first = ToSet[ineq]; ans = ({#1, InequalityToSetv4[And[a] /. vble -> #1, DeleteCases[vbles, vble]]} & ) /@ first; ans = DeleteCases[ans ,{aa_Integer, d_[False]}]; ans = Flatten[ans/.{aa_Integer,x_List}\[RuleDelayed] If[VectorQ[x],(Prepend[{#},aa]&/@ x),Prepend[#,aa]&/@ x],1]; ans ]; InequalityToSetv4[a_, {x_}] := ToSet[a]; ToSet[HoldPattern[Or[a__]]] := Module[{temp = (ToSet[#1, vbles] & ) /@ {a}}, If[MemberQ[temp, $Failed], $Failed, Union @@ temp]]; ToSet[HoldPattern[And[a__]]] := Module[{temp = ToSet /@ {a}, infints = Cases[temp, _Inequality]}, Select[Intersection @@ Cases[temp, List], good[#1, infints] & ]]; ToSet[HoldPattern[Inequality[a_, Less | LessEqual, c_]]] := $Failed; ToSet[HoldPattern[Equal[a_, b_]]] := If[IntegerQ[b], {b}, {}]; good[n_, infints_] := And @@ (good[n, #1] & ) /@ infints; good[n_, Inequality[a_, b_, c_]] := Inequality @@ {a, b, n}; ToSet[HoldPattern[Inequality[a_, LessEqual, b_, LessEqual, c_]]] := Select[Range[Ceiling[a], c], IntegerQ]; ToSet[HoldPattern[Inequality[a_, Less, b_, Less, c_]]] := DeleteCases[Select[Range[Ceiling[a], c], IntegerQ], a | c]; ToSet[HoldPattern[Inequality[a_, LessEqual, b_, Less, c_]]] := DeleteCases[Select[Range[Ceiling[a], c], IntegerQ], c]; ToSet[HoldPattern[Inequality[a_, Less, b_, LessEqual, c_]]] := DeleteCases[Select[Range[Ceiling[a], c], IntegerQ], a]; InequalityToSetv3[ineq_, v_] := Module[{is1,c,f}, ineq1=ineq /. { HoldPattern[LessEqual[x_, y_, z_]] :> Inequality[x, LessEqual, y, LessEqual,z], HoldPattern[ Less[x_, y_, z_] ] :> Inequality[x, Less, y, Less,z]}; is1 = List @@ Apply[List, If[Depth[ineq1] == 2, {ineq1}, ineq1], 1]; If[Union[Length /@ is1] != {5}, Message[LinearDiophantineSolve::infsol];Return[$Failed]]; easy = Select[is1, NumericQ[#1[[1]]] && NumericQ[#1[[-1]]] & ]; hard = Complement[is1, easy]; easy = easy /. {(a_Integer) | (a_Rational), z_, x_, y_, (b_Integer) | (b_Rational)} :> Prepend[{Ceiling[a], Floor[b]} + Which[ z === y === LessEqual, {0, 0}, z === Less && y === LessEqual && IntegerQ[a], {1, 0}, z === Less && y === LessEqual && !IntegerQ[a], {0, 0}, z === LessEqual && y === Less && IntegerQ[b], {0, -1}, z === LessEqual && y === Less && !IntegerQ[b], {0, 0}, z === y === Less && IntegerQ[a] && !IntegerQ[b], {1, 0}, z === y === Less && !IntegerQ[a] && IntegerQ[b], {0, -1}, z === y === Less && IntegerQ[a] && IntegerQ[b], {1, -1}, z === y === Less && IntegerQ[a] && !IntegerQ[b], {0, 0}], x]; hard = hard /. {a_, z_, x_, y_, b_} :> ({c, f} = {Ceiling[a], Floor[b]}; Prepend[Which[ z === y === LessEqual, {c, f}, z === Less && y === LessEqual, {c+Hold[If[iQ[a], 1, 0]], f}, z === LessEqual && y === Less, {c, f + Hold[If[iQ[b], -1, 0]]}, z === y === Less, {c + Hold[If[iQ[a], 1, 0]], f + Hold[If[iQ[b], -1, 0]]}], x]); tablearg = Prepend[ Join[easy, hard], If[Length[v] == 1, v[[1]], v]]; Flatten[Table @@ ReleaseHold[tablearg], Length[v] - 1]]; iQ[x_?NumericQ] := IntegerQ[x] Clear[Conductor]; (*Conductor[aa_List] := Module[{d = FoldList[GCD,0,aa], hf = HoldForm[Conductor[aa]]}, If[d[[-1]] != 1, Return[Infinity]]; If[LinearDiophantineSolve[ {aa[[1]],aa[[2]]}/d[[3]],aa[[3]]/d[[-1]], SolutionType->NonNegative] == {}, Message[Conductor::cantdo]; Return[hf], 1 + Sum[(d[[k]]/d[[k + 1]] - 1)*a[[k]], {k,3}]]] /; Length[aa] == 3 && And @@ (NonNegative /@ aa);*) Conductor[{x_Integer?NonNegative, y_Integer?NonNegative}] := If[GCD[x, y] == 1, (x - 1)(y - 1), Infinity]; LinearCongruenceSolve[{a_, b_}, m_] := Module[{e = ExtendedGCD[a, m]}, Sort[Flatten@ List[Mod[Which[ Mod[b, e[[1]]] != 0, {}, e[[1]] == 1, b*PowerMod[a, -1, m], True, (b*e[[2,1]] + m*Range[e[[1]]])/e[[1]]], m]]]]; ReducedResidues[n_] := Select[Range[n], GCD[#1, n] == 1 & ]; Global`SuperStar[Subscript[\[DoubleStruckCapitalZ], n_]] := ReducedResidues[n]; PowerMod[a_Integer, n_Integer, Infinity] := a^n; (* xx Fixed in version 3.5 *) LogStar[a_, n_] := Block[{$MaxExtraPrecision = 10000}, If[n < a, 0, LogStar[a, Log[a, n]] + 1]]; Tower[a_, n_] := If[n == 0, 1, a^Tower[a, n - 1]]; (* (a_) \[UpArrow] (n_) := Tower[a, n]; UpArrow[a_,n_] := Tower[a,n]; *) TowerMod[a_, n_, Infinity] := a \[UpArrow] n; TowerMod[a_, n_, 1] := 0; TowerMod[a_, 0, m_] := 1 /; m > 1; TowerMod[a_, 2, m_] := PowerMod[a, a, m]; TowerModPrimePower[_, 0, {p_, _}] := 1 /; p > 1; TowerModPrimePower[_, _, 1] := 0; TowerModPrimePower[a_, 2, {p_, e_}] := PowerMod[a, a, p^e]; TowerModPrimePower[a_, n_, {p_, e_}] := Module[{g, pp = p^e}, g = GCD[a, pp]; If[g > 1 && n - 1 > LogStar[a, Log[g, p^e]], 0, Mod[If[g == 1, 1, g^Tower[a,n-1]]* PowerMod[a/g, TowerMod[a, n - 1, CarmichaelLambda[p, e]], pp], pp]]]; TowerMod[b_, n_, m_] := Block[{$RecursionLimit=1000}, Module[{fi = FactorInteger[m],c}, ChineseRemainder[(TowerModPrimePower[b, Min[n, Ceiling[Log[2, #[[1]]^#[[2]]]]], #] &) /@ fi, Power @@ Transpose[fi]] ]]; EuclidSequence[1] = {2}; EuclidSequence[n_] := Flatten[ {EuclidSequence[n - 1], FactorInteger[Times @@ EuclidSequence[n - 1] + 1][[1,1]]}]; Options[RandomPrime] = {Certified -> True}; RandomPrime[d_, con_List:{}, opts___Rule] := Module[{n=Infinity}, cQ = (Certified /. {opts} /. Options[RandomPrime]); If[!cQ, While[Log[10., n] > d, n = NextPrime[Random[Integer,10^(d-{1,0})], con /.{} -> Sequence @@ {}]]; n, CertifiedPrime[d, con]]]; RandomPrimePrimeQ[d_, con_List:{}] := Module[{n=Infinity}, While[Log[10., n] > d, n=NextPrime[Random[Integer,10^(d-{1,0})], con /.{} -> Sequence @@ {}]]; n]; Eratosthenes[n_] := Module[{L = Range[2, n], i = 1}, While[L[[i]]^2 <= n, If[L[[i]] != 0, Do[L[[k]] = 0, {k, i^2 + 2*i, n - 1, L[[i]]}]]; i++]; Rest[Union[L]]]; Options[EratosthenesTable] ={FontSize->9, SurvivorCountOnly->False, Background->GrayLevel[0.9]}; EratosthenesTable[n_, opts___] := Module[{L = Range[2, n], ll, i=1}, ll = {L}; {scQ,fs,bg} = {SurvivorCountOnly, FontSize, Background} /. {opts} /. Options[EratosthenesTable]; If[scQ, Return[SieveOfE[n]]]; While[L[[i]]^2 <= n,If[L[[i]] != 0, Do[L[[k]] = 0, {k, i^2 + 2*i, n - 1, L[[i]]}]; AppendTo[ll, L]]; i++]; DisplayForm[StyleBox[GridBox[Transpose[ll] /. 0 -> "", ColumnLines -> 1, ColumnAlignments -> Right, RowSpacings -> Flatten[Table[{0.32, 0.32, 0.32, 1.2, 0.32}, {n/5}]], RowLines -> Flatten[Table[{0, 0, 0, 1, 0}, {n/5}]], GridFrame -> 1.5], FontSize -> fs, FontFamily -> Times, Background -> bg]]]; SieveOfE[x_] := Module[{ans,ans1}, ans=Range[2,x]; ans1={x-1}; Do[p=Prime[i]; ans=DeleteCases[ans, z_/;z!=p && Mod[z,p]==0]; ans1={ans1, Length[ans]}, {i, PrimePi[Sqrt[x]]}]; Flatten[ans1] ]; Options[PrattCertificate] = {Tree -> False, PrimeQTest -> True, AssumedPrimeBound -> 2, Background->GrayLevel[0.9], HighlightFirstLevel -> False, ShowWitnesses -> False, Show2s -> True, FontSize -> 12, FullPrattCertificate -> True}; PrattCertificate[2, level_Integer:0, opts___] := 2; PrattCertificate[p_Integer, level_Integer:0, opts___Rule] := Module[{w = 2, primes, treeQ, skipQ, specBound, show2Q}, {treeQ, skipQ, specBound, witQ, show2Q, fs, hfl, fullQ,bg} = {Tree, !PrimeQTest, AssumedPrimeBound, ShowWitnesses, Show2s, FontSize, HighlightFirstLevel, FullPrattCertificate, Background} /. {opts} /. Options[PrattCertificate]; input = HoldForm[PrattCertificate[p, opts]]; If[specBound > 2 && witQ, Message[PrattCertificate::nowitn]; Return[input]]; If[treeQ && witQ, Message[PrattCertificate::notree]; Return[input]]; witF = If[witQ && level == 0, {#1, Flatten[#1 //. {a_, b_, c_} :> {b, c}]} & , Identity]; If[!fullQ, specBound = Max[specBound, 5]]; If[p <= specBound, Return[p]]; If[!skipQ && !PrimeQ[p], Message[PrattCertificate::noprim, p]; Return[input]]; {primes, exps} = Transpose[FactorInteger[p - 1]]; If[!fullQ, prod = 1; i = Length[primes] + 1; While[prod^2 < p - 1, i--; prod *= primes[[i]]^exps[[i]]]; primes = Drop[primes, i - 1]]; If[level == 0, elim2sfcn[d_] := If[show2Q, d, DeleteCases[d, _List? (Union[Characters[#1[[1]]]] == {"\"", " ", "2"} & ), Infinity]]]; While[If[fullQ, PowerMod[w, (p - 1)/2, p] != p - 1, PowerMod[w, p - 1, p] != 1] || (temp = False; Do[If[If[fullQ, If[primes[[i]] == 2, False, PowerMod[w, (p - 1)/primes[[i]], p] == 1], GCD[PowerMod[w, (p - 1)/primes[[i]], p] - 1, p] != 1], temp = True; Break[]], {i, Length[primes]}]; temp), w++]; out = elim2sfcn[witF[If[treeQ && level == 0, CertificateTreeForm[#1, 0, 0, FontSize -> fs, Background->bg] &, Identity][ {p, w, PrattCertificate[primes, 1, PrimeQTest -> False, opts]} //. {2, x__} :> {If[treeQ && level == 0, List, Identity][2], x}]]]; If[!treeQ && !show2Q, out= DeleteCases[ out //. {2, x___}:> {x} , {}, Infinity]]; out //. If[hfl, {(s_String)?(coun[#1] == {weight[out]} & )} :> {StyleBox[s, FontWeight -> "Plain", FontColor -> RGBColor[1, 0, 0]]}, {}]] /; p > 2; indentFcn[ind_, x_] := StringJoin[Table[" ", {ind}], ToString[x]]; Options[CertificateTreeForm] = {FontSize->12, Background->GrayLevel[.9]}; CertificateTreeForm[m_Integer, indent_, level_, opts___] := If[level == 0, TableForm, Identity][ {indentFcn[indent, m]}]; CertificateTreeForm[{m_Integer}, indent_, level_, opts___] := If[level == 0, TableForm, Identity][{indentFcn[indent, m]}]; CertificateTreeForm[{{m_Integer}, other__}, indent_, level_, opts___] := Prepend[Join @@ (CertificateTreeForm[#1, indent, level + 1] & ) /@ {other}, indentFcn[indent, m]]; CertificateTreeForm[{p_, w_Integer, other_}, opts___Rule] := CertificateTreeForm[{p,w,other}, 0,0,opts]; CertificateTreeForm[{p_, w_Integer, other_}, indent_, level_, opts___Rule] := If[level == 0, {fs,bg} = {FontSize, Background} /. {opts} /. Options[CertificateTreeForm]; DisplayForm[StyleBox[GridBox[ List /@ ToString /@ MakeBoxes /@ #1, ColumnAlignments -> Left, GridFrame -> 1, System`GridFrameMargins -> {{1, 1}, {1.5, 1.5}}, System`RowSpacings -> (0.6)], FontFamily -> "Courier", FontSize -> fs, FontWeight -> "Plain", Background -> bg]] &, Identity][ Prepend[Flatten[(CertificateTreeForm[#, indent + 1 + If[IntegerQ[p], Length[IntegerDigits[p]],1], level + 1] & ) /@ other, 1], StringJoin[indentFcn[indent, p], " \[Rule] ", ToString[w]]]]; coun[s_] := First[Rest[Position[Characters[s], _?(# != " " &)]]]; weight[d_] := First[Rest[Flatten[(d /. x_String :> coun[x])[[1,1,1]]]]]; li = LogIntegral; QuadraticResidueQ[b_, p_?PrimeQ] := JacobiSymbol[b,p]==1; LegendreSymbol[b_,p_?PrimeQ] := (qr = QuadraticResidueQ[b,p]; Which[qr,1,Mod[b,p] != 0 && !qr, -1, True, 0]); FermatNumber[n_] := 2^2^n+1; (* SqrtModPrime is based on Shanks's method *) SqrtModPrime[a_, 2] := Mod[a, 2]; SqrtModPrime[a_, p_] := {} /; JacobiSymbol[a, p] == -1; SqrtModPrime[a_, p_] := 0 /; Mod[a, p] == 0; SqrtModPrime[a_, p_] := Min[PowerMod[a, (p+1)/4, p]*{-1, 1} + {p, 0}] /; Mod[p, 4] == 3; SqrtModPrime[a_, p_] := Block[{k = (p - 1)/2, i, L, n, r, s = 1, c, b}, While[EvenQ[k], k /= 2; s++]; {c, r, n} = PowerMod[{Nonresidue[p], a, a}, {k, (k+1)/2,k}, p]; While[n != 1, L = s; b = n; Do[ If[b == 1, b = c; s = i - 1, b = Mod[b b, p]], {i, L}]; {c, r, n} = Mod[b {b, r, b n}, p]]; Min[r, p - r]] /; Mod[p, 4] == 1; (*SqrtModPrimeLehmer[n_, p_] := Module[{s, h}, s = Which[Mod[n, p] == 0, 0, p == 2, Mod[n, 2], !QuadraticResidueQ[n, p], {}, Mod[p, 4] == 3, PowerMod[n, (p + 1)/4, p], True, h = 1; While[ QuadraticResidueQ[h^2 - 4*n, p], h++]; Mod[(p + 1)/2*(MatrixPowerMod[{{0, 1}, {-n, h}}, (p + 1)/2, p] .{2, h})[[1]], p]]; If[IntegerQ[s], Min[s, p - s], s]];*) (* Exponent = 1 case reduces to previous routine *) SqrtModPrimePower[a_, p_, 1] := Module[ {root = SqrtModPrime[a, p]}, If[root==={}, {}, Union[Mod[{p-root, root}, p]]]]; (* Non-coprime case, any prime, reduces to coprime case *) SqrtModPrimePower[a_, p_, r_] := Module[{s = 1, j}, If[Mod[a, p^r] == 0, Return[Range[0, p^r-1, p^Ceiling[r/2]]]]; While[Mod[a, p^s] == 0, s++]; s--; If[OddQ[s], {}, Flatten@Table[ (SqrtModPrimePower[a/p^s, p, r-s] + (j-1) p^(r-s)) p^(s/2), {j, p^(s/2)}]]] /; Mod[a, p] == 0 && r > 1; (* Coprime case, odd prime *) SqrtModPrimePower[a_, p_, r_] := Module[{j, root = SqrtModPrime[a, p]}, If[!NumberQ[root], {}, Do[root = Mod[root - PowerMod[2 root, -1, p]*(root^2-a), p^j], {j, 2, r}]; Sort[{-1,1}*root + {p^r, 0}]]] /; p != 2 && Mod[a, p] != 0 && r > 1; (* Coprime case, prime is 2, no solution *) SqrtModPrimePower[a_, 2, r_] := {} /; (r==2 && Mod[a, 4] == 3) || (3 <= r && OddQ[a] && Mod[a, 8] != 1); (* Coprime case modulo 4 *) SqrtModPrimePower[a_, 2, 2] := {1, 3} /; Mod[a, 4] == 1; (* Coprime case, prime is 2, exponent is at least 3, solutions exist *) SqrtModPrimePower[a_, 2, r_] := Module[{root = Nest[Mod[(#^3 + (2-a)#)/2, 2^r]&, 1, r-3]}, Union[ {root, 2^r-root}, Mod[{root, 2^r-root} + 2^(r-1), 2^r]]] /; OddQ[a] && 3 <= r; SqrtModAll[a_, n_] := Module[{factors = Transpose@FactorInteger[n]}, rootlist = SqrtModPrimePower[a, factors[[1]], factors[[2]]]; If[MemberQ[rootlist, {}], Return[{}]]; PrependTo[rootlist, ChineseRemainder[List[##], factors[[1]]^factors[[2]]]&]; Union@Flatten[Outer @@ rootlist, Length[rootlist] - 2]] /; n > 1; SqrtModAll[a_, 1] := {0}; PairwiseCoprimeQ[m_List] := Apply[LCM,m]==Apply[ Times,m]; PreviousPrime[n_] := Module[ {i=Floor[n-1]}, While[!PrimeQ[i], i--];i] ; NextPrime[n_, {i_,m_}] := Module[ {k=Mod[n,m], j=n, ii = Mod[i,m]}, j += ii - k + If[ii < k, m, 0]; While[!PrimeQ[j], j += m];j]; NextPrime[n_] := NextPrime[N[n]] /; NumericQ[n]; L[1] := Log; L[k_] := L[k] = x /. First[NDSolve[{Derivative[1][x][t] == L[k - 1][t - 1]/t, x[k] == 0}, x, {t, k, 15.5}, PrecisionGoal -> 18, WorkingPrecision -> 25]]; KnuthTrabbPardo[1, a_ /; a <= 12] := 1 + Sum[(-1)^k*L[k][a], {k, 1, Floor[a]}]; KnuthTrabbPardo[2, a_ /; a <= 12] := 1 - Sum[(-1)^k*(k - 1)*L[k - 1][a],{k, 2, Floor[a]}]; KnuthTrabbPardo[x_, a_ /; a > 12] := (Message[KnuthTrabbPardo::toobig]; HoldForm[KnuthTrabbPardo[x,a]]); StringToInteger[message_String] := (ToCharacterCode[message] - 31) . (100^Range[StringLength[message] - 1, 0, -1]); IntegerToString[codedmessage_] := FromCharacterCode[ IntegerDigits[codedmessage, 100] + 31]; RSAEncode[message_String, e_Integer, n_Integer] := PowerMod[StringToInteger[message], e, n] /; 2 StringLength[message] <= Length[IntegerDigits[n]]-1; RSADecode[m_, k_Integer, n_Integer] := IntegerToString[PowerMod[m, k, n]] /; m <= n; RSAEncode[message_, e_Integer, n_Integer] := Module[ {size = Quotient[Length[IntegerDigits[n]] - 2,2]}, Map[RSAEncode[#, e, n]&, Apply[StringJoin, Partition[Join[Characters[message], Table[" ", {Mod[-StringLength[message],size]}]], size], {1}]]] /; 2 StringLength[message] > Length[IntegerDigits[n]]-1; RSADecode[m_List, k_Integer, n_Integer] := StringJoin @@ (RSADecode[#, k, n]& /@ m); Nonresidue[p_?OddQ]:=Module[{n=0}, While[JacobiSymbol[Prime[++n],p]==1]; Prime[n]]; (* algorithm that follows is due to H. Wilf *) DivisorSigma1[k_, {}] := 1; DivisorSigma1[k_, {{p_, e_}, x___}] := (1 + (-1)^e*p^(k*(e + 1)))/(1 + p^k)* DivisorSigma1[k, {x}] + p^k*(1 - p^(k*(2*(e - 2 + Mod[e, 2])/2 + 2)))/ (1 - p^(2*k))*DivisorSigma[k, Times @@ Apply[Power, {x}, {1}]] /; k != 0; DivisorSigma1[0, {{p_, e_}, x___}] := Module[{c, d}, c = Quotient[e, 2]; d = (e - 1 + Mod[e, 2] - 1)/2; If[Mod[p, 4] == 3, c + 1, 2*c + 1]*DivisorSigma1[0, {x}] + (d + 1)*(DivisorSigma[0, Times @@ Apply[Power, {x}, {1}]] - DivisorSigma1[0, {x}])]; DivisorSigma1[k_, n_Integer] := ( fi = FactorInteger[n]; n1 = Times @@ Apply[Power, Cases[fi, {p_ /; Mod[p, 4] == 1, _}], {1}]; n3 = Cases[fi, {p_ /; Mod[p, 4] == 3, _}]; DivisorSigma[k, n1]*DivisorSigma1[k, n3]) /; n > 10^7; DivisorSigma1[k_,n_]:= (Plus @@ (Select[Divisors[n], Mod[#, 4] == 1 &]^k))/; n <= 10^7; DivisorSigma3[k_, n_] := DivisorSigma[k, NumberTheory`NumberTheoryFunctions`Private`OddPart[n]] - DivisorSigma1[k, n]; (* Following code overrides what is in the NumberTheoryFunctions package. Delete if package changes? *) SumOfSquaresR[6, n_] := 16*n^2*(DivisorSigma1[-2, n] - DivisorSigma3[-2, n]) - 4*(DivisorSigma1[2, n] - DivisorSigma3[2, n]); AlternatingDivisorSigma[k_, {p_, e_}] := -(1-p^(k(e+1)))/(1-p^k) /; OddQ[p]; AlternatingDivisorSigma[k_, {2, e_}] := (1-2^(k(e+1)))/(1-2^k) - 2; AlternatingDivisorSigma[k_, n_Integer] := (fi = FactorInteger[n]; -(-1)^Length[fi] * Times @@ (AlternatingDivisorSigma[k, #]&) /@ fi); SumOfSquaresR[8, n_] := 16 (-1)^n AlternatingDivisorSigma[3, n]; SumOfSquaresR[4, n_Integer] := If[EvenQ[n], 24, 8]*DivisorSigma[1, NumberTheory`NumberTheoryFunctions`Private`OddPart[n]]; SumOfSquaresR[2, n_Integer?Positive] := SumOfSquaresR[2, n]= If[MemberQ[fi = FactorInteger[n], {p_ /; Mod[p, 4] == 3, _?OddQ}], 0, 4*Times @@ Cases[fi, {p_ /; Mod[p, 4] == 1, a_} :> a + 1]]; Options[GaussianPrimePlot] = {PrimeColor -> RGBColor[0, 0, 0], UnitsColor -> GrayLevel[0.6], CompositeColor -> GrayLevel[1], DiskColor -> GrayLevel[1]}; boxes[{x_, y_}] := Rectangle /@ {{x, y}, {y, x}, {-x, y}, {y, -x}, {x, -y}, {-y, x}, {-x, -y}, {-y, -x}}; GaussianPrimePlot[n_, opts___] := ({pc, uc, cc, dc} = {PrimeColor, UnitsColor, CompositeColor, DiskColor} /. {opts} /. Options[GaussianPrimePlot]; Show[Graphics[{{dc, Disk[{0, 0}, n + 0.5]}, Circle[{0, 0}, n + 0.5], {{uc, boxes[{1, 0}]}, {pc, (boxes[{Re[#], Im[#]}] & ) /@ Select[Flatten[Table[a + b*I, {a, 0, n}, {b, 0, a}]], Abs[#] <= n && GaussianPrimeQ[#]& ]}, If[dc === cc, {}, {cc, (boxes[{Re[#1], Im[#1]}] & ) /@ DeleteCases[Select[Flatten[Table[ a + b*I, {a, 0, n}, {b, 0, a}]], Abs[#] <= n && !PrimeQ[#, GaussianIntegers -> True] & ], 1]}]} /. {x_?NumericQ, y_} :> {x, y} - 0.5}], FilterOptions[Graphics, opts], AspectRatio -> Automatic]); AbsoluteDifferences[x_List] := Abs[Apply[Subtract, Partition[x, 2, 1], {1}]]; Options[IteratedPrimeGaps] = {FontSize -> 10, Background-> RGBColor[.7,1,.7], HighlightColor -> Hue[0]}; IteratedPrimeGaps[n_, iters_, opts___] := ({hc,bc,fs} = {HighlightColor, Background, FontSize} /. {opts} /. Options[IteratedPrimeGaps]; If[iters > n, Message[IteratedPrimeGaps::bigitr]; Return[HoldForm[IteratedPrimeGaps[n, iters]]]]; DisplayForm[StyleBox[GridBox[ Transpose[Prepend[Transpose[ (Take[#1, n] & ) /@ NestList[AbsoluteDifferences, Prime[Range[2*n]], iters]] /. {1 -> StyleBox["1", FontWeight-> "Bold", FontSize -> fs, FontColor-> hc], m_Integer :> ToString[m]}, Join[{" Primes"}, Table[StringJoin["Differences ", ToString[i]], {i, iters}]]] ], GridFrame -> 2, System`GridFrameMargins -> {{1, 1}, {2, 2}}, System`ColumnWidths -> {7, -0.8 + Length[IntegerDigits[n]]}, ColumnAlignments -> {Left, Center}, RowLines -> {1, 0}, RowSpacings -> {2}, ColumnLines -> {3, 1}], Background -> bc, FontSize -> fs, FontFamily -> "Times"]]); If[!v4Q, IntegerExponent[n_, b_Integer] := Module[{m=n,r=0}, While[Mod[m,b] == 0,m/=b;r++];r] /; b>1 && n>0; IntegerExponent[0, b_Integer] := Infinity]; Options[MillerRabinPrimeQ] = {RandomBases->True}; MillerRabinPrimeQ[n_Integer, d_Integer:1, opts___Rule] := Module[{bSeq, m, r}, If[1= 1; SRKp=37;SRKn=13;SRKa=SRKn-2; SRKW=Table[Mod[i+SRKa,SRKp],{i,SRKn}]; SRKV=Table[Mod[i (i+SRKa),SRKp],{i,SRKn}]; ToMod37Integer[s_String] := ( c=First[ToCharacterCode[s]] - 55; Which[c>9, c, -8 True}; SRKCorrect[s_String] := Module[{}, s1=S1[s]; s2=S2[s]; sde=Select[Range[13], Mod[s1 SRKV[[#]] - s2 SRKW[[#]],37]==0&]; te = Select[Range[12], Mod[s1 (SRKV[[#+1]]-SRKV[[#]]) - s2(SRKW[[#+1]]-SRKW[[#]]),37]==0&]; Which[ s1 ==s2 == 0, s, Length[sde]==1, StringReplacePart[ s, FromMod37Integer[ Mod[-s1 PowerMod[SRKW[[sde[[1]]]],-1,37] + ToMod37Integer[ StringTake[s, sde]],37]], sde[[{1,1}]]], Length[te]==1,StringReplacePart[ s, StringReverse[StringTake[s, {te[[1]], te[[1]]+1}]], {te[[1]], te[[1]]+1}]]]; FirstSPSPWitness[n_, start_:2] := Module[{b = start}, If[PrimeQ[n], Message[FirstSPSPWitness::noncom,n]; Return[HoldForm[FirstSPSPWitness[n,start]]]]; While[StrongPseudoprimeQ[b, n],b++]; b]; If[!v4Q, FromContinuedFraction[x_?VectorQ] := Normal[ContinuedFractionForm[x]]; FromContinuedFraction[x_List] := Normal[ PeriodicCFForm[Drop[x,-1],x[[-1]]]] /; Depth[x] > 2]; cfSqrt[d_] := Module[{P = 0, Q = 1, a, a0 = Floor[Sqrt[d]], periodic = {}}, a = a0; While[a !=2 a0, P = a*Q - P; Q = Quotient[d - P^2, Q]; a = Quotient[P + a0, Q]; periodic = {periodic, a}]; {a0, Flatten[periodic]}] /; !IntegerQ[Sqrt[d]] (* If[!v4Q, (* the following xx calls use the CF4 *) (*ContinuedFraction[x_Rational] := ContinuedFraction[x][[1]];*) ContinuedFraction[Sqrt[x_Integer]] := cfSqrt[x] /; !IntegerQ[Sqrt[x]]; ContinuedFraction[x_Integer] :=ContinuedFraction[x][[1]]; ContinuedFraction[x_?QuadraticIrrationalQ] := ( i = InputForm[ContinuedFraction[x]]; Append[i[[1,1]], i[[1,2]]]); ContinuedFraction[x_, n_] := ContinuedFraction[x,n][[1]]; Attributes[ContinuedFraction] = Listable]; *) Options[HarmonicRearrangement] = {MaxTerms -> 70, Background -> RGBColor[1,1,.6], ShowSeries -> True, ShowPartialSums -> True, ShowSignChanges -> True, FontSize -> 12}; (* g and collapse are auxiliary functions used to turn a list of series terms into a list of partial sums at the sign changes. *) HRg[s_List, a_] := If[Sign[a] == Sign[Last[Last[s]]], Append[Drop[s, -1], Append[Last[s], a]], Append[s, {a}]]; collapse[series_] := Rest[FoldList[Plus, 0, (Plus @@ #) & /@ Fold[HRg, {{First[series]}}, Rest[series]]]]; HarmonicRearrangement[target_, opts___] := Module[ {maximum, outputSeries, partials, signChanges, currentSum = 0., j = 1, iEven = 2, iOdd = 1, counters = {}, makeNewSeries, jOld = 1, newSeries = {}, NTarget = N[target]}, {maximum, outputSeries, partials, signChanges,bg,fs} = {MaxTerms, ShowSeries, ShowPartialSums, ShowSignChanges, Background, FontSize} /. {opts} /. Options[HarmonicRearrangement]; cw= 45 12./fs; (* makeNewSeries[] is a function that turns a list of sign changes, which is what counters will be, into a list of terms of the rearranged series. *) makeNewSeries[clist_] := Module[ {j = {1, 2}, i = 2}, If[newSeries != {}, newSeries, newSeries = 1 / (Flatten[ clist /. n_Integer :> Range[j[[i = 3-i]], (j[[i]] += 2 n) - 2, 2]] /. x_?EvenQ :> -x) ]]; While[j <= maximum, (* The next two While loops do the adding until an overshoot and subtracting to an undershoot *) While[currentSum <= NTarget, j++; currentSum += 1/iOdd; iOdd += 2]; AppendTo[counters, j - jOld]; jOld = j; While[currentSum > NTarget, j++; currentSum -= 1/iEven; iEven += 2]; AppendTo[counters, j - jOld]; jOld = j]; (* The following code generates output according to the settings of the output options. *) Identity & /@ (sss=Select[{ If[outputSeries, DisplayForm@MakeBoxes[#, StandardForm]&@ ( DisplayForm[RowBox[Append[ Prepend[ makeNewSeries[counters] /. Rational[x_,y_] :> RowBox[{Sign[x]/.{1->"+",-1->"-"},FractionBox[Abs[x],y]}], ToString[StringForm[ "The rearranged series for `` is: ",target]]],"\[CenterDot]\[CenterDot]\[CenterDot]"]]]) ], If[partials, StringForm[ "The partial sums just before the sign changes: `` \n", Append[ collapse[makeNewSeries[counters]//N], "\[Ellipsis]"]]], If[signChanges, StringForm[ "Numbers of terms that are positive, negative, positive,...: ``", Append[counters, "\[Ellipsis]"]]]}, # =!= Null &]); DisplayForm[StyleBox[GridBox[List /@ MakeBoxes /@ sss, System`ColumnWidths -> cw, ColumnAlignments->Left, GridFrame->True, RowLines-> 1, RowSpacings->{3,2}], Background -> bg, FontFamily -> Times, FontSize->fs]]]; PellSolve[d_Integer] := Module[ {p, q=1, q0=0, P=0, Q=1, a, a0=Floor[Sqrt[d]]}, a = a0; While[a != 2 a0, P = a*Q - P; Q = Quotient[d - P^2, Q]; a = Quotient[P + a0, Q]; {q, q0} = {a*q + q0, q}]; p = Sqrt[If[Mod[p = d*q0^2 - 1, 4] <= 1, p, p + 2]]; {{p, q0}, {d*q0, p}}] /; !IntegerQ[Sqrt[d]]; (* old: slow PellSolve[d_] := Module[{cf, con, n, first}, cf = ContinuedFraction[Sqrt[d]]; n = Length[cf[[-1]]]; con = Convergents[ FromContinuedFraction[Join[Flatten[cf], cf[[-1]]]]]; first = ({Numerator[#1], Denominator[#1]} & )[con[[n]]]; {first, {d, 1}*Reverse[first]}]; *) BirthDate[n_Integer] := Select[SumOfSquaresRepresentations[2, n], 1900 < #1[[2]] < 1990 && 1 <= Quotient[#1[[1]], 100] <= 12 && 1 <= Mod[#1[[1]], 100] <= 31 & ] /. {monthday_Integer, year_} :> StringForm["`` ``, ``", {"January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November", "December"}[[Quotient[monthday, 100]]], Mod[monthday, 100], year] /. {ans_} :> ans; BLSPrimeProof[n_?PrimeQ] := Module[{b = 2, bvals, P = oddprimes = Rest[First /@ FactorInteger[n - 1]]}, If[P == {}, While[PowerMod[b, (n - 1)/2, n] != n - 1, b++]; StringForm["n\[ThinSpace]=\[ThinSpace]`` is prime; the prime divisors of n-1 and their witnesses are ``", n, Transpose[{{2, b}}]], While[P != {}, If[PowerMod[b, (n - 1)/2, n] == n - 1, P1 = Select[P, PowerMod[b, (n - 1)/(2*#1), n] != n - 1 & ]; Scan[Function[k, bvals[k] = b], P1]; P = Complement[P, P1]]; b++]; StringForm[ "n\[ThinSpace]=\[ThinSpace]`` is prime; the odd prime divisors of n-1 and their witnesses are ``", n, Transpose[{oddprimes, bvals /@ oddprimes}]]]]; SFTF[x_] := StandardForm[TraditionalForm[x]]; FromCFSquared[a_] := Module[{x}, periodicsoln = x /. Solve[x == FromContinuedFraction[Append[a[[-1]], x]], x]; RootReduce[FullSimplify[Expand[ FromContinuedFraction[Append[Drop[a, -1], periodicsoln[[2]]]]]^2]]]; (*CFSquareRoot[n_] :=(nn= CFSpecial[0,1,n]; Append[nn[[1]], nn[[2]]]); *) LucasU[{P_Integer, Q_Integer}, mod_:Infinity][n_Integer] := First[LucasUPair[{P,Q}, mod][n]]; LucasV[{P_Integer, Q_Integer}, mod_:Infinity][n_Integer] := Mod[LucasUPair[{P,Q},mod][n].{-P,2},mod]; LucasU[{P_Integer, Q_Integer}, mod_:Infinity][n_List] := (LucasU[{P,Q}, mod][#])& /@ n; LucasV[{P_Integer, Q_Integer}, mod_:Infinity][n_List] := (LucasV[{P,Q}, mod][#])& /@ n; Lucas[{P_Integer, Q_Integer}, mod_:Infinity][n_List] := (Lucas[{P,Q}, mod][#])& /@ n; LucasUPair[{P_,Q_}, mod_Integer][n_Integer] := Fold[Mod[ If[#2==0,{#[[1]](2 #[[2]]-P #[[1]]), #[[2]]^2 - Q #[[1]]^2}, {#[[2]]^2 - Q #[[1]]^2,#[[2]](P#[[2]]-2 Q #[[1]] ) }], mod]&,{0,1},IntegerDigits[n,2]]; LucasUPair[{P_,Q_}][n_Integer] := LucasUPair[{P,Q}, Infinity][n]; LucasUPair[{P_,Q_}, Infinity][n_Integer] := Fold[ If[#2==0,{#[[1]](2 #[[2]]-P #[[1]]), #[[2]]^2 - Q #[[1]]^2}, {#[[2]]^2 - Q #[[1]]^2,#[[2]](P#[[2]]-2 Q #[[1]] ) }] &,{0,1},IntegerDigits[n,2]]; Lucas[{P_,Q_}, mod_:Infinity][n_Integer] := Mod[{{1,0},{-P,2}}. LucasUPair[{P,Q},mod][n], mod]; Clear[LegendrePhi,PrimePiMod4]; Attributes[LegendrePhi] = Listable; breakpoint = 20000; ppbk = PrimePi[breakpoint]; LegendrePhi[n_Integer, 0] := n; LegendrePhi[n_Integer,a_Integer]:= (If[Length[primeset] != ppbk, primeset = Prime[Range[ppbk]]]; If[n < primeset[[a+1]], 1]; LegendrePhi[n,a] = n - Apply[Plus, LegendrePhi[Quotient[n, Take[primeset,a]], Range[0, a-1]]]); LegendrePhi[_,n_Integer,_] := 0 /; n <= 0; LegendrePhi[1, n_Integer,1] := 1 + Quotient[n-1, 4]/; n >= 1; LegendrePhi[3, n_Integer,1] := Quotient[n+1, 4] /; n >= 1; LegendrePhi[1, n_Integer, a_Integer] := (If[Length[primeset] != ppbk, primeset = Prime[Range[ppbk]]]; LegendrePhi[1,n,a] = LegendrePhi[1,n,1] - Apply[Plus, MapIndexed[LegendrePhi[Mod[#, 4], Quotient[n,#],#2[[1]]]&, Take[primeset,{2, a}]]]); LegendrePhi[3, n_Integer, a_Integer] := LegendrePhi[n, a] - LegendrePhi[1,n,a]; Sieved\[CapitalPhi][x_, z_] := LegendrePhi[x, PrimePi[z]- If[PrimeQ[z], 1, 0]]; Sieved\[CapitalPhi][x_, z_] := (u = Log[z]/Log[x]; 1+primesininterval[x^u, x] + doublets[x,u]) /; z >= x^(1/3); Sieved\[CapitalPhi][x_, z_] := (u = Log[z]/Log[x]; 1+primesininterval[z,x] + doublets[x,u] + triplets[x,u])/; x^(1/3) > z >= x^(1/4); (* the following allow for speedups for cases between the fourth root and square root, based on the fast PrimePi *) doublets[x_, u_] := Module[{sum=0}, Do[sum += PrimePi[x/Prime[i]]-i+1,{i,PrimePi[x^u]+If[PrimeQ[x^u],0,1], PrimePi[x^0.5]}]; sum]; triplets[x_,u_] := Module[{sum=0,q},Do[q=Prime[i]; sum+=doublets[x/q,Log[q]/Log[x/q]],{i, PrimePi[x^u]+If[PrimeQ[x^u],0,1],PrimePi[x^(1/3)]}]; sum]; primesininterval[a_, b_] := PrimePi[b] - PrimePi[a] + If[PrimeQ[a], 1, 0]; Subscript[\[Pi], j_][x_] := PrimePiMod4[j][x]; (* PrimePiMod4[0][_?NumericQ] := 0; PrimePiMod4[2][x_?NumericQ] := If[x >= 2, 1, 0]; PrimePiMod4[1][x_?NumericQ] := 0 /; x < 5; PrimePiMod4[3][x_?NumericQ] := 0 /; x < 3; PrimePiMod4[1][x_?NumericQ] := 1 /; 5 <= x < 13; PrimePiMod4[1][n_?NumericQ ] := ( If[Length[primeset] != ppbk, primeset = Prime[Range[ppbk]]]; If[Length[primesmod] != ppbk, primesmod= FoldList[Plus, 0, Mod[Rest[primeset], 4] /. 3->0]]; primesmod[[PrimePi[n]]]) /; n < breakpoint; PrimePiMod4[1][x_] := PrimePiMod4[1][Floor[x]] /; !IntegerQ[x]&&NumericQ[x]; PrimePiMod4[1][x_Integer] := Module[{m, fcn, p, z = x^(1/3)}, m = PrimePi[z]; PrimePiMod4[1][z] + LegendrePhi[1,x,m] - 1 - Plus @@ ((p = Prime[#]; fcn = PrimePiMod4[Mod[p, 4]]; fcn[Quotient[x,p]] - fcn[p - If[# == m + 1, (p - z)/2, 1]]) & ) /@ Range[m + 1, PrimePi[Sqrt[x]]]]; *) PrimePiMod4[0][_?NumericQ] := 0; PrimePiMod4[2][x_?NumericQ] := If[x >= 2, 1, 0]; PrimePiMod4[1][x_?NumericQ] := 0 /; x < 5; PrimePiMod4[3][x_?NumericQ] := 0 /; x < 3; PrimePiMod4[1][x_?NumericQ] := 1 /; Inequality[5, LessEqual, x, Less, 13]; PrimePiMod4[1][n_?NumericQ] := (If[Length[primeset] != ppbk, primeset = Prime[Range[ppbk]]]; If[Length[primesmod] != ppbk, primesmod = FoldList[Plus, 0, Mod[Rest[primeset], 4] /. 3 -> 0]]; primesmod[[PrimePi[n]]]) /; n < breakpoint; PrimePiMod4[1][x_?NumericQ] := PrimePiMod4[1][Floor[x]] /; !IntegerQ[x]; PrimePiMod4[1][x_Integer] := Module[{m, c, pp, qq, p, z = x^(1/3), zz = PrimePi[x^(1/2)]}, m = PrimePi[z]; pp = NextPrime[z]; c = PrimePiMod4[1][pp]; qq = Quotient[x, pp]; c - If[Mod[pp, 4] == 1, 1 + PrimePiMod4[1][qq] - PrimePiMod4[1][(pp + z)/2.], PrimePi[qq] - PrimePiMod4[1][qq] - PrimePi[(pp + z)/2.] + PrimePiMod4[1][(pp + z)/2.; z]] + LegendrePhi[1, x, m] - zz + m + Plus @@ ((p = Prime[#]; qq = Quotient[x, p]; If[Mod[p, 4] == 1, c++; c - PrimePiMod4[1][qq], PrimePiMod4[1][qq] + #1 - Pi[qq] - c]) & ) /@ Range[m + 2, zz]]; PrimePiMod4[3][x_?NumericQ] := PrimePi[x] - PrimePiMod4[1][x] - 1; PrimePi[R_?NumericQ, GaussianIntegers->True] := If[R>Sqrt[2],4,0]+ 4 PrimePiMod4[3][R] +8 PrimePiMod4[1][R^2]; \[Pi][R_?NumericQ, GaussianIntegers->True] := If[R>Sqrt[2],4,0]+ 4 PrimePiMod4[3][R] + 8 PrimePiMod4[1][R^2]; PrimeNumberRace[a_Integer, b_Integer, step_Integer:1, opts___Rule] := With[{z = PrimePi[a]}, temp = Block[{$DisplayFunction = Identity}, ListPlot[data = FoldList[Plus, PrimePiMod4[3][a] - PrimePiMod4[1][a], Mod[Prime[Range[z + 1, PrimePi[b]]], 4] - 2]; data[[Range[1, Length[data], step]]],Evaluate[FilterOptions[ListPlot,opts]], PlotJoined -> True, Frame -> True, Axes -> None, FrameTicks -> {Automatic, Automatic, None, None}, GridLines -> {{}, {0}}]]; ft = FullOptions[temp, FrameTicks]; Show[temp, FrameTicks -> {ft[[1]] /. {(xx_)?(NumericQ[#1] && #1 >= 0 & ), (yy_)?NumericQ, zz__List} :> {xx, Prime[Round[Pi[a] + step*xx]], zz}, ft[[2]], {}, {}}]]; Options[LegendrePiFull] = {FontSize -> 9, Background->GrayLevel[0.9]}; LegendrePiFull[x_, opts___] := Block[ {$RecursionLimit = Infinity}, Ldata = {}; If[Length[primeset] != 1229, primeset = Eratosthenes[10000]]; {bg, fs} = {Background,FontSize} /. {opts} /. Options[LegendrePiFull]; DisplayForm[StyleBox[GridBox[GridBoxFix[ Prepend[Transpose[{rr = Reverse[Drop[FixedPointList[Floor[Sqrt[#1]] & , x], -2]]; col1 = Prepend[rr, 1], col2 = Prepend[Prepend[Drop[rr, -1], 1], 1], col3 = Prepend[FoldList[(AppendTo[Ldata, ll = LegendrePhi[#2, #1]]; ll - 1 + #1) & , 0, Drop[rr, -1]], 0], col4 = Append[Prepend[Ldata, 1], LegendrePhi[x, Last[col3]]], col6 = col3 + col4 - 1}], {"Square roots,\[MediumSpace]m", "\[LeftFloor]\!\(\@m\)\[RightFloor]", "\[Pi](\[LeftFloor]\!\(\@m\)\[RightFloor])", "\[Phi](m,\[MediumSpace]\[Pi](\[LeftFloor]\!\(\@m\)\[RightFloor]))", "\[Pi](m)\[MediumSpace]=\[MediumSpace]\[Pi](\[LeftFloor]\!\(\@m\)\[RightFloor])\[MediumSpace]+\[Phi](m,\[MediumSpace]\[Pi](\[LeftFloor]\!\(\@m\)\[RightFloor]))\[MediumSpace]-\[MediumSpace]1"}]], ColumnAlignments -> Right, ColumnLines -> 1, RowLines -> {2, 1}, GridFrame -> 1], Background -> bg, FontSize -> fs, FontFamily -> "Times"]]] /; x < 10^8 + 1; Lucas\[CapitalPhi][{P_, Q_}, {p_, e_}] := p^(e - 1)* (p - JacobiSymbol[P^2 - 4*Q, p]) /; GCD[Q, p] == 1; Lucas\[CapitalPhi][{P_, Q_}, n_?OddQ] := Times @@ (Lucas\[CapitalPhi][{P,Q}, #] &) /@ FactorInteger[n] /; GCD[Q, n] == 1; Lucas\[CapitalLambda][{P_, Q_}, (n_)?OddQ] := LCM @@ (Lucas\[CapitalPhi][{P, Q}, #] & ) /@ FactorInteger[n] /; GCD[Q, n] == 1; LucasRank[{_, _}, 1] = 1; LucasRank[{P_, Q_}, n_?OddQ] := Select[ Divisors[Lucas\[CapitalLambda][{P, Q}, n]], LucasU[{P,Q}, n][#1] == 0 &,1][[1]]; EulerPseudoprimeTest[a_, n_?OddQ] := Module[ {r = PowerMod[a,(n-1)/2,n]}, MemberQ[{1, n-1}, r] && r == JacobiSymbol[a,n]]; EulerPseudoprimeQ[a_,n_?OddQ] := !PrimeQ[n] && EulerPseudoprimeTest[a,n]; Options[LucasCertificate] = {Tree->False, FontSize->12, Background->GrayLevel[0.9]}; Attributes[LucasCertificateBasic] = Listable; LucasCertificateBasic[2] = 2; LucasCertificateBasic[n_Integer /; PrimeQ[n] && n > 2] := Module[{Q = 2, primes = Rest[First /@ FactorInteger[n + 1]]}, While[JacobiSymbol[1 - 4 Q, n] != -1 || LucasV[{1, Q}, n][(n + 1)/2] != 0 || Or @@ Map[LucasV[{1, Q}, n][(n + 1)/(2 #)] == 0 &, primes], Q++]; {n, Q, LucasCertificateBasic[primes]}]; LucasCertificate[p_Integer?PrimeQ, opts___] := ({fs,tr,rs}={FontSize,Tree,RowSpacings}/. {opts}/. Options[LucasCertificate]; If[tr, CertificateTreeForm[LucasCertificateBasic[p],0,0,opts, Sequence @@ Options[LucasCertificate]], LucasCertificateBasic[p]]); QuadraticPseudoprimeTest[{P_, Q_}, n_] := Module[{d = P^2 - 4*Q}, Which[1 < (g = GCD[n, 2*Q*d]) < n, False, g == n, PrimeQ[n], g == 1, \[CurlyEpsilon] = JacobiSymbol[d, n]; Lucas[{P,Q}, n][n - \[CurlyEpsilon]] == {0, Mod[2*Q^((1 - \[CurlyEpsilon])/2), n]}]]; QuadraticPseudoprimeQ[{P_, Q_}, n_Integer] := !PrimeQ[n] && QuadraticPseudoprimeTest[{P, Q}, n]; QuadraticPseudoprimeQ[par_, n_List] := QuadraticPseudoprimeQ[par, #] /@ n QuadraticPseudoprimeTest[False, _] := False; QuadraticPseudoprimeTest[True, _] := True; LucasPseudoprimeTest[False, _] := False; LucasPseudoprimeTest[True, _] := True; LucasVPseudoprimeTest[False, _] := False; LucasVPseudoprimeTest[True, _] := True; LucasPseudoprimeTest[{P_, Q_}, n_] := Module[{d = P^2 - 4*Q}, Which[1 < (g = GCD[n, 2*Q*d]) < n, False, g == n, PrimeQ[n], g == 1, 0 == LucasU[{P,Q}, n][n - JacobiSymbol[d, n]]]]; LucasPseudoprimeQ[pars_, n_] := !PrimeQ[n] && LucasPseudoprimeTest[pars, n]; LucasPseudoprimeQ[pars_, n_List] := LucasPseudoprimeTest[pars, #]& /@ n; LucasVPseudoprimeTest[{P_, Q_}, n_] := Module[{d = P^2 - 4*Q}, Which[1 < (g = GCD[n, 2*Q*d]) < n, False, g == n, PrimeQ[n], g == 1, eps= JacobiSymbol[d, n]; Mod[2 PowerMod[Q,(1-eps)/2,n],n] == LucasV[{P,Q}, n][n - eps]]]; LucasVPseudoprimeQ[pars_, n_] := !PrimeQ[n] && LucasVPseudoprimeTest[pars, n]; LucasVPseudoprimeQ[pars_, n_List] := LucasVPseudoprimeTest[pars, #]& /@ n; Options[UVSequence] = {FontSize -> 9, UseGridBox -> False,Background-> GrayLevel[0.9]}; GridBoxFix[x_] := Map[MakeBoxes, x, {2}]; UVSequence[{P_, Q_}, n_?OddQ, opts___] := Module[{e,m,pows,fs,bg,gQ}, {gQ, fs, bg}= {UseGridBox, FontSize, Background} /. {opts} /. Options[UVSequence]; e = IntegerExponent[m =n - JacobiSymbol[P^2-4Q, n], 2]; m = Quotient[m, 2^e]; uvseq = Transpose[(Drop[#1, -1] & ) /@ NestList[Mod[{#1[[1]]*#1[[2]], #1[[2]]^2 - 2*#1[[3]], #1[[3]]^2}, n] & , Append[Lucas[{P,Q}, n][m], PowerMod[Q, m, n]], e]]; If[gQ, pows = "m"*"2"^Range[0, e]; DisplayForm[StyleBox[GridBox[ GridBoxFix[ Prepend[Transpose[Prepend[Transpose[uvseq /. {n-1->-1,n-P->-P,n-2->-2}], {"U", "V"}]], Prepend[Append[Drop[pows, -1], StringForm["\!\(2\^``\) m = ``", e, m 2^e]], "Powers"]]], GridFrame -> 2, ColumnLines -> {2, 1}, ColumnAlignments -> Right, RowLines -> 1], FontSize -> fs, Background -> bg]],uvseq]]; MethodA[n_ /; OddQ[n] && !IntegerQ[Sqrt[n]], start_Integer:5] := Module[{d=start}, While[(j=JacobiSymbol[d,n]) == 1, d=-d - 2 Sign[d]]; If[j == -1, {1, (1-d)/4}, If[GCD[d,n] < n, False, PrimeQ[n]]]]; MethodAPositive[n_ /; OddQ[n] && !IntegerQ[Sqrt[n]], start_Integer:5] := Module[{d=start}, While[(j=JacobiSymbol[d,n]) == -1, d = -d - 2 Sign[d]]; If[j == 1, {1, (1-d)/4}, If[GCD[d, n] < n, False, PrimeQ[n]]]]; MethodAStar[n_ /;OddQ[n] && !IntegerQ[Sqrt[n]], new_:{5,5}] := Module[{d=5}, While[(j=JacobiSymbol[d,n]) == 1, d=-d-2Sign[d]]; If[j == -1, {1,(1-d)/4} /.{1,-1}->new, If[GCD[d,n] < n, False, PrimeQ[n]]] ]; MethodAStarPositive[n_ /;OddQ[n] && !IntegerQ[Sqrt[n]], new_:{5,5}] := Module[{d=5}, While[(j=JacobiSymbol[d,n]) == -1, d=-d-2Sign[d]]; If[j == 1, {1,(1-d)/4} /.{1,-1}->new, If[GCD[d,n] < n, False, PrimeQ[n]]] ] Options[PseudoprimeSearch] = {FermatBase -> 2, TraceInterval -> None, PrintImmediately -> False, NumberSought -> Infinity, LucasParameters -> Automatic (* MethodA, MethodAStar, StartAtMinusSeven *) , Type -> Fermat (* or: Strong, Lucas, Quadratic, StrongQuadratic, Fibonacci, QuadraticFibonacci *) }; PseudoprimeSearch[n0_:1, n1_Integer, opts___Rule] := Module[{ans={}, j=0,ti, lp, type, nmax, prQ, fb, hf = HoldForm[PseudoprimeSearch[n0, n1, opts]]}, If[n1 >= 2^31, Message[PseudoprimeSearch::toobig]; Return[hf]]; {ti, lp, type, nmax, prQ, fb} = {TraceInterval, LucasParameters, Type, NumberSought, PrintImmediately, FermatBase} /. {opts} /. Options[PseudoprimeSearch]; ans = {}; ti = ti /. Automatic -> Round[(n1-n0+1)/10]; met = type /. {Fermat -> (PseudoprimeQ[fb, #]&), Fibonacci -> (LucasPseudoprimeQ[{1,-1}, #] &), QuadraticFibonacci -> (QuadraticPseudoprimeQ[{1,-1}, #] &), Strong -> (StrongPseudoprimeQ[fb, #] &), Lucas -> Which[ListQ[lp], LucasPseudoprimeQ[lp, #]&, lp===MethodA, LucasPseudoprimeQ[MethodA[#],#]&, lp===MethodAStar || lp===Automatic, LucasPseudoprimeQ[MethodAStar[#],#]&, lp===StartAtMinusSeven, LucasPseudoprimeQ[MethodA[#,-7],#]&], LucasV -> Which[ListQ[lp], LucasVPseudoprimeQ[lp, #]&, lp===MethodA, LucasVPseudoprimeQ[MethodA[#],#]&, lp===MethodAStar || lp===Automatic, LucasVPseudoprimeQ[MethodAStar[#],#]&, lp===StartAtMinusSeven, LucasVPseudoprimeQ[MethodA[#,-7],#]&], Quadratic -> Which[ListQ[lp], QuadraticPseudoprimeQ[lp, #]&, lp===MethodA, QuadraticPseudoprimeQ[MethodA[#],#]&, lp===MethodAStar || lp===Automatic, QuadraticPseudoprimeQ[MethodAStar[#],#]&, lp===StartAtMinusSeven, QuadraticPseudoprimeQ[MethodA[#,-7],#]& ], StrongQuadratic -> Which[ ListQ[lp], StrongQuadraticPseudoprimeQ[lp, #]&, lp===MethodA, StrongQuadraticPseudoprimeQ[MethodA[#],#]&, lp===MethodAStar || lp===Automatic, StrongQuadraticPseudoprimeQ[MethodAStar[#],#]&, lp===StartAtMinusSeven, StrongQuadraticPseudoprimeQ[MethodA[#,-7],#]& ]}; Do[count = i; (* in case of interrupt *) If[IntegerQ[ti] && Mod[i - 1, ti] == 0, ans = Flatten[ans]; Print[StringForm["Tracing: ``", i - 1]]]; If[met[i], ans = {ans,i}; If[prQ, Print[i]]; j++; If[j == nmax, Break[]]], {i, n0 + 1 - Mod[n0,2], n1}]; Flatten[ans]]; StrongQuadraticPseudoprimeTest[{P_, Q_}, n_?OddQ] := Module[{d = P^2 - 4*Q, g, u, v}, g = GCD[Q*d, n]; Which[1 < g < n, False, g == n, PrimeQ[n], True, {u, v} = UVSequence[{P, Q}, n]; v[[-1]] == Mod[2*Q^((1 - JacobiSymbol[d, n])/2), n] && (u[[1]] == 0 || MemberQ[v, 0])]]; StrongQuadraticPseudoprimeQ[{P_, Q_}, n_?OddQ] := !PrimeQ[n] && StrongQuadraticPseudoprimeTest[{P, Q}, n]; StrongQuadraticPseudoprimeQ[True, _] := True; StrongQuadraticPseudoprimeQ[False, _] := False; QuadraticCongruenceSolve[{a_Integer, b_Integer, c_Integer}, 2] := LinearCongruenceSolve[{a + b, -c}, 2]; QuadraticCongruenceSolve[{a_Integer, b_Integer, c_Integer}, p_?PrimeQ] := LinearCongruenceSolve[{b, -c}, 2] /; p > 2 && Mod[a, p] == 0; QuadraticCongruenceSolve[{a_Integer, b_Integer, c_Integer}, p_?PrimeQ] := (ainv = PowerMod[a, -1, p]; fourinv = PowerMod[4, -1, p]; Sort[Mod[SqrtModAll[ fourinv ainv^2 b^2 - ainv c, p] - (p + 1)/2*ainv*b, p]]) /; p > 2 && Mod[a, p] != 0; QuadraticCongruenceSolve[{a_Integer, b_Integer, c_Integer}, {p_, 1}] := QuadraticCongruenceSolve[{a, b, c}, p]; QuadraticCongruenceSolve[{a_Integer, b_Integer, c_Integer}, {p_, e_ /; e > 1}] := Module[{prev = QuadraticCongruenceSolve[{a, b, c}, {p, e - 1}]}, If[prev == {}, {}, Sort[Flatten[(#1 + LinearCongruenceSolve[{2 a # + b, -(a #^2 + b # + c)/p^(e - 1)}, p] p^(e - 1) & ) /@ prev]]]]; QuadraticCongruenceSolve[{a_Integer, b_Integer, c_Integer}, n_] := Module[{fi = FactorInteger[n]}, ppow = (QuadraticCongruenceSolve[{a, b, c}, #] & ) /@ fi; If[MemberQ[ppow, {}], {}, Sort[(ChineseRemainder[#, Apply[Power, fi, {1}]] & ) /@ Distribute[ppow, List]]]] /; !PrimeQ[n]; LucasParameter[n_?OddQ] := Module[{P = 3}, While[GCD[P, n] != 1 || JacobiSymbol[P^2 - 4, n] == 1, P++]; P]; PocklingtonPrimeWitness[n_?OddQ, q_Integer] := PocklingtonPrimeWitness[n, {q}]; PocklingtonPrimeWitness[n_?OddQ, primes_] := Module[{b = 2, nn = n - 1}, Scan[nn /= #^IntegerExponent[nn, #]&, primes]; If[n - 1 < nn Sqrt[n], Return[False]]; While[PowerMod[b, n - 1, n] != 1 || Or @@ (GCD[n, PowerMod[b, (n - 1)/#1, n] - 1] != 1 & ) /@ primes, b++]; b]; PocklingtonPrimeWitness[n_?OddQ] := PocklingtonPrimeWitness[n, First /@ FactorInteger[n - 1]]; primeProd = Apply[Times, Prime[Range[3,100]]]; QuickPrimeQ[n_] := GCD[primeProd,n]==1 && PowerMod[2,n-1,n]==1; Options[CertifiedPrime] = {ShowCertificate -> False}; CertifiedPrime[d_Integer, {}, opts___Rule] := CertifiedPrime[d, opts]; CertifiedPrime[d_Integer, opts___Rule] := If[ShowCertificate /. {opts} /. Options[CertifiedPrime],Identity, First][ PrimeAndCertificate[d]]; PrimeAndCertificate[d_Integer] := Module[{p, r, q, q1 = PrimeAndCertificate[Quotient[d + 1, 2] + 1]}, q = q1[[1]]; r = Random[Integer, {10.^(d - 1), 10.^d}/(q)]; r = r-Mod[r,6] + If[EvenQ[r], 0, Mod[q, 3] /. 1 -> 4]; While[r += 6; p = 1 + q*r; !QuickPrimeQ[p] || GCD[p, PowerMod[2, r, p] - 1] != 1]; Prepend[q1, p]] /; d >= 8; PrimeAndCertificate[d_ /; d < 8] := {RandomPrimePrimeQ[d]} CertifiedPrime[d_Integer, {a_,m_}, opts___Rule] := If[ShowCertificate /. {opts} /. Options[CertifiedPrime],Identity, First][ PrimeAndCertificate[d,{a,m}]]; PrimeAndCertificate[d_Integer, {a_, m_}] := Module[{res=ChineseRemainder[{#, a}, {6,m}] & /@ {1,-1}, lcm = LCM[6,m], p, r, q, q1 = PrimeAndCertificate[Quotient[d + 1, 2] + 1]}, q = q1[[1]]; r = Random[Integer, {10.^(d - 1), 10.^d}/(q)]; r = r- Mod[r,lcm]+ (res[[If[EvenQ[r],1,2]]]-1)*PowerMod[q,-1,lcm]; While[r += lcm; p = 1 + q*r; !QuickPrimeQ[p] || GCD[p, PowerMod[2, r, p] - 1] != 1]; Prepend[q1, p]] /; d >= 8; PrimeAndCertificate[d_ /; d < 8, {a_,m_}] := {RandomPrimePrimeQ[d, {a,m}]}; GaussianNorm[a_?GaussianIntegerQ] := Abs[a]^2; GaussianIntegerQ[a_Integer] := True; GaussianIntegerQ[Complex[a_Integer, b_Integer]] := True; GaussianIntegerQ[a_] := False; GaussianPrimeQ[p_?GaussianIntegerQ] := PrimeQ[p, GaussianIntegers -> True]; (*GaussianMod[a_?GaussianIntegerQ, p_?GaussianIntegerQ] := Module[{pp = GaussianNorm[p]}, Mod[Re[a] - Im[a]*Re[p]*PowerMod[Im[p], -1, pp], pp]] /; !IntegerQ[p] && !MatchQ[p, Complex[0, _Integer]] && GCD[Re[p], Im[p]] == 1; GaussianMod[a_?GaussianIntegerQ, p_?GaussianIntegerQ] := Mod[a, p]; *) GaussianMod[u_?GaussianIntegerQ, z_?GaussianIntegerQ] := Module[{d,q, a = Re[z], b=Im[z], x=Re[u], y=Im[u]}, {d, e} = ExtendedGCD[a,b]; q = Quotient[y,d]; Mod[x - a q e[[2]] + b q e[[1]], (a^2+b^2)/d] +I (y - q d)]; QuadraticMod[a_Integer, n_Integer] := Mod[a,n]; QuadraticMod[a_?QuadraticIrrationalQ, n_Integer] := Module[{d = Cases[{RootReduce[a]}, Sqrt[_], Infinity]}, If[d!={}, d=d[[1]] ]; Expand[ a /. d -> x /. Rational[e_,f_] :> Mod[e*PowerMod[f,-1,n],n]]/.i_Integer :> Mod[i, n] /.x -> d]; QuadraticPowerMod[a_?QuadraticIrrationalQ, n_Integer, m_Integer] := Module[{aa = RootReduce[a]}, Fold[QuadraticMod[Expand[QuadraticMod[Expand[#1*#1], m]* If[#2 == 1, aa, 1]], m]&, QuadraticMod[aa,m], Rest[IntegerDigits[n, 2]]]]; GaussianPrimesInExpandedFirstOctant[rad_, k_] := Select[Flatten[Table[{a, b}, {a, -Ceiling[k], rad}, {b, -Ceiling[k], Min[a + Floor[k*Sqrt[2]], Sqrt[rad^2 - a^2]]}], 1], GaussianPrimeQ[#1[[1]] + #1[[2]]*I] & ]; vectors[k_] := vectors[k] = Select[Flatten[Table[{i, j}, {i, -Floor[k], Floor[k]}, {j, -Floor[k], Floor[k]}], 1], #1 != {0, 0} && #1[[1]]^2 + #1[[2]]^2 <= k^2 && Mod[#1[[1]] - #1[[2]], 2] == 0 & ]; neighbors[pt_, k_] := Select[(pt + #1 & ) /@ vectors[k], GaussianPrimeQ[#1[[1]] + #1[[2]]*I] & ]; neighbors[{1, 1}, k_] := Select[Join[({2, 1} + #1 & ) /@ vectors[k], {{2, 1}, {1, 2}, {-1, 1}, {1, -1}}], (#1[[1]] - 1)^2 + (#1[[2]] - 1)^2 <= k && GaussianPrimeQ[#1[[1]] + #1[[2]]*I] & ]; symmetrize[{a_Integer, b_Integer}] := {{a, b}, {-a, b}, {a, -b}, {-a, -b}, {b, a}, {-b, a}, {b, -a}, {-b, -a}}; symmetrize[{p_, q_}] := Transpose[{symmetrize[p], symmetrize[q]}]; symmetrize[pts_] := Flatten[symmetrize /@ pts, 1] /; Depth[pts] == 3; symmetrize[{}] = {}; edges[Gp_, d_] := (Sort[{Gp, #1}] & ) /@ neighbors[Gp, d]; fix[pairs_] := Scan[Which[mark[#1[[1]]], mark[#1[[2]]] = True, mark[#1[[2]]], mark[#1[[1]]] = True] &, pairs]; Options[ReachableGaussianPrimes] = {ReachablesOnly -> False, UseLines -> True, ReachableStyle -> PointSize[0.01], PrintTheCounts -> False, UnreachableStyle -> PointSize[0.01], DiskBackground -> RGBColor[1, 1, 1]}; ReachableGaussianPrimes[d_, rad_, opts___] := ({bg, reachsty, unreachsty, lineQ, roQ, prQ} = {DiskBackground, ReachableStyle, UnreachableStyle, UseLines, ReachablesOnly, PrintTheCounts} /. {opts} /. Options[ReachableGaussianPrimes]; If[ !lineQ, roQ = True]; If[ !ListQ[unreachsty], unreachsty = {unreachsty}]; If[ !ListQ[reachsty], reachsty = {reachsty}]; roads = Union[Flatten[(edges[#, d]&) /@ GaussianPrimesInExpandedFirstOctant[rad, d], 1]]; Clear[mark]; mark[_] := False; mark[{1, 1}] = True; unmarkedPairs = roads; FixedPoint[(fix[unmarkedPairs]; Length[unmarkedPairs = Select[unmarkedPairs,!mark[#[[1]]] || !mark[#[[2]]]&]])&, Length[unmarkedPairs]]; firstOctant = Select[roads, mark[#[[1]]]&]; If[!lineQ, firstOctant = Union[Cases[firstOctant, {x_Integer, y_} /; 0 <= y <= x, Infinity]]]; roadNetwork = Union[Flatten[symmetrize /@ firstOctant, 1]]; unreachablePrimes = If[roQ, {}, Union[symmetrize[Union[Select[Flatten[roads, 1], (!mark[#1]) && #1[[1]]^2 + #1[[2]]^2 <= rad^2 & ]]]]]; If[!lineQ && prQ, Print[StringForm[ "There are `` primes in the first octant and `` \ in the entire component.", Length[firstOctant], Length[roadNetwork]]]]; Show[Graphics[{{bg, Disk[{0, 0}, rad]}, Circle[{0, 0}, rad], Join[unreachsty, Point /@ unreachablePrimes], {Thickness[0.0004], Join[reachsty, If[lineQ, Line, Point] /@ roadNetwork]}}], FilterOptions[Graphics, opts], AspectRatio -> Automatic, PlotRange -> All, Frame -> False]); QuarticResidueSymbol[\[Alpha]_?GaussianIntegerQ, p_] := (GaussianMod[PowerMod[\[Alpha], (GaussianNorm[p] - 1)/4, p], p] /. {If[IntegerQ[p], p - 1, GaussianNorm[p] - 1] -> -1, GaussianMod[I, p] -> I, GaussianMod[-I, p] -> -I}) /; GaussianIntegerQ[p] && GaussianNorm[p] != 2 Options[CFFactor] = {TraceIntervalFactor -> None, FactorBaseSize->Automatic, ShowData -> False, TargetMargin -> 30}; testnumber = Apply[Times,Prime[Range[PrimePi[10000]]]]; CFFactor[n_Integer, opts___] := Module[{ tif, k = 1, X, Y, Z, fac=False, i=1, target, g, a0, a, p0 = 1, p1,fbs, j=2}, {sdQ, tif, tm,fbs} = {ShowData, TraceIntervalFactor, TargetMargin, FactorBaseSize} /. {opts} /. Options[CFFactor]; If[1 < (g = GCD[n, testnumber]) < n, Return[{g, n/g}]]; While[!IntegerQ[r = n^(1/Prime[j])] && j <= 3, j++]; If[IntegerQ[r], Return[{r, n/r}]]; If[PrimeQ[n], Return[{n}]]; X = p1 = a0 = Floor[Sqrt[n]]; Y = n - a0^2; factorbase = fb[n, fbs /. Automatic->Floor[10*1.1^Log[10, n]]]; If[IntegerQ[factorbase], If[sdQ, Print["A prime in the factor base worked."]]; Return[Sort[{factorbase, n/factorbase}]]]; prod = Times @@ factorbase; pvals = rows = yvals = {}; ycount = count = 0; target = Length[factorbase] + tm; If[sdQ, Print[StringForm[ "Beginning the search for `` factorable Y-values...", target]]]; While[count < target, ycount++; k = 1 - k; a = Quotient[X + a0, Y]; {p1, p0} = {Mod[a*p1 + p0, n], p1}; X = a*Y - X; Y = Mod[(-1)^k*p1^2, n]; If[IntegerQ[Z == Sqrt[Y]] && 1 < (g = GCD[n, p1 - Z]) < n, If[sdQ, Print[StringForm["`` y-values have been examined and the p-y pair `` has a square y that succeeds", ycount, {p1, Y}]]]; Return[Sort[{g, n/g}]]]; If[PowerMod[prod, 2^Floor[Log[2., Y]], Y] == 0, count++; pvals = {pvals, p1}; yvals = {yvals, Y}; If[NumericQ[tif] && Mod[count, Max[1, Round[target*tif]]] == 0, Print[StringForm["``% done", NumberForm[Round[100*count/target], 2]]]]; rows = {rows, k, exponentsBinary[Y, factorbase]}]]; If[sdQ, Print[StringForm[ "It took `` Y-values to find the `` good values needed.", ycount,target]]]; rows = Transpose[Partition[Flatten[rows], Length[factorbase] + 1]]; {pvals,yvals} = Flatten /@ {pvals,yvals}; mat = NullSpace[rows, Modulus -> 2]; lm = Length[mat]; coeffs = (Flatten[Position[#1, 1]] & ) /@ mat; If[sdQ, badcount = goodcount = 0; Do[XX = Mod[Times @@ pvals[[coeffs[[i]]]], n]; YY = Mod[Sqrt[Times @@ yvals[[coeffs[[i]]]]], n]; If[XX != YY && XX != n - YY, g = GCD[n, XX - YY]; If[1 < g < n, goodcount++], badcount++], {i, Length[coeffs]}]; If[badcount==Length[coeffs], Print[StringForm["All `` of the square combinations had X\[ThinSpace]=\[ThinSpace]±Y.", badcount]], Print[StringForm["`` of the `` square combinations had X\[ThinSpace]=\[ThinSpace]±Y;\n`` of the remaining `` failed to split n.", badcount, Length[coeffs], Length[coeffs] - badcount - goodcount, Length[coeffs] - badcount]]]]; While[!fac && i <= Length[coeffs], XX = Mod[Times @@ pvals[[coeffs[[i]]]], n]; YY = Mod[Sqrt[Times @@ yvals[[coeffs[[i]]]]], n]; If[XX != YY && XX != -YY, g = GCD[n, XX - YY]; If[1 < g < n, fac = True]]; i++]; If[fac, Sort[{g, n/g}], $Failed]]; exponentsBinary[y_, helpers_] := Module[{ p = DeleteCases[{#, IntegerExponent[y, #]}& /@ helpers, {_, 0}], prs, ans = helpers}, prs = First /@ p; MapIndexed[(ans[[#2[[1]]]] = If[!MemberQ[prs, #1], 0, Cases[p, {#1, e_} :> Mod[e, 2]][[1]]])&, helpers]]; fb[n_Integer, len_Integer] := Module[{smallQ = False, fac, j}, If[EvenQ[n],Return[2]]; u = Select[Prime[Range[2*len]], # == 2 || (j = JacobiSymbol[n, #1]; If[j == 0, smallQ = True; fac = GCD[#1, n]]; j == 1)&]; If[smallQ, fac, u]]; PerrinPseudoprimeQ[n_Integer] := (!PrimeQ[n] && LinearRecurrence[{1,1,0},{3,0,2},n,n]==0); Options[AdmissibleQ] = {GaussianIntegers->False}; AdmissibleQ[c_, opts___] := Module[{primes,gi,l=Length[c],i,p=0}, gi=GaussianIntegers/.{opts}/.Options[AdmissibleQ]; If[!FreeQ[c,Complex],gi=True]; If[gi,primes= Select[Flatten[ Table[a+b I,{a,0,Floor[Sqrt[l]]},{b,0,Floor[Sqrt[l-a^2]]}]], GaussianPrimeQ]; And@@(Length[Union[Mod[c,#]]] != GaussianNorm[#]&) /@ primes, i=1; ans=True; While[(p=Prime[i]) <= l && ans, i++; ans = (Length[Union[Mod[c,p]]] != p)]; ans]]; ZetaGaussian[s_] := (Zeta[s]*(Zeta[s, 1/4] - Zeta[s, 3/4]))/4^s; distancesq[a_, b_] := Abs[b - a]^2; maxnorm[set_] := Max[Apply[distancesq, DeleteCases[Distribute[{set, set}, List], {x_, x_}], {1}]]; NextNorm[L_] := Module[{t = L + 1,tt}, While[!(IntegerQ[tt = Sqrt[t]] && Mod[tt, 4] == 3 && PrimeQ[tt] || Mod[t, 4] == 1 && PrimeQ[t]), t++]; t]; Gprime[n_] := FactorInteger[n, GaussianIntegers -> True][[-1,1]]; Options[CriticalPrimeInterval] = {GaussianIntegers -> False}; CriticalPrimeInterval[A_, opts___] := ( imagQ = GaussianIntegers /. {opts} /. Options[CriticalPrimeInterval]; If[!FreeQ[A, Complex], imagQ = True]; If[imagQ, Join[{1 + I}, Select[Range[Sqrt[maxnorm[A]]], Mod[#, 4] == 3 && PrimeQ[#]&], Gprime /@ Select[Range[maxnorm[A]], Mod[#,4]==1 && PrimeQ[#] &]], Select[Range[Sqrt[maxnorm[A]]], PrimeQ]]); PrimeZetaGaussian[1, ___] := Infinity; PrimeZeta[1, ___] := Infinity; Options[PrimeZeta] = {GaussianIntegers -> False}; PrimeZeta[s_, prec_Integer:6, opts___Rule] := (imagQ = GaussianIntegers /. {opts} /. Options[HardyLittlewood]; If[ !FreeQ[A, Complex], imagQ = True]; Plus @@ Table[MoebiusMu[k]* Log[If[imagQ, ZetaGaussian, Zeta][k*s]]/k, {k, IndexBoundPrimeZeta[s, Ceiling[Log[10, 1/2*10^(prec + 1)*2^s]]]}]); Options[HardyLittlewood] := {GaussianIntegers -> False}; HardyLittlewood[A_List, prec_Integer:6, opts___Rule] := Module[{n = Length[A], primess, imagQ}, imagQ = GaussianIntegers /. {opts} /. Options[HardyLittlewood]; If[!FreeQ[A, Complex], imagQ = True]; primes = CriticalPrimeInterval[A, opts]; max = If[imagQ, Max[GaussianNorm /@ primes], Max[primes]]; If[imagQ, primes = Join[{1 + I}, Select[Range[Sqrt[max]], Mod[#1, 4] == 3 && PrimeQ[#1] & ], Flatten[Gprime /@ Select[Range[max], Mod[#1, 4] == 1 && PrimeQ[#1] & ] /. Complex[a_, b_] :> {a+b*I, b+a*I}]]]; If[ !AdmissibleQ[A, opts], Return[0]]; If[imagQ, primes3mod4 = Select[primes, IntegerQ]; prime2 = If[MemberQ[primes, 1 + I], {1 + I}, {}]; primesimag = Cases[DeleteCases[primes, 1 + I], _Complex]; factor2 = If[prime2 == {1 + I}, 2^(n - 1), 1]; factor3mod4 = Times @@ ((#^2 - Length[Union[Mod[A, #]]])/#^2 * (#^2/(#^2 - 1))^n & ) /@ primes3mod4; factorimag = Times @@ ((gn = GaussianNorm[#1]; (gn - Length[Union[Mod[A, #1]]])/gn*(gn/(gn - 1))^n) & ) /@ primesimag; rat = Times @@ {factor2, factor3mod4, factorimag}; HLFunction[rat, n, NextNorm[max], GaussianIntegers -> True], rat = Times @@ ((# - Length[Union[Mod[A, #]]])/# * (#/(# - 1))^n & ) /@ primes; HLFunction[rat, n, NextPrime[max]]]]; Options[HLFunction] = {GaussianIntegers -> False}; Format[HoldPattern[HLFunction[rat_, n_Integer, k_, opts___]]] := (imagQ = GaussianIntegers /. {opts} /. Options[HLFunction]; DisplayForm[RowBox[Map[MakeBoxes[#]&, {FactorForm[FactorInteger[Numerator[rat]]]/ If[IntegerQ[rat], 1, FactorForm[FactorInteger[Denominator[rat]]]], EulerProduct[n, k, GaussianIntegers -> imagQ]}, {1}]]]); HLFunction /: N[HLFunction[rat_, n_, k_,opts___Rule], prec_:6] := (N[rat* EulerProduct[n, k, opts], prec]) /. -1. -> 1; HLFunction /: Round[HLFunction[rat_, n_, k_,opts___Rule], prec_:6] := Round[(N[rat* EulerProduct[n, k, opts], prec])]; IntegrateHold /: N[IntegrateHold[f_, {x_, a_, b_}], prec_:6] := Module[{temp = Integrate[f /. x -> ToExpression[x], ToExpression[x]]}, temp /. x -> ToExpression[b]]; Clear[ConstellationEstimate, IntegralFunction,ConstellationEstimateFunction]; Options[ConstellationEstimate] = Options[ConstellationEstimateFunction] = {GaussianIntegers -> False, EvaluateIntegral -> True, FormatOutputNicely->True}; fixConEst[expr_] := expr /. xx_Integer /; xx != -1 :> ExponentForm[xx]; Options[IntegralFunction] = {GuassianIntegers -> False}; IntegralFunction[A_, t_, opts___] := IntegralFunction[A, t,opts] =( imagQ = GaussianIntegers /. {opts} /. Options[IntegralFuntion]; If[!FreeQ[A, Complex], imagQ = True]; If[imagQ, Integrate[t/Log[t]^Length[A], t] /. ExpIntegralEi[2*Log[t]] -> LogIntegral[t^2], Integrate[1/Log[t]^Length[A], t]]); ConstellationEstimateFunction[{a_}, opts___Rule] := Module[{t}, {intQ, imagQ, intformatQ} = {EvaluateIntegral, GaussianIntegers, FormatOutputNicely} /. {opts} /. Options[ConstellationEstimateFunction]; functemp = If[intformatQ,Function[z, SFTF[If[!FreeQ[z, _Real],Rationalize, Identity][z] /.eform]], Identity]; If[!FreeQ[{a}, Complex], imagQ = True]; If[imagQ, Function[x, Evaluate[4/Pi^3* (If[intQ, IntegralFunction[{a},x,GaussianIntegers -> True], IntegrateHold[ToString[t]/Log[ToString[t]], {"t", 0, x}]])]], Function[x, Evaluate[ (If[intQ, IntegralFunction[{a},x,GaussianIntegers->False], IntegrateHold[1/Logg[ToString["t"]], {"t", 0, x}]])]]]]; ConstellationEstimate[A_, opts___Rule][x_] := Module[{t}, Print[hell]; {intQ, imagQ, intformatQ} = {EvaluateIntegral, GaussianIntegers, FormatOutputNicely} /. {opts} /. Options[ConstellationEstimate]; If[!FreeQ[A, Complex], imagQ = True]; functemp = If[intformatQ && !NumericQ[x],Function[z, SFTF[If[!FreeQ[z, _Real], Rationalize, Identity][z] /.eform]], Identity]; functemp[If[imagQ && !NumericQ[x], fixspecialcase, Identity][ConstellationEstimateFunction[A,opts][x]]]]; ConstellationEstimateFunction[A_, opts___Rule] := Module[{t}, {intQ, imagQ, intformatQ} = {EvaluateIntegral, GaussianIntegers, FormatOutputNicely} /. {opts} /. Options[ConstellationEstimateFunction]; functemp = If[intformatQ, $Post = Function[z, z /. {1. -> 1, -1. -> -1}]; Function[z, SFTF[If[!FreeQ[z, _Real],Rationalize, Identity][z] /.eform]], Clear[$Post];Identity]; If[!FreeQ[A, Complex], imagQ = True]; If[imagQ, Function[x, Evaluate[HardyLittlewood[A, GaussianIntegers -> True]* Pi/4 * (2/Pi)^Length[A]* (If[intQ, IntegralFunction[A,x,GaussianIntegers -> True], IntegrateHold["tt"/Log["tt"]^Length[A], {"t", 0, x}]])]], Function[x, Evaluate[HardyLittlewood[A]* (If[intQ, IntegralFunction[A,x,GaussianIntegers->False], IntegrateHold[1/Log["t"]^Length[A], {"t", 0, x}]])]]]] /; Length[A] > 1; Format[HoldPattern[IntegrateHold[f_, {x_, a_, b_}]]] := DisplayForm[RowBox[{SubsuperscriptBox["\[Integral]", (MakeBoxes[#] &)[a], (MakeBoxes[#] &)[b]], RowBox[{(MakeBoxes[#] &)[f], RowBox[{"\[DifferentialD]", ToString[x]}]}]}]]; Options[EulerProduct] = {GaussianIntegers -> False}; EulerProduct /: N[EulerProduct[len_, L_, opts___], prec_:6] := (imagQ = GaussianIntegers /. {opts} /. Options[EulerProduct]; If[imagQ, EulerProductNGaussian, EulerProductNRational][Round[len], Round[L], Round[prec]]); Format[EulerProduct[n_, k_, opts___]] := (imagQ = GaussianIntegers /. {opts} /. Options[EulerProduct]; base = If[imagQ, "N(p)", "p"]; DisplayForm[RowBox[{UnderscriptBox["\[Product]", RowBox[{base, "\[GreaterEqual]\[ThinSpace]", ToString[k]}]], FractionBox[RowBox[{RowBox[{"(", RowBox[{base, "-", ToString[n]}], ")"}], " ", If[n - 1 > 1, SuperscriptBox[base, ToString[n - 1]], base]}], SuperscriptBox[RowBox[{"(", RowBox[{base, "-", "1"}], ")"}], ToString[n]]]}]]); EulerProductNRational[len_, L_, prec_] := Module[{term1}, sum1 = 0; cp = Drop[Prime[Range[PrimePi[L]]], -1]; bd = IndexBoundHLProductRational[prec, len, Max[cp]]; Do[smallprimecorrex = -Plus @@ (1/cp^j); factor = (len - len^j)/j; term1 = Re[N[factor*(PrimeZeta[j, If[j == 2, Max[prec, 17], Max[prec, prec + 1 + Ceiling[Log[10, Abs[smallprimecorrex/(term1/factor)]]] ]]] + smallprimecorrex), Max[prec, 17]]]; sum1 += term1; If[traceQ && Mod[j, traceI] == 0, Print[{j, term1, sum1, Abs[term1/sum1]}]], {j, 2, bd}]; N[E^sum1, prec]]; EulerProductNGaussian[len_, L_, prec_] := Module[{term1}, cp = Select[Join[{1 + I}, Select[Range[Sqrt[L]], Mod[#1, 4] == 3 && PrimeQ[#1] &], Gprime /@ Select[Range[L], Mod[#1, 4] == 1 && PrimeQ[#1] & ]], GaussianNorm[#1] < L & ]; cp1 = GaussianNorm /@ Cases[DeleteCases[cp, 1 + I], _Complex]; cp3 = GaussianNorm /@ Select[cp, IntegerQ]; sum1 = 0; bd = IndexBoundHLProductGaussian[prec, len, Max[GaussianNorm /@ cp]]; If[traceQ, Print["The number of terms in the main sum is ", bd]]; Do[If[j > 2 && traceQ && Mod[j, traceI] == 0, Print[{j, term1, sum1, {Precision[term1], Precision[sum1]}, Abs[term1/sum1]}]]; factor = (len - len^j)/j; smallprimecorrex = -(Plus @@ (2/cp1^j) + Plus @@ (1/cp3^j) + If[MemberQ[cp, 1 + I], 1/2^j, 0]); term1 = Re[N[factor*(PrimeZeta[j, If[j == 2, Max[prec, 17], Max[prec, prec + 1 + Ceiling[Log[10, Abs[smallprimecorrex/(term1/factor)]]] ]], GaussianIntegers -> True] + smallprimecorrex), Max[17, prec]]]; sum1 += term1, {j, 2, bd}]; N[E^sum1, prec]] eform = {x_Integer /; Abs[x] > 10009 :> fixConEst[x], x_Real /; N[x - Round[x]] == 0 :> Round[x]}; fixspecialcase[z_] := ToString[4/\[Pi]^3, StandardForm]* \[Pi]^3 /4 z; IndexBoundHLProductGaussian[prec_, len_, L_] := (M = NextNorm[L]; Floor[ProductLog[(4*10^(prec + 1)*M* (M + 1)*Log[M/len])/(M - len)]/Log[M/len]]); IndexBoundPrimeZeta[s_, acc_] := Floor[ProductLog[(2^(1 + s)*10^(1 + acc)*s*(1 + s)* Log[2])/((1 - 2^s)*(1 - s))]/(s*Log[2])]; IndexBoundHLProductRational[prec_, len_, L_] := If[L == 2, Max[Ceiling[2*ProductLog[10^(prec/2)*Sqrt[len/(4 - len)]* Log[4/len]]/Log[4/len]], Floor[ProductLog[3/2*10^(prec+1)*Log[3/2]]/Log[3/2]]], M = NextPrime[L]; Ceiling[2*ProductLog[1/2*10^(prec/2)* Sqrt[len*(M - 1)/(M - 1 - len)]*Log[(M - 1)/len]]/ Log[(M - 1)/len]]]; (* from the 3.0 CF package *) Module[{iter, cf}, If[!v4Q, Clear[ContinuedFraction]; ContinuedFraction[x_Rational] := Identity[Floor[ cf[x] ] ]; ContinuedFraction[x_Real?MachineNumberQ, n_Integer?Positive] := Identity[ Floor[ cf[ SetPrecision[x, $MachinePrecision], n-1 ] ] ]; ContinuedFraction[x_Real, n_Integer?Positive] := Identity[ Floor[ cf[ x, n-1 ] ] ]; ContinuedFraction[x_Rational, n_Integer?Positive] := Identity[ Floor[ cf[ x, n-1 ] ] ]; (* We need to dynamically adjust the amount of precision to obtain the required order. This is done using an increment which varies according to the difference between the precision and the number of terms obtained (plus some fixed increment). *) ContinuedFraction[x_?NumericQ, n_Integer?Positive] := Module[{i, len, nm1, nx, prec, precinc, res}, len = 0; out = {}; precinc = 10; prec = n + precinc; For[ i=0; nm1 = n-1, (* Start *) len < n && i++ < n, (* Test *) prec += n + precinc - len, (* Increment *) nx = N[x, prec]; res = cf[nx, nm1]; len = Length[res]; ]; If[len < n, Message[ContinuedFraction::incomp, n] ]; Identity[ Floor[res] ] ]; iter[num_]:= If[ #==0, {}, (* Stopping criterion; something not NumberQ. *) Divide[1,#] ]& @ FractionalPart[num]; (* Remove terminating empty lists. *) cf[num_]:= If[ Last[#]==={}, Drop[#,-1], # ]& @ FixedPointList[ iter, num, SameTest->(!NumberQ[#2]&) ]; cf[num_,0]:= {num}; cf[num_,ord_]:= If[ Last[#]==={}, Drop[#,-1], # ]& @ FixedPointList[ iter, num, ord, SameTest->(!NumberQ[#2]&) ]; ]] (* end extract/revision from the 3.0 package *) (* Now my modification to quadratic irrationals *) If[!v4Q, ContinuedFraction[d_?QuadraticIrrationalQ] := ( hf = HoldForm[ContinuedFraction[d]]; If[d < 0, Message[ContinuedFraction::badqua]; Return[hf]]; CFSpecial[#[[1]], #[[3]], #[[2]]]& @ (QuadraticNormalForm[d]) ); ContinuedFraction[n_Integer] := ContinuedFractionForm[{n}]; Format[ContinuedFractionForm[{a_Integer}]]:= a]; bad[m_, n_] := !Positive[m] || !Positive[n] || OddQ[m*n] || GCD[m, n] != 1; PPT[m_Integer, n_Integer] := Module[{mm, nn}, If[bad[m, n], Message[PPT::badarg]; Return[HoldForm[PPT[m, n]]]]; mm = Max[m, n]; nn = Min[m, n]; If[OddQ[m*n] || GCD[m, n] != 1, {}, Sort[{mm^2 - nn^2, 2*mm*nn, mm^2 + nn^2}]]] Options[PythagoreanTriples] = {UseGridBox -> True, ShowSquares -> True, FontSize -> 12, Assumptions -> None, HeadingTextColor -> RGBColor[1, 0, 0], Background -> RGBColor[0.6, 1, 1]}; boldfont[s_, c_, sz_] := StyleBox[s, FontWeight -> Bold, FontSize -> sz, FontColor -> c, FontFamily -> "Times"]; PythagoreanTriples[start_Integer:1, howmany_Integer, opts___Rule] := Module[{gridQ, r, ass}, r = Ceiling[ ((howmany + start)/0.252)^(1/1.94)]; {sqQ, gridQ, ass, bg, htc, sz} = {ShowSquares, UseGridBox, Assumptions, Background, HeadingTextColor, FontSize} /. {opts} /. Options[PythagoreanTriples]; data = Sort[Select[Select[Flatten[Table[Flatten[ {m, n, p = If[ !bad[m, n], PPT[m, n], {}], If[sqQ, p^2, {}]}], {n, r}, {m, n + 1, r}], 1], Length[#1] > 2 & ], ass /. {None -> True, Global`a -> #1[[3]], Global`b -> #1[[4]], Global`c -> #1[[5]]} & ]]; data = Take[data, {If[ass === None, start, 1; start], If[ass === None, start + howmany - 1, Min[Length[data], start + howmany - 1]]}]; If[!gridQ, (#1[[{3, 4, 5}]] & ) /@ data, DisplayForm[StyleBox[GridBox[Prepend[data, Flatten[(boldfont[#1, htc, sz + 2] & ) /@ Join[{"m", "n", "a", "b", "c"}, If[sqQ, (SuperscriptBox[#1, "2"] & ) /@ {"a", "b", "c"}, {}]]]], GridFrame -> 2, ColumnAlignments -> Right, RowAlignments -> Baseline, ColumnLines -> {0, 1, 0, 0, 1, 0}, RowLines -> Flatten[Prepend[Array[{0, 0, 0, 0, 1} & , Ceiling[Length[data]/5.]], 2]], RowSpacings -> {2, 1}], Background -> bg, FontSize -> sz]]]]; (* the following is from Greenberg's paper in Mathematica in Edu. and Research, vol 8 No 3-4. It has not all been declared yet. *) Dijkstra[e_, p_, b_] := Module[{x}, parent = untraversed = Range[b]; dist = p[[1]]; dist[[1]] = 0; Scan[(parent[[#1]] = 1) & , e[[1]]]; While[untraversed != {}, x = First[untraversed]; Scan[If[dist[[#1]] < dist[[x]], x = #1] &, untraversed]; untraversed = Complement[untraversed, {x}]; Scan[If[dist[[#1]] > dist[[x]] + p[[x,#1]], dist[[#1]] = dist[[x]] + p[[x,#1]]; parent[[#1]] = x] & , e[[x]]]]; {parent, dist}]; DiophSolution[L_, parent_, dist_, b_] := Module[{y = Mod[L, b] + 1, x, path, r, s}, If[L >= dist[[y]], node[t_] := parent[[t]]; path = Drop[FixedPointList[node, y], -1]; r = dist[[path]]; s = Drop[r, -1] - Rest[r]; x = Flatten[{(L - dist[[y]])/b, (Count[s, #1] & ) /@ c}], {}]]; Frobenius[dist_, b_] := Max[dist] - b; f[x_] := Select[x - b, #1 > 0 & ]; FrobeniusNonRepresentables[aa_] := (a = Sort[aa]; b = a[[1]]; If[GCD @@ a != 1, Message[FrobeniusNonRepresentables::gcd]; Return[HoldForm[FrobeniusNonRepresentables[aa]]]]; c = Rest[a]; e = Reverse /@ Table[Mod[c+x,b],{x,0,b-1}]+1; p = Array[Infinity & , {b, b}]; Do[p[[i,i]] = 0; p[[i,e[[i]]]] = Reverse[c], {i, b}]; dist = Dijkstra[e, p, b][[-1]]; DeleteCases[Sort[Flatten[FixedPointList[ f, Rest[dist] - b]]], Infinity]); FrobeniusSolve[aa_List, (L_Integer)?NonNegative] := (a=Sort[aa]; If[GCD @@ aa != 1, Message[FrobeniusSolve::gcd]; Return[HoldForm[FrobeniusSolve[aa,L]]]]; perm = Position[aa, #][[1,1]] & /@ a; b = a[[1]]; c = Rest[a]; e = Reverse /@ Table[Mod[c + x, b], {x, 0, b - 1}] + 1; p = Array[Infinity & , {b, b}]; Do[p[[i,i]] = 0; p[[i,e[[i]]]] = Reverse[c], {i, b}]; {parent, dist} = Dijkstra[e, p, b]; ans=DiophSolution[L, parent, dist, b]; If[ans != {},ans[[perm]],{}]); Conductor[aa_] := (a=Sort[aa];b = a[[1]]; c = Rest[a]; e = Reverse /@ Table[Mod[c + x, b], {x, 0, b - 1}] + 1; p = Array[Infinity & , {b, b}]; Do[p[[i,i]] = 0; p[[i, e[[i]]]] = Reverse[c], {i, b}]; Frobenius[Dijkstra[e, p, b][[-1]], b] + 1); Clear[RootsModPrime]; RootsModPrime[f_, x_, p_] := (x /. Solve[f == 0 && Modulus == p, x]) /; Exponent[f, x] == 1; RootsModPrime[1, x_, p_] := {}; RootsModPrime[poly_, x_, 2] := Select[{0,1}, EvenQ[poly /. x -> #] &]; RootsModPrime[ff_, x_, p_] := Module[ {c = -1, h, hf=HoldForm[ModularRoots[f, x, p]], f = PolynomialMod[ff, p], g }, f = PolynomialMod[f, p]; If[f==0, Return[Range[0, p-1]]]; f= PolynomialMod[f * PowerMod[Coefficient[f, x^Exponent[f, x]], -1, p], p]; If[Exponent[f, x] == 0, Return[{}]]; g=PolynomialGCD[f, PolynomialPowerMod[x, p, {f, p}] - x, Modulus -> p]; If[Exponent[g, x] < 2, Return[RootsModPrime[g, x, p]]]; While[!( 1<= Exponent[h, x] < Exponent[g, x]) && c < 10, c++; h = PolynomialGCD[g, PolynomialPowerMod1[x - c, (p - 1)/2, {g, p}] - 1, Modulus -> p]]; If[c == 10, Message[RootsModPrime::nosplt];Return[hf]]; Union[RootsModPrime[h, x, p], RootsModPrime[PolynomialQuotient[ g, h, x, Modulus -> p], x, p]]]; PolynomialPowerMod1[p_, n_, {g_, pr_}] := PolynomialPowerMod[p, n, {g, pr}] /; g =!= 1; PolynomialPowerMod1[_, _, {1, _}] := 0; RootsModPrime[f_, x_, p_]:= QuadraticCongruenceSolve[Reverse@CoefficientList[f,x],p] /; Exponent[f, x]==2; RootLift[poly_, x_, a_List, p_, k_] := Module[ {deriv=D[poly,x], sols = a,s, newsols}, Do[newsols = {}; s = Transpose[{poly/p^j/. x -> sols,deriv/.x-> sols}]; Do[Which[ Mod[s[[i,1]],p] != 0 && Mod[s[[i,2]], p] !=0 , AppendTo[newsols, Mod[sols[[i]]+ Mod[-s[[i,1]] * PowerMod[s[[i,2]],-1,p], p^(j+1)] p^j , p^(j+1)]], Mod[s[[i,1]], p] == 0 && Mod[s[[i,2]], p] != 0 , newsols=Append[newsols, sols[[i]]], Mod[s[[i,1]], p] == 0 && Mod[s[[i,2]], p] == 0 , newsols = Join[newsols, Table[Mod[sols[[i]]+m p^j, p^(j+1)], {m,0,p-1}]], True, {}], {i, Length[s]}]; sols = newsols, {j, k-1} ]; Sort[sols]]; RootsModPrimePower[poly_, x_, p_, k_] := RootLift[poly, x, RootsModPrime[poly, x, p], p, k]; RootsMod[poly_, x_, n_] := Module[{fac=FactorInteger[n]}, Sort[ChineseRemainder[#, #[[1]]^#[[2]]& /@ fac]& /@ Distribute[List @@ RootLift[poly,x, RootsModPrime[poly,x,#[[1]]],#[[1]],#[[2]]]& /@ fac,List]]]; buch[1][u_] := 1/u; buch[2][u_] := (1 + Log[-1 + u])/u; buch[3][u_] := (12 + Pi^2 + 12 Log[-1 + u] + 12 Log[-2 + u] Log[-1 + u] + 12 PolyLog[2, 2 - u])/(12 u); buch[4][u_] := (-(Pi^2*(-1 + Log[128])) + Re[6*Log[2]*(Log[4] - Log[3 - u])*Log[3 - u] + 6*Log[2]*Log[-2 + u]*Log[4*(-2 + u)] + 3*(4 + Pi^2 - 4*(-1 + Log[6 - 2*u] - 2*Log[-3 + u])* Log[-2 + u])*Log[-1 + u] + 3*(4 - 2*Log[2]^3 + Log[2]*Log[9/4]*Log[16] - 4*Log[4/3]*PolyLog[2, 3/4] - 2*Log[81/16]* PolyLog[2, 3/2] + Log[16]*PolyLog[2, 4]) + 12*Log[(2*(-2 + u))/(-3 + u)]*PolyLog[2, 1 + (-3 + u)^(-1)] + 12*Log[(-1 + u)/(2*(-2 + u))]* PolyLog[2, 1 + (-2 + u)^(-1)] + 12*(1 + Log[2*(-2 + u)])*PolyLog[2, 2 - u] + 12*(Log[-2 + u] - Log[(2*(-2 + u))/(-3 + u)])* PolyLog[2, (3 - u)/2] - 12*Log[(-1 + u)/(2*(-2 + u))]* PolyLog[2, -2 + u] - 12*Log[(2*(-2 + u))/(-3 + u)]* PolyLog[2, (2*(-2 + u))/(-3 + u)] + 12*(Log[-2 + u] + Log[(-1 + u)/(2*(-2 + u))])* PolyLog[2, (-1 + u)/2] - 12*Log[(-1 + u)/(2*(-2 + u))]* PolyLog[2, (-1 + u)/(2*(-2 + u))] + 12*PolyLog[3, -1/2] - 12*PolyLog[3, 3/4] + 24*PolyLog[3, 3/2] - 9*PolyLog[3, 4] - 12*PolyLog[3, 1 + (-3 + u)^(-1)] - 12*PolyLog[3, (3 - u)/2] + 12*PolyLog[3, (2*(-2 + u))/(-3 + u)] - 3*PolyLog[3, (-2 + u)^2] - 12*PolyLog[3, (-1 + u)/2] + 12*PolyLog[3, (-1 + u)/(2*(-2 + u))] - 12*PolyLog[3, (-1 + u)/(-2 + u)] + (21*Zeta[3])/2])/ (12*u); Buchstab\[Omega][u_] := buch[Floor[u] - If[u == Floor[u] && u > 1, 1, 0]][u] /; 1 <= u <= 5; buchto10[i_][u_] := buch[i][u] /; 1 <= i <= 3; buchinitialize := Module[{x, u}, Do[buchto10[i][u_] := Evaluate[x[u] /. First[NDSolve[{ x'[u] ==(buchto10[i-1][u-1] - x[u])/u, x[i] == buchto10[i-1][i]}, {x[u]}, {u, i, i + 1}, WorkingPrecision -> 30, PrecisionGoal -> 20]]], {i, 4, 10}]]; buchdoneQ = False; Buchstab\[Omega]To10[u_] := (If[ !buchdoneQ, buchinitialize; buchdoneQ = True]; buchto10[Floor[u] - If[u == Floor[u] && u > 1, 1, 0]][u]); SetAttributes[{MiddleRemainder}, Listable]; SquareAbove[n_] := Unfactor[FactorInteger[n] /. {x_Integer, y_} :> {x, y + Mod[y,2]}]; midremain[z_, _, x_] := x /; x <= z; midremain[z_, m_, x_] := midremain[z, x, Mod[m, x]]; MiddleRemainder[n_, x_, f_] := midremain[ Floor[Sqrt[n/f]], n, x]; PrimitiveReps[n_, f_, g_] := Module[{temp}, ans = Complement[(If[IntegerQ[temp = Sqrt[(n - f*#1^2)/g]], {#1, temp}, {}] & ) /@ MiddleRemainder[n, Select[SqrtModAll[-g*PowerMod[f, -1, n], n], #1 <= n/2 &], f], {{}}]; If[f == g, Union[ans, ans /. {x_Integer, y_} :> {y, x}], ans]] /; f + g + 1 <= n; PrimitiveReps[n_, f_, g_] := Module[{ans = {}, ans1 = {}, ans2 = {}}, If[IntegerQ[s = Sqrt[n/g]], ans = {{0, s}}]; If[IntegerQ[s = Sqrt[n/f]], ans1 = {{s, 0}}]; If[n == f + g, ans2 = {{1, 1}}]; Join[ans, ans1, ans2]] /; n <= f + g; QuadraticRepresentation[n_?Negative, _, _] := {} QuadraticRepresentation[0, _, _] := {{0,0}} QuadraticRepresentation[n_Integer?Positive, f_Integer?Positive, g_Integer?Positive] := ( d = GCD[f, g]; If[Mod[n, d] != 0, {}, QuadraticRepresentation[ n/d, f/d, g/d]]) /; GCD[f, g] != 1; QuadraticRepresentation[n_Integer?Positive, f_Integer?Positive, g_Integer?Positive] := Module[{d=GCD[f,n],d1}, d1 = SquareAbove[d]; (# * {1, Sqrt[d1]} &) /@ QuadraticRepresentation[ n/d, f/d, g d1/d] ] /; GCD[f, n] != 1; QuadraticRepresentation[n_Integer?Positive, f_Integer?Positive, g_Integer?Positive] := Module[{ ans = {}, ans1 = {}, ans2 = {}}, If[IntegerQ[s = Sqrt[n/g]], ans = {{0, s}}]; If[IntegerQ[s = Sqrt[n/f]], ans1 = {{s, 0}}]; If[n == f + g, ans2 = {{1, 1}}]; Union[ans, ans1, ans2]] /; n <= f + g; QuadraticRepresentation[n_Integer?Positive, f_Integer?Positive, g_Integer?Positive] := Union @@ (Sqrt[#]*PrimitiveReps[n/#, f, g] & ) /@ Select[Divisors[n], IntegerQ[Sqrt[#1]] & ] /; f + g + 1 <= n; Options[PartialFactor] = {TimeLimit->10}; PartialFactor[n_Integer, u_, opts___] := Module[{U,F,tm,ansnew,ans,z,p}, {tm}={TimeLimit} /. {opts} /.Options[PartialFactor]; ans=Drop[FactorInteger[n, FactorComplete\[Rule]False],-1]; F=Unfactor[ans]; U=n/F; If[F>n^u, Return[ans]]; While[F a+1]]]; SumOfSquaresR1[3,n_] := 0 /; Mod[FourFreePart1[n], 8] == 7; SumOfSquaresR1[3,n_] := (SumOfSquaresR1[3, n/4]) /; Mod[n,4] == 0; SumOfSquaresR1[3,n_] := (Module[{n1=OddPartTrue[Sqrt[n]],ss}, ss=Select[FactorInteger[n1],Mod[#[[1]],4]==3&]; If[ss=={},6 n1,{q,a}=Transpose[ss]; 6 n1/Times@@(q^a) (Apply[Times,q^a+2 (q^a-1)/(q-1)])]])/; IntegerQ[Sqrt[n]]; rprim[n_?OddQ]:=(24 Sum[JacobiSymbol[s,n],{s,n/4}])/;Mod[n,4]==1; rprim[n_?OddQ]:=(8 Sum[JacobiSymbol[s,n],{s,n/2}])/;Mod[n,8]==3; rprim[n_?OddQ]:=0/;Mod[n,8]==7; (*odd square-free case;because rprim sums to n/4, we restrict this case to small n*) SumOfSquaresR1[3, n_?OddQ] := (Plus@@(rprim/@(n/Select[Divisors[n],IntegerQ[Sqrt[#]]&])))/; n<500 && MoebiusMu[n] != 0; SumOfSquaresR1[3, n_?EvenQ]:=Module[{s=0}, Do[s += SumOfSquaresR1[2,n-k^2], {k, Sqrt[N[n]]}]; 2 s + SumOfSquaresR1[2, n]]; SumOfSquaresR1[4,n_Integer?Positive]:= If[EvenQ[n],24,8]*DivisorSigma[1,OddPartTrue[n]]; SumOfSquaresR1[5,n_?((IntegerQ[Sqrt[#]] && (#>1) && !IntegerQ[Log[2,#]])&)]:= Module[{nn,n1,p,b},nn=Sqrt[n]; n1=OddPartTrue[nn]; {p,b}=Transpose@FactorInteger[n1]; 10 (8^(1+Log[2,nn/n1])-1)/7*Times@@((p^(3 b+3)-p^(3 b+1)+p-1)/(p^3-1))]; SumOfSquaresR1[6,n_Integer?Positive]:= 16*n^2*(DivisorSigma1[-2,n]-DivisorSigma3[-2,n])-4*(DivisorSigma1[2,n]- DivisorSigma3[2,n]); SumOfSquaresR1[8,n_Integer?Positive]:=16 (-1)^n AlternatingDivisorSigma[3,n]; SumOfSquaresR1[d_, n_Integer?Positive]:=Module[{ s = 0}, Do[s += SumOfSquaresR1[d-1,n-k^2], {k, Sqrt[N[n]]}]; 2 s+SumOfSquaresR1[d-1, n]]; freqs[list_List] := {Length[#],First[#]}&/@Split[Sort[list]]; SumOfSquaresRPrimitive[2, n_Integer] := Which[ n == 1 || n == 2, 4, Mod[n, 4] == 0, 0, True, dd = DeleteCases[ Mod[First /@ FactorInteger[n], 4], 2]; If[MemberQ[dd, 3], 0, 4 2^Length[dd]]]; SumOfSquaresRPrimitive[3, n_] := If[Mod[n,4]==0, 0, SumOfSquaresR[3, n] / P[n]]; (* by Grosswald (4.9 and 4.10) *) SumOfSquaresRPrimitive[k_Integer, n_Integer]:= ( If[MoebiusMu[n] != 0, Return[SumOfSquaresR[k,n]]]; Plus @@ (((Multinomial @@ (First /@ freqs[#]))*(2^Plus @@ Sign[#])&) /@ Select[SumOfSquaresRepresentations[k,n], GCD @@ # == 1 &])); MyKroneckerSymbol[1|-1, 0] = 1; MyKroneckerSymbol[_Integer, 0] := 0; MyKroneckerSymbol[a_Integer, -1] := Sign[a + 0.1]; MyKroneckerSymbol[a_, b_?OddQ]:= If[b > 0, JacobiSymbol[a, b], Sign[a] JacobiSymbol[a, -b]] MyKroneckerSymbol[a_Integer?EvenQ, b_Integer?EvenQ] := 0; (*KroneckerSymbol[a_Integer,b_Integer?EvenQ]:= Module[{bb=Abs[b]/2,v=1},While[EvenQ[bb],bb/=2;v=-v]; If[v\[Equal]-1,1,{1,0,-1,0,-1,0,1,0}[[Mod[a,8,1]]]]*If[b<0&&a<0,-1,1]* JacobiSymbol[a,bb]];*) MyKroneckerSymbol[a_Integer, b_Integer?EvenQ] := ( Off[CompiledFunction::"cfsa"]; kk = kscom[a,b]; On[CompiledFunction::"cfsa"]; kk); kscom = Compile[{{a,_Integer,0},{b,_Integer,0}}, Module[{bb = Quotient[Abs[b], 2], v}, v = 1; While[EvenQ[bb], bb = Quotient[bb, 2]; v = -v]; If[v == -1, 1, {1,0,-1,0,-1,0,1,0}[[Mod[a,8] /. 0->8]]* If[b<0 && a<0, -1, 1] * JacobiSymbol[a, bb]]]]; (* old code KStable = {1,0,-1,0,-1,0,1,0}; MyKroneckerSymbol[a_Integer,b_Integer?OddQ]:= JacobiSymbol[a,b] /; NonNegative[b]; MyKroneckerSymbol[a_Integer,b_Integer?OddQ]:= Sign[a+.1] JacobiSymbol[a,-b] /; Negative[b]; MyKroneckerSymbol[a_Integer,0]:= If[Abs[a]==1,1,0]; MyKroneckerSymbol[a_Integer?EvenQ,b_Integer?EvenQ]:=0; MyKroneckerSymbol[a_Integer,b_Integer]:= Module[{k,v=IntegerExponent[b,2],bb,aa=a}, bb=Abs[b/2^v]; k=If[EvenQ[v],1,KStable[[Mod[a,8] /. 0->1]]]; If[b<0 && a<0,k = -k]; While[aa !=0, v=IntegerExponent[aa,2]; aa/= 2^v; If[OddQ[v], k *= KStable[[Mod[bb,8,1]]]]; k*=If[OddQ[(aa-1)/2] && OddQ[(bb-1)/2],-1,1]; {aa,bb} ={Mod[bb,Abs[aa]],Abs[aa] }]; If[bb>1, 0, k]]; *) q[n_] := Module[{nn, a}, a = IntegerExponent[n, 4]; nn = Mod[n/4^a, 8]; Which[nn == 7, 0, nn == 3, 2^(-a), True, 3*2^(-a - 1)]]; P[n_] := Module[{p, b}, Times @@ ( ({p, b} = #; 1 + Sum[p^(-j), {j, 1, b - 1}] + 1/(p^b*(1 - JacobiSymbol[-(n/p^(2*b)), p]/p))) &) /@ (DeleteCases[FactorInteger[n], {2, _} | {_, 1}] /. {x_Integer, y_Integer} :> {x, Quotient[y, 2]})]; w[d_] := Which[d < -4, 2, d == -4, 4, d == -3, 6]; MyFundamentalDiscriminantQ[d_] := d != 1 && (Mod[d, 4] == 1 && MoebiusMu[d] != 0 || Mod[d, 4] == 0 && MoebiusMu[d/4] != 0 && Mod[d/4, 4] >= 2); FundamentalDiscriminantRep[d_] := Module[{fi = FactorInteger[d]}, {Times @@ Cases[fi, {p_, (e_)?OddQ} :> p], Cases[fi, {p_, _?(#1 > 1 & )} :> p]}] /; Mod[d, 4] == 1; FundamentalDiscriminantRep[d_] := Module[{dd, fi, oddsquarefree}, fi = FactorInteger[d/2^IntegerExponent[d, 2]]; oddsquarefree = {Times @@ Cases[fi, {p_, e_?OddQ} :> p], Cases[fi, {p_, _?(#1 > 1 & )} :> p]}; If[Mod[oddsquarefree[[1]], 4] == 1, {oddsquarefree[[1]], Prepend[oddsquarefree[[2]], 2]}, {4*oddsquarefree[[1]], If[Mod[d, 8] == 0, Prepend[oddsquarefree[[2]], 2], oddsquarefree[[2]]]}]] /; Mod[d, 4] == 0; (* we use r3 to get h, via formulas in Grosswald *) ClassNumber[-4 | -3] = 1; ClassNumber[n_Integer] := SumOfSquaresR[3,-n]/(24P[-n]) /; n < 0 && Mod[n, 8] == 5; ClassNumber[n_Integer] := SumOfSquaresR[3,-n]/(16 q[-n] P[-n]) /; n < 0 && Mod[n,4] == 0 && q[-n] != 0; (* Call Vardi's package code in the next case NO: better to use Cohen's erfc formula ClassNumber[n_Integer] := Length[ClassList[n]] /; n < 0 && Mod[n, 8] == 1 && MoebiusMu[n] != 0; *) ClassNumber[d_] := Module[{rr = Range[Sqrt[-N[d]/(2 Pi) Log[-d]]]}, Round[KroneckerSymbol[d, rr].( Sqrt[-N[d]] * E^(Pi rr^2/d) / (Pi rr) + Erfc[rr * Sqrt[-Pi/N[d]]])]] /; d < 0 && Mod[d, 8] == 1 && MoebiusMu[d] != 0; (* Use Cohen p. 233 in the nonfundamental case *) ClassNumber[n_Integer] := Module[{f, plis}, {f, plis} = FundamentalDiscriminantRep[n]; w[n] ClassNumber[f]/w[f] Sqrt[n/f]* Times @@ (1 - (If[#==2,KroneckerSymbol, JacobiSymbol][f, #] & ) /@ plis/plis)] /; n < 0 && Mod[n, 4] <= 1 && !FundamentalDiscriminantQ[n]; VisualSieve[] := Module[{ar, f, g, radius, pbound, tsty, smsty, bg, xtop, vals, prod, strikers, nonstrikers}, ar = 0.35; radius = 1.25; diskcolor = RGBColor[0.5, 1, 0.5]; bg = RGBColor[0.6, 0.6, 1]; tsty = {FontSize -> 14, FontFamily -> Times}; smsty = {FontSize -> 9, FontFamily -> "Times"}; f[x_] := {x^2, x}; g[x_] := {x^2, -x}; pbound = 11; xtop = 121; vals = Range[2, xtop]; prod = 1; strikers = Range[2, pbound]; nonstrikers = {}; Do[p = Prime[i]; stemp = Show[Block[{$DisplayFunction = Identity}, {Plot[{-Sqrt[x], Sqrt[x]}, {x, 0, xtop*1.05}, PlotStyle -> {{RGBColor[1, 1, 0], Dashing[{0.001, 0.005}]}}], Graphics[{(Line[{f[p], g[#1]}] & ) /@ Select[Range[2, xtop/p], GCD[#1, prod] == 1 & ], ({Text[#1, {#1^2, #1}, TextStyle -> tsty], Text[#1, {#1^2, -#1}, TextStyle -> tsty]} & ) /@ strikers, {GrayLevel[0.33], ({Text[#1, {#1^2, #1}, TextStyle -> tsty], Text[#1, {#1^2, -#1}, TextStyle -> tsty]} & ) /@ nonstrikers}, ({{diskcolor, Disk[{#1, 0}, radius*{1, ar}]}, {RGBColor[1, 0, 0], Text[#1, {#1, 0}, TextStyle -> smsty]}} & ) /@ Select[vals, #1 < p^2 & ], {{(Text[#1, {#1, 0}, TextStyle -> smsty] & ) /@ Select[vals, #1 >= p^2 & ]}}}]}], PlotRange -> {{-3, xtop*1.05}, 1.15*{-Sqrt[xtop], Sqrt[xtop]}}, AspectRatio -> ar, Axes -> False, FrameTicks -> None, Background -> bg, Frame -> True]; prod *= p; nonstrikers = Join[nonstrikers, Select[strikers, GCD[#1, prod] != 1 & ]]; strikers = Select[strikers, GCD[#1, prod] == 1 & ]; vals = Select[vals, #1 == p || Mod[#1, p] != 0 & ], {i, 5}]]; mult[array_, a_, m_] := Mod[RotateRight[array] - a*array, m]; square[array_, m_] := Mod[Plus @@ Mod[Partition[Append[ListConvolve[array, array, {1, -1}, 0], 0], Length[array]], m], m]; ArrayPolynomialPowerMod[a_, power_, r_, m_] := Fold[If[#2 == 1, mult[square[#1, m], a, m], square[#1, m]] & , Join[{-a, 1}, Table[0, {r - 2}]], Rest[IntegerDigits[power, 2]]]; AKSPolynomial[p_, a_, r_] := Module[ {iden = ArrayPolynomialPowerMod[a, p, r, p]}, iden[[{0, Mod[p, r]} + 1]] = Mod[iden[[{0, Mod[p, r]} + 1]] + {a, -1}, p]; iden]; AKSBPPolynomial[n_] := Module[{r = 2, ans = True}, While[If[Mod[n, r] == 0, ans = False; Break[]]; Mod[n^2 - 1, r] == 0, r++]; AKSPolynomial[n, 1, r]] /; n > 6; Protect[ClassNumber, CarmichaelLambda, Pi, UpArrow, PrimePi, UpArrow, Normal, Solve, PowerMod, PrimitiveRoot, NextPrime, ChineseRemainderTheorem,ChineseRemainder, ContinuedFraction, ContinuedFractionForm, QuadraticRepresentation]; End[]; EndPackage[]; (* PowerQMod[a_, n_, m_]:=Module[{x}, {}!=RootsModPrime[x^n-a, x,m]] *) (* fourthpower[0, _] := True; fourthpower[m_, p_] := (s=SqrtMod[m,p]; If[!IntegerQ[s], Return[False]]; LegendreSymbol[s,p]==1|| LegendreSymbol[-s,p]==1) ModPowerQ[a_,r_,p_] := Module[{d=GCD[r,p-1]}, Mod[a,p]==0 || d==1 ||1==PowerMod[a, (p-1)/d,p] ]; *)