home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / NUMBERTH.PAK / NUMBERTH.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  10.7 KB  |  330 lines
  1.  
  2. (*:Version: Mathematica 2.0 *)
  3.  
  4. (*:Name: NumberTheory`NumberTheoryFunctions` *)
  5.  
  6. (*:Title: Number Theory Functions *)
  7.  
  8. (*:Author: Ilan Vardi
  9.  
  10. *)
  11.  
  12. (*:Keywords:
  13.     number theory, Chinese Remainder Theorem, 
  14.         class number, primitive roots
  15. *)
  16.  
  17. (*:Requirements: none. *)
  18.  
  19. (*:Warnings: SqrtMod needs to have d a residue for every prime 
  20.              factor of n. QuadraticRepresentation needs to have -d
  21.              a quadratic residue only for primes appearing to odd powers.
  22.  
  23. *)
  24.  
  25. (*:Limitations: SquareFreeQ, SqrtMod, and QuadraticRepresentation 
  26.                 depend on obtaining the factorization of n. 
  27.  
  28.                 SqrtMod[d, n] needs to have odd n.
  29.        
  30.                 For a discussion of when QuadraticRepresentation returns
  31.                 and answer, see Cox's book.
  32.  
  33.                 PrimitiveRoot does not need to find the factorization 
  34.                 of n, but of p-1 if n is a power of p, or twice 
  35.                 a prime power of p, p an odd prime.
  36.  
  37.                 ClassList and so ClassNumber has only been implemented
  38.                 for negative discriminants. The implementation uses
  39.                 the simplest algorithm and is slow for large inputs.
  40. *)
  41.  
  42. (*:Sources: 
  43.  
  44.   G.H. Hardy and E.M. Wright, An Introduction to the Theory Of
  45.   Numbers, Oxford University Press, 1988.
  46.  
  47.   D.E. Knuth, Seminumerical Algorithms, Addison Wesley, 1981.
  48.  
  49.   D.A. Buell, Binary Quadratic Forms, Springer Verlag, 1989.
  50.  
  51.   Stan Wagon, Mathematica in Action, Freeman, 1991.
  52.  
  53.   Stan Wagon, The Euclidean Algorithm Strikes Again,
  54.   American Math. Monthly # 97, (1990), 125--124. 
  55.  
  56.   K. Hardy, J.B. Muskat, and K.S. Williams, A determinisitic
  57.   algorithm for solving n = fu^2 + gv^2 in coprime integers
  58.   u and v, Math. Comp. #55 (1990) 327-343.
  59.  
  60.   D.A. Cox, Primes of the Form x^2 + n y^2, Wiley, 1989.  
  61.  
  62. *)
  63.  
  64. (*:Summary: This package implements some standard functions from
  65. number theory.
  66.  
  67. *)
  68.  
  69.  
  70. BeginPackage["NumberTheory`NumberTheoryFunctions`"]
  71.  
  72. Needs["NumberTheory`FactorIntegerECM`"]
  73.  
  74. SquareFreeQ::usage = "SquareFreeQ[n] is True if n is a square free
  75. integer, False otherwise."
  76.  
  77.  
  78. NextPrime::usage = "NextPrime[n] gives the smallest prime p such
  79. that p > n."
  80.  
  81. ChineseRemainderTheorem::usage = "ChineseRemainderTheorem[list1, list2]
  82. gives the minimal non-negative integer solution of Mod[r, list2] == list1."
  83.  
  84. SqrtMod::usage = "SqrtMod[d, n] is the square root of d (mod n), where
  85. n is an odd number."
  86.  
  87. QuadraticRepresentation::usage = "QuadraticRepresentation[d, n] gives
  88. {x, y} where x^2 + d y^2 = n, for odd n and d positive, whenever
  89. such a representation exists. "
  90.  
  91. PrimitiveRoot::usage ="PrimitiveRoot[n] gives a cyclic generator of
  92. the multiplicative group (mod n), when n is a prime power or twice
  93. a prime power."
  94.  
  95. ClassList::usage = "ClassList[d] gives a list of inequivalent quadratic
  96. forms of discriminant d, where d < 0 is a square free integer of the
  97. form 4 n + 1, and  a quadratic form a x^2 + b x y + c y^2 is represented
  98. as {a, b, c}."
  99.  
  100. ClassNumber::usage = "ClassNumber[d] gives the number of inequivalent
  101. quadratic forms of discriminant d, where d < 0 is a square free integer
  102. of the form 4 n + 1."
  103.  
  104.  
  105. Begin["`Private`"]
  106.  
  107.  
  108. (* SquareFreeQ checks for a square factor of d. This may require
  109.    finding the complete factorization of d. The second argument
  110.    ntd is there so as to avoid doing trial division on the factors
  111.    of d as well. *)
  112.  
  113.  
  114. SquareFreeQ[d_Integer, ntd___]:=    
  115. Block[{i = 1, gcd, product = If[ntd === {},Times @@ Prime[Range[1229]]], 
  116.        q = d, flag = False, pr}, 
  117.       If[q < 10^5, Return[MoebiusMu[q] != 0]]; 
  118.       If[ntd === {}, gcd = GCD[q, product]; 
  119.                     q = Quotient[q, gcd];
  120.                     If[GCD[q, gcd] > 1, Return[False]]]; 
  121.       If[q == 1 || PrimeQ[q], Return[True]];
  122.       If[GCD[PowerMod[2, q, q] -2 , q] > 1,  (* Prime Power *)
  123.          Return[MoebiusMu[q] != 0]];  
  124.       If[Select[q^(1/Select[Range[Log[10000, q]], PrimeQ]), IntegerQ] != {},
  125.          Return[False]];                   (* Perfect Power *)
  126.       pr = FactorIntegerECM[q];
  127.       If[GCD[pr, q/pr] > 1, False,
  128.          SquareFreeQ[Min[pr, q/pr], 1] && SquareFreeQ[Max[pr, q/pr], 1]
  129.       ]] 
  130.  
  131. NextPrime[n_]:=    
  132.  Block[{i = n + 1 + Mod[n, 2]},
  133.         If[n < 10^20,   
  134.            While[Not[PrimeQ[i]], i += 2],
  135.             prod = Apply[Times, Prime[Range[Log[n]]]];
  136.             While[True,
  137.                   If[GCD[ i, prod] == 1, If[PrimeQ[i], Break[]]];
  138.                   i += 2]
  139.            ];  
  140.        i]
  141.  
  142.  
  143.  
  144. (* Written by Stan Wagon, see Mathematica in Action.  *)
  145.  
  146. SetAttributes[PowerMod, Listable]
  147.  
  148. ChineseRemainderTheorem[list1_, list2_]:=
  149.  Block[{m = Times @@ list2, mm},
  150.         mm = m/list2;
  151.         Mod[Plus @@ (list1 mm PowerMod[mm, -1, list2]), m]
  152.        ] 
  153.  
  154.  
  155. SqrtMod::nres := "Warning: n has a prime factor p for which d is a 
  156. nonresidue mod p."
  157.  
  158. HalfMod[n_, p_]:= Quotient[If[EvenQ[n], n, p - n], 2]
  159.  
  160. SqrtMod[d_Integer, n_Integer]:= Mod[Sqrt[d], n] /; IntegerQ[Sqrt[d]]
  161.  
  162. SqrtMod[d_Integer, q_Integer] := 
  163.        If[PrimeQ[q], 
  164.  
  165. (* Shank's modular Sqrt algorithm, Knuth, Vol. 2, page 626. Note that
  166.    finding a non residue below takes < 2 (Log [n])^2 time by GRH, and
  167.    this avoids using Random[Integer,{1,p}]. *)
  168.  
  169.        Block[{e, p = q, q2,z,n,v,w,k,t,a,count=0},
  170.               
  171.              e = 1; qq = (p-1)/2;
  172.              While[EvenQ[qq], e++; qq /= 2];
  173.              If[ e == 1, 
  174.                  z = -1,    
  175.                  n = 2;                                
  176.                  While[JacobiSymbol[n,p] != -1, n++];  
  177.                  z = PowerMod[n,qq,p]
  178.                ];
  179.                v = PowerMod[d,(qq+1)/2,p]; 
  180.                w = Mod[PowerMod[d,-1, p] Mod[ v v ,p], p]; 
  181.                While[w != 1, count++; 
  182.                      k = 1; t = Mod[w w,p];
  183. (*  The following loop is very slow for primes of the form 
  184.     h 2^n + 1, large  *)
  185.                      While[t != 1, k++; t = Mod[t t,p]]; 
  186.                      a = PowerMod[z,2^(e-k-1),p];
  187.                      e = k;
  188.                      v = Mod[v a, p];
  189.                      z = Mod[a a,p];
  190.                      w = Mod[w z,p];
  191.                     ]; 
  192.                 Return[v]
  193.               ],    
  194.  
  195.  
  196. (* When q is a prime power p^k = q, SqrtMod[d, q] is computed using 
  197.    a form of Newton's iteration 
  198.  
  199.       x -> (x + d/x)/2
  200.  
  201.    This takes a p^k solution to a p^(2k) solution but requires a
  202.    modular division at each step. 
  203.  
  204.    One can compute this without any modular inversions (except modulo p) 
  205.    by first computing 
  206.  
  207.       x = 1/SqrtMod[d, p] (mod p)
  208.  
  209.    and then iterating 
  210.  
  211.       x -> (3 x - d x^3)/2
  212.  
  213.    which does not require modular inversion. The last step is 
  214.    to multiply by d to get the answer. 
  215.  
  216.    The reason this works is due to the p-adic identity 
  217.  
  218.       Sqrt[d] = d x /Sqrt[1 - (1 - d x^2)]
  219.  
  220.       = d x Sum[Binomial[2 k, k]/4^k (1 - d x^2)^k, {k, 0, Infinity}]
  221.  
  222.    whenever  x^2 = d (mod p^k), some k > 0.
  223.  
  224.    This method is quite efficient in practice and a direct implementation
  225.    computed the square root of -1 mod 5^100000 in 750 seconds on a DEC3100
  226.    (it took twice as long to check the answer). 
  227.  
  228. *)
  229.    
  230.  
  231.         Block[{fi = FactorInteger[q]}, 
  232.                If[Select[First /@ fi, JacobiSymbol[d, #] == -1&] != {},
  233.                   Message[SqrtMod::nres]; Return[]];
  234.           ChineseRemainderTheorem @@  
  235.           Transpose[Map[Function[x, 
  236.               Block[{p, pn, qn, qr, dd, sqrtm, sqrtmin, log},
  237.                   {p, pn} = x;
  238.                    qn = p^pn; 
  239.                    If[pn > 1, log = Ceiling[Log[2, N[pn]]]];
  240.                   sqrtm = SqrtMod[d, p];
  241.                   If[pn == 1, {sqrtm, p}, 
  242.   {If[Mod[d + 1, qn] == 0,  (* d = -1 is special *)
  243.       Mod[Mod[Quotient[#[[1]] (#[[1]]^2 + 3), 2], #[[2]]]& @
  244.       Nest[{Mod[Quotient[#[[1]] (#[[1]]^2 + 3), 2], #[[2]]], #[[2]]^2}&,
  245.        {sqrtm, p^2}, log-1], qn], 
  246. (* else *)    
  247.       Mod[d (HalfMod[Mod[#[[1]] (3 - d #[[1]]^2), #[[2]]], #[[2]]]& @
  248.       Nest[{HalfMod[Mod[#[[1]] (3 - d #[[1]]^2), #[[2]]], #[[2]]], #[[2]]^2}&,
  249.       {PowerMod[sqrtm, -1, p], p^2}, log-1]), qn]],
  250.     qn}]]],
  251.         fi]
  252.            ]]] /; OddQ[q] && JacobiSymbol[d, q] == 1 
  253.  
  254.  
  255. (* Uses generalization of Cornacchia's algorithm, see Stan Wagon's article 
  256.    and Hardy, Muskat, and Williams. See Cox's book for the theory of when
  257.    a representation x^2 + d y^2 = n exists.    *)
  258.  
  259. QuadraticRepresentation::norep := "Warning: n cannot be represented as
  260. x^2 + d y^2 since it splits into factors not in the principal ideal
  261. class of Q(Sqrt[-d])."
  262.  
  263. QuadraticRepresentation[d_, n_]:= 
  264. Block[{x = SqrtMod[-d, n], y = n},
  265.                If[2 x > n, x = n - x];
  266.                While[x^2 >= n, {x, y} = {Mod[y, x], x}];
  267.                y = Sqrt[(n - x^2)/d];
  268.                If[!IntegerQ[y], 
  269.                   Message[QuadraticRepresentation::norep],
  270.                   {x, y}]
  271.              ]  /; OddQ[n] && JacobiSymbol[-d, n] == 1 && d > 0
  272.  
  273.  
  274. (* Written with Ferrell S. Wheeler  *)
  275.  
  276. PrimitiveRoot[2] = 1
  277.  
  278. PrimitiveRoot[4] = 3
  279.  
  280. PrimitiveRoot[n_Integer]:= 
  281. If[OddQ[#], #, n/2 + #]& @ PrimitiveRoot[n/2] /; 
  282.    Mod[n, 4] == 2 && n > 2
  283.  
  284. PrimitiveRoot[p_Integer] :=  
  285.                Block[{g = 2, flist = (p-1)/First /@ FactorInteger[p-1]},
  286.                      While[JacobiSymbol[g, p] == 1 ||
  287.                            Scan[If[PowerMod[g, #, p] == 1, 
  288.                                    Return[True]]&, flist],
  289.                            g++];
  290.                       g]  /; PrimeQ[p] && p > 2
  291.  
  292. PrimitiveRoot[n_Integer]:= 
  293.  Block[{q, fi = FactorInteger[n], g},
  294.         If[Length[fi] > 1, Return[]];
  295.         q = fi[[1,1]];
  296.         g = PrimitiveRoot[q];
  297.         If[PowerMod[g, q - 1, q^2] == 1, g += q];
  298.         g]  /;  GCD[PowerMod[2, n, n]-2, n] > 1
  299.  
  300. (* ClassList only works for negative d, and is inefficient for large d. *)
  301.  
  302. ClassList[d_Integer] := 
  303.      Block[{a,b,c,list},
  304.            a = 1; list = {};
  305.            While[a <= N[Sqrt[-d/3]],
  306.                  b = 1 - a;
  307.                  While[b <= a, 
  308.                        c= (b^2 - d)/(4 a);
  309.                        If[Mod[c,1] == 0 && c >= a && GCD[a,b,c] == 1 &&
  310.                             Not[a == c && b < 0],
  311.                           AppendTo[list,{a,b,c}]
  312.                          ];
  313.                         b++
  314.                        ];
  315.                    a++
  316.                   ];
  317.              Return[list]
  318.             ]     /;           d < 0 && Mod[d,4] == 1  && SquareFreeQ[d] 
  319.  
  320. ClassNumber[d_Integer] := Length[ClassList[d]];
  321.  
  322.  
  323. End[]   (* NumberTheory`NumberTheoryFunctions`Private` *)
  324.  
  325. Protect[ClassNumber, ClassList, SquareFreeQ, ChineseRemainderTheorem,
  326.         NextPrime, SqrtMod, PrimitiveRoot, QuadraticRepresentation]
  327.  
  328.  
  329. EndPackage[] (* NumberTheory`NumberTheoryFunctions` *)
  330.