(*:Mathematica Version: 4.0 *)

(*:Package Version: 1.4 *)

(*:Copyright: Copyright 1992-1999, Wolfram Research, Inc. *)

(*:Name: NumberTheory`NumberTheoryFunctions` *)

(*:Title: Number Theory Functions *)

(*:Author: Ilan Vardi *)

(*:History:
Original package by Ilan Vardi, Wolfram Research, Inc., 1992.
Support added for cases where PrimeQ returns True for composites,
	ECM, Wolfram Research, Inc., 1993.
SumOfSquaresR and SumOfSquaresRepresentations by Stan Wagon,
	Macalester College, 1995.
PrimitiveRoot by Mark Sofroniou,
	Wolfram Research, Inc., Sept. 1996.
Enhancements to QuadraticRepresentation, ECM,
	Wolfram Research, Inc., March 1997.
ChineseRemainderTheorem by Stan Wagon, Macalester College, and
        Daniel Lichtblau, Wolfram Research, Inc., Jan. 1998.  Code is
        based in part on earlier code by Mike McGeachie and Craig Ortner,
        students at Macalester College.
Renamed ChineseRemainderTheorem to be ChineseRemainder, retaining support
	for older symbol, ECM, WRI, Feb. 1998.
*)

(*: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.

For SumOfSquaresR[d, n], the cases d = 2, 4, 6, 8 can handle large
integer values of n, so long as n can be factored. The
other cases use recursion, calling (d-1, m) where m takes
on Sqrt[n] many values. Thus only modest size values of
m and d can be used.
*)

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

SqrtMod[d, n] needs to have an odd n.
       
For a discussion of when QuadraticRepresentation returns
an 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 also ClassNumber have only been implemented
for negative discriminants. The implementation uses
the simplest algorithm and is slow for large inputs.

For SumOfSquaresR[d, n], d = 2, 4, 6, or 8, some simple formulas are called,
so n can be quite large (its factorization will be used).  For other
values of d a recursion is used that calls r[d-1, n] approximately Sqrt[n]
times, so n must not be too large.
*)

(*:Sources: 

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

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

Emil Grosswald, Representations of Integers as
  Sums of Squares, Springer, New York, 1985.

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

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.E. Knuth, Seminumerical Algorithms, Addison-Wesley, 1981.

Stan Wagon, Mathematica in Action, Freeman, 1991.

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

Stan Wagon, Complex Factoring, Mathematica in Education and Research, 1996,
  to appear.
*)

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


BeginPackage["NumberTheory`NumberTheoryFunctions`",
	"NumberTheory`FactorIntegerECM`"]

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

ChineseRemainderTheorem::usage =
"ChineseRemainderTheorem[list1, list2] gives the minimal
nonnegative integer solution of Mod[r, list2] == list1.
The solution is unique modulo the LCM of list2.
ChineseRemainderTheorem[list1, list2] is an obsolete function,
superceded by ChineseRemainder[list1, list2].";

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

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

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

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

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

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

SumOfSquaresR::usage = "SumOfSquaresR[d, n] gives the number of
representations of an integer n as a sum of d squares; standard notation for
this is r[d, n]."

SumOfSquaresRepresentations::usage = "SumOfSquaresRepresentations[d, n] gives
the list of all representations of the positive integer n as a sum of d
squares, ignoring order and signs."


Unprotect[ChineseRemainder, ChineseRemainderTheorem,
	ClassNumber, ClassList,
	NextPrime, PrimitiveRoot, QuadraticRepresentation, SquareFreeQ, SqrtMod]

Begin["`Private`"]


(* SquareFreeQ checks for a square factor of d. This may require
   finding the complete factorization of d.
*)


SquareFreeQ[d_Integer]:=    
Block[{q = d, gcd, product = Times @@ Prime[Range[1229]], pr},
      If[q < 10^5, Return[MoebiusMu[q] != 0]];
      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[-3] := -2
NextPrime[-2] := 2 
NextPrime[-1] := 2 
NextPrime[0] := 2 
NextPrime[1] := 2

NextPrime[r_Real] := NextPrime[Floor[r]] 

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


(* Chinese Remainder code is by Stan Wagon, Macalester College, and
Daniel Lichtblau, Wolfram Research, Inc. It is based in part on
earlier code by Mike McGeachie and Craig Ortner (students at
Macalester College). January 1998. *)



ChineseRemainder[a_List, m_List] := If[
  Max[m] < 3000000 && LCM @@ m == Times @@ m,
    CRTStandard[a, m],
    Catch[First[CRTlist[a, m, 1, Length[a]]]]] /;
        Length[a] == Length[m]

ChineseRemainderTheorem[a_List, m_List] :=
  (
  If[
  Max[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}]]


SqrtMod::nres := "Warning: `1` has a prime factor `2` for which `3` is a 
			nonresidue mod `2`."
SqrtMod::fail :=
"Warning: $IterationLimit exceeded; PrimeQ[`1`] may be incorrect."

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, qq, z, n, v, w, k, t, a, itr = 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,
                     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,
		       itr++;
		       If[itr > $IterationLimit,
                         Message[SqrtMod::fail, p];
			 Return[$Failed]
		       ];
		       k++;
		       t = Mod[t t, p];
		     ];
		     If[e-k-1 < 0, Return[$Failed]];
                     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]
       ], (* end prime q Block *)


(* 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], mmm},
	  (* q is not Prime *)
          If[Select[First /@ fi,
		 (mmm = #; JacobiSymbol[d, mmm] == -1)&] =!= {},
                  Message[SqrtMod::nres, q, mmm, d]; Return[$Failed]];
          ChineseRemainderTheorem @@  
            Transpose[
	     Map[Function[x, 
              Block[{p, pn, qn, sqrtm, 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}
		    ] (* end If pn == 1 *)
              ] (* end Block *) ], (* end Function *)
              fi] (* end Map *)
            ] (* end Transpose *)
       ] (* end composite q Block *)
   ] /; OddQ[q] && q != 1 && JacobiSymbol[d, q] == 1 
	(* Require that q != 1 because 1 is neither prime nor composite. *)


(* 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: `1` splits into factors not in the principal ideal class of
Q(Sqrt[-`2`]), so a representation of `1` as x^2 + `2` y^2 may be missed."

QuadraticRepresentation[d_, n_]:= 
Block[{x, y, result, root1, quadrep, roots, scan},
  (
	result
  ) /; ( root1 = SqrtMod[-d, n];
         {x, y} = {root1, n};
         (* first try to find representation using root provided by SqrtMod *)
	 quadrep = quadraticRepresentation[d, n, x, y];
         If[quadrep === $Failed, 
	    If[FreeQ[roots = Roots[z^2 == -d, z, Modulus -> n], Roots] &&
			Head[roots] === Or,
		roots = Union[  Apply[List, Map[Sort[{#[[2]], n-#[[2]]}]&,
			 roots]]  ];
		roots = Select[Map[#[[1]]&, roots],
			 (# =!= root1 && # =!= (n-root1))&];
	        scan = Scan[ (
		              {x, y} = {#, n};
	        	      (* try to find representation using a root
					provided by Roots *)
			      quadrep = quadraticRepresentation[d, n, x, y];
			      If[quadrep =!= $Failed, Return[OK]]
			     )&, roots];
		If[scan =!= OK,
                   Message[QuadraticRepresentation::norep, n, d]; False,
	           result = quadrep;  True
		],
		(* ! (FreeQ[roots, Roots] && Head[roots] === Or) *)
		Message[QuadraticRepresentation::norep, n, d]; False
	    ],
	    (* quadrep =!= $Failed *)
	    result = quadrep;  True
	 ] (* end If quadrep === $Failed *)
       )
]  /; n > 0 && OddQ[n] && d > 0 && JacobiSymbol[-d, n] == 1 

quadraticRepresentation[d_, n_, x_, y_] :=
  Module[{xx = x, yy = y},
         If[2 xx > n, xx = n - xx];
         While[xx^2 >= n, {xx, yy} = {Mod[yy, xx], xx}];
         yy = Sqrt[(n - xx^2)/d];
	 If[IntegerQ[yy],
	    {xx, yy},
	    $Failed
	 ]
  ]


(* The following method suggested by Dan Lichtblau. *)
QuadraticRepresentation[d_, n_]:=
Module[{g, sqrt, quadrep, x, y, xx, yy},
  (
  {xx, yy}
  ) /; (g = GCD[d, n];  IntegerQ[g] && g =!= 1) &&
       (sqrt = Sqrt[d/g];
	IntegerQ[sqrt]) &&
       FreeQ[quadrep = QuadraticRepresentation[g, n/g], 
		QuadraticRepresentation] && 
       ({x, y} = quadrep;  xx = g y;
	IntegerQ[yy = x/sqrt])
] /; n > 0 && OddQ[n] && d > 0


(* Original version by Ilan Vardi and Ferrell S. Wheeler failed to
   return unevaluated for composites. *)

PrimitiveRoot[p_Integer?Positive] := 
  Module[{res},
    res /; NumberQ[ res = pr[p] ]
  ];

pr[2] = 1;
pr[3] = 2;
pr[4] = 3;

(* filter out case of m=2*2^a *)
pr[n_]:=$Failed/;Mod[n,4]==0

(* handle case m=2*p^a, p is odd prime *)
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;

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

ClassList::bad =
"Argument `1` is not a negative square-free integer of the form 4 n + 1.";

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] 

ClassList[n_] := Null/;(Message[ClassList::bad, n]; False)

ClassNumber[d_Integer] :=
	With[{ans = ClassList[d]},
		Length[ans] /; ListQ[ans]
	]


(* Written by Stan Wagon, see Complex Factoring,
	Mathematica in Education and Research, 1996  *)

SetAttributes[{SumOfSquaresR, SumOfSquaresRepresentations}, Listable]

SumOfSquaresRepresentations[1, n_?Positive] :=
 If[IntegerQ[Sqrt[n]], {Sqrt[n]}, {}]

SumOfSquaresRepresentations[2, n_?Positive] :=
 SumOfSquaresRepresentations[2, n] = Union[
  Map[Sort[{Re[#], Im[#]}] &, Select[
    Divisors[n, GaussianIntegers -> True],
    Abs[#]^2 == n &]]]

SumOfSquaresRepresentations[d_?Positive, n_?Positive] :=
 (SumOfSquaresRepresentations[d, n] =

  Union[Map[Sort, Join[
   Append[#, 0] & /@ SumOfSquaresRepresentations[d - 1, n],
   If[IntegerQ[Sqrt[n]],
    {Append[Array[0&, d-1], Sqrt[n]]}, {}],
   Flatten[Table[Append[#, z]& /@
   SumOfSquaresRepresentations[d - 1, n - z^2],
                  {z, Sqrt[n - 1]}], 1]] ]])

SumOfSquaresRepresentations[d_?Positive, n_?Negative] := {}

SumOfSquaresRepresentations[d_?Positive, 0] := {Array[0&, d]}

FourFreePart[n_] :=
  (m = n; While[Mod[m, 4] == 0, m /= 4]; m)

SumOddDivs[n_] := Plus @@ Select[Divisors[n], OddQ]

OddPart[n_] := Module[{m = n}, While[EvenQ[m], m /= 2]; m]

(* DivisorSigma utilities, that consider only primes 1 mod 4 or
   3 mod 4 but are otherwise like DivisorSigma; also
   AlternatingDivisorSigma[k,n] gives the sum of (-1)^d
   d^k over the divisors d of n; all from Stan Wagon *)
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}]]

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])

DivisorSigma3[k_, n_] := DivisorSigma[k, OddPart[n]] - DivisorSigma1[k, 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 definitions *)

SumOfSquaresR[_, _?Negative] := 0;
SumOfSquaresR[_, 0] := 1;
SumOfSquaresR[0, _] := 0;

(* Following definition for SumOfSquaresR[2, _] is due to
   Stan Wagon. *)
SumOfSquaresR[2, n_Integer?Positive] :=
    SumOfSquaresR[2, n] =
        If[MemberQ[fi = FactorInteger[n],
                {p_ /; Mod[p,4] == 3, _?OddQ}],
            0,
            4 Apply[Times, Cases[fi,{p_ /; Mod[p,4] == 1, a_} :> a+1 ]]
        ]

SumOfSquaresR[3, n_] := (SumOfSquaresR[3, n] = 0)/; Mod[FourFree[n], 8] == 7;
SumOfSquaresR[3, n_] := (SumOfSquaresR[3, n] = SumOfSquaresR[3, n/4])/;
       Mod[n, 4] == 0;

SumOfSquaresR[3, n_] := (SumOfSquaresR[3, n] = 
Module[{n1 = OddPart[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] := (rprim[n] =
  24 Sum[JacobiSymbol[s,n], {s, n/4}]) /; Mod[n, 4] == 1

rprim[n_?OddQ] := (rprim[n] =
  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
restrcit this case to small n *)

SumOfSquaresR[3, n_?OddQ] := (SumOfSquaresR[3,n] =
  Plus @@ (rprim /@ (n / Select[Divisors[n], IntegerQ[Sqrt[#]]&]))) /;
     SquareFreeQ[n] && n < 500

SumOfSquaresR[3, n_?EvenQ] := SumOfSquaresR[3, n] =
    SumOfSquaresR[2, n] + 2 *
        Apply[Plus, SumOfSquaresR[2, n - Range[Sqrt[N[n]]]^2]]

SumOfSquaresR[4, n_Integer?Positive] := SumOfSquaresR[4,n] =
  If[EvenQ[n], 24, 8] * DivisorSigma[1, OddPart[n]]

SumOfSquaresR[5, n_?((IntegerQ[Sqrt[#]] && (# > 1) &&
                        !MatchQ[FactorInteger[#], {{2,_}}])&)] :=
  SumOfSquares[5, n] =
    Module[{nn, n1, p, b},
        nn = Sqrt[n];
        n1 = OddPart[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))
    ]

SumOfSquaresR[6, n_Integer?Positive] := SumOfSquares[6, n] =
  16*n^2*(DivisorSigma1[-2, n] - DivisorSigma3[-2, n]) -
    4*(DivisorSigma1[2, n] - DivisorSigma3[2, n])

SumOfSquaresR[8, n_Integer?Positive] := SumOfSquaresR[8,n] =
  16 (-1)^n AlternatingDivisorSigma[3, n]

SumOfSquaresR[d_, n_Integer?Positive] :=
    SumOfSquaresR[d, n] =
      SumOfSquaresR[d-1, n] + 2 *
      Apply[Plus, SumOfSquaresR[d-1, n - Range[Sqrt[N[n]]]^2]]

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

SetAttributes[ {PrimitiveRoot}, Listable];

SetAttributes[ {ChineseRemainder, ChineseRemainderTheorem,
	ClassNumber, ClassList, NextPrime, PrimitiveRoot,
	QuadraticRepresentation, SquareFreeQ, SqrtMod},
ReadProtected ]

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


EndPackage[] (* NumberTheory`NumberTheoryFunctions` *)

(* :Examples:

SumOfSquaresRepresentations[2, 100]

SumOfSquaresRepresentations[3, 100]

Apply[Plus, (%^2), 2]	(* check *)

SumOfSquaresRepresentations[5, 100]

SumOfSquaresRepresentations[{1, 2, 3, 4, 5}, 64]

SumOfSquaresR[5, 100]

SumOfSquaresR[2, 10^30]

SumOfSquaresR[4, 10^30]

SumOfSquaresR[6, 10^30]

SumOfSquaresR[8, 10^30]

SumOfSquaresR[Range[10], 100]

(* The asymptotic average value of r2 is Pi *)
Sum[N[SumOfSquaresR[2, k]], {k, 200}] / 200

(* A curious alternating sum equals - Pi Log[2] *)
Sum[N[(-1)^n SumOfSquaresR[2, n] / n], {n, 1, 200}]

- Pi Log[2] //N

QuadraticRepresentation[14,65]
	{3, 2}

QuadraticRepresentation[12,21]
	{3, 1}

QuadraticRepresentation[2, 9]
	{1, 2}

QuadraticRepresentation[3, 49]
	{1, 4}

QuadraticRepresentation[15, 991]
	{16, 7}

QuadraticRepresentation[2, 81]
	{7, 4}

QuadraticRepresentation[49, 277]
	{9, 2}

QuadraticRepresentation[48, 147]
	unevaluated


*)

