Chapter Eight:

Some Algorithms of Number Theory

8.1 The Ancient and Modern Euclidean Algorithm

GCD[10^500+134270000, 10^500+12466625460000000]//Timing

gcd[a_, b_] := gcd[b, Mod[a, b]]
gcd[a_, 0]  := a

gcd[55492163484645888, 33787461185758215936]

FullGCD[a_, b_] := FullGCD[a, b, {a}, {}]
FullGCD[a_, b_, r_, q_] :=
	FullGCD[b, Mod[a,b], Append[r, b], Append[q, Quotient[a,b]]]
FullGCD[a_, 0, r_, q_]  := {r, q}

FullGCD[33787461185715936, 92163484645888]

Map[Length[Digits[#]] &, First[%]]

FullGCDRound[a_, b_] := FullGCDRound[a, b, {a}, {}]
FullGCDRound[a_, b_, r_, q_] :=
	FullGCDRound[b, Min[Mod[a,b], b-Mod[a,b]],
		    Append[r, b], Append[q, Round[a/b]]]
FullGCDRound[a_, 0, r_, q_]  := {r, q}

FullGCDRound[33787461185715936, 92163484645888]

Table[(Plus @@ EulerPhi /@ Range[200 n]) /
      (200 n + Binomial[200 n, 2]),            {n, 10}]//N

8.2 The Extended Euclidean Algorithm

ExtendedGCD[Prime[1000], 1000000]

FullExtendedGCD[a_, b_] := FullExtendedGCD[a, b, {a}, {}]
FullExtendedGCD[a_, b_, r_, q_] := FullExtendedGCD[
	b, Mod[a,b], Append[r, b], Append[q, Floor[a/b]]]
FullExtendedGCD[a_, 0, r_, q_]  :=
	Block[{s = {1, 0}, t = {0, 1}},
	Scan[(AppendTo[t, t[[-2]] - t[[-1]] #];
		  AppendTo[s, s[[-2]] - s[[-1]] #]) &, q];
	{r, q, s, t}]

FullExtendedGCD[1000000, 7919]

%[[3]]*1000000 + %[[4]]*7919

FullExtendedGCD[73, 27]//TableForm

Q[{}] := 1 
Q[{q_}] := q
Q[qlist_] := First[qlist] Q[Rest[qlist]] + Q[Drop[qlist,2]] /;
										Length[qlist] > 1

Q[{2,1,2,2,1,2}]

Q[{126, 3, 1, 1, 2, 3, 1, 1, 56}]

Table[Q[Drop[{126, 3, 1, 1, 2, 3, 1, 1, 56}, n]],{n, 9, 0, -1}]

Expand[Q[{a,b,c,d,e}]]

Binet[n_] := (GoldenRatio^n - (-1/GoldenRatio)^n )/Sqrt[5]
Attributes[Binet] = Listable
Binet[Range[0, 10]]//N

Note: I know that following will slow down all instances of Power, and so should be avoided. But I could not find a simple way to get 0 by using replacement rules.

Unprotect[Power]
GoldenRatio^2 = 1 + GoldenRatio
Protect[Power];

Simplify[Binet[n]+ Binet[n+1] - Binet[n+2]]

8.3 Linear Diophantine Equations: The Problem of Shielding Variable Names

ExtendedGCD[28,91][[2]] 147 / GCD[28,91]

Simplify[ 28 x + 91 y /. {x -> -63 + 13 m , y -> 21 - 4 m}]

LinearDiophantine[a_ x_ + b_ y_ == c_] := Block[{s, t, d},
	{d, {s, t}} = ExtendedGCD[a,b];
	If[Mod[c, d] == 0,
		{ {x -> c s / d, y -> c t / d}, {b/d, -a/d} },
		{}]]

LinearDiophantine[28 x + 91 y == 147]

LinearDiophantine[28 t + 91 u == 147]

Linear Diophantine Package

Remove[LinearDiophantine]  
		(* avoids rebooting to use following package *)

BeginPackage["LinearDiophantine`"];

LinearDiophantine::usage =
 "LinearDiophantine[(a x + b y) == c] returns a solution to the equation in integers, together with a pair such that the sum of the particular solution with all integer multiples of the pair yield the general solution. If no solution exists, {} is returned.";

Begin["LinearDiophantine`Private`"];

LinearDiophantine[a_ x_ + b_ y_ == c_] := Block[{s, t, d},
	{d, {s, t}} = ExtendedGCD[a,b];
	If[Mod[c, d] == 0,
		{ {x -> c s / d, y -> c t / d}, {b/d, -a/d} },
		{}]]

End[];
EndPackage[];

LinearDiophantine[28 t + 91 u == 147]

8.4 The Chinese Remainder Theorem

SetAttributes[PowerMod, Listable]
Unprotect[PowerMod];
PowerMod[a_, 1, n_] := Mod[a, n]; 
Protect[PowerMod];

(* the bug in the exponent-1 case that the preceding addresses
   is gone in version 2.0 *)

CRT[alist_, mlist_] := Block[{m = Times @@ mlist, mm},
	mm = m/mlist;
	Mod[Apply[Plus, alist*mm*PowerMod[mm, -1, mlist]], m]]

CRT[{0,1,2,3}, {4,9,25,49}]

Mod[14752, {4,9,25,49}]

n = 5
CRT[Range[0, -n + 1, -1], Prime[Range[1, n]]^2]

8.5 Continued Fractions

Solve[3 + 1/x == Pi]

%//N

<<:NumberTheory:ContinuedFractions.m

ContinuedFraction[N[Pi,50], 12]

Normal[%]

ContinuedFraction[N[Sqrt[2], 50], 12]

ContinuedFraction[N[GoldenRatio, 50], 12]

ContinuedFraction[N[E, 50], 12]

Normal[ContinuedFractionForm[{3,7,15,1,292,1,1,1,2,1,3,"1"}]]

quotients[a_, b_] := (AppendTo[qlist, Floor[a/b]];
				      quotients[b, Mod[a, b]])
quotients[a_, 0] := 0

ContinuedFraction[r_Rational] := (qlist = {};
	quotients[Numerator[r], Denominator[r]];
	ContinuedFractionForm[qlist])

CFFormatted[r_Rational] :=
	Block[{qlist = ContinuedFraction[r][[1]]},
		Normal@ContinuedFractionForm[
			Rest[Append[qlist, ToString[Last[qlist]]]]]]
	
CFFormatted[x_Real, n_Integer?Positive] := 
	Block[{qlist = ContinuedFraction[x, n][[1]]},
		Normal[ContinuedFractionForm[
			Rest[Append[qlist, ToString[Last[qlist]]]]]]]

ContinuedFraction[43/99]

CFFormatted[43/99]

ContinuedFraction[1000000/7919]

Normal[%]

convergents[q_Rational] :=
	Block[{qlist = ContinuedFraction[q][[1]]},
	Map[Normal[ContinuedFractionForm[Take[qlist, #]]] &,
				Range[Length[qlist]]]]

convergents[1000000/7919]

convergents[r_Real, n_Integer] :=
	Block[{qlist = ContinuedFraction[r, n][[1]]},
	Map[Normal[ContinuedFractionForm[Take[qlist, #]]] &,
				Range[Length[qlist]]]]

convergents[N[Pi,100], 10]

PeriodicCF::usage = "PeriodicCF[initial:{}, periodic] returns the quadratic irrational represented by an eventually periodic simple continued fraction; the default is the case of no nonperiodic part"

QuadConjugate[a_Integer] = a
QuadConjugate[a_:0 + b_:1 * Sqrt[c_Integer]] := a - b Sqrt[c]

PeriodicCF[initial_List:{}, periodic_List] := 
     Block[{raw, persoln},
     persoln = Solve[x == 
          Normal[ContinuedFractionForm[Append[periodic,x]]]];
     persoln = If[N[persoln[[1,1,2]]] < 0, persoln[[2,1,2]],
                                           persoln[[1,1,2]]];
     raw = Together[
      Normal[ContinuedFractionForm[Append[initial, persoln]]]];   
  Simplify[
    Expand[Numerator[raw]   QuadConjugate[Denominator[raw]]]/
    Expand[Denominator[raw] QuadConjugate[Denominator[raw]]] ]]

PeriodicCF[{1}]

PeriodicCF[{1,2,3,4}]

PeriodicCF[{2}]

PeriodicCF[{2,2},{1}]

PeriodicCF[{1},{2,3}]

PeriodicCF[{10},{100,1000}]

8.6 Egyptian Fractions

EgyptianFraction[0]  := {}
EgyptianFraction[q_] := Prepend[
		EgyptianFraction[q - 1/Ceiling[1/q]], Ceiling[1/q]]
Attributes[EgyptianFraction] = Listable

EgyptianFraction[4/23]

EgyptianFraction[Range[0,1,1/23]]

Remainders[0]  := {}
Remainders[q_] := Prepend[
		Remainders[q - 1/Ceiling[1/q]], q - 1/Ceiling[1/q]]

Remainders[19/23]

Table[EgyptianFraction[{n/131, (65 - n)/131}], {n, 1, 32}]

EgyptianFractionOdd[0, old_]  := {}
EgyptianFractionOdd[q_, old_:1] := Block[
	{m = Max[old+1, Ceiling[1/q]]},
	Prepend[While[EvenQ[m], m++];
	EgyptianFractionOdd[q - 1/m, m], m]]
Attributes[EgyptianFractionOdd] = Listable

EgyptianFractionOdd[Range[1,4]/5]

Ceiling[Log[10, EgyptianFractionOdd[4/13]//N]]

RemaindersOdd[0, old_]  := {}
RemaindersOdd[q_, old_:1] := Block[
	{m = Max[old+1, Ceiling[1/q]]},
	Prepend[While[EvenQ[m], m++];
	RemaindersOdd[q - 1/m, m], q - 1/m]]

Map[Numerator, RemaindersOdd[999/1001]]

8.7 Prime Certificates

version 2.0 notes:

much of this is in the :NumberTheory:PrimeQ.m package distributed with version 2.0.

uses Listability of PowerMod; not built-in; 
see section on Chinese Remainder Theorem

PrattCertify::usage = "PrattCertify[p] returns a Pratt certificate of p's primality. If p is not a prime number the program does not halt."

SetAttributes[PrattCertify, Listable]

PrattCertify[2] := 2

PrattCertify[p_Integer] := Block[{x = 2, primes},
   primes = First[Transpose[FactorInteger[p - 1]]];
   While[
	Apply[Or, PowerMod[x, (p-1)/#, p] == 1 & /@ primes], x++];
   Return[{p, x, PrattCertify[primes]}]] /; p > 2

PrattCertify[7919]

PrattTree::usage = "PrattTree[p] returns a Pratt certificate of p's primality in the form of a tree. If p is not a prime number, the program will probably halt with a NOT PRIME message."

Attributes[PrattTree] = Listable
 
space[n_] := StringJoin @@ Table["     ", {n}];

PrattTree[2, level_] := Print[space[level], 2];

PrattTree[p_Integer, level_:0] := Block[
   {x = 2, primes = First[Transpose[FactorInteger[p-1]]]},
   While[Apply[Or, PowerMod[x, (p-1)/#, p] == 1 & /@ primes],
                                                          x++];
	If[level == 0 && PowerMod[x, p-1, p] != 1,
		Return["NOT PRIME"],
			Print[space[level], p, " ±±± ", x];
			PrattTree[primes, level+1];]]   /; p > 2

PrattTree[7919]

Bonus on Pratt Certification of Primes

The following modification to PrattCertify runs 10% faster on large inputs. Idea due to Ferrell Wheeler. The point is that a quadratic residue cannot be a witness (primitive root) because of Euler's Criterion.

SetAttributes[PrattCertify1, Listable]

PrattCertify1[2] := 2

PrattCertify1[p_Integer] := Block[{x = 2, primes},
   primes = First[Transpose[FactorInteger[p - 1]]];
   While[JacobiSymbol[x, p] == 1 ||
	Apply[Or, PowerMod[x, (p-1)/#, p] == 1 & /@ primes], x++];
   Return[{p, x, PrattCertify1[primes]}]] /; p > 2

version 2.0 notes: Several useful number theory functions will be distributed with version 2.0. These include a prime certifying package by Ilan Vardi which uses Pratt's method for small primes and a state-of-the-art elliptic curve method for larger primes.