(* :Mathematica Version: 6.0 *) (* :Copyright: Copyright 1992-2007, Wolfram Research, Inc. *) (* :Name: NumberTheory`NumberTheoryFunctions` *) (* :Title: Number Theory Functions *) (* :Summary: This package is obsolete; the functionality is loaded by default. *) (* :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. - 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. - Reimplemented PrimitiveRoot as a kernel function - Mark Sofroniou, Wolfram Research, Inc., April. 1999. - Moved ChineseRemainder code and parts of SqrtMod code to StartUp/Algebra/ RNumberTheory.m, so that they can be used by Reduce - Adam Strzebonski, Wolfram Research, Inc., March 2000. - KroneckerSymbol by David Terr, Wolfram Research, Inc., Jan. 2001. - ClassNumber modified by David Terr, Wolfram Research, Inc., Jan. 2001. - Moved SumOfSquaresRepresentations code to StartUp/Algebra/RNumberTheory.m, so that it can be used by Reduce - Adam Strzebonski, Wolfram Research, Inc., January 2001. - SqrtMod modified by David Terr, Wolfram Research, Inc., Feb. 2001. - SqrtModList by David Terr, Wolfram Research, Inc., Feb. 2001. - FundamentalDiscriminantQ by David Terr, Wolfram Research, Inc., June 2001. - PreviousPrime by David Terr, Wolfram Research, Inc., Sept. 2001. - Random[Prime,{a,b}] by David Terr, Wolfram Research, Inc., Oct. 2001. - WhichRootOfUnity[a] by David Terr, Wolfram Research, Inc., Dec. 2001. - AliquotSequence by David Terr, Wolfram Research, Inc., Jan. 2002. - AliquotCycle by David Terr, Wolfram Research, Inc., Jan. 2002. - LeastPrimeFactor by David Terr, Wolfram Research, Inc., Feb. 2002. - PrimePowerQ by David Terr, Wolfram Research, Inc., Feb. 2002. - I removed the restriction that the discriminant be fundamental in ClassList. David Terr, Wolfram Research, Inc., Aug. 2002. - PrimeFactorList by David Terr, Wolfram Research, Inc., Aug. 2002. - OrderedSumOfSquaresRepresentations (former SumOfSquaresRepresentations with order of terms reversed) by David Terr, Wolfram Research, Inc., Jan. 2003 - SumOfSquaresRepresentations modified by David Terr, Wolfram Research, Inc., Jan. 2003 - Removed PrimitiveRoot, which is now in numbertheory.mc - David Terr. Wolfram Research, Inc., Feb. 2003. - Modified ClassNumber to work with positive as well as negative discriminants David Terr, Wolfram Research, Inc., Mar. 2003. - Sped up ClassNumber for positive discriminants by a factor of 3 by using my own regulator. David Terr, Wolfram Research, Inc., Mar. 2003. - The NumberTheoryFunctions package is obsolete; the functionality is loaded by default. Most functions are modified to call kernel functions. Charles Pooh, Wolfram Research Inc., August 2006. *) (* :Keywords: number theory, Chinese Remainder Theorem, class number, primitive roots, Fundamental discriminant, Aliquot sequence *) (* :Requirements: None. *) (* :Warnings: 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. AliquotSequence, AliquotCycle, and SumOfFactors compute various integer sequences based on sums of divisors. The definitions were taken from David Moews' website at http://xraysgi.ims.uconn.edu:8080/amicable.html *) (* :Limitations: SquareFreeQ, SqrtMod, SqrtModList, and QuadraticRepresentation depend on obtaining the factorization of n. For a discussion of when QuadraticRepresentation returns an answer, see Cox's book. 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. QuadraticRepresentation 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. *) (* :Sources: D.A. Buell, Binary Quadratic Forms, Springer-Verlag, 1989. D.A. Cox, Primes of the Form x^2 + n y^2, Wiley, 1989. H. Cohen, A Course in Computational Algebraic Number Theory, Springer-Verlag, 1993. 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. *) Message[General::obspkg, "NumberTheory`NumberTheoryFunctions`"] BeginPackage["NumberTheory`NumberTheoryFunctions`"] AliquotCycle::usage = "AliquotCycle[n] gives the repeating cycle at the end of the Aliquot \ sequence of n." AliquotSequence::usage = "AliquotSequence[n] gives the Aliquot sequence {s_0, s_1, s_2, ..., s_k} \ of n, where s_0 = n and s_(i+1) is the sum of factors of s_i other than s_i. \ The sequence either terminates with s_k = 0, with s_k = s_j for some jAutomatic, TermIncrement->0, ShowProgress->False, MaxIterations->Infinity, MaxTerms->Infinity}; Options[AliquotSequence] = {SumOfFactorsType->Automatic, TermIncrement->0, ShowProgress->False, MaxIterations->Infinity, MaxTerms->Infinity}; Options[SumOfFactors] = {SumOfFactorsType->Automatic}; Begin["`Private`"] issueObsoleteFunMessage[fun_, context_] := (Message[fun::obspkgfn, fun, context]; ) (* Aliquot and related divisor sum sequences *) as[n_, opts___?OptionQ] := Module[{type, increment, printQ, maxlen, maxsize, i, l, m, res = {n}, done = False}, If[ n < 0 || !IntegerQ[n], Message[AliquotSequence::badarg, n]; Return[$Failed]; ]; {type, increment, printQ, maxlen, maxsize} = {SumOfFactorsType, TermIncrement, ShowProgress, MaxIterations, MaxTerms} /. Flatten[{opts,Options[AliquotSequence]}]; If [ !MemberQ[ {Automatic,Unitary,Biunitary,Exponential,ModifiedExponential,Infinitary}, type], Message[SumOfFactors::type, type]; Return[$Failed] ]; If[ !IntegerQ[increment], Message[AliquotSequence::increment, increment]; Return[$Failed] ]; If[ !MemberQ[{True,False},printQ], Message[AliquotSequence::progress, printQ]; Return[$Failed] ]; If[ (!IntegerQ[maxlen] || maxlen < 1) && maxlen =!= Infinity, Message[AliquotSequence::maxit, maxlen]; Return[$Failed] ]; If[ (!IntegerQ[maxsize] || maxsize < 1) && maxsize =!= Infinity, Message[AliquotSequence::maxterms, maxsize]; Return[$Failed] ]; If[ printQ, Print[{0,n}]//TableForm ]; For[i = 1, i <= maxlen && !done, i++, l = If[i > 1, res[[i]], n]; m = SumOfFactors[l,SumOfFactorsType->type] + increment; If[ printQ, Print[{i,m}]//TableForm ]; If[m <= 0 || MemberQ[res, m] || m > maxsize, done = True]; res = Append[res, m]; ]; res ]; AliquotSequence[n_, opts___?OptionQ]:= (issueObsoleteFunMessage[AliquotSequence, "NumberTheory`NumberTheoryFunctions`"]; With[{res = as[n,opts]}, res /; res =!= $Failed ]); NormalizeList[l_List] := Module[{min = Min[l], res = l, leadingTerm = l[[1]]}, While[leadingTerm != min, res = Append[res, leadingTerm]; res = Delete[res, 1]; leadingTerm = res[[1]]; ]; res ]; ac[n_, opts___?OptionQ] := Module[{as = as[n, opts], increment, maxlen, maxsize, k, len, t}, If[ as === $Failed, Return[$Failed]]; {increment, maxlen, maxsize} = {TermIncrement, MaxIterations, MaxTerms} /. Flatten[{opts,Options[AliquotCycle]}]; len = Length[as]; t = as[[len]]; If[ len == maxlen + 1 || t > maxsize, Message[AliquotCycle::open]; Return[$Failed]; ]; If[t == 0 && increment == 0, Return[{0}]]; For[k = len - 1, as[[k]] != t, k--]; NormalizeList[Take[as, k - len]] ]; AliquotCycle[n_Integer, opts___?OptionQ]:= (issueObsoleteFunMessage[AliquotCycle, "NumberTheory`NumberTheoryFunctions`"]; With[{res = ac[n,opts]}, res /; res =!= $Failed ]); (* ClassList *) ClassList[d_Integer /; d < 0 ] := (issueObsoleteFunMessage[ClassList, "NumberTheory`NumberTheoryFunctions`"]; 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] ]) ClassList[n_] := (Message[ClassNumber::baddisc, n]; Null /; False) (* ClassNumber *) ClassNumber[d_Integer /; FundamentalDiscriminantQ[d]] := (issueObsoleteFunMessage[ClassNumber, "NumberTheory`NumberTheoryFunctions`"]; With[{res = NumberTheory`ClassNumber[d]}, res /; FreeQ[res, $Failed | NumberTheory`ClassNumber] ]); ClassNumber[n_] := (Message[ClassNumber::baddisc, n]; Null /; False) (* FundamentalDiscriminantQ *) FundamentalDiscriminantQ[d_Integer]:= (issueObsoleteFunMessage[FundamentalDiscriminantQ, "NumberTheory`NumberTheoryFunctions`"]; With[{res = NumberTheory`FundamentalDiscriminantQ[d]}, res /; FreeQ[res, $Failed | NumberTheory`FundamentalDiscriminantQ] ]) FundamentalDiscriminantQ[___] := False (* LeastPrimeFactor *) lpf[n_] := Module[{i, p, found = False}, If[n <= 1 || ! IntegerQ[n], Message[LeastPrimeFactor::badarg, n]; Return[$Failed]; ]; If[PrimeQ[n], Return[n] ]; For[i = 1, !found, i++, p = Prime[i]; If[Mod[n, p] == 0, found = True; ]; ]; p ]; LeastPrimeFactor[n_] := (issueObsoleteFunMessage[LeastPrimeFactor, "NumberTheory`NumberTheoryFunctions`"]; With[{res = lpf[n]}, res /; res =!= $Failed ]); LeastPrimeFactor[a_, b__] := (Message[LeastPrimeFactor::nargs, Length[{b}] + 1]; Null /; False) LeastPrimeFactor[] := (Message[LeastPrimeFactor::nargs, 0]; Null /; False) (* OrderedSumOfSquaresRepresentations *) OrderedSumOfSquaresRepresentations[d_Integer?Positive, n_Integer] := (issueObsoleteFunMessage[OrderedSumOfSquaresRepresentations, "NumberTheory`NumberTheoryFunctions`"]; With[{res = Sort[#, Greater] & /@ PowersRepresentations[n, d, 2]}, res /; FreeQ[res, $Failed | PowersRepresentations] ]) (* PreviousPrime *) PreviousPrime[x_] := (issueObsoleteFunMessage[PreviousPrime, "NumberTheory`NumberTheoryFunctions`"]; With[{res = -NextPrime[-x]}, res /; FreeQ[res, $Failed | NextPrime] ]) (* PrimeFactorList *) PrimeFactorList[0 | 1] := {} PrimeFactorList[n:(_Integer | _Rational)] := (issueObsoleteFunMessage[PrimeFactorList, "NumberTheory`NumberTheoryFunctions`"]; FactorInteger[n][[All, 1]]) PrimeFactorList[any_] := (Message[PrimeFactorList::badarg, any]; Null /; False) (* QuadraticRepresentation *) qr[d_, n_] := Block[{roots, res, rep, x, y}, If[JacobiSymbol[-d, n] == 1, roots = SqrtModList[-d, n]; res = Catch[Scan[ If[(rep = quadraticRepresentation[d, n, #, n]) =!= $Failed, Throw[rep]]&, roots];$Failed]; If[res =!= $Failed, Return[res]] ]; (* try reduce method *) res = Reduce[{x^2 + d y^2 == n, Element[{x,y},Integers], x >= 0, y >= 0}]; If[res === False, Return[{}]]; If[Head[res] === Reduce, Message[QuadraticRepresentation::norep, n, d]; $Failed, First[{x,y} /. Reduce`ReduceToRules[res]] ] ] 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 ] ] QuadraticRepresentation[d_, n_] /; (n > 0 && OddQ[n] && d > 0) := (issueObsoleteFunMessage[QuadraticRepresentation, "NumberTheory`NumberTheoryFunctions`"]; Module[{res = qr[d, n]}, res /; res =!= $Failed ]) (* RandomPrime *) Random[Prime, range_List] := With[{res = RandomPrime[range]}, res /; FreeQ[res, $Failed | RandomPrime] ] Random[Prime, n_ /; IntegerQ[n] && Positive[n]] := With[{res = RandomPrime[{1, n}]}, res /; FreeQ[res, $Failed | RandomPrime] ] (* SqrtMod and SqrtModList *) sm[n_, m_] := Module[{sml, res}, If[ !IntegerQ[n], Message[PowerMod::arg1, n]; Return[$Failed]; ]; If[ !IntegerQ[m], Message[PowerMod::arg2, m]; Return[$Failed]; ]; If[ m<=0, Message[PowerMod::arg2, m]; Return[$Failed]; ]; sml = Internal`SqrtModList[n, m]; If[ sml == {}, Message[PowerMod::root, 2, n, m]; res = $Failed, res = Internal`SqrtMod[n,m]; ]; res ]; SqrtMod[n_, m_] := (issueObsoleteFunMessage[SqrtMod, "NumberTheory`NumberTheoryFunctions`"]; With[{res = sm[n,m]}, res /; FreeQ[res, $Failed | Internal`SqrtMod] ]) SqrtModList[n_Integer, m_Integer] := (issueObsoleteFunMessage[SqrtModList, "NumberTheory`NumberTheoryFunctions`"]; With[{res = Internal`SqrtModList[n,m]}, res /; FreeQ[res, $Failed | Internal`SqrtModList] ]) (* SumOfFactors *) sof[n_Integer, opts___?OptionQ]:= Module[{type, fi, p, d, e, q, i, len, mult=1, j, res=1}, If[ n == 0, Return[0] ]; If[n < 0 || !(IntegerQ[n] && n >= 0), Message[SumOfFactors::badarg, n]; Return[$Failed]; ]; type = SumOfFactorsType /. Flatten[{opts, Options[SumOfFactors]}]; If[ type == Automatic, Return[ DivisorSigma[1, n] - n ] ]; If [ !MemberQ[ {Unitary,Biunitary,Exponential,ModifiedExponential,Infinitary}, type], Message[SumOfFactors::type, type]; Return[$Failed] ]; fi = FactorInteger[n]; len = Length[fi]; For[ i=1, i<=len, i++, p = fi[[i,1]]; e = fi[[i,2]]; q = p^e; res *= Switch[type, Unitary, q+1, Biunitary, (p*q-1)/(p-1) - If[EvenQ[e],p^(e/2),0], Exponential, 1+p + If[e>1,p^e + Sum[p^d*If[Mod[e,d]==0,1,0],{d,2,e-1}], 0], ModifiedExponential, 1+p^e + Sum[p^d*If[Mod[e,d]==0&&Mod[e+1,d+1]==0,1,0],{d,1,e-1}], Infinitary, 1+p^e + Sum[p^d*If[BitAnd[e,d]==d,1,0],{d,1,e-1}] ]; ]; res - n ] SumOfFactors[n_Integer, opts___?OptionQ]:= (issueObsoleteFunMessage[SumOfFactors, "NumberTheory`NumberTheoryFunctions`"]; With[{res = sof[n, opts]}, res /; FreeQ[res, $Failed] ]) (* SumOfSquaresR *) SumOfSquaresR[d_, n_] := (issueObsoleteFunMessage[SumOfSquaresR, "NumberTheory`NumberTheoryFunctions`"]; With[{res = SquaresR[d, n]}, res /; FreeQ[res, $Failed | SquaresR] ]) (* SumOfSquaresRepresentations *) SignedPermutations[a_List] := Module[{bd, perm0, perm, perms = Permutations[a], i, alen, bdlen, plen, res = {}}, alen = Length[a]; plen = Length[perms]; For[i = 1, i <= plen, i++, perm0 = perms[[i]]; res = Append[res, perm0]; For[j = 1, j < 2^alen, j++, bd = IntegerDigits[j, 2]; bdlen = Length[bd]; bd = Flatten[Prepend[bd,Table[0,{k,alen-bdlen}]]]; perm = perm0; For[k = 1, k <= alen, k++, If[bd[[k]] == 1, perm[[k]] = -perm[[k]]; ]; ]; res = Append[res, perm]; ]; ]; Sort[Union[res]] ] SumOfSquaresRepresentations[d_Integer?Positive, n_Integer] := (issueObsoleteFunMessage[SumOfSquaresRepresentations, "NumberTheory`NumberTheoryFunctions`"]; Module[{len,preres=OrderedSumOfSquaresRepresentations[d,n],sp,splen,res={}}, len = Length[preres]; For[i=1,i<=len,i++, sp = SignedPermutations[preres[[i]]]; splen = Length[sp]; For[j=1, j<=splen, j++, res = Append[res,sp[[j]]]; ]; ]; Sort[Union[res]] ]) (* WhichRootOfUnity *) wrou[x_] := Module[{n,d,k,x1=x,y}, If[!System`Private`AlgebraicNumberQ[x], x1 = Simplify[x]; If[!System`Private`AlgebraicNumberQ[x1], Message[RootOfUnityQ::nalg, x]; Return[$Failed]; ]; ]; If[RootOfUnityQ[x1], If[x1 == 1, Return[{1, 1}]]; y = ToRadicals[RootReduce[x1]]; k = FullSimplify[Log[y]/(2*Pi*I)]; d = Denominator[k]; n = Mod[Numerator[k],d]; Return[{d, n}], Message[WhichRootOfUnity::nonroot, x]; Return[$Failed]; ]; ] WhichRootOfUnity[x_] := (issueObsoleteFunMessage[WhichRootOfUnity, "NumberTheory`NumberTheoryFunctions`"]; With[{res = wrou[x]}, res /; FreeQ[res, $Failed] ]) End[] SetAttributes[ { SumOfSquaresR, OrderedSumOfSquaresRepresentations, SumOfSquaresRepresentations, PreviousPrime }, { Listable, ReadProtected, Protected } ] SetAttributes[ { AliquotCycle, AliquotSequence, Automatic, Biunitary, ClassList, ClassNumber, Exponential, FundamentalDiscriminantQ, Infinitary, LeastPrimeFactor, MaxIterations, MaxTerms, ModifiedExponential, PrimeFactorList, QuadraticRepresentation, Random, ShowProgress, SqrtMod, SqrtModList, SumOfFactors, SumOfFactorsType, TermIncrement, Unitary, WhichRootOfUnity }, { ReadProtected, Protected } ] EndPackage[]