1/* 2 * factorial2 - implementation of different factorial related functions 3 * 4 * Copyright (C) 2013,2021 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 factorial toomcook specialfunctions; 36 37 38/* 39 Factorize a factorial and put the result in a 2-column matrix with pi(n) rows 40 mat[ primes , exponent ] 41 Result can be restricted to start at a prime different from 2 with the second 42 argument "start". That arguments gets taken at face value if it prime and 43 smaller than n, otherwise the next larger prime is taken if that prime is 44 smaller than n. 45*/ 46 47define __CZ__factor_factorial(n,start){ 48 local prime prime_list k pix stop; 49 50 51 if(!isint(n)) return 52 newerror("__CZ__factor_factorial(n,start): n is not integer"); 53 if(n < 0) return newerror("__CZ__factor_factorial(n,start): n < 0"); 54 if(n == 1) return newerror("__CZ__factor_factorial(n,start): n == 1"); 55 56 if(start){ 57 if(!isint(start) && start < 0 && start > n) 58 return newerror("__CZ__factor_factorial(n,start): value of " 59 "parameter 'start' out of range"); 60 if(start == n && isprime(n)){ 61 prime_list = mat[1 , 2]; 62 prime_list[0,0] = n; 63 prime_list[0,1] = 1; 64 } 65 else if(!isprime(start) && nextprime(start) >n) 66 return newerror("__CZ__factor_factorial(n,start): value of parameter " 67 "'start' out of range"); 68 else{ 69 if(!isprime(start)) prime = nextprime(start); 70 else prime = start; 71 } 72 } 73 else 74 prime = 2; 75 76 pix = pix(n); 77 if(start){ 78 pix -= pix(prime) -1; 79 } 80 prime_list = mat[pix , 2]; 81 82 k = 0; 83 84 do { 85 prime_list[k ,0] = prime; 86 prime_list[k++,1] = __CZ__prime_divisors(n,prime); 87 prime = nextprime(prime); 88 }while(prime <= n); 89 90 return prime_list; 91} 92 93/* 94 95 subtracts exponents of n_1! from exponents of n_2! with n_1<=n_2 96 97 Does not check for size or consecutiveness of the primes or a carry 98*/ 99 100define __CZ__subtract_factored_factorials(matrix_2n,matrix_n){ 101 local k ret len1,len2,tmp count p e; 102 len1 = size(matrix_n)/2; 103 len2 = size(matrix_2n)/2; 104 if(len2<len1){ 105 106 swap(len1,len2); 107 tmp = matrix_n; 108 matrix_n = matrix_2n; 109 matrix_2n = tmp; 110 } 111 tmp = mat[len1,2]; 112 k = 0; 113 114 for(;k<len1;k++){ 115 p = matrix_2n[k,0]; 116 e = matrix_2n[k,1] - matrix_n[k,1]; 117 if(e!=0){ 118 tmp[count ,0] = p; 119 tmp[count++,1] = e; 120 } 121 } 122 ret = mat[count + (len2-len1),2]; 123 for(k=0;k<count;k++){ 124 ret[k,0] = tmp[k,0]; 125 ret[k,1] = tmp[k,1]; 126 } 127 128 free(tmp); 129 for(k=len1;k<len2;k++){ 130 ret[count,0] = matrix_2n[k,0]; 131 ret[count++,1] = matrix_2n[k,1]; 132 } 133 return ret; 134} 135 136/* 137 138 adds exponents of n_1! to exponents of n_2! with n_1<=n_2 139 140 Does not check for size or consecutiveness of the primes or a carry 141*/ 142 143define __CZ__add_factored_factorials(matrix_2n,matrix_n){ 144 local k ret len1,len2,tmp; 145 len1 = size(matrix_n)/2; 146 len2 = size(matrix_2n)/2; 147 if(len2<len1){ 148 swap(len1,len2); 149 tmp = matrix_n; 150 matrix_n = matrix_2n; 151 matrix_2n = tmp; 152 } 153 ret = mat[len2,2]; 154 k = 0; 155 for(;k<len1;k++){ 156 ret[k,0] = matrix_2n[k,0]; 157 ret[k,1] = matrix_2n[k,1] + matrix_n[k,1]; 158 } 159 for(;k<len2;k++){ 160 ret[k,0] = matrix_2n[k,0]; 161 ret[k,1] = matrix_2n[k,1]; 162 } 163 return ret; 164} 165 166/* 167 Does not check if all exponents are positive 168 169 170 timings 171 this comb comb-this rel. k/n 172; benchmark_binomial(10,13) 173n=2^13 k=2^10 0.064004 0.016001 + 0.76923076923076923077 174n=2^13 k=2^11 0.064004 0.048003 + 0.84615384615384615385 175n=2^13 k=2^12 0.068004 0.124008 - 0.92307692307692307692 176; benchmark_binomial(10,15) 177n=2^15 k=2^10 0.216014 0.024001 + 0.66666666666666666667 178n=2^15 k=2^11 0.220014 0.064004 + 0.73333333333333333333 179n=2^15 k=2^12 0.228014 0.212014 + 0.8 180n=2^15 k=2^13 0.216013 0.664042 - 0.86666666666666666667 181n=2^15 k=2^14 0.240015 1.868117 - 0.93333333333333333333 182; benchmark_binomial(11,15) 183n=2^15 k=2^11 0.216014 0.068004 + 0.73333333333333333333 184n=2^15 k=2^12 0.236015 0.212013 + 0.8 185n=2^15 k=2^13 0.216013 0.656041 - 0.86666666666666666667 186n=2^15 k=2^14 0.244016 1.872117 - 0.93333333333333333333 187; benchmark_binomial(11,18) 188n=2^18 k=2^11 1.652103 0.100006 + 0.61111111111111111111 189n=2^18 k=2^12 1.608101 0.336021 + 0.66666666666666666667 190n=2^18 k=2^13 1.700106 1.140071 + 0.72222222222222222222 191n=2^18 k=2^14 1.756109 3.924245 - 0.77777777777777777778 192n=2^18 k=2^15 2.036127 13.156822 - 0.83333333333333333333 193n=2^18 k=2^16 2.172135 41.974624 - 0.88888888888888888889 194n=2^18 k=2^17 2.528158 121.523594 - 0.94444444444444444444 195; benchmark_binomial(15,25) 196n=2^25 k=2^15 303.790985 38.266392 + 0.6 197; benchmark_binomial(17,25) 198n=2^25 k=2^17 319.127944 529.025062 - 0.68 199*/ 200 201define benchmark_binomial(s,limit){ 202 local ret k A B T1 T2 start end N K; 203 N = 2^(limit); 204 for(k=s;k<limit;k++){ 205 K = 2^k; 206 start=usertime();A=binomial(N,K);end=usertime(); 207 T1 = end-start; 208 start=usertime();B=comb(N,K);end=usertime(); 209 T2 = end-start; 210 print "n=2^"limit,"k=2^"k," ",T1," ",T2,T1<T2?"-":"+"," "k/limit; 211 if(A!=B){ 212 print "false"; 213 break; 214 } 215 } 216} 217 218define __CZ__multiply_factored_factorial(matrix,stop){ 219 local prime result shift prime_list k k1 k2 expo_list pix count start; 220 local hb flag; 221 222 result = 1; 223 shift = 0; 224 225 226 if(!ismat(matrix)) 227 return newerror("__CZ__multiply_factored_factorial(matrix): " 228 "argument matrix not a matrix "); 229 if(!matrix[0,0]) 230 return 231 newerror("__CZ__multiply_factored_factorial(matrix): " 232 "matrix[0,0] is null/0"); 233 234 if(!isnull(stop)) 235 pix = stop; 236 else 237 pix = size(matrix)/2-1; 238 239 if(matrix[0,0] == 2 && matrix[0,1] > 0){ 240 shift = matrix[0,1]; 241 if(pix-1 == 0) 242 return 2^matrix[0,1]; 243 } 244 245 /* 246 This is a more general way to do the multiplication, so any optimization 247 must have been done by the caller. 248 */ 249 k = 0; 250 /* 251 The size of the largest exponent in bits is calculated dynamically. 252 Can be done more elegantly and saves one run over the whole array if done 253 inside the main loop. 254 */ 255 hb =0; 256 for(k=0;k<pix;k++){ 257 k1=highbit(matrix[k,1]); 258 if(hb < k1)hb=k1; 259 } 260 261 k2 = pix; 262 start = 0; 263 if(shift) start++; 264 265 for(k1=hb;k1>=0;k1--){ 266 /* 267 the cut-off for T-C-4 ist still too low, using T-C-3 here 268 TODO: check cutoffs 269 */ 270 result = toomcook3square(result); 271 272 for(k=start; k<=k2; k++) { 273 if((matrix[k,1] & (1 << k1)) != 0) { 274 result *= matrix[k,0]; 275 } 276 } 277 } 278 279 result <<= shift; 280 return result; 281} 282 283/* 284 Compute binomial coefficients n!/(k!(n-k)!) 285 286 One of the rare cases where a formula once meant to ease manual computation 287 is actually the (asymptotically) fastest way to do it (in July 2013) for 288 the extreme case binomial(2N,N) but for a high price, the memory 289 needed is pi(N)--theoretically. 290*/ 291define binomial(n,k){ 292 local ret factored_n factored_k factored_nk denom num quot K prime_list prime; 293 local pix diff; 294 295 if(!isint(n) || !isint(k)) 296 return newerror("binomial(n,k): input is not integer"); 297 if(n<0 || k<0) 298 return newerror("binomial(n,k): input is not >= 0"); ; 299 if(n<k ) return 0; 300 if(n==k) return 1; 301 if(k==0) return 1; 302 if(k==1) return n; 303 if(n-k==1) return n; 304 /* 305 cut-off depends on real size of n,k and size of n/k 306 The current cut-off is to small for large n, e.g.: 307 for 2n=2^23, k=n-n/2 the quotient is q=2n/k=0.25. Empirical tests showed 308 that 2n=2^23 and k=2^16 with q=0.0078125 are still faster than the 309 builtin function. 310 311 The symmetry (n,k) = (n,n-k) is of not much advantage here. One way 312 might be to get closer to k=n/2 if k<n-k but only if the difference 313 is small and n very large. 314 */ 315 if(n<2e4 && !isdefined("test8900")) return comb(n,k); 316 if(n<2e4 && k< n-n/2 && !isdefined("test8900")) return comb(n,k); 317 /* 318 This should be done in parallel to save some memory, e.g. no temporary 319 arrays are needed, all can be done inline. 320 The theoretical memory needed is pi(k). 321 Which is still a lot. 322 */ 323 324 prime = 2; 325 pix = pix(n); 326 prime_list = mat[pix , 2]; 327 K = 0; 328 do { 329 prime_list[K ,0] = prime; 330 diff = __CZ__prime_divisors(n,prime)- 331 ( __CZ__prime_divisors(n-k,prime)+__CZ__prime_divisors(k,prime)); 332 if(diff != 0) 333 prime_list[K++,1] = diff; 334 prime = nextprime(prime); 335 }while(prime <= k); 336 337 do { 338 prime_list[K ,0] = prime; 339 diff = __CZ__prime_divisors(n,prime)-__CZ__prime_divisors(n-k,prime); 340 if(diff != 0) 341 prime_list[K++,1] = diff; 342 prime = nextprime(prime); 343 }while(prime <= n-k); 344 345 do { 346 prime_list[K ,0] = prime; 347 prime_list[K++,1] = __CZ__prime_divisors(n,prime); 348 prime = nextprime(prime); 349 }while(prime <= n); 350 ##print K,pix(k),pix(n-k),pix(n); 351 ##factored_k = __CZ__factor_factorial(k,1); 352 ##factored_nk = __CZ__factor_factorial(n-k,1); 353 354 ##denom = __CZ__add_factored_factorials(factored_k,factored_nk); 355 ##free(factored_k,factored_nk); 356 ##num = __CZ__factor_factorial(n,1); 357 ##quot = __CZ__subtract_factored_factorials( num , denom ); 358 ##free(num,denom); 359 360 ret = __CZ__multiply_factored_factorial(`prime_list,K-1); 361 362 return ret; 363} 364 365/* 366 Compute large catalan numbers C(n) = binomial(2n,n)/(n+1) with 367 cut-off: (n>5e4) 368 Needs a lot of memory. 369*/ 370define bigcatalan(n){ 371 if(!isint(n) )return newerror("bigcatalan(n): n is not integer"); 372 if( n<0) return newerror("bigcatalan(n): n < 0"); 373 if( n<5e4 && !isdefined("test8900") ) return catalan(n); 374 return binomial(2*n,n)/(n+1); 375} 376 377/* 378 df(-111) = -1/3472059605858239446587523014902616804783337112829102414124928 379 7753332469144201839599609375 380 381 df(-3+1i) = 0.12532538977287649201-0.0502372106177184607i 382 df(2n + 1) = (2*n)!/(n!*2^n) 383*/ 384define __CZ__double_factorial(n){ 385 local n1 n2 diff prime pix K prime_list k; 386 prime = 3; 387 pix = pix(2*n)+1; 388 prime_list = mat[pix , 2]; 389 K = 0; 390 do { 391 prime_list[K ,0] = prime; 392 diff = __CZ__prime_divisors(2*n,prime)-( __CZ__prime_divisors(n,prime)); 393 if(diff != 0) 394 prime_list[K++,1] = diff; 395 prime = nextprime(prime); 396 }while(prime <= n); 397 do { 398 prime_list[K ,0] = prime; 399 prime_list[K++,1] = __CZ__prime_divisors(2*n,prime); 400 prime = nextprime(prime); 401 }while(prime <= 2*n); 402 return __CZ__multiply_factored_factorial(prime_list,K); 403/* 404 n1=__CZ__factor_factorial(2*n,1); 405 n1[0,1] = n1[0,1]-n; 406 n2=__CZ__factor_factorial(n,1); 407 diff=__CZ__subtract_factored_factorials( n1 , n2 ); 408 return __CZ__multiply_factored_factorial(diff); 409*/ 410} 411 412##1, 1, 3, 15, 105, 945, 10395, 135135, 2027025, 34459425, 654729075, 413##13749310575, 316234143225, 7905853580625, 213458046676875, 414##6190283353629375, 191898783962510625, 6332659870762850625, 415##221643095476699771875, 8200794532637891559375 416 417## 1, 2, 8, 48, 384, 3840, 46080, 645120, 10321920, 185794560, 418##3715891200, 81749606400, 1961990553600, 51011754393600, 419##1428329123020800, 42849873690624000, 1371195958099968000, 420##46620662575398912000, 1678343852714360832000, 63777066403145711616000 421define doublefactorial(n){ 422 local n1 n2 diff eps ret; 423 if(!isint(n) ){ 424 /* 425 Probably one of the not-so-good ideas. See result of 426 http://www.wolframalpha.com/input/?i=doublefactorial%28a%2Bbi%29 427 */ 428 eps=epsilon(epsilon()*1e-2); 429 ret = 2^(n/2-1/4 * cos(pi()* n)+1/4) * pi()^(1/4 * 430 cos(pi()* n)-1/4)* gamma(n/2+1); 431 epsilon(eps); 432 return ret; 433 } 434 if(n==2) return 2; 435 if(n==3) return 3; 436 switch(n){ 437 case -1: 438 case 0 : return 1;break; 439 case 2 : return 2;break; 440 case 3 : return 3;break; 441 case 4 : return 8;break; 442 default: break; 443 } 444 if(isodd(n)){ 445 /* 446 TODO: find reasonable cutoff 447 df(2n + 1) = (2*n)!/(n!*2^n) 448 */ 449 if(n>0){ 450 n = (n+1)//2; 451 return __CZ__double_factorial(n); 452 } 453 else{ 454 if(n == -3 ) return -1; 455 n = ((-n)-1)/2; 456 return ((-1)^-n)/__CZ__double_factorial(n); 457 } 458 } 459 else{ 460 /* 461 I'm undecided here. The formula for complex n is valid for the negative 462 integers, too. 463 */ 464 n = n>>1; 465 if(n>0){ 466 if(!isdefined("test8900")) 467 return factorial(n)<<n; 468 else 469 return n!<<n; 470 } 471 else 472 return newerror("doublefactorial(n): even(n) < 0"); 473 } 474} 475 476/* 477 Algorithm 3.17, 478 Donald Kreher and Douglas Simpson, 479 Combinatorial Algorithms, 480 CRC Press, 1998, page 89. 481*/ 482static __CZ__stirling1; 483static __CZ__stirling1_n = -1; 484static __CZ__stirling1_m = -1; 485 486define stirling1(n,m){ 487 local i j k; 488 if(n<0)return newerror("stirling1(n,m): n <= 0"); 489 if(m<0)return newerror("stirling1(n,m): m < 0"); 490 if(n<m) return 0; 491 if(n==m) return 1; 492 if(m==0 || n==0) return 0; 493 /* We always use the list */ 494 /* 495 if(m=1){ 496 if(iseven(n)) return -factorial(n-1); 497 else return factorial(n-1); 498 } 499 if(m == n-1){ 500 if(iseven(n)) return -binomial(n,2); 501 else return -binomial(n,2); 502 } 503 */ 504 if(__CZ__stirling1_n >= n && __CZ__stirling1_m >= m){ 505 return __CZ__stirling1[n,m]; 506 } 507 else{ 508 __CZ__stirling1 = mat[n+1,m+1]; 509 __CZ__stirling1[0,0] = 1; 510 for(i=1;i<=n;i++) 511 __CZ__stirling1[i,0] = 0; 512 for(i=1;i<=n;i++){ 513 for(j=1;j<=m;j++){ 514 if(j<=i){ 515 __CZ__stirling1[i, j] = __CZ__stirling1[i - 1, j - 1] - (i - 1)\ 516 * __CZ__stirling1[i - 1, j]; 517 } 518 else{ 519 __CZ__stirling1[i, j] = 0; 520 } 521 } 522 } 523 __CZ__stirling1_n = n; 524 __CZ__stirling1_m = m; 525 return __CZ__stirling1[n,m]; 526 } 527} 528 529define stirling2(n,m){ 530 local k sum; 531 if(n<0)return newerror("stirling2(n,m): n < 0"); 532 if(m<0)return newerror("stirling2(n,m): m < 0"); 533 if(n<m) return 0; 534 if(n==0 && n!=m) return 0; 535 if(n==m) return 1; 536 if(m==0 )return 0; 537 if(m==1) return 1; 538 if(m==2) return 2^(n-1)-1; 539 /* 540 There are different methods to speed up alternating sums. 541 This one doesn't. 542 */ 543 if(isdefined("test8900")){ 544 for(k=0;k<=m;k++){ 545 sum += (-1)^(m-k)*comb(m,k)*k^n; 546 } 547 return sum/(m!); 548 } 549 else{ 550 for(k=0;k<=m;k++){ 551 sum += (-1)^(m-k)*binomial(m,k)*k^n; 552 } 553 return sum/factorial(m); 554 } 555} 556 557static __CZ__stirling2; 558static __CZ__stirling2_n = -1; 559static __CZ__stirling2_m = -1; 560define stirling2caching(n,m){ 561 local nm i j ; 562 if(n<0)return newerror("stirling2iter(n,m): n < 0"); 563 if(m<0)return newerror("stirling2iter(n,m): m < 0"); 564 /* no shortcuts here */ 565 566 if(n<m) return 0; 567 if(n==0 && n!=m) return 0; 568 if(n==m) return 1; 569 if(m==0 )return 0; 570 if(m==1) return 1; 571 if(m==2) return 2^(n-1)-1; 572 573 nm = n-m; 574 if(__CZ__stirling2_n >= n && __CZ__stirling2_m >= m){ 575 return __CZ__stirling2[n,m]; 576 } 577 else{ 578 __CZ__stirling2 = mat[n+1,m+1]; 579 __CZ__stirling2[0,0] = 1; 580 for(i=1;i<=n;i++){ 581 __CZ__stirling2[i,0] = 0; 582 for(j=1;j<=m;j++){ 583 if(j<=i){ 584 __CZ__stirling2[i, j] = __CZ__stirling2[i -1, j -1] + (j )\ 585 * __CZ__stirling2[i - 1, j]; 586 } 587 else{ 588 __CZ__stirling2[i, j] = 0; 589 } 590 } 591 } 592 } 593 __CZ__stirling2_n = (n); 594 __CZ__stirling2_m = (m); 595 return __CZ__stirling2[n,m]; 596} 597 598define bell(n){ 599 local sum s2list k A; 600 601 if(!isint(n)) return newerror("bell(n): n is not integer"); 602 if(n < 0) return newerror("bell(n): n is not positive"); 603 /* place some more shortcuts here?*/ 604 if(n<=15){ 605 mat A[16] = { 606 1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975, 678570, 607 4213597, 27644437, 190899322, 1382958545 608 }; 609 return A[n]; 610 } 611 /* Start by generating the list of stirling numbers of the second kind */ 612 s2list = stirling2caching(n,n//2); 613 if(iserror(s2list)) 614 return newerror("bell(n): could not build stirling num. list"); 615 sum = 0; 616 for(k=1;k<=n;k++){ 617 sum += stirling2caching(n,k); 618 } 619 return sum; 620} 621 622define subfactorialrecursive(n){ 623 if(n==0) return 1; 624 if(n==1) return 0; 625 if(n==2) return 1; 626 return n * subfactorialrecursive(n-1) + (-1)^n; 627} 628 629/* This is, quite amusingly, faster than the very same algorithm in 630 PARI/GP + GMP*/ 631define subfactorialiterative(n){ 632 local k temp1 temp2 ret; 633 if(n==0) return 1; 634 if(n==1) return 0; 635 if(n==2) return 1; 636 temp1 = 0; 637 ret = 1; 638 for(k=3;k<=n;k++){ 639 temp2 = temp1; 640 temp1 = ret; 641 ret = (k-1) *(temp1 + temp2); 642 } 643 return ret; 644} 645 646define subfactorial(n){ 647 local epsilon eps ret lnfact; 648 if(!isint(n))return newerror("subfactorial(n): n is not integer."); 649 if(n < 0)return newerror("subfactorial(n): n < 0"); 650 return subfactorialiterative(n); 651} 652 653define risingfactorial(x,n){ 654 local num denom quot ret; 655 if(n == 1) return x; 656 if(x==0) return newerror("risingfactorial(x,n): x == 0"); 657 if(!isint(x) || !isint(n)){ 658 return gamma(x+n)/gamma(x); 659 } 660 if(x<1)return newerror("risingfactorial(x,n): integer x and x < 1"); 661 if(x+n < 1)return newerror("risingfactorial(x,n): integer x+n and x+n < 1"); 662 if(x<9000&&n<9000){ 663 return (x+n-1)!/(x-1)!; 664 } 665 else{ 666 num = __CZ__factor_factorial(x+n-1,1); 667 denom = __CZ__factor_factorial(x-1,1); 668 quot = __CZ__subtract_factored_factorials( num , denom ); 669 free(num,denom); 670 ret = __CZ__multiply_factored_factorial(quot); 671 return ret; 672 } 673} 674 675define fallingfactorial(x,n){ 676 local num denom quot ret; 677 if(n == 0) return 1; 678 679 if(!isint(x) || !isint(n)){ 680 if(x == n) return gamma(x+1); 681 return gamma(x+1)/gamma(x-n+1); 682 } 683 else{ 684 if(x<0 || x-n < 0) 685 return newerror("fallingfactorial(x,n): integer x<0 or x-n < 0"); 686 if(x == n) return factorial(x); 687 if(x<9000&&n<9000){ 688 return (x)!/(x-n)!; 689 } 690 else{ 691 num = __CZ__factor_factorial(x,1); 692 denom = __CZ__factor_factorial(x-n,1); 693 quot = __CZ__subtract_factored_factorials( num , denom ); 694 free(num,denom); 695 ret = __CZ__multiply_factored_factorial(quot); 696 return ret; 697 } 698 } 699} 700 701 702/* 703 * restore internal function from resource debugging 704 * report important interface functions 705 */ 706config("resource_debug", resource_debug_level),; 707if (config("resource_debug") & 3) { 708 print "binomial(n,k)"; 709 print "bigcatalan(n)"; 710 print "doublefactorial(n)"; 711 print "subfactorial(n)"; 712 print "stirling1(n,m)"; 713 print "stirling2(n,m)"; 714 print "stirling2caching(n,m)"; 715 print "bell(n)"; 716 print "subfactorial(n)"; 717 print "risingfactorial(x,n)"; 718 print "fallingfactorial(x,n)"; 719} 720