/* * factorial - implementation of different algorithms for the factorial * * Copyright (C) 2013 Christoph Zurnieden * * Calc is open software; you can redistribute it and/or modify it under * the terms of the version 2.1 of the GNU Lesser General Public License * as published by the Free Software Foundation. * * Calc is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General * Public License for more details. * * A copy of version 2.1 of the GNU Lesser General Public License is * distributed with calc under the filename COPYING-LGPL. You should have * received a copy with calc; if not, write to Free Software Foundation, Inc. * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. * * Under source code control: 2013/08/11 01:31:28 * File existed as early as: 2013 */ /* * hide internal function from resource debugging */ static resource_debug_level; resource_debug_level = config("resource_debug", 0); /* get dependencies */ read -once toomcook; /* A simple list to keep things...uhm...simple?*/ static __CZ__primelist = list(); /* Helper for primorial: fill list with primes in range a,b */ define __CZ__fill_prime_list(a,b) { local k; k=a; if(isprime(k))k--; while(1){ k = nextprime(k); if(k > b) break; append(__CZ__primelist,k ); } } /* Helper for factorial: how often prime p divides the factorial of n */ define __CZ__prime_divisors(n,p) { local q,m; q = n; m = 0; if (p > n) return 0; if (p > n/2) return 1; while (q >= p) { q = q//p; m += q; } return m; } /* Wrapper. Please set cut-offs to own taste and hardware. */ define factorial(n){ local prime result shift prime_list k k1 k2 expo_list pix cut primorial; result = 1; prime = 2; if(!isint(n)) { return newerror("factorial(n): n is not an integer"); ## or gamma(n)? } if(n < 0) return newerror("factorial(n): n < 0"); if(n < 9000 && !isdefined("test8900")) { ## builtin is implemented with splitting but only with ## Toom-Cook 2 (by Karatsuba (the father)) return n!; } shift = __CZ__prime_divisors(n,prime); prime = 3; cut = n//2; pix = pix(cut); prime_list = mat[pix]; expo_list = mat[pix]; k = 0; /* Peter Borwein's algorithm @Article{journals/jal/Borwein85, author = {Borwein, Peter B.}, title = {On the Complexity of Calculating Factorials.}, journal = {J. Algorithms}, year = {1985}, number = {3}, url = {http://dblp.uni-trier.de/db/journals/jal/jal6.html#Borwein85} */ do { prime_list[k] = prime; expo_list[k++] = __CZ__prime_divisors(n,prime); prime = nextprime(prime); }while(prime <= cut); /* size of the largest exponent in bits */ k1 = highbit(expo_list[0]); k2 = size(prime_list)-1; for(;k1>=0;k1--){ /* the cut-off for T-C-4 ist still to low, using T-C-3 here TODO: check cutoffs */ result = toomcook3square(result); /* almost all time is spend in this loop, so cutting of the upper half of the primes makes sense */ for(k=0; k<=k2; k++) { if((expo_list[k] & (1 << k1)) != 0) { result *= prime_list[k]; } } } primorial = primorial( cut, n); result *= primorial; result <<= shift; return result; } /* Helper for primorial: do the product with binary splitting TODO: do it without the intermediate list */ define __CZ__primorial__lowlevel( a, b ,p) { local c; if( b == a) return p ; if( b-a > 1){ c= (b + a) >> 1; return __CZ__primorial__lowlevel( a , c , __CZ__primelist[a] ) * __CZ__primorial__lowlevel( c+1 , b , __CZ__primelist[b] ) ; } return __CZ__primelist[a] * __CZ__primelist[b]; } /* Primorial, Product of consecutive primes in range a,b Originally meant to do primorials with a start different from 2, but found out that this is faster at about a=1,b>=10^5 than the builtin function pfact(). With the moderately small list a=1,b=10^6 (78498 primes) it is 3 times faster. A quick look-up showed what was already guessed: pfact() does it linearly. (BTW: what is the time complexity of the primorial with the naive algorithm?) */ define primorial(a,b) { local C1 C2; if(!isint(a)) return newerror("primorial(a,b): a is not an integer"); else if(!isint(b)) return newerror("primorial(a,b): b is not an integer"); else if(a < 0) return newerror("primorial(a,b): a < 0"); else if( b < 2 ) return newerror("primorial(a,b): b < 2"); else if( b < a) return newerror("primorial(a,b): b < a"); else{ /* last prime < 2^32 is also max. prime for nextprime()*/ if(b >= 4294967291) return newerror("primorial(a,b): max. prime exceeded"); if(b == 2) return 2; /* Can be extended by way of pfact(b)/pfact(floor(a-1/2)) for small a */ if(a<=2 && b < 10^5) return pfact(b); /* TODO: use pix() and a simple array (mat[])instead*/ __CZ__primelist = list(); __CZ__fill_prime_list(a,b); C1 = size(__CZ__primelist)-1; return __CZ__primorial__lowlevel( 0, C1,1) } } /* * restore internal function from resource debugging * report important interface functions */ config("resource_debug", resource_debug_level),; if (config("resource_debug") & 3) { print "factorial(n)"; print "primorial(a, b)"; }