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