/cln-1.3.2/src/integer/misc/combin/cl_I_factorial.cc

  1// factorial().
  2
  3// General includes.
  4#include "base/cl_sysdep.h"
  5
  6// Specification.
  7#include "cln/integer.h"
  8
  9
 10// Implementation.
 11
 12#include "integer/cl_I.h"
 13#include "integer/misc/combin/cl_I_combin.h"
 14
 15namespace cln {
 16
 17  // Methode:
 18  // n <= 10 -> Ergebnis (Fixnum) aus Tabelle
 19  // Sonst:
 20  //   Zweierpotenzen extra am Schluß durch einen Shift um
 21  //   ord2(n!) = sum(k>=1, floor(n/2^k) ) = n - logcount(n)  Bits.
 22  //   Für k>=1 wird jede ungerade Zahl m im Intervall n/2^k < m <= n/2^(k-1)
 23  //   genau k mal gebraucht (als ungerader Anteil von m*2^0,...,m*2^(k-1) ).
 24  //   Zur Bestimmung des Produkts aller ungeraden Zahlen in einem Intervall
 25  //   a < m <= b verwenden wir eine rekursive Funktion, die nach Divide-and-
 26  //   Conquer das Produkt über die Intervalle a < m <= c und c < m <= b
 27  //   (c := floor((a+b)/2)) bestimmt und beide zusammenmultipliziert. Dies
 28  //   vermeidet, daß oft große Zahlen mit ganz kleinen Zahlen multipliziert
 29  //   werden.
 30
 31const cl_I factorial (uintL n) // assume n >= 0 small
 32{
 33	static uintV const fakul_table [] = {
 34        1,
 35        1UL,
 36        1UL*2,
 37        #if (cl_value_len>=4)
 38        1UL*2*3,
 39        #if (cl_value_len>=6)
 40        1UL*2*3*4,
 41        #if (cl_value_len>=8)
 42        1UL*2*3*4*5,
 43        #if (cl_value_len>=11)
 44        1UL*2*3*4*5*6,
 45        #if (cl_value_len>=14)
 46        1UL*2*3*4*5*6*7,
 47        #if (cl_value_len>=17)
 48        1UL*2*3*4*5*6*7*8,
 49        #if (cl_value_len>=20)
 50        1UL*2*3*4*5*6*7*8*9,
 51        #if (cl_value_len>=23)
 52        1UL*2*3*4*5*6*7*8*9*10,
 53        #if (cl_value_len>=27)
 54        1UL*2*3*4*5*6*7*8*9*10*11,
 55        #if (cl_value_len>=30)
 56        1UL*2*3*4*5*6*7*8*9*10*11*12,
 57        #if (cl_value_len>=34)
 58        1UL*2*3*4*5*6*7*8*9*10*11*12*13,
 59        #if (cl_value_len>=38)
 60        1UL*2*3*4*5*6*7*8*9*10*11*12*13*14,
 61        #if (cl_value_len>=42)
 62        1UL*2*3*4*5*6*7*8*9*10*11*12*13*14*15,
 63        #if (cl_value_len>=46)
 64        1UL*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16,
 65        #if (cl_value_len>=50)
 66        1UL*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17,
 67        #if (cl_value_len>=54)
 68        1UL*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18,
 69        #if (cl_value_len>=58)
 70        1UL*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19,
 71        #if (cl_value_len>=63)
 72        1UL*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19*20,
 73        #if (cl_value_len>=67)
 74        ...
 75        #endif
 76        #endif
 77        #endif
 78        #endif
 79        #endif
 80        #endif
 81        #endif
 82        #endif
 83        #endif
 84        #endif
 85        #endif
 86        #endif
 87        #endif
 88        #endif
 89        #endif
 90        #endif
 91        #endif
 92        #endif
 93        #endif
 94	};
 95
 96      if (n < sizeof(fakul_table)/sizeof(cl_I))
 97        { return UV_to_I(fakul_table[n]); }
 98        else
 99        { var cl_I prod = 1; // bisheriges Produkt := 1
100          var uintL k = 1;
101          var uintL A = n;
102          var uintL B = n; // obere Intervallgrenze floor(n/2^(k-1))
103          loop
104            { // 'A' enthält floor(n/2^(k-1)).
105              A = A >> 1; // untere Grenze floor(n/2^k)
106              // 'A' enthält floor(n/2^k).
107              // Bilde Teilprodukt prod(A < i <= B & oddp(i), i)
108              //       = prod(floor((A-1)/2) < i <= floor((B-1)/2), 2*i+1)
109              // wobei B = floor(n/2^(k-1)), A = floor(n/2^k) = floor(B/2).
110              { var uintL b = floor(B-1,2);
111                if (b==0) break; // B=2 oder B=1 -> Produkt fertig
112                var uintL a = floor(A-1,2);
113                prod = expt_pos(cl_I_prod_ungerade(a,b),k) * prod; // aufmultiplizieren
114              }
115              k = k+1;
116              B = A;
117            }
118          return prod << (n - logcount(n));
119        }
120}
121// Bit complexity (N := n): O(log(N)^2*M(N)).
122
123}  // namespace cln
124