drygdryg/calc
1
0
Fork 0
mirror of https://github.com/lcn2/calc.git synced 2025-08-16 01:03:29 +03:00
Code Issues Packages Projects Releases Wiki Activity
Files
db77e29a232ecb63ee7c85b5369b69d1653470d5
calc/cal/factorial2.cal
Landon Curt Noll db77e29a23 convert ASCII TABs to ASCII SPACEs 
Converted all ASCII tabs to ASCII spaces using a 8 character
tab stop, for all files, except for all Makefiles (plus rpm.mk).
The `git diff -w` reports no changes.
2024-07-11 22:03:52 -07:00

720 lines
20 KiB
Plaintext
Raw Blame History

/*
* factorial2 - implementation of different factorial related functions
*
* Copyright (C) 2013,2021 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 factorial toomcook specialfunctions;
/*
Factorize a factorial and put the result in a 2-column matrix with pi(n) rows
mat[ primes , exponent ]
Result can be restricted to start at a prime different from 2 with the second
argument "start". That arguments gets taken at face value if it prime and
smaller than n, otherwise the next larger prime is taken if that prime is
smaller than n.
*/
define __CZ__factor_factorial(n,start){
local prime prime_list k pix stop;
if(!isint(n)) return
newerror("__CZ__factor_factorial(n,start): n is not integer");
if(n < 0) return newerror("__CZ__factor_factorial(n,start): n < 0");
if(n == 1) return newerror("__CZ__factor_factorial(n,start): n == 1");
if(start){
if(!isint(start) && start < 0 && start > n)
return newerror("__CZ__factor_factorial(n,start): value of "
"parameter 'start' out of range");
if(start == n && isprime(n)){
prime_list = mat[1 , 2];
prime_list[0,0] = n;
prime_list[0,1] = 1;
}
else if(!isprime(start) && nextprime(start) >n)
return newerror("__CZ__factor_factorial(n,start): value of parameter "
"'start' out of range");
else{
if(!isprime(start)) prime = nextprime(start);
else prime = start;
}
}
else
prime = 2;
pix = pix(n);
if(start){
pix -= pix(prime) -1;
}
prime_list = mat[pix , 2];
k = 0;
do {
prime_list[k ,0] = prime;
prime_list[k++,1] = __CZ__prime_divisors(n,prime);
prime = nextprime(prime);
}while(prime <= n);
return prime_list;
}
/*
subtracts exponents of n_1! from exponents of n_2! with n_1<=n_2
Does not check for size or consecutiveness of the primes or a carry
*/
define __CZ__subtract_factored_factorials(matrix_2n,matrix_n){
local k ret len1,len2,tmp count p e;
len1 = size(matrix_n)/2;
len2 = size(matrix_2n)/2;
if(len2<len1){
swap(len1,len2);
tmp = matrix_n;
matrix_n = matrix_2n;
matrix_2n = tmp;
}
tmp = mat[len1,2];
k = 0;
for(;k<len1;k++){
p = matrix_2n[k,0];
e = matrix_2n[k,1] - matrix_n[k,1];
if(e!=0){
tmp[count ,0] = p;
tmp[count++,1] = e;
}
}
ret = mat[count + (len2-len1),2];
for(k=0;k<count;k++){
ret[k,0] = tmp[k,0];
ret[k,1] = tmp[k,1];
}
free(tmp);
for(k=len1;k<len2;k++){
ret[count,0] = matrix_2n[k,0];
ret[count++,1] = matrix_2n[k,1];
}
return ret;
}
/*
adds exponents of n_1! to exponents of n_2! with n_1<=n_2
Does not check for size or consecutiveness of the primes or a carry
*/
define __CZ__add_factored_factorials(matrix_2n,matrix_n){
local k ret len1,len2,tmp;
len1 = size(matrix_n)/2;
len2 = size(matrix_2n)/2;
if(len2<len1){
swap(len1,len2);
tmp = matrix_n;
matrix_n = matrix_2n;
matrix_2n = tmp;
}
ret = mat[len2,2];
k = 0;
for(;k<len1;k++){
ret[k,0] = matrix_2n[k,0];
ret[k,1] = matrix_2n[k,1] + matrix_n[k,1];
}
for(;k<len2;k++){
ret[k,0] = matrix_2n[k,0];
ret[k,1] = matrix_2n[k,1];
}
return ret;
}
/*
Does not check if all exponents are positive
timings
this comb comb-this rel. k/n
; benchmark_binomial(10,13)
n=2^13 k=2^10 0.064004 0.016001 + 0.76923076923076923077
n=2^13 k=2^11 0.064004 0.048003 + 0.84615384615384615385
n=2^13 k=2^12 0.068004 0.124008 - 0.92307692307692307692
; benchmark_binomial(10,15)
n=2^15 k=2^10 0.216014 0.024001 + 0.66666666666666666667
n=2^15 k=2^11 0.220014 0.064004 + 0.73333333333333333333
n=2^15 k=2^12 0.228014 0.212014 + 0.8
n=2^15 k=2^13 0.216013 0.664042 - 0.86666666666666666667
n=2^15 k=2^14 0.240015 1.868117 - 0.93333333333333333333
; benchmark_binomial(11,15)
n=2^15 k=2^11 0.216014 0.068004 + 0.73333333333333333333
n=2^15 k=2^12 0.236015 0.212013 + 0.8
n=2^15 k=2^13 0.216013 0.656041 - 0.86666666666666666667
n=2^15 k=2^14 0.244016 1.872117 - 0.93333333333333333333
; benchmark_binomial(11,18)
n=2^18 k=2^11 1.652103 0.100006 + 0.61111111111111111111
n=2^18 k=2^12 1.608101 0.336021 + 0.66666666666666666667
n=2^18 k=2^13 1.700106 1.140071 + 0.72222222222222222222
n=2^18 k=2^14 1.756109 3.924245 - 0.77777777777777777778
n=2^18 k=2^15 2.036127 13.156822 - 0.83333333333333333333
n=2^18 k=2^16 2.172135 41.974624 - 0.88888888888888888889
n=2^18 k=2^17 2.528158 121.523594 - 0.94444444444444444444
; benchmark_binomial(15,25)
n=2^25 k=2^15 303.790985 38.266392 + 0.6
; benchmark_binomial(17,25)
n=2^25 k=2^17 319.127944 529.025062 - 0.68
*/
define benchmark_binomial(s,limit){
local ret k A B T1 T2 start end N K;
N = 2^(limit);
for(k=s;k<limit;k++){
K = 2^k;
start=usertime();A=binomial(N,K);end=usertime();
T1 = end-start;
start=usertime();B=comb(N,K);end=usertime();
T2 = end-start;
print "n=2^"limit,"k=2^"k," ",T1," ",T2,T1<T2?"-":"+"," "k/limit;
if(A!=B){
print "false";
break;
}
}
}
define __CZ__multiply_factored_factorial(matrix,stop){
local prime result shift prime_list k k1 k2 expo_list pix count start;
local hb flag;
result = 1;
shift = 0;
if(!ismat(matrix))
return newerror("__CZ__multiply_factored_factorial(matrix): "
"argument matrix not a matrix ");
if(!matrix[0,0])
return
newerror("__CZ__multiply_factored_factorial(matrix): "
"matrix[0,0] is null/0");
if(!isnull(stop))
pix = stop;
else
pix = size(matrix)/2-1;
if(matrix[0,0] == 2 && matrix[0,1] > 0){
shift = matrix[0,1];
if(pix-1 == 0)
return 2^matrix[0,1];
}
/*
This is a more general way to do the multiplication, so any optimization
must have been done by the caller.
*/
k = 0;
/*
The size of the largest exponent in bits is calculated dynamically.
Can be done more elegantly and saves one run over the whole array if done
inside the main loop.
*/
hb =0;
for(k=0;k<pix;k++){
k1=highbit(matrix[k,1]);
if(hb < k1)hb=k1;
}
k2 = pix;
start = 0;
if(shift) start++;
for(k1=hb;k1>=0;k1--){
/*
the cut-off for T-C-4 ist still too low, using T-C-3 here
TODO: check cutoffs
*/
result = toomcook3square(result);
for(k=start; k<=k2; k++) {
if((matrix[k,1] & (1 << k1)) != 0) {
result *= matrix[k,0];
}
}
}
result <<= shift;
return result;
}
/*
Compute binomial coefficients n!/(k!(n-k)!)
One of the rare cases where a formula once meant to ease manual computation
is actually the (asymptotically) fastest way to do it (in July 2013) for
the extreme case binomial(2N,N) but for a high price, the memory
needed is pi(N)--theoretically.
*/
define binomial(n,k){
local ret factored_n factored_k factored_nk denom num quot K prime_list prime;
local pix diff;
if(!isint(n) || !isint(k))
return newerror("binomial(n,k): input is not integer");
if(n<0 || k<0)
return newerror("binomial(n,k): input is not >= 0"); ;
if(n<k ) return 0;
if(n==k) return 1;
if(k==0) return 1;
if(k==1) return n;
if(n-k==1) return n;
/*
cut-off depends on real size of n,k and size of n/k
The current cut-off is to small for large n, e.g.:
for 2n=2^23, k=n-n/2 the quotient is q=2n/k=0.25. Empirical tests showed
that 2n=2^23 and k=2^16 with q=0.0078125 are still faster than the
builtin function.
The symmetry (n,k) = (n,n-k) is of not much advantage here. One way
might be to get closer to k=n/2 if k<n-k but only if the difference
is small and n very large.
*/
if(n<2e4 && !isdefined("test8900")) return comb(n,k);
if(n<2e4 && k< n-n/2 && !isdefined("test8900")) return comb(n,k);
/*
This should be done in parallel to save some memory, e.g. no temporary
arrays are needed, all can be done inline.
The theoretical memory needed is pi(k).
Which is still a lot.
*/
prime = 2;
pix = pix(n);
prime_list = mat[pix , 2];
K = 0;
do {
prime_list[K ,0] = prime;
diff = __CZ__prime_divisors(n,prime)-
( __CZ__prime_divisors(n-k,prime)+__CZ__prime_divisors(k,prime));
if(diff != 0)
prime_list[K++,1] = diff;
prime = nextprime(prime);
}while(prime <= k);
do {
prime_list[K ,0] = prime;
diff = __CZ__prime_divisors(n,prime)-__CZ__prime_divisors(n-k,prime);
if(diff != 0)
prime_list[K++,1] = diff;
prime = nextprime(prime);
}while(prime <= n-k);
do {
prime_list[K ,0] = prime;
prime_list[K++,1] = __CZ__prime_divisors(n,prime);
prime = nextprime(prime);
}while(prime <= n);
##print K,pix(k),pix(n-k),pix(n);
##factored_k = __CZ__factor_factorial(k,1);
##factored_nk = __CZ__factor_factorial(n-k,1);
##denom = __CZ__add_factored_factorials(factored_k,factored_nk);
##free(factored_k,factored_nk);
##num = __CZ__factor_factorial(n,1);
##quot = __CZ__subtract_factored_factorials( num , denom );
##free(num,denom);
ret = __CZ__multiply_factored_factorial(`prime_list,K-1);
return ret;
}
/*
Compute large catalan numbers C(n) = binomial(2n,n)/(n+1) with
cut-off: (n>5e4)
Needs a lot of memory.
*/
define bigcatalan(n){
if(!isint(n) )return newerror("bigcatalan(n): n is not integer");
if( n<0) return newerror("bigcatalan(n): n < 0");
if( n<5e4 && !isdefined("test8900") ) return catalan(n);
return binomial(2*n,n)/(n+1);
}
/*
df(-111) = -1/3472059605858239446587523014902616804783337112829102414124928
7753332469144201839599609375
df(-3+1i) = 0.12532538977287649201-0.0502372106177184607i
df(2n + 1) = (2*n)!/(n!*2^n)
*/
define __CZ__double_factorial(n){
local n1 n2 diff prime pix K prime_list k;
prime = 3;
pix = pix(2*n)+1;
prime_list = mat[pix , 2];
K = 0;
do {
prime_list[K ,0] = prime;
diff = __CZ__prime_divisors(2*n,prime)-( __CZ__prime_divisors(n,prime));
if(diff != 0)
prime_list[K++,1] = diff;
prime = nextprime(prime);
}while(prime <= n);
do {
prime_list[K ,0] = prime;
prime_list[K++,1] = __CZ__prime_divisors(2*n,prime);
prime = nextprime(prime);
}while(prime <= 2*n);
return __CZ__multiply_factored_factorial(prime_list,K);
/*
n1=__CZ__factor_factorial(2*n,1);
n1[0,1] = n1[0,1]-n;
n2=__CZ__factor_factorial(n,1);
diff=__CZ__subtract_factored_factorials( n1 , n2 );
return __CZ__multiply_factored_factorial(diff);
*/
}
##1, 1, 3, 15, 105, 945, 10395, 135135, 2027025, 34459425, 654729075,
##13749310575, 316234143225, 7905853580625, 213458046676875,
##6190283353629375, 191898783962510625, 6332659870762850625,
##221643095476699771875, 8200794532637891559375
## 1, 2, 8, 48, 384, 3840, 46080, 645120, 10321920, 185794560,
##3715891200, 81749606400, 1961990553600, 51011754393600,
##1428329123020800, 42849873690624000, 1371195958099968000,
##46620662575398912000, 1678343852714360832000, 63777066403145711616000
define doublefactorial(n){
local n1 n2 diff eps ret;
if(!isint(n) ){
/*
Probably one of the not-so-good ideas. See result of
http://www.wolframalpha.com/input/?i=doublefactorial%28a%2Bbi%29
*/
eps=epsilon(epsilon()*1e-2);
ret = 2^(n/2-1/4 * cos(pi()* n)+1/4) * pi()^(1/4 *
cos(pi()* n)-1/4)* gamma(n/2+1);
epsilon(eps);
return ret;
}
if(n==2) return 2;
if(n==3) return 3;
switch(n){
case -1:
case 0 : return 1;break;
case 2 : return 2;break;
case 3 : return 3;break;
case 4 : return 8;break;
default: break;
}
if(isodd(n)){
/*
TODO: find reasonable cutoff
df(2n + 1) = (2*n)!/(n!*2^n)
*/
if(n>0){
n = (n+1)//2;
return __CZ__double_factorial(n);
}
else{
if(n == -3 ) return -1;
n = ((-n)-1)/2;
return ((-1)^-n)/__CZ__double_factorial(n);
}
}
else{
/*
I'm undecided here. The formula for complex n is valid for the negative
integers, too.
*/
n = n>>1;
if(n>0){
if(!isdefined("test8900"))
return factorial(n)<<n;
else
return n!<<n;
}
else
return newerror("doublefactorial(n): even(n) < 0");
}
}
/*
Algorithm 3.17,
Donald Kreher and Douglas Simpson,
Combinatorial Algorithms,
CRC Press, 1998, page 89.
*/
static __CZ__stirling1;
static __CZ__stirling1_n = -1;
static __CZ__stirling1_m = -1;
define stirling1(n,m){
local i j k;
if(n<0)return newerror("stirling1(n,m): n <= 0");
if(m<0)return newerror("stirling1(n,m): m < 0");
if(n<m) return 0;
if(n==m) return 1;
if(m==0 || n==0) return 0;
/* We always use the list */
/*
if(m=1){
if(iseven(n)) return -factorial(n-1);
else return factorial(n-1);
}
if(m == n-1){
if(iseven(n)) return -binomial(n,2);
else return -binomial(n,2);
}
*/
if(__CZ__stirling1_n >= n && __CZ__stirling1_m >= m){
return __CZ__stirling1[n,m];
}
else{
__CZ__stirling1 = mat[n+1,m+1];
__CZ__stirling1[0,0] = 1;
for(i=1;i<=n;i++)
__CZ__stirling1[i,0] = 0;
for(i=1;i<=n;i++){
for(j=1;j<=m;j++){
if(j<=i){
__CZ__stirling1[i, j] = __CZ__stirling1[i - 1, j - 1] - (i - 1)\
* __CZ__stirling1[i - 1, j];
}
else{
__CZ__stirling1[i, j] = 0;
}
}
}
__CZ__stirling1_n = n;
__CZ__stirling1_m = m;
return __CZ__stirling1[n,m];
}
}
define stirling2(n,m){
local k sum;
if(n<0)return newerror("stirling2(n,m): n < 0");
if(m<0)return newerror("stirling2(n,m): m < 0");
if(n<m) return 0;
if(n==0 && n!=m) return 0;
if(n==m) return 1;
if(m==0 )return 0;
if(m==1) return 1;
if(m==2) return 2^(n-1)-1;
/*
There are different methods to speed up alternating sums.
This one doesn't.
*/
if(isdefined("test8900")){
for(k=0;k<=m;k++){
sum += (-1)^(m-k)*comb(m,k)*k^n;
}
return sum/(m!);
}
else{
for(k=0;k<=m;k++){
sum += (-1)^(m-k)*binomial(m,k)*k^n;
}
return sum/factorial(m);
}
}
static __CZ__stirling2;
static __CZ__stirling2_n = -1;
static __CZ__stirling2_m = -1;
define stirling2caching(n,m){
local nm i j ;
if(n<0)return newerror("stirling2iter(n,m): n < 0");
if(m<0)return newerror("stirling2iter(n,m): m < 0");
/* no shortcuts here */
if(n<m) return 0;
if(n==0 && n!=m) return 0;
if(n==m) return 1;
if(m==0 )return 0;
if(m==1) return 1;
if(m==2) return 2^(n-1)-1;
nm = n-m;
if(__CZ__stirling2_n >= n && __CZ__stirling2_m >= m){
return __CZ__stirling2[n,m];
}
else{
__CZ__stirling2 = mat[n+1,m+1];
__CZ__stirling2[0,0] = 1;
for(i=1;i<=n;i++){
__CZ__stirling2[i,0] = 0;
for(j=1;j<=m;j++){
if(j<=i){
__CZ__stirling2[i, j] = __CZ__stirling2[i -1, j -1] + (j )\
* __CZ__stirling2[i - 1, j];
}
else{
__CZ__stirling2[i, j] = 0;
}
}
}
}
__CZ__stirling2_n = (n);
__CZ__stirling2_m = (m);
return __CZ__stirling2[n,m];
}
define bell(n){
local sum s2list k A;
if(!isint(n)) return newerror("bell(n): n is not integer");
if(n < 0) return newerror("bell(n): n is not positive");
/* place some more shortcuts here?*/
if(n<=15){
mat A[16] = {
1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975, 678570,
4213597, 27644437, 190899322, 1382958545
};
return A[n];
}
/* Start by generating the list of stirling numbers of the second kind */
s2list = stirling2caching(n,n//2);
if(iserror(s2list))
return newerror("bell(n): could not build stirling num. list");
sum = 0;
for(k=1;k<=n;k++){
sum += stirling2caching(n,k);
}
return sum;
}
define subfactorialrecursive(n){
if(n==0) return 1;
if(n==1) return 0;
if(n==2) return 1;
return n * subfactorialrecursive(n-1) + (-1)^n;
}
/* This is, quite amusingly, faster than the very same algorithm in
PARI/GP + GMP*/
define subfactorialiterative(n){
local k temp1 temp2 ret;
if(n==0) return 1;
if(n==1) return 0;
if(n==2) return 1;
temp1 = 0;
ret = 1;
for(k=3;k<=n;k++){
temp2 = temp1;
temp1 = ret;
ret = (k-1) *(temp1 + temp2);
}
return ret;
}
define subfactorial(n){
local epsilon eps ret lnfact;
if(!isint(n))return newerror("subfactorial(n): n is not integer.");
if(n < 0)return newerror("subfactorial(n): n < 0");
return subfactorialiterative(n);
}
define risingfactorial(x,n){
local num denom quot ret;
if(n == 1) return x;
if(x==0) return newerror("risingfactorial(x,n): x == 0");
if(!isint(x) || !isint(n)){
return gamma(x+n)/gamma(x);
}
if(x<1)return newerror("risingfactorial(x,n): integer x and x < 1");
if(x+n < 1)return newerror("risingfactorial(x,n): integer x+n and x+n < 1");
if(x<9000&&n<9000){
return (x+n-1)!/(x-1)!;
}
else{
num = __CZ__factor_factorial(x+n-1,1);
denom = __CZ__factor_factorial(x-1,1);
quot = __CZ__subtract_factored_factorials( num , denom );
free(num,denom);
ret = __CZ__multiply_factored_factorial(quot);
return ret;
}
}
define fallingfactorial(x,n){
local num denom quot ret;
if(n == 0) return 1;
if(!isint(x) || !isint(n)){
if(x == n) return gamma(x+1);
return gamma(x+1)/gamma(x-n+1);
}
else{
if(x<0 || x-n < 0)
return newerror("fallingfactorial(x,n): integer x<0 or x-n < 0");
if(x == n) return factorial(x);
if(x<9000&&n<9000){
return (x)!/(x-n)!;
}
else{
num = __CZ__factor_factorial(x,1);
denom = __CZ__factor_factorial(x-n,1);
quot = __CZ__subtract_factored_factorials( num , denom );
free(num,denom);
ret = __CZ__multiply_factored_factorial(quot);
return ret;
}
}
}
/*
* restore internal function from resource debugging
* report important interface functions
*/
config("resource_debug", resource_debug_level),;
if (config("resource_debug") & 3) {
print "binomial(n,k)";
print "bigcatalan(n)";
print "doublefactorial(n)";
print "subfactorial(n)";
print "stirling1(n,m)";
print "stirling2(n,m)";
print "stirling2caching(n,m)";
print "bell(n)";
print "subfactorial(n)";
print "risingfactorial(x,n)";
print "fallingfactorial(x,n)";
}
Reference in New Issue View Git Blame Copy Permalink