
(*:Version: Mathematica 2.0 *)

(*:Name: NumberTheory`NumberTheoryFunctions` *)

(*:Title: Number Theory Functions *)

(*:Author: Ilan Vardi

*)

(*:Keywords:
	number theory, Chinese Remainder Theorem, 
        class number, primitive roots
*)

(*:Requirements: none. *)

(*:Warnings: SqrtMod needs to have d a residue for every prime 
             factor of n. QuadraticRepresentation needs to have -d
             a quadratic residue only for primes appearing to odd powers.

*)

(*:Limitations: SquareFreeQ, SqrtMod, and QuadraticRepresentation 
                depend on obtaining the factorization of n. 

                SqrtMod[d, n] needs to have odd n.
       
                For a discussion of when QuadraticRepresentation returns
                and answer, see Cox's book.

                PrimitiveRoot does not need to find the factorization 
                of n, but of p-1 if n is a power of p, or twice 
                a prime power of p, p an odd prime.

                ClassList and so ClassNumber has only been implemented
                for negative discriminants. The implementation uses
                the simplest algorithm and is slow for large inputs.
*)

(*:Sources: 

  G.H. Hardy and E.M. Wright, An Introduction to the Theory Of
  Numbers, Oxford University Press, 1988.

  D.E. Knuth, Seminumerical Algorithms, Addison Wesley, 1981.

  D.A. Buell, Binary Quadratic Forms, Springer Verlag, 1989.

  Stan Wagon, Mathematica in Action, Freeman, 1991.

  Stan Wagon, The Euclidean Algorithm Strikes Again,
  American Math. Monthly # 97, (1990), 125--124. 

  K. Hardy, J.B. Muskat, and K.S. Williams, A determinisitic
  algorithm for solving n = fu^2 + gv^2 in coprime integers
  u and v, Math. Comp. #55 (1990) 327-343.

  D.A. Cox, Primes of the Form x^2 + n y^2, Wiley, 1989.  

*)

(*:Summary: This package implements some standard functions from
number theory.

*)


BeginPackage["NumberTheory`NumberTheoryFunctions`"]

Needs["NumberTheory`FactorIntegerECM`"]

SquareFreeQ::usage = "SquareFreeQ[n] is True if n is a square free
integer, False otherwise."


NextPrime::usage = "NextPrime[n] gives the smallest prime p such
that p > n."

ChineseRemainderTheorem::usage = "ChineseRemainderTheorem[list1, list2]
gives the minimal non-negative integer solution of Mod[r, list2] == list1."

SqrtMod::usage = "SqrtMod[d, n] is the square root of d (mod n), where
n is an odd number."

QuadraticRepresentation::usage = "QuadraticRepresentation[d, n] gives
{x, y} where x^2 + d y^2 = n, for odd n and d positive, whenever
such a representation exists. "

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."

ClassList::usage = "ClassList[d] gives a list of inequivalent quadratic
forms of discriminant d, where d < 0 is a square free integer of the
form 4 n + 1, and  a quadratic form a x^2 + b x y + c y^2 is represented
as {a, b, c}."

ClassNumber::usage = "ClassNumber[d] gives the number of inequivalent
quadratic forms of discriminant d, where d < 0 is a square free integer
of the form 4 n + 1."


Begin["`Private`"]


(* SquareFreeQ checks for a square factor of d. This may require
   finding the complete factorization of d. The second argument
   ntd is there so as to avoid doing trial division on the factors
   of d as well. *)


SquareFreeQ[d_Integer, ntd___]:=    
Block[{i = 1, gcd, product = If[ntd === {},Times @@ Prime[Range[1229]]], 
       q = d, flag = False, pr}, 
      If[q < 10^5, Return[MoebiusMu[q] != 0]]; 
      If[ntd === {}, gcd = GCD[q, product]; 
                    q = Quotient[q, gcd];
                    If[GCD[q, gcd] > 1, Return[False]]]; 
      If[q == 1 || PrimeQ[q], Return[True]];
      If[GCD[PowerMod[2, q, q] -2 , q] > 1,  (* Prime Power *)
         Return[MoebiusMu[q] != 0]];  
      If[Select[q^(1/Select[Range[Log[10000, q]], PrimeQ]), IntegerQ] != {},
         Return[False]];                   (* Perfect Power *)
      pr = FactorIntegerECM[q];
      If[GCD[pr, q/pr] > 1, False,
         SquareFreeQ[Min[pr, q/pr], 1] && SquareFreeQ[Max[pr, q/pr], 1]
      ]] 

NextPrime[n_]:=    
 Block[{i = n + 1 + Mod[n, 2]},
        If[n < 10^20,   
           While[Not[PrimeQ[i]], i += 2],
            prod = Apply[Times, Prime[Range[Log[n]]]];
            While[True,
                  If[GCD[ i, prod] == 1, If[PrimeQ[i], Break[]]];
                  i += 2]
           ];  
       i]
 


(* Written by Stan Wagon, see Mathematica in Action.  *)

SetAttributes[PowerMod, Listable]

ChineseRemainderTheorem[list1_, list2_]:=
 Block[{m = Times @@ list2, mm},
        mm = m/list2;
        Mod[Plus @@ (list1 mm PowerMod[mm, -1, list2]), m]
       ] 


SqrtMod::nres := "Warning: n has a prime factor p for which d is a 
nonresidue mod p."

HalfMod[n_, p_]:= Quotient[If[EvenQ[n], n, p - n], 2]

SqrtMod[d_Integer, n_Integer]:= Mod[Sqrt[d], n] /; IntegerQ[Sqrt[d]]

SqrtMod[d_Integer, q_Integer] := 
       If[PrimeQ[q], 

(* Shank's modular Sqrt algorithm, Knuth, Vol. 2, page 626. Note that
   finding a non residue below takes < 2 (Log [n])^2 time by GRH, and
   this avoids using Random[Integer,{1,p}]. *)

       Block[{e, p = q, q2,z,n,v,w,k,t,a,count=0},
              
             e = 1; qq = (p-1)/2;
             While[EvenQ[qq], e++; qq /= 2];
             If[ e == 1, 
                 z = -1,    
                 n = 2;                                
                 While[JacobiSymbol[n,p] != -1, n++];  
                 z = PowerMod[n,qq,p]
               ];
               v = PowerMod[d,(qq+1)/2,p]; 
               w = Mod[PowerMod[d,-1, p] Mod[ v v ,p], p]; 
               While[w != 1, count++; 
                     k = 1; t = Mod[w w,p];
(*  The following loop is very slow for primes of the form 
    h 2^n + 1, large  *)
                     While[t != 1, k++; t = Mod[t t,p]]; 
                     a = PowerMod[z,2^(e-k-1),p];
                     e = k;
                     v = Mod[v a, p];
                     z = Mod[a a,p];
                     w = Mod[w z,p];
                    ]; 
                Return[v]
              ],    


(* When q is a prime power p^k = q, SqrtMod[d, q] is computed using 
   a form of Newton's iteration 

      x -> (x + d/x)/2

   This takes a p^k solution to a p^(2k) solution but requires a
   modular division at each step. 

   One can compute this without any modular inversions (except modulo p) 
   by first computing 

      x = 1/SqrtMod[d, p] (mod p)

   and then iterating 

      x -> (3 x - d x^3)/2

   which does not require modular inversion. The last step is 
   to multiply by d to get the answer. 

   The reason this works is due to the p-adic identity 

      Sqrt[d] = d x /Sqrt[1 - (1 - d x^2)]

      = d x Sum[Binomial[2 k, k]/4^k (1 - d x^2)^k, {k, 0, Infinity}]

   whenever  x^2 = d (mod p^k), some k > 0.

   This method is quite efficient in practice and a direct implementation
   computed the square root of -1 mod 5^100000 in 750 seconds on a DEC3100
   (it took twice as long to check the answer). 

*)
   

        Block[{fi = FactorInteger[q]}, 
               If[Select[First /@ fi, JacobiSymbol[d, #] == -1&] != {},
                  Message[SqrtMod::nres]; Return[]];
          ChineseRemainderTheorem @@  
          Transpose[Map[Function[x, 
              Block[{p, pn, qn, qr, dd, sqrtm, sqrtmin, log},
                  {p, pn} = x;
                   qn = p^pn; 
                   If[pn > 1, log = Ceiling[Log[2, N[pn]]]];
                  sqrtm = SqrtMod[d, p];
                  If[pn == 1, {sqrtm, p}, 
  {If[Mod[d + 1, qn] == 0,  (* d = -1 is special *)
      Mod[Mod[Quotient[#[[1]] (#[[1]]^2 + 3), 2], #[[2]]]& @
      Nest[{Mod[Quotient[#[[1]] (#[[1]]^2 + 3), 2], #[[2]]], #[[2]]^2}&,
       {sqrtm, p^2}, log-1], qn], 
(* else *)	
      Mod[d (HalfMod[Mod[#[[1]] (3 - d #[[1]]^2), #[[2]]], #[[2]]]& @
      Nest[{HalfMod[Mod[#[[1]] (3 - d #[[1]]^2), #[[2]]], #[[2]]], #[[2]]^2}&,
      {PowerMod[sqrtm, -1, p], p^2}, log-1]), qn]],
    qn}]]],
        fi]
           ]]] /; OddQ[q] && JacobiSymbol[d, q] == 1 


(* Uses generalization of Cornacchia's algorithm, see Stan Wagon's article 
   and Hardy, Muskat, and Williams. See Cox's book for the theory of when
   a representation x^2 + d y^2 = n exists.    *)

QuadraticRepresentation::norep := "Warning: n cannot be represented as
x^2 + d y^2 since it splits into factors not in the principal ideal
class of Q(Sqrt[-d])."

QuadraticRepresentation[d_, n_]:= 
Block[{x = SqrtMod[-d, n], y = n},
               If[2 x > n, x = n - x];
               While[x^2 >= n, {x, y} = {Mod[y, x], x}];
               y = Sqrt[(n - x^2)/d];
               If[!IntegerQ[y], 
                  Message[QuadraticRepresentation::norep],
                  {x, y}]
             ]  /; OddQ[n] && JacobiSymbol[-d, n] == 1 && d > 0


(* Written with Ferrell S. Wheeler  *)

PrimitiveRoot[2] = 1

PrimitiveRoot[4] = 3

PrimitiveRoot[n_Integer]:= 
If[OddQ[#], #, n/2 + #]& @ PrimitiveRoot[n/2] /; 
   Mod[n, 4] == 2 && n > 2

PrimitiveRoot[p_Integer] :=  
               Block[{g = 2, flist = (p-1)/First /@ FactorInteger[p-1]},
                     While[JacobiSymbol[g, p] == 1 ||
                           Scan[If[PowerMod[g, #, p] == 1, 
                                   Return[True]]&, flist],
                           g++];
                      g]  /; PrimeQ[p] && p > 2

PrimitiveRoot[n_Integer]:= 
 Block[{q, fi = FactorInteger[n], g},
        If[Length[fi] > 1, Return[]];
        q = fi[[1,1]];
        g = PrimitiveRoot[q];
        If[PowerMod[g, q - 1, q^2] == 1, g += q];
        g]  /;  GCD[PowerMod[2, n, n]-2, n] > 1

(* ClassList only works for negative d, and is inefficient for large d. *)

ClassList[d_Integer] := 
     Block[{a,b,c,list},
           a = 1; list = {};
           While[a <= N[Sqrt[-d/3]],
                 b = 1 - a;
                 While[b <= a, 
                       c= (b^2 - d)/(4 a);
                       If[Mod[c,1] == 0 && c >= a && GCD[a,b,c] == 1 &&
                            Not[a == c && b < 0],
                          AppendTo[list,{a,b,c}]
                         ];
                        b++
                       ];
                   a++
                  ];
             Return[list]
            ]     /;           d < 0 && Mod[d,4] == 1  && SquareFreeQ[d] 

ClassNumber[d_Integer] := Length[ClassList[d]];


End[]   (* NumberTheory`NumberTheoryFunctions`Private` *)

Protect[ClassNumber, ClassList, SquareFreeQ, ChineseRemainderTheorem,
        NextPrime, SqrtMod, PrimitiveRoot, QuadraticRepresentation]


EndPackage[] (* NumberTheory`NumberTheoryFunctions` *)
