home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / NUMBERTH.PAK / BINOMIAL.M next >
Encoding:
Text File  |  1992-07-29  |  6.2 KB  |  213 lines
  1. (*:Version: Mathematica 2.0 *)
  2.  
  3. (*:Name: NumberTheory`Binomial` *)
  4.  
  5. (*:Title: Binomial coefficients *)
  6.  
  7. (*:Author: Ilan Vardi
  8.  
  9. *)
  10.  
  11. (*:Keywords: Binomial coefficients, discrete mathematics, factorials *)
  12.  
  13. (*:Requirements: none. *)
  14.  
  15. (*:Warnings: none *)
  16.  
  17. (* Sources: P. Goetgheluck, ``Prime Divisors of Binomial Coefficients,''
  18.             Mathematics of Computation #51 (1988), p. 325-329.
  19.  
  20.             I. Vardi, ``Computational Recreations in Mathematica,''
  21.             Addison Wesley 1991.
  22. *)
  23.  
  24. (*:Limitations:  FastBinomial is more efficient than the built-in
  25.                  Binomial only for arguments with hundreds of digits.
  26.    
  27.                  FastFactorial is more efficient than the built in
  28.                  Factorial only for numbers of thousands of digits.
  29.  
  30.                  BinomialMod requires a prime modulus
  31. *)
  32.  
  33. (*:Discussion: This package implements fast evaluation of binomial
  34.                coefficients. The evaluation of the standard binomial
  35.                coefficients is done by using their prime factorization
  36.                (see the references above). This reduces the number of 
  37.                multiplications by a factor of log. 
  38.                In the case of Binomial[2 n, n] one has an even better
  39.                algorithm due to its explicit factorization. 
  40.                These algorithms require the function PrimePi[x]. 
  41.  
  42.                The fast evaluation of Mod[Binomial[n, k], p] is 
  43.                done via a well known theorem of E. Lucas (see 
  44.                Computational Recreations, page 70) which reduces 
  45.                the computation to n, k less than p. This is done 
  46.                as before, but with all products in the computation
  47.                being taken mod p. 
  48. *)
  49.  
  50. BeginPackage["NumberTheory`Binomial`"]
  51.  
  52. FastBinomial::usage = "FastBinomial[n, k] returns Binomial[n, k]
  53. efficiently for values of n and k with hundreds of digits."
  54.  
  55. FastFactorial::usage = "FastFactorial[n] return n! efficiently
  56. for n having thousands of digits."
  57.  
  58. BinomialMod::usage = "BinomialMod[n, k, p] returns Mod[Binomial[n, k], p] 
  59. efficiently, where p is a prime number."
  60.  
  61. Begin["`Private`"]  
  62.  
  63. FastBinomial[n_, k_]:= Binomial[n, k] /; k < 500 
  64.  
  65. FastBinomial[n_, k_]:= FastBinomial[n, n - k] /; 2 k > n 
  66.  
  67. FastBinomial[n_, k_]:=
  68. (Times @@ 
  69. (#^(Quotient[n, #] - Quotient[k, #] - 
  70.     Quotient[n - k, #])& /@
  71.  Prime[Range[PrimePi[N[Sqrt[n]]] + 1, 
  72.              PrimePi[n]]])) *
  73.  Times @@ (#^(Plus @@ 
  74.              (Block[{r = #^Range[Log[n]/Log[#]]},
  75.  Quotient[n, r] - Quotient[k, r] - 
  76.  Quotient[n - k, r]]))& /@
  77.  Prime[Range[PrimePi[N[Sqrt[n]]]]])  /; n != 2 k
  78.  
  79. FastBinomial[n_, k_]:= 
  80. Times @@ Prime[Range[PrimePi[k] + 1, 
  81.                      PrimePi[n]]] *
  82. Times @@ Prime[Range[PrimePi[k/2] + 1, 
  83.                      PrimePi[2 k/3]]] *
  84. (Times @@ 
  85.   (#^(Quotient[2 k, #] - 2 Quotient[k , #]) & /@
  86.    Prime[Range[PrimePi[N[Sqrt[2k]]] + 1, 
  87.                PrimePi[2 k/5]]])) *
  88. Times @@ 
  89.    (#^(Plus @@ (Block[{r = #^Range[Log[2 k]/Log[#]]},
  90.        Quotient[2 k, r] - 2 Quotient[k, r]]))& /@
  91.        Prime[Range[PrimePi[N[Sqrt[2 k]]]]]) /; 
  92.                       n == 2 k
  93.  
  94. FastFactorial[0]= 1
  95.  
  96. FastFactorial[n_]:= 
  97. If[OddQ[n], 
  98.    n FastFactorial[n-1], 
  99.    FastBinomial[n, n/2] FastFactorial[n/2]^2]
  100.  
  101. (* Modular calculations: ProductMod gives the product of a list of 
  102.    numbers mod p keeping intermediate results small *)
  103.  
  104. ProductMod[list_, m_]:= Fold[Mod[#1 #2, m]&, 1, list]
  105.  
  106. BinomialMod[n_, k_, p_]:= Mod[Binomial[n, k], p] /; n <= 100
  107.  
  108. Attributes[BinomialMod]= {Listable}
  109.  
  110. BinomialMod[n_, k_, p_]:= 0 /; n < k
  111.  
  112. BinomialMod[n_, k_, p_]:= BinomialMod[n, n - k, p] /; 2 k > n
  113.  
  114. t[n_, k_, p_]:= Block[{tt = 
  115.                       {BinomialMod[n, k , p]//Timing,
  116.                        Mod[Binomial[n, k], p]//Timing}},
  117.                   {tt[[1, 1]], tt[[2, 1]], tt[[1, 2]] == tt[[2,2]]}]
  118.  
  119. (* If n < p then use above algorithms with modular multiplications *)
  120.  
  121. BinomialMod[n_, k_, p_]:= 
  122. Mod[ProductMod[
  123. (#^(Quotient[n, #] - Quotient[k , #] - 
  124.     Quotient[n - k, #])& /@
  125.  Prime[Range[PrimePi[N[Sqrt[n]]] + 1, 
  126.              PrimePi[n]]]), p] *
  127. ProductMod[(#^(Plus @@ 
  128.              (Block[{r = #^Range[Log[n]/Log[#]]},
  129.  Quotient[n, r] - Quotient[k, r] - 
  130.  Quotient[n - k, r]]))& /@
  131.  Prime[Range[PrimePi[N[Sqrt[n]]]]]), p], 
  132.  p]  /; n < p && n != 2 k
  133.  
  134. BinomialMod[n_, k_, p_]:= 
  135. Mod[
  136. ProductMod[Prime[Range[PrimePi[k] + 1, PrimePi[n]]], p] *
  137. ProductMod[Prime[Range[PrimePi[k/2] + 1, PrimePi[2 k/3]]], p] *
  138. ProductMod[
  139.   (#^(Quotient[2 k, #] - 2 Quotient[k , #]) & /@
  140.    Prime[Range[PrimePi[N[Sqrt[2k]]] + 1, 
  141.                PrimePi[2 k/5]]]), p] *
  142. ProductMod[
  143.    (#^(Plus @@ (Block[{r = #^Range[Log[2 k]/Log[#]]},
  144.        Quotient[2 k, r] - 2 Quotient[k, r]]))& /@
  145.        Prime[Range[PrimePi[N[Sqrt[2 k]]]]]), p], 
  146.     p]   /;  n < p && n == 2 k
  147.  
  148.  
  149. (* For n >= p use Lucas' theorem *)
  150.  
  151. BinomialMod[n_, k_, p_]:= 
  152. Block[{kd = IntegerDigits[k, p], nkd = IntegerDigits[n - k, p], 
  153.        nd, min},
  154.        min = Min[Length[kd], Length[nkd]];
  155.        If[Select[Take[nkd, -min] + Take[kd, -min], # >= p&] != {},
  156.           Return[0]];
  157.        nd = IntegerDigits[n, p]; 
  158.        min = Min[Length[nd], Length[kd]]; 
  159.        ProductMod[BinomialMod[Take[nd, -min], Take[kd, -min], p], p]]
  160.        
  161. End[] (* NumberTheory`Binomial` *)   
  162.  
  163. Protect[FastBinomial, FastFactorial, BinomialMod]
  164.  
  165. EndPackage[]   (* NumberTheory`Binomial` *)   
  166.  
  167. (* Tests: 
  168.  
  169. Test[FastBinomial[1000, 501] == Binomial[1000, 501], 
  170.      True
  171. ]
  172.  
  173. Test[FastBinomial[2000, 1000] == Binomial[2000, 1000],
  174.      True
  175. ]
  176.  
  177. Test[FastFactorial[1000] == 1000!, 
  178.      True
  179. ]
  180.  
  181. Test[FastFactorial[5000] == 5000!,
  182.      True
  183. ]                (* Crossover point seems to be 2,500 digits *)
  184.  
  185. Test[BinomialMod[200, 100, 19],
  186.      6
  187. ]
  188.  
  189. Test[BinomialMod[200, 99, 19], 
  190.      5
  191. ]
  192.  
  193. Test[BinomialMod[200, 100, 97],
  194.      40
  195. ]
  196.  
  197. Test[BinomialMod[2 10^4, 10^4, 65537],
  198.      4351
  199. ]
  200.  
  201. Test[BinomialMod[2 10^5, 10^5, Prime[10^8]],
  202.      1194079563
  203. ]                     (* 45 seconds on a DEC3100 *)
  204.  
  205. Test[BinomialMod[2 10^6, 10^6, 10^64 + 57],
  206. 8325550413494218558772064233832167627413775156731422716387786753 
  207. ]                     (*  10 minutes on a DEC3100    *)
  208.  
  209.  
  210. *)
  211.  
  212.  
  213.