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