1/*
2 * factorial - implementation of different algorithms for the factorial
3 *
4 * Copyright (C) 2013 Christoph Zurnieden
5 *
6 * Calc is open software; you can redistribute it and/or modify it under
7 * the terms of the version 2.1 of the GNU Lesser General Public License
8 * as published by the Free Software Foundation.
9 *
10 * Calc is distributed in the hope that it will be useful, but WITHOUT
11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12 * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
13 * Public License for more details.
14 *
15 * A copy of version 2.1 of the GNU Lesser General Public License is
16 * distributed with calc under the filename COPYING-LGPL.  You should have
17 * received a copy with calc; if not, write to Free Software Foundation, Inc.
18 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
19 *
20 * Under source code control:	2013/08/11 01:31:28
21 * File existed as early as:	2013
22 */
23
24
25/*
26 * hide internal function from resource debugging
27 */
28static resource_debug_level;
29resource_debug_level = config("resource_debug", 0);
30
31
32/*
33  get dependencies
34*/
35read -once toomcook;
36
37
38/* A simple list to keep things...uhm...simple?*/
39static __CZ__primelist = list();
40
41/* Helper for primorial: fill list with primes in range a,b */
42define __CZ__fill_prime_list(a,b)
43{
44  local k;
45  k=a;
46  if(isprime(k))k--;
47  while(1){
48    k = nextprime(k);
49    if(k > b) break;
50    append(__CZ__primelist,k );
51  }
52}
53
54/* Helper for factorial: how often prime p divides the factorial of n */
55define __CZ__prime_divisors(n,p)
56{
57    local q,m;
58    q = n;
59    m = 0;
60    if (p > n) return 0;
61    if (p > n/2) return 1;
62    while (q >= p) {
63        q  = q//p;
64        m += q;
65    }
66    return m;
67}
68
69/*
70  Wrapper. Please set cut-offs to own taste and hardware.
71*/
72define factorial(n){
73  local prime result shift prime_list k k1 k2 expo_list pix cut primorial;
74
75  result = 1;
76  prime  = 2;
77
78  if(!isint(n))  {
79    return newerror("factorial(n): n is not an integer"); ## or gamma(n)?
80  }
81  if(n < 0)      return newerror("factorial(n): n < 0");
82  if(n < 9000 && !isdefined("test8900"))   {
83    ## builtin is implemented with splitting but only with
84    ## Toom-Cook 2 (by Karatsuba (the father))
85    return n!;
86  }
87
88  shift = __CZ__prime_divisors(n,prime);
89  prime = 3;
90  cut = n//2;
91  pix   = pix(cut);
92  prime_list = mat[pix];
93  expo_list  = mat[pix];
94
95  k = 0;
96/*
97  Peter Borwein's algorithm
98
99  @Article{journals/jal/Borwein85,
100  author = {Borwein, Peter B.},
101  title  = {On the Complexity of Calculating Factorials.},
102  journal = {J. Algorithms},
103  year   = {1985},
104  number = {3},
105  url    = {http://dblp.uni-trier.de/db/journals/jal/jal6.html#Borwein85}
106*/
107
108  do {
109    prime_list[k]  = prime;
110    expo_list[k++] = __CZ__prime_divisors(n,prime);
111    prime          = nextprime(prime);
112  }while(prime <= cut);
113
114  /* size of the largest exponent in bits  */
115  k1 = highbit(expo_list[0]);
116  k2 = size(prime_list)-1;
117
118  for(;k1>=0;k1--){
119    /*
120       the cut-off for T-C-4 ist still to low, using T-C-3 here
121       TODO: check cutoffs
122     */
123    result = toomcook3square(result);
124    /*
125      almost all time is spend in this loop, so cutting of the
126      upper half of the primes makes sense
127    */
128    for(k=0; k<=k2; k++) {
129      if((expo_list[k] & (1 << k1)) != 0) {
130        result *= prime_list[k];
131      }
132    }
133
134  }
135  primorial = primorial( cut, n);
136  result *= primorial;
137  result <<= shift;
138  return result;
139}
140
141/*
142  Helper for primorial: do the product with binary splitting
143  TODO: do it without the intermediate list
144*/
145define __CZ__primorial__lowlevel( a, b ,p)
146{
147  local  c;
148  if( b == a) return p ;
149  if( b-a > 1){
150    c= (b + a) >> 1;
151    return  __CZ__primorial__lowlevel( a   , c , __CZ__primelist[a] )
152         *  __CZ__primorial__lowlevel( c+1 , b , __CZ__primelist[b] ) ;
153  }
154  return  __CZ__primelist[a] * __CZ__primelist[b];
155}
156
157/*
158   Primorial, Product of consecutive primes in range a,b
159
160  Originally meant to do primorials with a start different from 2, but
161  found out that this is faster at about a=1,b>=10^5 than the builtin
162  function pfact().  With the moderately small list a=1,b=10^6 (78498
163  primes) it is 3 times faster.  A quick look-up showed what was already
164  guessed: pfact() does it linearly.  (BTW: what is the time complexity
165  of the primorial with the naive algorithm?)
166*/
167define primorial(a,b)
168{
169  local C1 C2;
170  if(!isint(a)) return newerror("primorial(a,b): a is not an integer");
171  else if(!isint(b)) return newerror("primorial(a,b): b is not an integer");
172  else if(a < 0) return newerror("primorial(a,b): a < 0");
173  else if( b < 2 ) return newerror("primorial(a,b): b < 2");
174  else if( b < a) return newerror("primorial(a,b): b < a");
175  else{
176    /* last prime < 2^32 is also  max. prime for nextprime()*/
177    if(b >= 4294967291) return newerror("primorial(a,b): max. prime exceeded");
178    if(b ==  2) return 2;
179    /*
180      Can be extended by way of pfact(b)/pfact(floor(a-1/2)) for small a
181     */
182    if(a<=2 && b < 10^5) return pfact(b);
183    /* TODO: use pix() and a simple array (mat[])instead*/
184     __CZ__primelist = list();
185     __CZ__fill_prime_list(a,b);
186    C1 = size(__CZ__primelist)-1;
187    return __CZ__primorial__lowlevel( 0, C1,1)
188  }
189}
190
191
192/*
193 * restore internal function from resource debugging
194 * report important interface functions
195 */
196config("resource_debug", resource_debug_level),;
197if (config("resource_debug") & 3) {
198    print "factorial(n)";
199    print "primorial(a, b)";
200}
201