BeginPackage["Numbers`"] (* Copyright: Erich Kaltofen, 4/27/1990. Permission to use provided the copyright notice is not removed. *) NTrace::usage= "Type \"Numbers`NTrace[s__]:=Print[s]\" to trace this package." Numbers`QuotientRemainder::usage= "Numbers`QuotientRemainder[x_Integer, y_Integer] returns {q, r} such that x = q y + r and the absolute value of r is as small as possible." Numbers`Quotient::usage= "Numbers`Quotient[x_Integer, y_Integer] returns the quotient that leaves the absolutely smallest remainder." Numbers`Remainder::usage= "Numbers`Remainder[x_Integer, y_Integer] returns the absolutely smallest remainder." Begin["`quotients`"] QuotientRemainder[x_Integer, y_Integer]:= Block[{q = System`Quotient[x, y], r = System`Mod[x, y]}, If[2 r > Abs[y], Return[ List[q + Sign[y], r - Abs[y] ] ], Return[ List[q, r] ] ] ] Quotient[x_Integer, y_Integer]:=First[QuotientRemainder[x, y] ] Remainder[x_Integer, y_Integer]:=(QuotientRemainder[x, y])[[2]] End[] (* `quotients` *) PeanoMultiplication::usage= "PeanoMultiplication[m_Integer,n_Integer?NonNegative] realizes m * n by m+m+...+m, n times." BitList::usage= "BitList[n_Integer?Positive] returns {b0, b1,..., bk=1} such that n=Sum[bi 2^i, {i, 0, k}]." PeasantMultiplication::usage= "PeasantMultiplication[m_Integer,n_Integer?Positive] realizes m * n by the \"Russian peasant\" method, i.e., by binary exponentiation." BinaryExponentiation::usage= "BinaryExponentiation[m_, n_Integer?Positive] realizes m^n by binary exponentiation using the Horner evaluation of n in binary." Begin["`private`"] PeanoMultiplication[m_Integer,n_Integer?NonNegative]:= (* Peano definition of multiplication *) Block[{nn, Accu}, Accu = 0; For[nn = n, nn > 0, nn--, Accu = Accu + m ]; Return[Accu] ] BitList[n_Integer?Positive]:= (* Returns {b0, b1,..., bk=1} such that n=Sum[bi 2^i, {i, 0, k}] *) Block[{nn, list}, nn = n; list={}; While[ nn > 0, AppendTo[ list, Mod[nn, 2] ]; nn = System`Quotient[nn, 2] ]; Return[list] ] PeasantMultiplication[m_Integer,n_Integer?Positive]:= (* Russian peasant method for multiplication *) Block[{nn, bits, nu, Accu}, bits = BitList[n]; (* Print[StringForm["Bits of n = ``", bits] ]; *) Accu = m; (* Process digits of n from second highest to lowest *) For[nu = Length[bits]-1, nu >= 1, nu--, Accu = Accu + Accu; If[bits[[nu]]==1, Accu = Accu + m ]; (* Print[StringForm["Accumulator = ``", Accu] ] *) ]; Return[ Accu ] ] BinaryExponentiation[m_, n_Integer?Positive]:= Block[{nn, bits, nu, Accu}, bits = BitList[n]; (* Print[StringForm["Bits of n = ``", bits] ]; *) Accu = m; (* Process digits of n from second highest to lowest *) For[nu = Length[bits]-1, nu >= 1, nu--, Accu = Accu * Accu; If[bits[[nu]]==1, Accu = Accu * m ]; (* Print[StringForm["Accumulator = ``", Accu] ] *) ]; Return[ Accu ] ] End[] (* `private` *) GenericExponentiation::usage= "GenericExponentiation[m_, n_Integer?Positive, binop_:Times] realizes m o m o...o m, n times, by binary \"exponentiation\" using Horner evaluation of n, where o is an associative binary function passed in the binop argument." OptimalRadixDigitList::usage= "OptimalRadixDigitList[expo_Integer?Positive]" GenericThurberExponentiation::usage= "GenericThurberExponentiation[m_, optlist_, binop_:Times] computes m^e, where optlist == OptimalRadixDigitList[e] and binop is the binary associative operation used." ModularExponentiation::usage= "ModularExponentiation[a_Integer, b_Integer?Positive, m_Integer?Positive] returns a^b modulo m (the Mathematica PowerMod function), using GenericExponentiation." Begin["`exponentiation`"] GenericExponentiation[m_, n_Integer?Positive, binop_:Times]:= Block[{nn, bits, nu, Accu}, bits = BitList[n]; (* Print[StringForm["Bits of n = ``", bits] ]; *) Accu = m; (* Process digits of n from second highest to lowest *) For[nu = Length[bits]-1, nu >= 1, nu--, Accu = binop[Accu, Accu]; If[bits[[nu]]==1, Accu = binop[Accu, m] ]; (* Print[StringForm["Accumulator = ``", Accu] ] *) ]; Return[ Accu ] ] OptimalRadixDigitList[expo_Integer?Positive]:= (* Returns {rad, ind[k], ind[k-1],...,ind[0]} such that Thurber's scheme with radix r has a minimum number of multiplications; rad = {r, i} such that r = 2^i, 1<=i; inds = {b[s], j[s]} such that d[s] = 2^b[s] * (2 j[s] - 1), b[s]>0, or j[s] such that d[s] = 2 j[s] - 1, or 0 such that d[s] = 0; expo = Sum[d[s] * r^s, {s, 0, k}], the radix r representation of expo. *) Block[{n, quot, list, bestlist, i, r, d, b, ct, bestct, ind}, n = Ceiling[ N[ Log[2, expo] ] ]; bestct = 2 * n; For[i=1, (r=2^i) < n, i++, If[r==2, ct=0, ct = r / 2 (* to compute a^2, a^3, a^5,..., a^(r-1) *) ]; quot = expo; list={}; While[ quot > 0, d = Mod[quot, r]; (* next exponent digit *) If[d == 0, ind = 0, (* put 0 on index list *) (* Else: compute digit == 2^b * odd *) b = 0; While[ EvenQ[d], b = b+1; d = System`Quotient[d, 2] ]; (* now we have digit == 2^b * d *) If[quot == expo && (b > 1 || b==1 && d>1), (* cannot pre-multiply highest order digit; exceptions: digit odd or digit==2, whose powers are already computed *) ct = ct + 1 ]; If[b!=0, ind = List[b, (d+1)/2], ind = (d+1)/2 (* indices into odd powers array *) ]; ct = ct + 1 (* X = X * a^d while computing X = X^r *) ]; ct = ct + i; (* X = X^r *) PrependTo[ list, ind ]; quot = System`Quotient[quot, r]; ]; PrependTo[ list, List[r, i] ]; NTrace[StringForm["ct=``, list=``", ct, list] ]; If[ct < bestct, bestct = ct; bestlist = list] ]; Return[bestlist] ] GenericThurberExponentiation[m_, optlist_, binop_:Times]:= (* E. C. Thurber's exponentiation with fewer multiplications than the binary method [Duke Math. J., vol. 40, pp. 907-913 (1973).*) Block[{r, rho, diglist, mpowers, msquare, ind, b, i, nu, Accu}, r = optlist[[1,1]]; rho = optlist[[1,2]]; (* r = 2^rho *) diglist = Drop[optlist, 1]; (* Pre-compute m^3, m^5,..., m^(r-1) *) msquare = binop[m, m]; Array[mpowers, r/2]; mpowers[1] = m; For[i=2, i<=r/2, i++, mpowers[i] = binop[mpowers[i-1], msquare] ]; (* Take leading digit dk and initialize Accu to m^dk *) ind = First[diglist]; diglist = Drop[diglist, 1]; If[NumberQ[ind], Accu = mpowers[ind], (* leading digit odd *) (* else *) b = ind[[1]]; ind = ind[[2]]; d = 2^b * (2 ind - 1); (* reconstruct digit *) If[d == 2, Accu = msquare, (* else *) ind = d/2; Accu = binop[m, mpowers[ind] ] ] ]; (* Process digits of exponent from second highest to lowest *) For[nu = Length[diglist], nu >= 1, nu--, (* retrieve next digit *) ind = First[diglist]; diglist = Drop[diglist, 1]; If[NumberQ[ind], b = 0, (* else *) b = ind[[1]]; ind = ind[[2]] ]; For[i=1, i<=rho, i++, Accu = binop[Accu, Accu]; If[i==rho-b, (* adjust Accu by odd power such that 2^b will be accounted for by rest of loop *) If[ind != 0, Accu = binop[Accu, mpowers[ind] ] ] ] ]; NTrace[StringForm["Accumulator = ``", Accu] ] ]; Return[ Accu ] ] (* end GenericTurnerExponentiation *) ModularExponentiation[a_Integer, b_Integer?Positive, m_Integer?Positive]:= (* Computes a ^ b modulo m (the Mathematica PowerMod function) *) Block[{localmult}, localmult[x_, y_]:=Mod[x * y, m]; Return[GenericExponentiation[a, b, localmult] ] ] End[] (* `exponentiation` *) FindPrime::usage= "FindPrime[n_, pred_:OddQ] returns the first probabilistic prime p >= n such that pred[n] is true. " FermatPrimeQ::usage= "FermatPrimeQ[p_Integer?Positive] is a primality test for p based on Fermat's little theorem." Begin["`primes`"] FindPrime[n_, pred_:OddQ]:= Block[{m}, For[m = n, True, m++, If[pred[m], If[PrimeQ[m], Return[m] ] ] ] ] (* FindPrime *) FermatPrimeQ[p_Integer?Positive]:= (* Fermat primality test (fails on Carmichael numbers) *) Block[{isPrime, a, i, j, q}, isPrime = True; For[ i=1, i<=2 (* number of tries *), i++, a = Random[ Integer, {2, (p-1)/2} ]; If[ GCD[a, p] != 1, isPrime = False; Break[] ]; q = ModularExponentiation[ a, p-1, p]; If[ 1 != q, isPrime = False; Break[] ] ]; Return[ isPrime ] ] End[] (* `primes` *) CycleIndices::usage= "CycleIndices[x0_, next_, equal_] receives functions next[x_] and equal[x_,y_] and returns the lengths of the pre-period and period of the sequence {x0, next[x0], next[next[x0]], ...}. Equality is determined be the equal function argument." PollardRho::usage= "PollardRho[N_Integer?Positive] demonstrates the Pollard rho integer factoring algorithm. A factor is printed (twice), but cylce indices are returned." Begin["`factorization`"] CycleIndices[x0_, next_, equal_]:= (* Finds the preperiod length mu and the period length lambda of an eventually cycling function next with respect to the predicate equal. *) Block[{X, Y, Z, i, n, lambda}, (* R. Floyd's trick: First find n such that equal[next^(n)[x0], next^(2 n)[x0] ]. Then n is the first exponent in the period that is divisible by lambda. *) X = x0; Y = x0; For[ i = 0, True, i++, (* At this point X == next^(i)[x0], Y == next^(2 i)[x0], and Not[ equal[ next^(j)[x0], next^(2j)[x0] ] ] for all j < i. *) X = next[X]; Y = next[ next[Y] ]; NTrace[ StringForm["X = ``, Y = ``", X, Y] ]; If[equal[X, Y], n = i; Break[] ] ]; NTrace[StringForm["n = ``", n] ]; (* Search for shortest repetion in period *) Z = X; lambda = 0; For[ i = 0, lambda == 0, i++, (* At this point Z == next^(n+i)[x0] *) Z = next[Z]; If[equal[X, Z], lambda = i+1 ] ]; (* Search for the preperiod *) Y = x0; Z = x0; For[i = 0, i < lambda, i++, Z = next[Z] ]; For[mu = 0, Not[ equal[Y, Z] ], mu++, (* At this point we have Y == next^(mu)[x0] and Z = next^(mu+lambda)[x0]. *) Y = next[Y]; Z = next[Z] ]; Return[{mu, lambda}] ] PollardRho[N_Integer?Positive]:= (* Demonstration of the Pollard rho method for integer factoring *) Block[{randomf, gcd}, randomf[x_] := Mod[x^2 + 1, N]; gcd[x_, y_]:= Block[{g = GCD[x-y, N]}, If[g > 1, (* then *) Print[StringForm["GCD =``",g] ]; Return[True], (* else *) Return[False] ] ]; CycleIndices[3, randomf, gcd] ] End[] (* `factorization` *) EndPackage[] (* Numbers` *) Print["Package Numbers` loaded"]