1%
2%  Testing file for REDUCE Special Functions Package
3%
4%             Chris Cannam, ZIB Berlin
5%             October 1992 -> Feb 1993
6%        (only some of the time, of course)
7%
8%  Corrections and comments to neun@sc.zib-berlin.de
9%
10
11
12on savesfs;	% just in case it's off for some reason
13
14off bfspace;	% to provide more similarity between runs
15		% with old & new bigfloats
16
17let {sinh (~x) => (exp(x) - exp (-x))/2 };
18		% this will improve some results
19
20
21% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
22
23%     1. Bernoulli numbers
24
25% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
26
27
28off rounded;
29
30procedure do!*one!*bern(x);
31   write ("Bernoulli ", x, " is ", bernoulli x);
32
33do!*one!*bern(1);
34do!*one!*bern(2);
35do!*one!*bern(3);
36do!*one!*bern(13);
37do!*one!*bern(14);
38do!*one!*bern(300);
39do!*one!*bern(-2);
40do!*one!*bern(0);
41
42for n := 2 step 2 until 100 do do!*one!*bern n;
43
44on rounded;
45precision 100;
46
47do!*one!*bern(1);
48do!*one!*bern(2);
49do!*one!*bern(3);
50do!*one!*bern(13);
51do!*one!*bern(14);
52do!*one!*bern(300);
53do!*one!*bern(-2);
54do!*one!*bern(0);
55do!*one!*bern(38);
56do!*one!*bern(400);
57
58
59% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
60
61%     2. Gamma function
62
63% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
64
65
66on rounded;
67off complex;
68precision 40;
69
70algebraic procedure wg(x);
71   write ("gamma (", x, ")  ==>  ", gamma x);
72
73algebraic procedure wp(x);
74   write ("-- precision ", x, ",  from ", precision(x));
75
76wg (1/2);
77wg (3/2);
78
79write ("sqrt(pi)/2  ==>  ", sqrt(pi)/2);
80
81wp(10);
82
83for x := 0 step 5 until 100 do
84   << wg (1 + x/1000);
85      wg (-1 - x/13);
86      wp (8+floor(x/4)) >>;
87
88wg(1/1000000003);
89
90off rounded;
91
92gamma(17/2);
93gamma(-17/2);
94gamma(4);
95gamma(0);
96gamma(-4);
97gamma(-17/3);
98
99p := gamma(x**2) * gamma(x-y**gamma(y)) - (1/(gamma(4*(x-y))));
100y := 1/4;
101p;
102x := 3;
103p;
104y := -3/8;
105p;
106
107on rounded, complex;
108precision 50;
109
110p;
111
112off rounded, complex;
113clear y;
114
115p;
116
117
118
119% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
120
121%     3. Beta function.  Not very interesting
122
123% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
124
125
126algebraic procedure do!*one!*beta(x,y);
127   write ("Beta (", x, ",", y, ") = ", beta(x,y));
128
129do!*one!*beta(0,1);
130do!*one!*beta(2,-3);
131do!*one!*beta(3,2);
132do!*one!*beta(a+b,(c+d)**(b-a));
133do!*one!*beta(-3,4);
134do!*one!*beta(-3,2);
135do!*one!*beta(-3,-7.5);
136do!*one!*beta((pi * 10), exp(5));
137
138on rounded;
139precision 30;
140
141do!*one!*beta(0,1);
142do!*one!*beta(2,-3);
143do!*one!*beta(3,2);
144do!*one!*beta(a+b,(c+d)**(b-a));
145do!*one!*beta(-3,4);
146do!*one!*beta(-3,2);
147do!*one!*beta(-3,-7.5);
148do!*one!*beta((pi * 10), exp(5));
149
150
151
152% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
153
154%     4. Pochhammer notation
155
156% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
157
158
159off rounded;
160
161pochhammer(4,5);
162pochhammer(-4,5);
163pochhammer(4,-5);
164pochhammer(-4,-5);
165pochhammer(17/2,12);
166pochhammer(-17/2,12);
167pochhammer(1/3,14)*pochhammer(2/3,15);
168
169q := pochhammer(1/5,11)*pochhammer(2/5,11)*pochhammer(3/5,11)*
170      pochhammer(1-1/5,11)*pochhammer(1,11)*pochhammer(6/5,11)*
171       pochhammer(70/50,11)*pochhammer(8/5,11)*pochhammer(9/5,11);
172
173on complex;
174
175pochhammer(a+b*i,c)*pochhammer(a-b*i,c);
176
177a := 2;
178b := 3;
179c := 5;
180
181pochhammer(a+b*i,c)*pochhammer(a-b*i,c);
182
183off complex;
184on rounded;
185
186pochhammer(1/5,11)*pochhammer(2/5,11)*pochhammer(3/5,11)*
187 pochhammer(1-1/5,11)*pochhammer(1,11)*pochhammer(6/5,11)*
188  pochhammer(70/50,11)*pochhammer(8/5,11)*pochhammer(9/5,11);
189
190q;
191
192pochhammer(pi,floor (pi**8));
193pochhammer(-pi,floor (pi**7));
194pochhammer(1.5,floor (pi**8));
195
196
197
198% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
199
200%     5. Digamma function
201
202% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
203
204
205procedure do!*one!*psi(x);
206   << precision (precision(0) + 4)$
207      write("Psi of ", x, " is ", psi(x) ) >> ;
208
209clear x, y;
210
211z := x * ((x+y)**2 + (x**y));
212
213off rounded;
214
215do!*one!*psi(3);
216do!*one!*psi(pi);
217do!*one!*psi(1.005);
218do!*one!*psi(1.995);
219do!*one!*psi(74);
220do!*one!*psi(-1/2);
221do!*one!*psi(-3);
222do!*one!*psi(z);
223
224on rounded;
225precision 100;
226
227do!*one!*psi(3);
228do!*one!*psi(pi);
229do!*one!*psi(1.005);
230do!*one!*psi(1.995);
231do!*one!*psi(74);
232do!*one!*psi(-1/2);
233do!*one!*psi(-3);
234do!*one!*psi(z);
235
236precision 15;
237
238x := 8/3;
239y := 7/1000;
240do!*one!*psi(z);
241
242off rounded;
243
244clear x, y;
245
246df(psi(z), x);
247df(df(psi(z), y),x);
248int(psi(z), z);
249
250on rounded;
251
252for k := 1 step 0.1 until 2 do do!*one!*psi(k);
253
254off rounded;
255
256% PSI_SIMP.TST  F.J.Wright, 2 July 1993
257
258on evallhseqp;
259factor psi;  on rat, intstr, div;  % for neater output
260% Do not try using "off mcd"!
261
262psi(x+m) - psi(x+m-1) = 1/(x+m-1);
263psi(x+2) - psi(x+1) + 2*psi(x) = 1/(x+1) + 2*psi(x);
264psi(x+2) + 3*psi(x) = 4*psi(x) + 1/x + 1/(x+1);
265
266psi(x + 1) = psi(x) + 1/x;
267psi(x + 3/2) = psi(x + 1/2) + 1/(x + 1/2);
268psi(x - 1/2) = psi(x + 1/2) - 1/(x - 1/2);
269psi((x + 3a)/a);  psi(x/y + 3);
270
271off rat, intstr, div;  on rational;
272
273psi(x+m) - psi(x+m-1) = 1/(x+m-1);
274psi(x+2) - psi(x+1) + 2*psi(x) = 1/(x+1) + 2*psi(x);
275psi(x+2) + 3*psi(x) = 4*psi(x) + 1/x + 1/(x+1);
276
277psi(x + 1) = psi(x) + 1/x;
278psi(x + 3/2) = psi(x + 1/2) + 1/(x + 1/2);
279psi(x - 1/2) = psi(x + 1/2) - 1/(x - 1/2);
280psi((x + 3a)/a);  psi(x/y + 3);
281
282off rational;
283
284% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
285
286%     6. Polygamma functions
287
288% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
289
290
291procedure do!*one!*pg(n,x);
292   write ("Polygamma (", n, ") of ", x, " is ", polygamma(n,x));
293
294off rounded;
295
296do!*one!*pg(1,1/2);
297do!*one!*pg(1,1);
298do!*one!*pg(1,3/2);
299do!*one!*pg(1,1.005);
300do!*one!*pg(1,1.995);
301do!*one!*pg(1,1e-10);
302do!*one!*pg(2,1.45);
303do!*one!*pg(3,1.99);
304do!*one!*pg(4,-8.2);
305do!*one!*pg(5,0);
306do!*one!*pg(6,-5);
307do!*one!*pg(7,200);
308
309on rounded;
310precision 100;
311
312do!*one!*pg(1,1/2);
313do!*one!*pg(1,1);
314do!*one!*pg(1,3/2);
315do!*one!*pg(1,1.005);
316do!*one!*pg(1,1.995);
317do!*one!*pg(1,1e-10);
318do!*one!*pg(2,1.45);
319do!*one!*pg(3,1.99);
320do!*one!*pg(4,-8.2);
321do!*one!*pg(5,0);
322do!*one!*pg(6,-5);
323do!*one!*pg(7,200);
324
325off rounded;
326clear x;
327
328
329% Polygamma differentiation has already
330% been tried a bit in the psi section
331
332df(int(int(int(polygamma(3,x),x),x),x),x);
333
334clear w, y, z;
335
336
337
338% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
339
340%     7. Zeta function
341
342% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
343
344
345procedure do!*one!*zeta(n);
346   write ("Zeta of ", n, " is ", zeta n);
347
348off rounded;
349clear x, y, z;
350
351z := x * ((x+y)**5 + (x**y));
352
353do!*one!*zeta(0);
354
355for k := 4 step 2 until 35 do
356   do!*one!*zeta(k);
357
358for k := 1 step 2 until 31 do
359   do!*one!*zeta(-k);
360
361do!*one!*zeta(-17/3);
362do!*one!*zeta(190);
363do!*one!*zeta(300);
364do!*one!*zeta(0);
365do!*one!*zeta(-44);
366
367on rounded;
368clear x, y;
369
370for k := 3 step 3 until 36 do <<
371   precision (31+k*3);
372   do!*one!*zeta(k) >>;
373
374precision 20;
375
376do!*one!*zeta(-17/3);
377do!*one!*zeta(z);
378
379y := 3;
380x := pi;
381
382do!*one!*zeta(z);
383do!*one!*zeta(190);
384do!*one!*zeta(300);
385do!*one!*zeta(0);
386do!*one!*zeta(-44);
387
388off rounded;
389
390
391
392% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
393
394%     8. Kummer functions
395
396% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
397
398
399off rounded;
400
401t!*kummer!*a := { 1, 2.4, -1397/10  }$
402t!*kummer!*b := { 0, 1, pi, -pi, 26 }$
403
404for each a in t!*kummer!*a do
405   for each b in t!*kummer!*a do
406      for each z in t!*kummer!*a do
407      	 << write "KummerM(", a, ",", b, ",", z, ") = ",
408      	       kummerm(a,b,z);
409      	    write "KummerU(", a, ",", b, ",", z, ") = ",
410      	       kummeru(a,b,z) >>;
411
412on rounded;
413precision 30;
414t!*k!*c := 7;
415
416%  To test each and every possible combination of
417%  three arguments from t!*kummer!*b would take too
418%  long, but we want the possibility of trying most
419%  special cases.  Compromise: test every seventh
420%  possibility.
421
422for each a in t!*kummer!*b do
423   for each b in t!*kummer!*b do
424      for each z in t!*kummer!*b do
425      	 << if t!*k!*c = 7
426      	    then << write "KummerM(", a, ",", b, ",", z, ") = ",
427      	       	       kummerm(a,b,z);
428      	       	    write "KummerU(", a, ",", b, ",", z, ") = ",
429      	       	       kummeru(a,b,z);
430      	       	    t!*k!*c := 0 >>;
431      	    t!*k!*c := t!*k!*c + 1 >>;
432
433off rounded;
434clear x, y, z, t!*k!*c;
435
436df(df(kummerM(a,b,z),z),z);
437df(kummerU(a,b,z),z);
438
439z := ((x^2 + y)^5) + (x^(x+y));
440
441df(df(kummerM(a,b,z),y),x);
442df(kummerU(a,b,z),x);
443
444
445
446% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
447
448%     9. Bessel functions
449
450% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
451
452
453%  Lengthy test of the Bessel functions.  This isn't even
454%  remotely exhaustive of the special cases -- though a
455%  real person with lots of time could no doubt come up
456%  with a better lot of tests than this automated rubbish.
457%  Again, compromise by only actually doing one in five or
458%  nine.  If you want a really thorough test, you can
459%  easily change this to get it; but it'll take hours to
460%  run.
461
462clear p, q;
463
464hankel1(p,q);
465r := df(ws,q);
466
467on complex;
468
469r;
470p:=3/8;
471r;
472q := pi;
473r;
474
475on rounded;
476
477r;
478
479off complex, rounded;
480
481df(df(besselj(pp,qq)+rr * hankel1(pp*2,qq) *
482      bessely(pp-qq,qq),qq),qq);
483
484% Possible values for real args
485t!*bes!*vr := { 1, pi, -pi, 26 }$
486
487% Possible values for real and imaginary parts of complex args
488t!*bes!*vc := { 0, 3, -41/2 }$
489
490array s!*bes(4)$
491
492s!*bes(1) := "BesselJ"$
493s!*bes(2) := "BesselY"$
494s!*bes(3) := "BesselI"$
495s!*bes(4) := "BesselK"$
496
497pre := 16;
498precision pre;
499preord := 10**pre;
500t!*b!*c := 3;
501
502algebraic procedure do!*one!*bessel(s,n,z);
503   (if s = 1 then besselj(n,z)
504      else if s = 2 then bessely(n,z)
505      else if s = 3 then besseli(n,z)
506      else besselk(n,z));
507
508algebraic procedure pr!*bessel(s,n,z,k);
509   << if t!*b!*c = k
510      then
511      	 << on rounded;
512      	    bes1 := do!*one!*bessel(s,n,z);
513      	    precision(pre+5);
514      	    bes2 := do!*one!*bessel(s,n,z);
515      	    if bes1 neq 0
516      	       then disc := floor abs(100*(bes2-bes1)*preord/bes1)
517      	       else disc := 0;
518      	    precision pre;
519      	    write s!*bes(s), "(", n, ",", z, ") = ", bes1;
520      	    if not numberp disc then
521      	       << precom := !*complex;
522      	       	  on complex;
523      	       	  disc := disc;
524      	       	  if not precom then off complex >>;
525      	    if disc neq 0
526      	       then write "   (discrepancy ", disc, "% of one s.f.)";
527      	    if numberp disc and disc > 200 then
528      	       << write "***** WARNING  Significant Inaccuracy.";
529      	       	  write "      Lower  precision result:";
530      	       	  write "      ", bes1;
531      	       	  write "      Higher precision result:";
532      	       	  precision(pre+5); write "      ", bes2; precision pre >>;
533      	    off rounded;
534      	    t!*b!*c := 0 >>;
535      t!*b!*c := t!*b!*c + 1 >>;
536
537%  About to begin Bessel test.  We have a list of possible
538%  values, and we test every Bessel, with every value on the
539%  list as both order and argument.  Every Bessel is computed
540%  twice, to different precisions (differing by 3), and any
541%  discrepancy is reported.  The value reported is the diff-
542%  erence between the two computed values, expressed as a
543%  percentage of the unit of the least significant displayed
544%  digit.  A discrepancy between 100% and 200% means that the
545%  last digit of the displayed value was found to differ at
546%  higher precision; values greater than 200% are cause for
547%  concern.  An ideal discrepancy would be between about 1%
548%  and 20%; if the value is found to be zero, no discrepancy
549%  is reported.
550
551off msg;
552
553for s := 1:4 do
554   << write(" ... Testing ", s!*bes(s), " for real domains ... ");
555      for each n in t!*bes!*vr do
556         for each z in t!*bes!*vr do
557      	    pr!*bessel(s, n, z, 5) >>;
558
559on complex;
560
561for s := 1:3 do
562   << write (" ... Testing  ", s!*bes(s), " for complex domains ... ");
563      for each nr in t!*bes!*vc do
564      	 for each ni in t!*bes!*vc do
565      	    for each zr in t!*bes!*vc do
566      	       for each zi in t!*bes!*vc do
567      	       	  pr!*bessel(s, nr+ni*i, zr+zi*i, 9) >>;
568
569off complex;
570
571on msg;
572
573write (" ...");
574write ("Bessel test complete.");
575
576
577% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
578
579%     10. Incomplete Gamma and Beta functions (regularized)
580
581% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
582
583
584igamma(3,0);
585
586igamma(1,1);
587
588igamma(1,4);
589
590on rounded;
591
592igamma(1,1);
593
594igamma(1,4);
595
596igamma(2,4);
597
598igamma(0.5,4);
599
600beta(1,1,1);
601
602ibeta(1,2,1);
603
604ibeta(1,4,1);
605
606ibeta(2,4,1);
607
608ibeta(2,4,0.5);
609
610ibeta(0.12,0.43,0.9);
611
612precision 50;
613
614ibeta(0.12,0.43,0.9);
615
616precision  20;
617
618ibeta(0.12,0.43,0.9);
619
620on complex;
621
622ibeta(1+i,1,1.5*i);
623
624off rounded,complex;
625
626ibeta(3,2,x);
627
628
629% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
630
631%     11. Dilogarithm, polylogartihm and Lerch_phi
632
633% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
634
635
636polylog(n,0);
637
638polylog(2,1);
639
640polylog(3,1);
641
642polylog(2,i);
643
644df(polylog(a,x),x);
645
646polylog(1,x);
647
648precision reset;
649on rounded;
650
651polylog(2,1/3);
652
653off rounded;
654
655
656lerch_phi(3,4,1);
657
658lerch_phi(4,0,3);
659
660lerch_phi(x,a,1);
661
662lerch_phi(1,x,1);
663
664df(lerch_phi(x,3,4),x);
665
666
667% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
668
669%     12. Constants
670
671% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
672
673on rounded;
674off complex;
675
676precision 50;
677
678euler_gamma;
679
680Khinchin;
681
682golden_ratio;
683
684Catalan;
685
686on complex;
687
688euler_gamma;
689
690Khinchin;
691
692golden_ratio;
693
694Catalan;
695
696
697end;
698
699
700