home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / NUMBERTH.PAK / POLYNOMI.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  9.1 KB  |  279 lines
  1. (*:Version: Mathematica 2.0 *)
  2.  
  3. (*:Name: NumberTheory`PolynomialMod` *)
  4.  
  5. (* :Title:  PolynomialMod *)
  6.  
  7. (* :Author: Ilan Vardi *)
  8.  
  9. (*:Keywords: Polynomial operations modulo a prime > 2^16.
  10.  
  11. *)
  12.  
  13.  
  14. (* :Source: D.E. Knuth, ``Seminumerical Algorithms,'' Second Edition,
  15.             Addison Wesley 1981.
  16.  
  17. *)
  18.  
  19. (* :Warning: This top level code is much slower than the built in 
  20.              functions for primes < 2^16. For p = 2, this code is
  21.              expecially inefficient.
  22.              
  23. *)
  24.  
  25.  
  26. (* :Limitation: FactorMod can factor degree 20 polynomials modulo hundred
  27.                 digit primes in reasonable time. 
  28.                 These algorithms are for polynomials in one varable only.
  29. *)  
  30.  
  31. (* :Discussion:  This package supplants the built-in polynomial functions
  32.                  PolynomialGCD, FactorList, and Factor to primes larger 
  33.                  than 2^16. As well it generalizes PolynomialQuotient and
  34.                  PolynomialRemainder to arithmetic modulo a prime, and 
  35.                  introduces the function PolynomialPowerMod, which computes
  36.                  powers of a polynomial modulo a prime and another polynomial
  37.                  efficiently.
  38.  
  39.     The factoring algorithm is described in the Knuth reference,
  40.     Section 4.6.2: One first uses the derivative to factor the polynomial
  41.     f(x) into squarefree factors. One computes the degree d factorization
  42.     of f by taking GCD[f(x), x^(p^d) - x], for d = 1,2,.... This is 
  43.     computed quickly by using the power mod algorithm, i.e., repeated
  44.     squaring. Finally, one uses the probabilistic factoring algorithm of
  45.     Cantor-Zassenhaus to factor these into irreducible degree d factors.
  46.  
  47. *)
  48.  
  49. BeginPackage["NumberTheory`PolynomialMod`",
  50.     "NumberTheory`NumberTheoryFunctions`"]  
  51.  
  52. PolynomialQuotientMod::usage := "PolynomialQuotientMod[f, g, x, p]
  53. computes PolynomialMod[PolynomialQuotient[f, g, x], p]."
  54.  
  55. PolynomialRemainderMod::usage := "PolynomialRemainderMod[f, g, x, p]
  56. computes PolynomialMod[PolynomialRemainder[f, g, x], p], i.e., 
  57. f modulo g and p, where g is a polynomial and p is a prime."
  58.  
  59. PolynomialGCDMod::usage := "PolynomialGCDMod[f, g, h,..., p] computes
  60. PolynomialGCD[f, g, h,..., Modulus -> p], when p > 2^16."
  61.  
  62. PolynomialPowerMod::usage := "PolynomialPowerMod[f, n, {g, p}]
  63. computes f^n modulo g and p, where g is a polynomial and p is a prime."
  64.  
  65. FactorMod::usage := "FactorMod[f, p] computes Factor[f, Modulus -> p], 
  66. where p > 2^16."
  67.  
  68. FactorListMod::usage := "FactorListMod[f, p] computes 
  69. FactorList[f, Modulus -> p] where p > 2^16."
  70.  
  71. Begin["`Private`"]
  72.                         
  73. PolynomialQuotientMod[f_, g_, x_, p_]:= 
  74. Block[{fp = PolynomialMod[f, p], gp = PolynomialMod[g, p]}, 
  75.        {fp, gp} = 
  76.        PolynomialMod[
  77.         PowerMod[Coefficient[gp, x, Exponent[gp, x]], -1, p] {fp, gp}, 
  78.         p];
  79.        PolynomialMod[PolynomialQuotient[fp, gp, x], p]
  80.       ]
  81.  
  82. PolynomialRemainderMod[f_, g_, x_, p_]:= 
  83.     PolynomialMod[PolynomialRemainder[PolynomialMod[f, p], 
  84.                   PolynomialMod[g, p], x], p]
  85.  
  86. PolynomialGCDMod[f_, p_]:= PolynomialMod[f, p]
  87.  
  88. PolynomialGCDMod[f_, g_, p_]:= 
  89.     polynomialgcdmod[PolynomialMod[f, p], PolynomialMod[g, p],
  90.                      If[Variables[f] != {}, First[Variables[f]],
  91.                                             First[Variables[g]]],
  92.                      p]
  93.  
  94. PolynomialGCDMod[f_, g__, p_]:= 
  95.      PolynomialGCDMod[f, PolynomialGCDMod[g, p], p]
  96.  
  97. polynomialgcdmod[f_, g_, x_, p_Integer]:= 
  98.    polynomialgcdmod[g, f, x, p] /; Exponent[f, x] < Exponent[g, x] 
  99.  
  100. polynomialgcdmod[f_, g_, x_, p_Integer]:= f  /;  Exponent[g,x] == -Infinity
  101.  
  102. polynomialgcdmod[f_, g_, x_, p_Integer] := 
  103.           Block[{fp = PolynomialMod[f, p], gp = PolynomialMod[g, p], 
  104.                  q, r, monic},
  105.                      monic = Coefficient[gp, x, Exponent[gp, x]];
  106.                      If[monic != 1, 
  107.                         gp = PolynomialMod[PowerMod[monic, -1, p] gp, p]];
  108.                      q = PolynomialMod[PolynomialQuotient[fp, gp, x], p];
  109.                      r = PolynomialMod[PolynomialRemainder[fp, gp, x], p];
  110.                      If[Exponent[r, x] <= 0,
  111.                         If[r == 0, Return[gp], Return[1]],
  112.                      Return[polynomialgcdmod[gp, r, x, p]]
  113.                      ]
  114.              ]   /; Exponent[g, x] >= 0  
  115.  
  116.  
  117. PolynomialPowerMod[f_, n_Integer, {g_, p_Integer}]:=
  118. PolynomialPowerMod[f, n, First[Variables[g]], {g, p}]
  119.  
  120. PolynomialPowerMod[f_, n_Integer, x_, {g_, p_Integer}]:=
  121. Fold[PolynomialRemainderMod[#1^2 #2, g, x, p] &,
  122.      1, 
  123.      Block[{prm = PolynomialRemainderMod[f, g, x, p]},
  124.             If[# == 0, 1, prm]& /@ IntegerDigits[n, 2]]] /; n > 0
  125.  
  126.  
  127. FactorMod[f_, p_]:= Times @@ (#[[1]]^#[[2]]& /@ FactorListMod[f, p])
  128.  
  129. FactorListMod[f_, p_]:= 
  130. Block[{flm = factorlistmod[f, p]},
  131.       {#, Count[flm, #]}& /@ Union[flm]]
  132.  
  133. factorlistmod[f_, p_]:= factorlistmod[f, First[Variables[f]], p]
  134.  
  135. factorlistmod[f_, x_, p_]:= 
  136. Block[{fp = PolynomialMod[f, p], monic, factor, df},
  137.        monic = Coefficient[fp, x, Exponent[fp, x]];
  138.        If[monic != 1, 
  139.           Return[
  140.           Join[{monic}, 
  141.                factorlistmod[
  142.                  PolynomialMod[PowerMod[monic, -1, p] fp, p], x, p]]
  143.                 ]];
  144.        df = PolynomialMod[D[fp, x], p]; 
  145.        If[df === 0, 
  146.           Return[Flatten[
  147.                   Table[factorlistmod[f /. x -> x^(1/p), x, p], {p}]]]];
  148.        factor = PolynomialGCDMod[fp, df, p];
  149.        If[Exponent[factor, x] == 0, 
  150.           factorsquarefree[fp, x, 1, x, p], 
  151.           Join[factorlistmod[PolynomialQuotientMod[fp, factor, x, p], x,p],
  152.                factorlistmod[factor, x, p]]]]
  153.  
  154. factorsquarefree[f_, x_, d_, power_, p_]:= {f} /; Exponent[f, x] == d
  155.  
  156. factorsquarefree[f_, x_, d_, xpd_, p_]:= 
  157. Block[{power = PolynomialPowerMod[xpd, p, x, {f, p}], factor, q}, 
  158.        factor = PolynomialGCDMod[f, power - x, p];
  159.        If[Exponent[factor, x] == 0, 
  160.           Return[factorsquarefree[f, x, d + 1, power, p]]];
  161.        If[Degree[factor, x] == Degree[f, x], 
  162.           Return[factordegree[f, x, d, p]]];
  163.        q = PolynomialQuotientMod[f, factor, x, p];
  164.        Join[factorsquarefree[q, x, d + 1,
  165.                              PolynomialRemainderMod[power, q, x, p],
  166.                              p], 
  167.             factordegree[factor, x, d, p]]]
  168.  
  169.      
  170. factordegree[f_, x_, d_, p_]:= {f} /; Exponent[f, x] == d
  171.  
  172. factordegree[f_, x_, d_, p_]:= 
  173.   Block[{exp = Exponent[f, x], fp = PolynomialMod[f, p], 
  174.          s, t, i, q, pd = p^d},
  175.   While[True, 
  176.         s = Plus @@ Table[Random[Integer,{0, p}] x^i, {i,0,exp-2}] + 
  177.             x^(exp - 1); 
  178.         t = PolynomialGCDMod[s, fp, p]; 
  179.        If[0 < Exponent[t,x] < exp, Break[],
  180.          t =  PolynomialGCDMod[fp, 
  181.                 PolynomialPowerMod[s, (pd-1)/2, x, {fp, p}] - 1, p];
  182.             If[ 0 < Exponent[t, x] < exp, Break[]]
  183.         ]]; 
  184.     q = PolynomialQuotientMod[fp, t, x, p];
  185.     t = If[Exponent[q, x] < Exponent[t, x], q, t];
  186.     Join[factordegree[t, x, d, p], 
  187.          factordegree[PolynomialQuotientMod[f, t, x, p], x, d, p]]
  188.    ]
  189.  
  190. Protect[PolynomialQuotientMod, PolynomialRemainderMod, PolynomialGCDMod,
  191.         PolynomialPowerMod, FactorMod, FactorListMod]
  192.  
  193. End[]
  194.  
  195.  
  196. EndPackage[]
  197.  
  198.  
  199. (*:Tests:    (* Timings on DEC3100 *)
  200.  
  201. Test[PolynomialQuotientMod[1 + 3 x^2, 1 + x, x, 3], 
  202.      0
  203. ]
  204.  
  205. Test[PolynomialQuotientMod[1 + 2 x + 4 x^2 + 2 x^5, 1 + x, x, 5]
  206.      2 x + 2 x^2 + 3 x^3 + 2 x^4
  207. ]
  208.  
  209. Test[PolynomialQuotientMod[(1 + x)^3 (1+ x^3), x^2 + 2, x, 11] == 
  210.      PolynomialMod[PolynomialQuotient[(1 + x)^3 (1+ x^3), x^2 + 2, x], 11],
  211.      True
  212. ]
  213.  
  214. Test[PolynomialRemainderMod[1 + x + x^2 + x^3, 2 + x, x, 13],
  215.      8
  216. ]
  217.  
  218. Test[PolynomialRemainderMod[(1 + x)^3 (1+ x^3), x^2 + 2, x, 11] == 
  219.      PolynomialMod[PolynomialRemainder[(1 + x)^3 (1+ x^3), x^2 + 2, x], 11],
  220.      True
  221. ]
  222.  
  223. Test[PolynomialGCDMod[(1 + 14 x) (1 + x + x^3), 1 + x, 13],
  224.      1 + x
  225. ]
  226.  
  227. Test[PolynomialGCDMod[x^3 (1 + x) (1 + x^2)^3, x^2 (1 + x)^2, x^5 (1 + x), 13],
  228.      x^2 + x^3
  229. ]
  230.  
  231. Test[PolynomialGCDMod[x^2 -1, x - 1, Prime[10^8]],
  232.      2038074742 + x
  233. ]
  234.  
  235. Test[PolynomialPowerMod[1 + x, 5, {x^2 + 1, 7}], 
  236.      3 + 3 x
  237. ]
  238.  
  239. Test[PolynomialPowerMod[1 + x, 5, {x^2 + 1, 7}] == 
  240.      PolynomialRemainderMod[(1 + x)^5, x^2 + 1, x, 7], 
  241.      True
  242. ]
  243.  
  244. Test[PolynomialPowerMod[1 + x, 10^4, {x^3 + x^2 + 1, Prime[10^9]}],
  245.      8076973906 + 12198694781 x + 13390038726 x^2  
  246. ]                              (* 1 second *)
  247.  
  248.  
  249. Test[FactorMod[1 + 2 x^3, 3], 
  250.      2 (2 + x)^3
  251. ]
  252.  
  253. Test[FactorMod[Expand[(x + 1)^9 (x^4 + x + 1)], 3],
  254.      (1 + x)^9 (2 + x) (2 + x + x^2 + x^3)
  255. ]
  256.  
  257. Test[FactorMod[Expand[(x + 1)^9 (x^4 + x + 1)], 3] == 
  258.      Factor[Expand[(x + 1)^9 (x^4 + x + 1)], Modulus -> 3],
  259.      True
  260. ]
  261.  
  262. Test[FactorMod[Expand[(x+1)^2 (x^2 + 5)^2 (x^2 + 13)], 10^64 + 57],
  263.      (1 + x)^2 (5 + x^2)^2 (13 + x^3)
  264. ]                                      (* 100 seconds *)
  265.  
  266. Test[FactorListMod[1 + 2 x^3, 3],
  267.      {{2, 1}, {2 + x, 3}}
  268. ]
  269.  
  270. Test[FactorListMod[Expand[(x + 1)^9 (x^4 + x + 1)], 3] == 
  271.      FactorList[Expand[(x + 1)^9 (x^4 + x + 1)], Modulus -> 3],
  272.      True
  273.  
  274. *)
  275.  
  276.  
  277.  
  278.  
  279.