Artifact 0e9c301aa72676085d2075a9ef6ee902d308d6ef7eb215c45fc0ee5269da90d5:
- File
r35/xlog/specfn.log
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 82827) [annotate] [blame] [check-ins using] [more...]
Codemist Standard Lisp 3.54 for DEC Alpha: May 23 1994 Dump file created: Mon May 23 10:39:11 1994 REDUCE 3.5, 15-Oct-93 ... Memory allocation: 6023424 bytes +++ About to read file ndotest.red % % Testing file for REDUCE Special Functions Package % % Chris Cannam, ZIB Berlin % October 1992 -> Feb 1993 % (only some of the time, of course) % % Corrections and comments to neun@sc.zib-berlin.de % on savesfs; % just in case it's off for some reason off bfspace; % to provide more similarity between runs % with old & new bigfloats let {sinh (~x) => (exp(x) - exp (-x))/2 }; % this will improve some results % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % 1. Bernoulli numbers % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= off rounded; procedure do!*one!*bern(x); write ("Bernoulli ", x, " is ", bernoulli x); do*one*bern do!*one!*bern(1); - 1 Bernoulli 1 is ------ 2 do!*one!*bern(2); 1 Bernoulli 2 is --- 6 do!*one!*bern(3); Bernoulli 3 is 0 do!*one!*bern(13); Bernoulli 13 is 0 do!*one!*bern(14); 7 Bernoulli 14 is --- 6 do!*one!*bern(300); Bernoulli 300 is ( - 1863878995204859011995045341848156066182191846635905937518715320655775958174360523134990756922303410810482600528769479642021001218415879006164302955370460829146434807964717737195356935144151583424833154250047747433575584999029126775186293388721514970183351129809976971603227633930434923843984829580311593372565398574762880028289167635570012415606941367995702212211519561707046505473575241 )/866054419230 do!*one!*bern(-2); Bernoulli -2 is bernoulli(-2) do!*one!*bern(0); Bernoulli 0 is 1 for n := 2 step 2 until 100 do do!*one!*bern n; 1 Bernoulli 2 is --- 6 - 1 Bernoulli 4 is ------ 30 1 Bernoulli 6 is ---- 42 - 1 Bernoulli 8 is ------ 30 5 Bernoulli 10 is ---- 66 - 691 Bernoulli 12 is -------- 2730 7 Bernoulli 14 is --- 6 - 3617 Bernoulli 16 is --------- 510 43867 Bernoulli 18 is ------- 798 - 174611 Bernoulli 20 is ----------- 330 854513 Bernoulli 22 is -------- 138 - 236364091 Bernoulli 24 is -------------- 2730 8553103 Bernoulli 26 is --------- 6 - 23749461029 Bernoulli 28 is ---------------- 870 8615841276005 Bernoulli 30 is --------------- 14322 - 7709321041217 Bernoulli 32 is ------------------ 510 2577687858367 Bernoulli 34 is --------------- 6 - 26315271553053477373 Bernoulli 36 is ------------------------- 1919190 2929993913841559 Bernoulli 38 is ------------------ 6 - 261082718496449122051 Bernoulli 40 is -------------------------- 13530 1520097643918070802691 Bernoulli 42 is ------------------------ 1806 - 27833269579301024235023 Bernoulli 44 is ---------------------------- 690 596451111593912163277961 Bernoulli 46 is -------------------------- 282 - 5609403368997817686249127547 Bernoulli 48 is --------------------------------- 46410 495057205241079648212477525 Bernoulli 50 is ----------------------------- 66 - 801165718135489957347924991853 Bernoulli 52 is ----------------------------------- 1590 29149963634884862421418123812691 Bernoulli 54 is ---------------------------------- 798 - 2479392929313226753685415739663229 Bernoulli 56 is --------------------------------------- 870 84483613348880041862046775994036021 Bernoulli 58 is ------------------------------------- 354 - 1215233140483755572040304994079820246041491 Bernoulli 60 is ------------------------------------------------ 56786730 12300585434086858541953039857403386151 Bernoulli 62 is ---------------------------------------- 6 - 106783830147866529886385444979142647942017 Bernoulli 64 is ----------------------------------------------- 510 1472600022126335654051619428551932342241899101 Bernoulli 66 is ------------------------------------------------ 64722 - 78773130858718728141909149208474606244347001 Bernoulli 68 is ------------------------------------------------- 30 1505381347333367003803076567377857208511438160235 Bernoulli 70 is --------------------------------------------------- 4686 Bernoulli 72 is ( - 5827954961669944110438277244641067365282488301844260429 )/140100870 34152417289221168014330073731472635186688307783087 Bernoulli 74 is ---------------------------------------------------- 6 Bernoulli 76 is ( - 24655088825935372707687196040585199904365267828865801 )/30 Bernoulli 78 is 414846365575400828295179035549542073492199375372400483487 /3318 Bernoulli 80 is ( - 4603784299479457646935574969019046849794257872751288919656867 )/230010 Bernoulli 82 is 1677014149185145836823154509786269900207736027570253414881613 /498 Bernoulli 84 is ( - 2024576195935290360231131160111731009989917391198090877281083932477 )/ 3404310 Bernoulli 86 is 660714619417678653573847847426261496277830686653388931761996983 /6 Bernoulli 88 is ( - 1311426488674017507995511424019311843345750275572028644296919890574047 )/61410 Bernoulli 90 is 1179057279021082799884123351249215083775254949669647116231545215727922535 /272118 Bernoulli 92 is ( - 1295585948207537527989427828538576749659341483719435143023316326829946247 )/1410 Bernoulli 94 is 1220813806579744469607301679413201203958508415202696621436215105284649447 /6 Bernoulli 96 is ( - 211600449597266513097597728109824233673043954389060234150638733420050668349987259 )/4501770 Bernoulli 98 is 67908260672905495624051117546403605607342195728504487509073961249992947058239 /6 Bernoulli 100 is ( - 94598037819122125295227433069493721872702841533066936133385696204311395415197247711 )/33330 on rounded; precision 100; 12 do!*one!*bern(1); Bernoulli 1 is - 0.5 do!*one!*bern(2); Bernoulli 2 is 0.1666666666666666666666666666666666666666666666666666 666666666666666666666666666666666666666666666667 do!*one!*bern(3); Bernoulli 3 is 0 do!*one!*bern(13); Bernoulli 13 is 0 do!*one!*bern(14); Bernoulli 14 is 1.166666666666666666666666666666666666666666666666666 666666666666666666666666666666666666666666666667 do!*one!*bern(300); Bernoulli 300 is - 2.15214997327998682971981737675608819857345657278 3806701849343925846796995824343474766557085976828531 e+375 do!*one!*bern(-2); Bernoulli -2 is bernoulli(-2) do!*one!*bern(0); Bernoulli 0 is 1 do!*one!*bern(38); Bernoulli 38 is 4.883323189735931666666666666666666666666666666666666 666666666666666666666666666666666666666666666667e+14 do!*one!*bern(400); Bernoulli 400 is - 6.84694485580645336061625858231088359767823009718 0625742414781150311357197834589950328990573681886442 e+549 % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % 2. Gamma function % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= on rounded; off complex; precision 40; 100 algebraic procedure wg(x); write ("gamma (", x, ") ==> ", gamma x); wg algebraic procedure wp(x); write ("-- precision ", x, ", from ", precision(x)); wp wg (1/2); gamma (0.5) ==> 1.772453850905516027298167483341145182798 wg (3/2); gamma (1.5) ==> 0.8862269254527580136490837416705725913988 write ("sqrt(pi)/2 ==> ", sqrt(pi)/2); sqrt(pi)/2 ==> 0.8862269254527580136490837416705725913988 wp(10); -- precision 10, from 40 for x := 0 step 5 until 100 do << wg (1 + x/1000); wg (-1 - x/13); wp (8+floor(x/4)) >>; gamma (1) ==> 1 gamma (-1) ==> infinity -- precision 8, from 10 gamma (1.005) ==> 0.99713854 gamma ( - 1.3846154) ==> 2.7320314 -- precision 9, from 8 gamma (1.01) ==> 0.994325851 gamma ( - 1.76923077) ==> 2.89933597 -- precision 10, from 9 gamma (1.015) ==> 0.9915612888 gamma ( - 2.153846154) ==> - 2.919307224 -- precision 11, from 10 gamma (1.02) ==> 0.98884420326 gamma ( - 2.5384615385) ==> - 0.91247160689 -- precision 13, from 11 gamma (1.025) ==> 0.9861739631483 gamma ( - 2.923076923077) ==> - 2.407817725014 -- precision 14, from 13 gamma (1.03) ==> 0.98354995055382 gamma ( - 3.3076923076923) ==> 0.42665848359037 -- precision 15, from 14 gamma (1.035) ==> 0.980971560550586 gamma ( - 3.69230769230769) ==> 0.250121998146955 -- precision 16, from 15 gamma (1.04) ==> 0.9784382009142447 gamma ( - 4.076923076923077) ==> - 0.4868211588210416 -- precision 18, from 16 gamma (1.045) ==> 0.975949291822951489 gamma ( - 4.46153846153846154) ==> - 0.064315864992688343 -- precision 19, from 18 gamma (1.05) ==> 0.9735042655627756432 gamma ( - 4.846153846153846154) ==> - 0.07308480130893001268 -- precision 20, from 19 gamma (1.055) ==> 0.9711025662416699039 gamma ( - 5.2307692307692307692) ==> 0.026504298643014994546 -- precision 21, from 20 gamma (1.06) ==> 0.968743649511638364209 gamma ( - 5.61538461538461538462) ==> 0.00947958151841097406813 -- precision 23, from 21 gamma (1.065) ==> 0.96642698229883993296188 gamma (-6) ==> infinity -- precision 24, from 23 gamma (1.07) ==> 0.964152042541366448869499 gamma ( - 6.38461538461538461538462) ==> - 0.00224562428754660955672785 -- precision 25, from 24 gamma (1.075) ==> 0.9619183189344448686422338 gamma ( - 6.769230769230769230769231) ==> - 0.0014913552210721080976799 -- precision 26, from 25 gamma (1.08) ==> 0.95972531068282223532653464 gamma ( - 7.1538461538461538461538462) ==> 0.0009821350464726696000 373901 -- precision 28, from 26 gamma (1.085) ==> 0.9575725272601010524841387298 gamma ( - 7.538461538461538461538461538) ==> 0.00020813675527744238 31456972591 -- precision 29, from 28 gamma (1.09) ==> 0.95545948817480124076971838245 gamma ( - 7.9230769230769230769230769231) ==> 0.0003837247035972194 5415080451977 -- precision 30, from 29 gamma (1.095) ==> 0.95338572274293305256571929565 gamma ( - 8.30769230769230769230769230769) ==> - 0.0000487302698231202992921248405302 -- precision 31, from 30 gamma (1.1) ==> 0.9513507698668731836292487177265 gamma ( - 8.692307692307692307692307692308) ==> - 0.00002092711689267207910068115762837 -- precision 33, from 31 wg(1/1000000003); gamma (0.000000000999999997000000008999999973) ==> 1.00000000242278 433608752313084681e+9 off rounded; gamma(17/2); 2027025*sqrt(pi) ------------------ 256 gamma(-17/2); - 512*pi ------------------- 34459425*sqrt(pi) gamma(4); 6 gamma(0); infinity gamma(-4); infinity gamma(-17/3); 6*pi ------------------------ 17 17*sqrt(3)*gamma(----) 3 p := gamma(x**2) * gamma(x-y**gamma(y)) - (1/(gamma(4*(x-y)))); 2 gamma(y) gamma(x )*gamma( - y + x)*gamma(4*x - 4*y) - 1 p := -------------------------------------------------------- gamma(4*x - 4*y) y := 1/4; 1 y := --- 4 p; 2*gamma(1/4) 2 2 *x - 1 gamma(x )*gamma(4*x - 1)*gamma(---------------------) - 1 2*gamma(1/4) 2 ----------------------------------------------------------- gamma(4*x - 1) x := 3; x := 3 p; 2*gamma(1/4) 3*2 - 1 146313216000*gamma(---------------------) - 1 2*gamma(1/4) 2 ----------------------------------------------- 3628800 y := -3/8; - 3 y := ------ 8 p; (128*(2490343877896875*sqrt(pi)*gamma(( (8*pi)/(3*gamma(3/8)*sin((3*pi)/8)) 3*( - 3) (8*pi)/(gamma(3/8)*sin((3*pi)/8)) - 2 )/ (8*pi)/(3*gamma(3/8)*sin((3*pi)/8)) ( - 3) ) - 64))/( 7905853580625*sqrt(pi)) on rounded, complex; *** Domain mode rounded changed to complex-rounded precision 50; 33 p; -0.00000000058461000084165968732153392127582134179078414159599 + 3.7721251013859508830301986850709684723938237902095e-60*i off rounded, complex; *** Domain mode complex-rounded changed to complex clear y; p; gamma(y) 40320*gamma( - y + 3)*gamma( - 4*y + 12) - 1 ------------------------------------------------------ gamma( - 4*y + 12) % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % 3. Beta function. Not very interesting % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= algebraic procedure do!*one!*beta(x,y); write ("Beta (", x, ",", y, ") = ", beta(x,y)); do*one*beta do!*one!*beta(0,1); Beta (0,1) = infinity do!*one!*beta(2,-3); Beta (2,-3) = beta(2,-3) do!*one!*beta(3,2); 1 Beta (3,2) = ---- 12 do!*one!*beta(a+b,(c+d)**(b-a)); b (c + d) Beta (a + b,----------) = a (c + d) b (c + d) gamma(a + b)*gamma(----------) a (c + d) --------------------------------------------- a a b (c + d) *a + (c + d) *b + (c + d) gamma(------------------------------------) a (c + d) do!*one!*beta(-3,4); Beta (-3,4) = infinity do!*one!*beta(-3,2); Beta (-3,2) = beta(-3,2) do!*one!*beta(-3,-7.5); - 15 Beta (-3,-------) = infinity 2 do!*one!*beta((pi * 10), exp(5)); 5 5 gamma(e )*gamma(10*pi) Beta (10*pi,e ) = ------------------------ 5 gamma(e + 10*pi) on rounded; precision 30; 50 do!*one!*beta(0,1); Beta (0,1) = infinity do!*one!*beta(2,-3); Beta (2,-3) = beta(2,-3) do!*one!*beta(3,2); Beta (3,2) = 0.0833333333333333333333333333333 do!*one!*beta(a+b,(c+d)**(b-a)); b (c + d) Beta (a + b,----------) = a (c + d) b (c + d) gamma(a + b)*gamma(----------) a (c + d) --------------------------------------------- a a b (c + d) *a + (c + d) *b + (c + d) gamma(------------------------------------) a (c + d) do!*one!*beta(-3,4); Beta (-3,4) = infinity do!*one!*beta(-3,2); Beta (-3,2) = beta(-3,2) do!*one!*beta(-3,-7.5); Beta (-3, - 7.5) = infinity do!*one!*beta((pi * 10), exp(5)); Beta (31.4159265358979323846264338328,148.413159102576603421115580041 ) = 3.26162024071771351768890966259e-37 % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % 4. Pochhammer notation % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= off rounded; pochhammer(4,5); 6720 pochhammer(-4,5); pochhammer(-4,5) pochhammer(4,-5); pochhammer(4,-5) pochhammer(-4,-5); - 1 ------- 15120 pochhammer(17/2,12); 157783444591397625 -------------------- 4096 pochhammer(-17/2,12); - 17 pochhammer(-------,12) 2 pochhammer(1/3,14)*pochhammer(2/3,15); 2 1 pochhammer(---,15)*pochhammer(---,14) 3 3 q := pochhammer(1/5,11)*pochhammer(2/5,11)*pochhammer(3/5,11)* pochhammer(1-1/5,11)*pochhammer(1,11)*pochhammer(6/5,11)* pochhammer(70/50,11)*pochhammer(8/5,11)*pochhammer(9/5,11); 9 8 q := 39916800*pochhammer(---,11)*pochhammer(---,11) 5 5 7 6 4 *pochhammer(---,11)*pochhammer(---,11)*pochhammer(---,11) 5 5 5 3 2 1 *pochhammer(---,11)*pochhammer(---,11)*pochhammer(---,11) 5 5 5 on complex; pochhammer(a+b*i,c)*pochhammer(a-b*i,c); pochhammer(a - i*b,c)*pochhammer(a + i*b,c) a := 2; a := 2 b := 3; b := 3 c := 5; c := 5 pochhammer(a+b*i,c)*pochhammer(a-b*i,c); pochhammer(2 - 3*i,5)*pochhammer(2 + 3*i,5) off complex; on rounded; pochhammer(1/5,11)*pochhammer(2/5,11)*pochhammer(3/5,11)* pochhammer(1-1/5,11)*pochhammer(1,11)*pochhammer(6/5,11)* pochhammer(70/50,11)*pochhammer(8/5,11)*pochhammer(9/5,11); 2.38581559706937212593793381562e+67 q; 2.38581559706937212593793381562e+67 pochhammer(pi,floor (pi**8)); 2.47253079057195612235973919163e+33625 pochhammer(-pi,floor (pi**7)); 5.91750008101140123889628058136e+9185 pochhammer(1.5,floor (pi**8)); 1.88808885937650373473836451368e+33619 % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % 5. Digamma function % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= procedure do!*one!*psi(x); << precision (precision(0) + 4)$ write("Psi of ", x, " is ", psi(x) ) >> ; do*one*psi clear x, y; z := x * ((x+y)**2 + (x**y)); y 2 2 z := x*(x + x + 2*x*y + y ) off rounded; do!*one!*psi(3); 1 2*log(2) + psi(---) + psi(1) + 3 2 Psi of 3 is ---------------------------------- 2 do!*one!*psi(pi); Psi of pi is psi(pi) do!*one!*psi(1.005); 201 1 Psi of ----- is psi(-----) + 200 200 200 do!*one!*psi(1.995); 199 199*psi(-----) + 200 399 200 Psi of ----- is ---------------------- 200 199 do!*one!*psi(74); Psi of 74 is (269426837164032294756360054754050*log(2) 1 + 134713418582016147378180027377025*psi(---) 2 + 2138308231460573767907619482175*psi(1) + 667084944417653637854891458725877)/ 136851726813476721146087646859200 do!*one!*psi(-1/2); - 1 1 Psi of ------ is psi(---) + 2 2 2 do!*one!*psi(-3); Psi of -3 is infinity do!*one!*psi(z); y 2 2 y 3 2 2 Psi of x*(x + x + 2*x*y + y ) is psi(x *x + x + 2*x *y + x*y ) on rounded; precision 100; 62 do!*one!*psi(3); Psi of 3 is 0.9227843350984671393934879099175975689578406640600764011 9423276511513227322233532906305293670825325048537 do!*one!*psi(pi); Psi of 3.141592653589793238462643383279502884197169399375105820974944 59230781640628620899862803482534211706798214809 is 0.977213307 942006733292069486406182343640834609994325638009523286531810592477714 131730207565436292873435576949 do!*one!*psi(1.005); Psi of 1.005 is - 0.569020911344382831688295690950028809267968870263 46909469924377713720944057710319900770454501151058409 60011306543 do!*one!*psi(1.995); Psi of 1.995 is 0.419554603028108628651620272287243954692129976503405 96127576800237211077024989654159358528162438486780461 247996589393 do!*one!*psi(74); Psi of 74 is 4.297293118804667639106350362560367137919643935000983615 96464857079142072638532815275329153012799334519059885201 047316254 do!*one!*psi(-1/2); Psi of - 0.5 is 0.03648997397857652055902366700124443280684039533956 5892952872746128345029282945897851326282715415875401 36559070905154605166846 do!*one!*psi(-3); Psi of -3 is infinity do!*one!*psi(z); y 2 2 y 3 2 2 Psi of x*(x + x + 2*x*y + y ) is psi(x *x + x + 2*x *y + x*y ) precision 15; 132 x := 8/3; x := 2.66666666666667 y := 7/1000; y := 0.007 do!*one!*psi(z); Psi of 21.74768766103287773 is 3.056340330052438423 off rounded; clear x, y; df(psi(z), x); y 3 2 2 polygamma(1,x *x + x + 2*x *y + x*y ) y y 2 2 *(x *y + x + 3*x + 4*x*y + y ) df(df(psi(z), y),x); 2*y y 3 2 2 x *log(x)*polygamma(2,x *x + x + 2*x *y + x*y )*x*y 2*y y 3 2 2 + x *log(x)*polygamma(2,x *x + x + 2*x *y + x*y )*x y y 3 2 2 3 + 3*x *log(x)*polygamma(2,x *x + x + 2*x *y + x*y )*x y y 3 2 2 2 + 4*x *log(x)*polygamma(2,x *x + x + 2*x *y + x*y )*x *y y y 3 2 2 2 + x *log(x)*polygamma(2,x *x + x + 2*x *y + x*y )*x*y y y 3 2 2 + x *log(x)*polygamma(1,x *x + x + 2*x *y + x*y )*y y y 3 2 2 + x *log(x)*polygamma(1,x *x + x + 2*x *y + x*y ) y y 3 2 2 2 + 2*x *polygamma(2,x *x + x + 2*x *y + x*y )*x *y y y 3 2 2 2 + 2*x *polygamma(2,x *x + x + 2*x *y + x*y )*x y y 3 2 2 2 + 2*x *polygamma(2,x *x + x + 2*x *y + x*y )*x*y y y 3 2 2 + 2*x *polygamma(2,x *x + x + 2*x *y + x*y )*x*y y y 3 2 2 + x *polygamma(1,x *x + x + 2*x *y + x*y ) y 3 2 2 4 + 6*polygamma(2,x *x + x + 2*x *y + x*y )*x y 3 2 2 3 + 14*polygamma(2,x *x + x + 2*x *y + x*y )*x *y y 3 2 2 2 2 + 10*polygamma(2,x *x + x + 2*x *y + x*y )*x *y y 3 2 2 3 + 2*polygamma(2,x *x + x + 2*x *y + x*y )*x*y y 3 2 2 + 4*polygamma(1,x *x + x + 2*x *y + x*y )*x y 3 2 2 + 2*polygamma(1,x *x + x + 2*x *y + x*y )*y int(psi(z), z); y 3 2 2 log(gamma(x *x + x + 2*x *y + x*y )) on rounded; for k := 1 step 0.1 until 2 do do!*one!*psi(k); Psi of 1 is - 0.57721566490153286060651 Psi of 1.09999999999999999999999999 is - 0.423754940411076795168216226 Psi of 1.19999999999999999999999999483 is - 0.2890398965921882955472079690017 Psi of 1.2999999999999999999999999948304368 is - 0.16919088886679965563116117990224309 Psi of 1.39999999999999999999999999483043679294 is - 0.0613845445851161457306754873481741603465 Psi of 1.499999999999999999999999994830436792938162 is 0.036489973978 57652055902366216872537099062413 Psi of 1.5999999999999999999999999948304367929381615219 is 0.12604745 277347625190600271342271679469086378102 Psi of 1.69999999999999999999999999483043679293816152191799 is 0.2085 4787487349395667996417990930337500109911872963 Psi of 1.799999999999999999999999994830436792938161521917991479 is 0. 2849914332938615406087023630333301945809176983421812009 Psi of 1.8999999999999999999999999948304367929381615219179914789724 is 0.35618416116405971922472708037210055174041765918986424485621 Psi of 1.999999999999999999999999994830436792938161521917991478972404 77 is 0.422784335098467139393487906583570145998489083710421411 058968004 off rounded; % PSI_SIMP.TST F.J.Wright, 2 July 1993 on evallhseqp; factor psi; on rat, intstr, div; % for neater output % Do not try using "off mcd"! psi(x+m) - psi(x+m-1) = 1/(x+m-1); 1 1 -----------=----------- m + x - 1 m + x - 1 psi(x+2) - psi(x+1) + 2*psi(x) = 1/(x+1) + 2*psi(x); 1 1 2*psi(x) + -------=2*psi(x) + ------- x + 1 x + 1 psi(x+2) + 3*psi(x) = 4*psi(x) + 1/x + 1/(x+1); -1 -1 x + 2 x + 2 4*psi(x) + ---------=4*psi(x) + --------- x + 1 x + 1 psi(x + 1) = psi(x) + 1/x; -1 -1 psi(x) + x =psi(x) + x psi(x + 3/2) = psi(x + 1/2) + 1/(x + 1/2); 1 2 1 2 psi(x + ---) + ---------=psi(x + ---) + --------- 2 2*x + 1 2 2*x + 1 psi(x - 1/2) = psi(x + 1/2) - 1/(x - 1/2); 1 2 1 2 psi(x + ---) - ---------=psi(x + ---) - --------- 2 2*x - 1 2 2*x - 1 psi((x + 3a)/a); -1 1 2*(3*x + 8*x + 12) psi(---*x) + ---------------------- 2 2 x + 6*x + 8 psi(x/y + 3); -1 2 -1 y*(3*x + 2*x *y + 6*y) psi(x*y ) + -------------------------- 2 2 x + 3*x*y + 2*y off rat, intstr, div; on rational; psi(x+m) - psi(x+m-1) = 1/(x+m-1); 1 1 -----------=----------- m + x - 1 m + x - 1 psi(x+2) - psi(x+1) + 2*psi(x) = 1/(x+1) + 2*psi(x); 2*psi(x)*x + 2*psi(x) + 1 2*psi(x)*x + 2*psi(x) + 1 ---------------------------=--------------------------- x + 1 x + 1 psi(x+2) + 3*psi(x) = 4*psi(x) + 1/x + 1/(x+1); 2 4*psi(x)*x + 4*psi(x)*x + 2*x + 1 ------------------------------------ 2 x + x 2 4*psi(x)*x + 4*psi(x)*x + 2*x + 1 =------------------------------------ 2 x + x psi(x + 1) = psi(x) + 1/x; psi(x)*x + 1 psi(x)*x + 1 --------------=-------------- x x psi(x + 3/2) = psi(x + 1/2) + 1/(x + 1/2); 1 1 1 psi(x + ---)*x + ---*psi(x + ---) + 1 2 2 2 --------------------------------------- 1 x + --- 2 1 1 1 psi(x + ---)*x + ---*psi(x + ---) + 1 2 2 2 =--------------------------------------- 1 x + --- 2 psi(x - 1/2) = psi(x + 1/2) - 1/(x - 1/2); 1 1 1 psi(x + ---)*x - ---*psi(x + ---) - 1 2 2 2 --------------------------------------- 1 x - --- 2 1 1 1 psi(x + ---)*x - ---*psi(x + ---) - 1 2 2 2 =--------------------------------------- 1 x - --- 2 psi((x + 3a)/a); 1 2 2 8 psi(---*x)*x*(x + 6*x + 8) + 6*(x + 4*x + ---) 2 3 -------------------------------------------------- 2 x*(x + 6*x + 8) psi(x/y + 3); x 2 2 2 2 2 psi(---)*x*(x + 3*x*y + 2*y ) + 3*y*(x + 2*x*y + ---*y ) y 3 ------------------------------------------------------------ 2 2 x*(x + 3*x*y + 2*y ) off rational; % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % 6. Polygamma functions % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= procedure do!*one!*pg(n,x); write ("Polygamma (", n, ") of ", x, " is ", polygamma(n,x)); do*one*pg off rounded; do!*one!*pg(1,1/2); 1 1 Polygamma (1) of --- is polygamma(1,---) 2 2 do!*one!*pg(1,1); 2 pi Polygamma (1) of 1 is ----- 6 do!*one!*pg(1,3/2); 2 3 pi - 8 Polygamma (1) of --- is --------- 2 2 do!*one!*pg(1,1.005); 201 201 Polygamma (1) of ----- is polygamma(1,-----) 200 200 do!*one!*pg(1,1.995); 399 399 Polygamma (1) of ----- is polygamma(1,-----) 200 200 do!*one!*pg(1,1e-10); 1 1 Polygamma (1) of ------------- is polygamma(1,-------------) 10000000000 10000000000 do!*one!*pg(2,1.45); 29 29 Polygamma (2) of ---- is polygamma(2,----) 20 20 do!*one!*pg(3,1.99); 199 199 Polygamma (3) of ----- is polygamma(3,-----) 100 100 do!*one!*pg(4,-8.2); - 41 - 41 Polygamma (4) of ------- is polygamma(4,-------) 5 5 do!*one!*pg(5,0); Polygamma (5) of 0 is infinity do!*one!*pg(6,-5); Polygamma (6) of -5 is infinity do!*one!*pg(7,200); Polygamma (7) of 200 is ( 17726903568079021516229662669344581999513453462651624712682433450930396569413638060869523209196416456149981841411350319061723276968944831274405326841106087088051862580482604049179034792880997036306450715491134920388121743532178874504889199026667348488846246287054278592284199023548572837863642592371308428279823642858774260678182876022899303888558378920013392172441409952343098219076353830605697314645020537255667718925997023946826863459069360054494011515187931095594830347819766125471751832158677864778297413490571816668099535649328289203181141151677544880676008798702466775519536091686976761627732998715124407667063278175501392084867643186474774933465396962161891279302009565617687762042880000000000000000000000 8 *pi - 168202274324609919911886603179553832485831631970776943913214724306758011473978941132413094748296521617074167117661053498983406945306976674698602419479427329202365001144489134708765120356232374484239386020849061429951991522474455058548094378661234695692419462486402225369012574020278308089090741397551205530205566242796380021997272735088025054393694714504643986352458446128585149055841912367918690170304028682697422603414064799513570602284100682485151815951427876744662169487167539572896985974962832514966524381217035170993953423670212502913136819535576467936946089479260579835796147751358821073694464012485544356514018523618843496912692835966265011413341612149624871304965950289830955139972163546512706990894178649281 )/ 33237944190148165342930617505021091249087725242471796336279562720494493567650571364130356017243280855281215952646281848240731144316771558639509987827073913290097242338404882592210690236651869443074595091545877975727728269122835389696667248175001278416586711788226772360532873169153574070994329860696203303024669330360201738771592892542936194791046960475025110323327643660643309160768163432385682464959413507354376972986244419900300368985755050102176271590977370804240306902162061485259534685297520996459307650294822156252686629342490542255964639659395396651267516497567125204099130171913081428051999372590858264375743646579065110159126830974640203000247619304053546148691267935533164553830400000000000000000000000 on rounded; precision 100; 63 do!*one!*pg(1,1/2); Polygamma (1) of 0.5 is 4.9348022005446793094172454999380755676568497 03620395313206674688110022411209602621500886701859276116 do!*one!*pg(1,1); Polygamma (1) of 1 is 1.644934066848226436472415166646025189218949901 206798437735558229370007470403200873833628900619758705 do!*one!*pg(1,3/2); Polygamma (1) of 1.5 is 0.9348022005446793094172454999380755676568497 036203953132066746881100224112096026215008867018592761159 do!*one!*pg(1,1.005); Polygamma (1) of 1.005 is 1.63299415675568097525358697861102698996066 5426629281017423768306472667920740136690588658095268582677 do!*one!*pg(1,1.995); Polygamma (1) of 1.995 is 0.64696082864058235124763992710782769524368 64784225077824256922477813603546774633500406650324443323374 do!*one!*pg(1,1e-10); Polygamma (1) of 0.0000000001 is 1.0000000000000000000164493406660781 50558729660065732639916420752057756424780180942170336642993019363e+20 do!*one!*pg(2,1.45); Polygamma (2) of 1.45 is - 0.903837403076257688232309538529867192860 090308111446658245367522559171482627166928803699348863802546 do!*one!*pg(3,1.99); Polygamma (3) of 1.99 is 0.502907132416865320180510990862570197186151 6895613080444480934082942919777908933506993776280887706183 do!*one!*pg(4,-8.2); Polygamma (4) of - 8.2 is 74935.512595774270527120292307813135522959 05493466277917617629718610711391593937121099755096024434196 do!*one!*pg(5,0); Polygamma (5) of 0 is infinity do!*one!*pg(6,-5); Polygamma (6) of -5 is infinity do!*one!*pg(7,200); Polygamma (7) of 200 is 5.7240937253925583738370296815447272453468073 68195956154743567847385926644129753317064645661270233713e-14 off rounded; clear x; % Polygamma differentiation has already % been tried a bit in the psi section df(int(int(int(polygamma(3,x),x),x),x),x); polygamma(1,x) clear w, y, z; % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % 7. Zeta function % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= procedure do!*one!*zeta(n); write ("Zeta of ", n, " is ", zeta n); do*one*zeta off rounded; clear x, y, z; z := x * ((x+y)**5 + (x**y)); y 5 4 3 2 2 3 4 5 z := x*(x + x + 5*x *y + 10*x *y + 10*x *y + 5*x*y + y ) do!*one!*zeta(0); - 1 Zeta of 0 is ------ 2 for k := 4 step 2 until 35 do do!*one!*zeta(k); 4 pi Zeta of 4 is ----- 90 6 pi Zeta of 6 is ----- 945 8 pi Zeta of 8 is ------ 9450 10 pi Zeta of 10 is ------- 93555 12 691*pi Zeta of 12 is ----------- 638512875 14 2*pi Zeta of 14 is ---------- 18243225 16 3617*pi Zeta of 16 is -------------- 325641566250 18 43867*pi Zeta of 18 is ---------------- 38979295480125 20 174611*pi Zeta of 20 is ------------------ 1531329465290625 22 155366*pi Zeta of 22 is ------------------- 13447856940643125 24 236364091*pi Zeta of 24 is ----------------------- 201919571963756521875 26 1315862*pi Zeta of 26 is ---------------------- 11094481976030578125 28 6785560294*pi Zeta of 28 is -------------------------- 564653660170076273671875 30 6892673020804*pi Zeta of 30 is ------------------------------ 5660878804669082674070015625 Zeta of 32 is zeta(32) Zeta of 34 is zeta(34) do!*one!*zeta(-17/3); - 17 - 17 Zeta of ------- is zeta(-------) 3 3 do!*one!*zeta(190); Zeta of 190 is zeta(190) do!*one!*zeta(300); Zeta of 300 is zeta(300) do!*one!*zeta(0); - 1 Zeta of 0 is ------ 2 do!*one!*zeta(-44); Zeta of -44 is 0 on rounded; clear x, y; for k := 3 step 3 until 36 do << precision (31+k*3); do!*one!*zeta(k) >>; Zeta of 3 is 1.202056903159594285399738161511449990765 Zeta of 6 is 1.017343061984449139714517929790920527901817490033 Zeta of 9 is 1.002008392826082214417852769232412060485605851394888756 549 Zeta of 12 is 1.00024608655330804829863799804773967096041608845800340 4533040952133 Zeta of 15 is 1.00003058823630702049355172851064506258762794870685817 7506569932893332267156 Zeta of 18 is 1.00000381729326499983985646164462193973045469721895333 1143174429987630039542650045638 Zeta of 21 is 1.00000047693298678780646311671960437304596644669478493 7600207487376596839087898159833876638564 Zeta of 24 is 1.00000005960818905125947961244020793580122750391883730 2795864246972321724495355468544848206832825003614 Zeta of 27 is 1.00000000745071178983542949198100417060411945471903188 2565829993239578352147606271570867900837100313523764933 952 Zeta of 30 is 1.00000000093132743241966818287176473502121981356795513 6816185008613360441960672940496363503624604027929086312 123388047291 Zeta of 33 is 1.00000000011641550172700519775929738354563095165224717 2763593256517739947029124624567548673934974376008810870 912845774213829513369 Zeta of 36 is 1.00000000001455192189104198423592963224531842098380889 4124038069139542218571745865030220152998942329578185363 084791339999779092891491916899 precision 20; 139 do!*one!*zeta(-17/3); Zeta of - 5.6666666666666666667 is - 0.0018766468228592287697 do!*one!*zeta(z); y 5 4 3 2 2 3 4 5 Zeta of x*(x + x + 5*x *y + 10*x *y + 10*x *y + 5*x*y + y ) is y 6 5 4 2 3 3 2 4 5 zeta(x *x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + x*y ) y := 3; y := 3 x := pi; x := 3.1415926535897932385 do!*one!*zeta(z); Zeta of 27548.203250393209469 is 1 do!*one!*zeta(190); Zeta of 190 is 1.0 do!*one!*zeta(300); Zeta of 300 is 1 do!*one!*zeta(0); Zeta of 0 is - 0.5 do!*one!*zeta(-44); Zeta of -44 is 0 off rounded; % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % 8. Kummer functions % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= off rounded; t!*kummer!*a := { 1, 2.4, -1397/10 }$ t!*kummer!*b := { 0, 1, pi, -pi, 26 }$ for each a in t!*kummer!*a do for each b in t!*kummer!*a do for each z in t!*kummer!*a do << write "KummerM(", a, ",", b, ",", z, ") = ", kummerm(a,b,z); write "KummerU(", a, ",", b, ",", z, ") = ", kummeru(a,b,z) >>; KummerM(1,1,1) = e KummerU(1,1,1) = kummeru(1,1,1) 12 2/5 2 KummerM(1,1,----) = e *e 5 12 12 KummerU(1,1,----) = kummeru(1,1,----) 5 5 - 1397 1 KummerM(1,1,---------) = ------------ 10 7/10 139 e *e - 1397 - 1397 KummerU(1,1,---------) = kummeru(1,1,---------) 10 10 12 12 KummerM(1,----,1) = kummerm(1,----,1) 5 5 12 KummerU(1,----,1) = 5 2 2*pi 12 12 2*gamma(---)*sin(------)*(gamma(----)*e - kummerm(1,----,1)) 5 5 5 5 -------------------------------------------------------------- 12 12*pi 5*gamma(----)*sin(-------) 5 5 12 12 12 12 KummerM(1,----,----) = kummerm(1,----,----) 5 5 5 5 12 12 2 2*pi KummerU(1,----,----) = (gamma(---)*sin(------) 5 5 5 5 2/5 2/5 12 2 4/5 2/5 12 12 *(5*e *5 *gamma(----)*e - 12*2 *3 *kummerm(1,----,----))) 5 5 5 4/5 2/5 12 12*pi /(30*2 *3 *gamma(----)*sin(-------)) 5 5 12 - 1397 12 - 1397 KummerM(1,----,---------) = kummerm(1,----,---------) 5 10 5 10 12 - 1397 2 2*pi KummerU(1,----,---------) = (2*gamma(---)*sin(------)*( 5 10 5 5 7/10 2/5 12 - 1397 139 - 1397*e *1397 *kummerm(1,----,---------)*e 5 10 2/5 12 7/10 2/5 12 - 10*10 *gamma(----)))/(6985*e *1397 *gamma(----) 5 5 12*pi 139 *sin(-------)*e ) 5 - 1397 - 1397 KummerM(1,---------,1) = kummerm(1,---------,1) 10 10 - 1397 KummerU(1,---------,1) = 10 1397 - 1397 1397*pi 1397*gamma(------)*kummerm(1,---------,1)*sin(---------) + 10*e*pi 10 10 10 -------------------------------------------------------------------- 1417 1397*pi 10*gamma(------)*sin(---------) 10 10 - 1397 12 - 1397 12 KummerM(1,---------,----) = kummerm(1,---------,----) 10 5 10 5 - 1397 12 KummerU(1,---------,----) = ( 10 5 48689401707698071401637346296148709517635305466315546784369062686038955387470672028214421172492426257970894859785608060010798787187515790195560655355904 2/5 2/5 7/10 2 *e *2 *3 *e *pi + 20045966895736519148343810333245571434480905551887283539400953780162950579324387945234775543212890625 7/10 1397 *5 *gamma(------) 10 - 1397 12 1397*pi *kummerm(1,---------,----)*sin(---------))/( 10 5 10 143492962746861268062589909328887411843098822848155215027923792270314606867032125592231750488281250 7/10 1417 1397*pi *5 *gamma(------)*sin(---------)) 10 10 - 1397 - 1397 - 1397 - 1397 KummerM(1,---------,---------) = kummerm(1,---------,---------) 10 10 10 10 - 1397 - 1397 KummerU(1,---------,---------) = (1397*( 10 10 152159982579315173600331019038628152238885433889572221061155552992026662790488679161680024820756008730757778143037042693923932846444113056473446432154278794079710763645820340853851602689287025084815217783532480956517472235723731181252401977856952544597761994384769996636245679210281354768553595796003820698711669598622602752713703151447293557592795912091076094345980353443096817710446904761297306235868715700448137787799653833690839289133 7/10 *( - 1397) *pi + 10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 7/10 7/10 1397 - 1397 - 1397 *e *10 *gamma(------)*kummerm(1,---------,---------) 10 10 10 1397*pi 139 *sin(---------)*e ))/( 10 100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 7/10 7/10 1417 1397*pi 139 *e *10 *gamma(------)*sin(---------)*e ) 10 10 12 12 KummerM(----,1,1) = kummerm(----,1,1) 5 5 12 12 KummerU(----,1,1) = kummeru(----,1,1) 5 5 12 12 12 12 KummerM(----,1,----) = kummerm(----,1,----) 5 5 5 5 12 12 12 12 KummerU(----,1,----) = kummeru(----,1,----) 5 5 5 5 12 - 1397 12 - 1397 KummerM(----,1,---------) = kummerm(----,1,---------) 5 10 5 10 12 - 1397 12 - 1397 KummerU(----,1,---------) = kummeru(----,1,---------) 5 10 5 10 12 12 KummerM(----,----,1) = e 5 5 12 12 KummerU(----,----,1) = 5 5 2 - 2 2*pi 2*gamma(---)*kummerm(1,------,1)*sin(------) + 5*e*pi 5 5 5 ------------------------------------------------------- 12 12*pi 5*gamma(----)*sin(-------) 5 5 12 12 12 2/5 2 KummerM(----,----,----) = e *e 5 5 5 12 12 12 2/5 4/5 2/5 2 KummerU(----,----,----) = (6*e *2 *3 *e *pi 5 5 5 2/5 2 - 2 12 2*pi 4/5 + 5 *gamma(---)*kummerm(1,------,----)*sin(------))/(6*2 5 5 5 5 2/5 12 12*pi *3 *gamma(----)*sin(-------)) 5 5 12 12 - 1397 1 KummerM(----,----,---------) = ------------ 5 5 10 7/10 139 e *e 12 12 - 1397 7/10 2/5 2 KummerU(----,----,---------) = ( - 4*e *10 *gamma(---) 5 5 10 5 - 2 - 1397 2*pi 139 2/5 *kummerm(1,------,---------)*sin(------)*e + 1397*1397 *pi)/( 5 10 5 7/10 2/5 12 12*pi 139 1397*e *1397 *gamma(----)*sin(-------)*e ) 5 5 12 - 1397 12 - 1397 KummerM(----,---------,1) = kummerm(----,---------,1) 5 10 5 10 12 - 1397 KummerU(----,---------,1) = ( 5 10 1431 1431 1417 1417 10*gamma(------)*kummerm(------,------,1)*pi + 1397*gamma(------) 10 10 10 10 1397 12 12 - 1397 *gamma(------)*gamma(----)*kummerm(----,---------,1) 10 5 5 10 1397*pi 1431 1417 12 *sin(---------))/(10*gamma(------)*gamma(------)*gamma(----) 10 10 10 5 1397*pi *sin(---------)) 10 12 - 1397 12 12 - 1397 12 KummerM(----,---------,----) = kummerm(----,---------,----) 5 10 5 5 10 5 12 - 1397 12 KummerU(----,---------,----) = ( 5 10 5 20045966895736519148343810333245571434480905551887283539400953780162950579324387945234775543212890625 7/10 1417 1397 *5 *gamma(------)*gamma(------) 10 10 12 12 - 1397 12 1397*pi *gamma(----)*kummerm(----,---------,----)*sin(---------) + 5 5 10 5 10 48689401707698071401637346296148709517635305466315546784369062686038955387470672028214421172492426257970894859785608060010798787187515790195560655355904 2/5 7/10 1431 1431 1417 12 *2 *3 *gamma(------)*kummerm(------,------,----) 10 10 10 5 *pi)/( 143492962746861268062589909328887411843098822848155215027923792270314606867032125592231750488281250 7/10 1431 1417 *5 *gamma(------)*gamma(------) 10 10 12 1397*pi *gamma(----)*sin(---------)) 5 10 12 - 1397 - 1397 12 - 1397 - 1397 KummerM(----,---------,---------) = kummerm(----,---------,---------) 5 10 10 5 10 10 12 - 1397 - 1397 KummerU(----,---------,---------) = (1397*( 5 10 10 152159982579315173600331019038628152238885433889572221061155552992026662790488679161680024820756008730757778143037042693923932846444113056473446432154278794079710763645820340853851602689287025084815217783532480956517472235723731181252401977856952544597761994384769996636245679210281354768553595796003820698711669598622602752713703151447293557592795912091076094345980353443096817710446904761297306235868715700448137787799653833690839289133 7/10 1431 *( - 1397) *gamma(------) 10 1431 1417 - 1397 *kummerm(------,------,---------)*pi + 10 10 10 10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 7/10 1417 1397 12 *10 *gamma(------)*gamma(------)*gamma(----) 10 10 5 12 - 1397 - 1397 1397*pi *kummerm(----,---------,---------)*sin(---------)))/( 5 10 10 10 100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 7/10 1431 1417 12 1397*pi *10 *gamma(------)*gamma(------)*gamma(----)*sin(---------)) 10 10 5 10 - 1397 - 1397 KummerM(---------,1,1) = kummerm(---------,1,1) 10 10 - 1397 - 1397 KummerU(---------,1,1) = kummeru(---------,1,1) 10 10 - 1397 12 - 1397 12 KummerM(---------,1,----) = kummerm(---------,1,----) 10 5 10 5 - 1397 12 - 1397 12 KummerU(---------,1,----) = kummeru(---------,1,----) 10 5 10 5 - 1397 - 1397 - 1397 - 1397 KummerM(---------,1,---------) = kummerm(---------,1,---------) 10 10 10 10 - 1397 - 1397 - 1397 - 1397 KummerU(---------,1,---------) = kummeru(---------,1,---------) 10 10 10 10 - 1397 12 - 1397 12 KummerM(---------,----,1) = kummerm(---------,----,1) 10 5 10 5 - 1397 12 KummerU(---------,----,1) = ( 10 5 1411 - 1397 12 1411*pi - 7055*gamma(------)*kummerm(---------,----,1)*sin(---------)*pi 10 10 5 10 1397 12 2 - 2794*gamma(------)*gamma(----)*gamma(---) 10 5 5 - 1411 - 2 1397*pi 2*pi *kummerm(---------,------,1)*sin(---------)*sin(------))/(50 10 5 10 5 12 12*pi *gamma(----)*sin(-------)*pi) 5 5 - 1397 12 12 - 1397 12 12 KummerM(---------,----,----) = kummerm(---------,----,----) 10 5 5 10 5 5 - 1397 12 12 2/5 1397 KummerU(---------,----,----) = ( - 1397*5 *gamma(------) 10 5 5 10 12 2 - 1411 - 2 12 *gamma(----)*gamma(---)*kummerm(---------,------,----) 5 5 10 5 5 1397*pi 2*pi 4/5 2/5 1411 *sin(---------)*sin(------) - 8466*2 *3 *gamma(------) 10 5 10 - 1397 12 12 1411*pi 4/5 2/5 *kummerm(---------,----,----)*sin(---------)*pi)/(60*2 *3 10 5 5 10 12 12*pi *gamma(----)*sin(-------)*pi) 5 5 - 1397 12 - 1397 - 1397 12 - 1397 KummerM(---------,----,---------) = kummerm(---------,----,---------) 10 5 10 10 5 10 - 1397 12 - 1397 2/5 1411 KummerU(---------,----,---------) = ( - 1411*1397 *gamma(------) 10 5 10 10 - 1397 12 - 1397 1411*pi 2/5 *kummerm(---------,----,---------)*sin(---------)*pi + 4*10 10 5 10 10 1397 12 2 *gamma(------)*gamma(----)*gamma(---) 10 5 5 - 1411 - 2 - 1397 1397*pi 2*pi *kummerm(---------,------,---------)*sin(---------)*sin(------))/( 10 5 10 10 5 2/5 12 12*pi 10*1397 *gamma(----)*sin(-------)*pi) 5 5 - 1397 - 1397 KummerM(---------,---------,1) = e 10 10 - 1397 - 1397 KummerU(---------,---------,1) = 10 10 1397 1417 1417 1397*gamma(------)*(gamma(------)*e - kummerm(1,------,1)) 10 10 10 ------------------------------------------------------------ 1417 10*gamma(------) 10 - 1397 - 1397 12 2/5 2 KummerM(---------,---------,----) = e *e 10 10 5 - 1397 - 1397 12 1397 KummerU(---------,---------,----) = (1397*gamma(------)*( 10 10 5 10 71746481373430634031294954664443705921549411424077607513961896135157303433516062796115875244140625 2/5 7/10 1417 2 *e *5 *gamma(------)*e - 10 24344700853849035700818673148074354758817652733157773392184531343019477693735336014107210586246213128985447429892804030005399393593757895097780327677952 2/5 7/10 1417 12 *2 *3 *kummerm(1,------,----)))/( 10 5 717464813734306340312949546644437059215494114240776075139618961351573034335160627961158752441406250 7/10 1417 *5 *gamma(------)) 10 - 1397 - 1397 - 1397 1 KummerM(---------,---------,---------) = ------------ 10 10 10 7/10 139 e *e - 1397 - 1397 - 1397 1397 KummerU(---------,---------,---------) = (1397*gamma(------)*( - 10 10 10 10 212567495663303297519662433596963528677722951143732392822434307529861247918312684788866994674596144196868616065822748643411734186482425939893404665719527475329355936813211016172830688956933974043486859243594875896254908713306052460209605563066162704803073506155523685300835213856763052611669373327017337516100202429275776045541043302571869099957135889191233303801334553760006254341494325951532336811508595833526048489556116405666102486918801 7/10 7/10 *e *( - 1397) 1417 - 1397 139 *kummerm(1,------,---------)*e + 10 10 100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 7/10 1417 *10 *gamma(------)))/( 10 1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 7/10 7/10 1417 139 *e *10 *gamma(------)*e ) 10 on rounded; precision 30; 20 t!*k!*c := 7; t*k*c := 7 % To test each and every possible combination of % three arguments from t!*kummer!*b would take too % long, but we want the possibility of trying most % special cases. Compromise: test every seventh % possibility. for each a in t!*kummer!*b do for each b in t!*kummer!*b do for each z in t!*kummer!*b do << if t!*k!*c = 7 then << write "KummerM(", a, ",", b, ",", z, ") = ", kummerm(a,b,z); write "KummerU(", a, ",", b, ",", z, ") = ", kummeru(a,b,z); t!*k!*c := 0 >>; t!*k!*c := t!*k!*c + 1 >>; KummerM(0,0,0) = 1 KummerU(0,0,0) = kummeru(0,0,0) KummerM(0,1,3.14159265358979323846264338328) = 1 KummerU(0,1,3.14159265358979323846264338328) = kummeru(0,1, 3.14159265358979323846264338328) KummerM(0,3.14159265358979323846264338328,26) = 1 KummerU(0,3.14159265358979323846264338328,26) = kummeru(0, 3.14159265358979323846264338328,26) KummerM(0,26,1) = 1 KummerU(0,26,1) = kummeru(0,26,1) KummerM(1,0, - 3.14159265358979323846264338328) = kummerm(1,0, - 3.14159265358979323846264338328) KummerU(1,0, - 3.14159265358979323846264338328) = kummeru(1,0, - 3.14159265358979323846264338328) KummerM(1,3.14159265358979323846264338328,0) = 1 KummerU(1,3.14159265358979323846264338328,0) = - 0.466942206924259859983394813234 KummerM(1, - 3.14159265358979323846264338328,3.1415926535897932384626 4338328) = 2692.89480079631357528203659153 KummerU(1, - 3.14159265358979323846264338328,3.1415926535897932384626 4338328) = 0.12955419429695280640964490566 KummerM(1,26,26) = 7.74565667206271943920803547594 KummerU(1,26,26) = kummeru(1,26,26) KummerM(3.14159265358979323846264338328,1,1) = 10.2259987795328570162 092950355 KummerU(3.14159265358979323846264338328,1,1) = kummeru( 3.14159265358979323846264338328,1,1) KummerM(3.14159265358979323846264338328,3.141592653589793238462643383 28, - 3.14159265358979323846264338328) = 0.04321391826377224977441773 71717 KummerU(3.14159265358979323846264338328,3.141592653589793238462643383 28, - 3.14159265358979323846264338328) = ( - 0.137891580772667438638357954721* ( - 3.14159265358979323846264338328) **2.14159265358979323846264338328 + 1.22938452238615186004844111454)/ 2.14159265358979323846264338328 ( - 3.14159265358979323846264338328) KummerM(3.14159265358979323846264338328,26,0) = 1 KummerU(3.14159265358979323846264338328,26,0) = kummeru( 3.14159265358979323846264338328,26,0) KummerM( - 3.14159265358979323846264338328,0,3.1415926535897932384626 4338328) = kummerm( - 3.14159265358979323846264338328,0, 3.14159265358979323846264338328) KummerU( - 3.14159265358979323846264338328,0,3.1415926535897932384626 4338328) = kummeru( - 3.14159265358979323846264338328,0, 3.14159265358979323846264338328) KummerM( - 3.14159265358979323846264338328,1,26) = 6.1852222656472280 0173513559462e+5 KummerU( - 3.14159265358979323846264338328,1,26) = kummeru( - 3.14159265358979323846264338328,1,26) KummerM( - 3.14159265358979323846264338328, - 3.14159265358979323846264338328,1) = 2.718281828459045235360287471 35 KummerU( - 3.14159265358979323846264338328, - 3.14159265358979323846264338328,1) = 19.24195644060284656613462373 48 KummerM( - 3.14159265358979323846264338328,26, - 3.14159265358979323846264338328) = 1.42892253084220157246185820464 KummerU( - 3.14159265358979323846264338328,26, - 3.14159265358979323846264338328) = kummeru( - 3.14159265358979323846264338328,26, - 3.14159265358979323846264338328) KummerM(26,1,0) = 1 KummerU(26,1,0) = kummeru(26,1,0) KummerM(26,3.14159265358979323846264338328,3.141592653589793238462643 38328) = 3.91638029828702661403357541917e+5 KummerU(26,3.14159265358979323846264338328,3.141592653589793238462643 38328) = 4.18995913372628050180833640475e-32 KummerM(26, - 3.14159265358979323846264338328,26) = 2.194719452653224 19268333614674e+34 KummerU(26, - 3.14159265358979323846264338328,26) = 0 off rounded; clear x, y, z, t!*k!*c; df(df(kummerM(a,b,z),z),z); kummerm(4,5,z) ---------------- 2 df(kummerU(a,b,z),z); - 2*kummeru(3,4,z) z := ((x^2 + y)^5) + (x^(x+y)); x + y 10 8 6 2 4 3 2 4 5 z := x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y df(df(kummerM(a,b,z),y),x); 2*x + 2*y (3*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y ) 2 2*x + 2*y *log(x) *x + 3*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y ) 2*x + 2*y *log(x)*x + 3*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y ) x + y *log(x)*y + 30*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y ) 10 x + y *log(x)*x + 15*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y ) 9 x + y *log(x)*x + 120*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y ) 8 x + y *log(x)*x *y + 60*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y ) 7 x + y *log(x)*x *y + 180*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y ) 6 2 x + y *log(x)*x *y + 90*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y ) 5 2 x + y *log(x)*x *y + 120*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y ) 4 3 x + y *log(x)*x *y + 60*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y ) 3 3 x + y *log(x)*x *y + 30*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y ) 2 4 x + y *log(x)*x *y + 15*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y ) 4 x + y *log(x)*x*y + 15*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 9 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y )*x + x + y 15*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 8 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y )*x *y x + y + 60*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 7 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y )*x *y x + y + 60*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 6 2 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y )*x *y x + y + 90*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 5 2 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y )*x *y x + y + 90*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 4 3 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y )*x *y x + y + 60*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 3 3 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y )*x *y x + y + 60*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 2 4 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y )*x *y x + y + 15*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 4 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y )*x*y x + y + 15*x *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y )*y + x + y 4*x *kummerm(3,4, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y ) 2 x + y *log(x) *x + 4*x *kummerm(3,4, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y ) x + y *log(x)*x + 4*x *kummerm(3,4, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y ) x + y *log(x)*y + 4*x *kummerm(3,4, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y ) + 150 *kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y ) 18 *x + 1200*kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 16 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y )*x *y x + y 10 8 6 2 4 3 + 4200*kummerm(4,5,x + x + 5*x *y + 10*x *y + 10*x *y 2 4 5 14 2 + 5*x *y + y )*x *y + 8400*kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 12 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y )*x 3 *y + 10500*kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 10 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y )*x 4 *y + 8400*kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 8 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y )*x *y x + y 10 8 6 2 4 3 + 4200*kummerm(4,5,x + x + 5*x *y + 10*x *y + 10*x *y 2 4 5 6 6 + 5*x *y + y )*x *y + 1200*kummerm(4,5, x + y 10 8 6 2 4 3 2 4 5 4 7 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y )*x *y x + y 10 8 6 2 4 3 + 150*kummerm(4,5,x + x + 5*x *y + 10*x *y + 10*x *y 2 4 5 2 8 + 5*x *y + y )*x *y + 160*kummerm(3,4, x + y 10 8 6 2 4 3 2 4 5 8 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y )*x + x + y 10 8 6 2 4 3 480*kummerm(3,4,x + x + 5*x *y + 10*x *y + 10*x *y 2 4 5 6 + 5*x *y + y )*x *y + 480*kummerm(3,4, x + y 10 8 6 2 4 3 2 4 5 4 2 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y )*x *y x + y 10 8 6 2 4 3 + 160*kummerm(3,4,x + x + 5*x *y + 10*x *y + 10*x *y 2 4 5 2 3 + 5*x *y + y )*x *y )/(6*x) df(kummerU(a,b,z),x); (2*kummeru(3,4, x + y 10 8 6 2 4 3 2 4 5 x + x + 5*x *y + 10*x *y + 10*x *y + 5*x *y + y x + y x + y x + y 10 8 )*( - x *log(x)*x - x *x - x *y - 10*x - 40*x *y 6 2 4 3 2 4 - 60*x *y - 40*x *y - 10*x *y ))/x % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % 9. Bessel functions % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % Lengthy test of the Bessel functions. This isn't even % remotely exhaustive of the special cases -- though a % real person with lots of time could no doubt come up % with a better lot of tests than this automated rubbish. % Again, compromise by only actually doing one in five or % nine. If you want a really thorough test, you can % easily change this to get it; but it'll take hours to % run. clear p, q; hankel1(p,q); hankel1(p,q) r := df(ws,q); hankel1(p - 1,q)*q - hankel1(p,q)*p r := ------------------------------------- q on complex; r; (besselj(p - 1,q)*q - besselj(p,q)*p + i*bessely(p - 1,q)*q - i*bessely(p,q)*p)/q p:=3/8; 3 p := --- 8 r; - 5 5*pi 3*pi ( - 8*i*besselj(------,q)*cos(------)*sin(------)*q 8 8 8 - 5 5*pi 3*pi + 8*besselj(------,q)*sin(------)*sin(------)*q 8 8 8 - 3 5*pi + 3*i*besselj(------,q)*sin(------) 8 8 5 3*pi + 8*i*besselj(---,q)*sin(------)*q 8 8 3 3*pi 5*pi - 3*i*besselj(---,q)*cos(------)*sin(------) 8 8 8 3 5*pi 3*pi 5*pi - 3*besselj(---,q)*sin(------)*sin(------))/(8*sin(------) 8 8 8 8 3*pi *sin(------)*q) 8 q := pi; q := pi r; - 5 5*pi 3*pi ( - 8*i*besselj(------,pi)*cos(------)*sin(------)*pi 8 8 8 - 5 5*pi 3*pi + 8*besselj(------,pi)*sin(------)*sin(------)*pi 8 8 8 - 3 5*pi + 3*i*besselj(------,pi)*sin(------) 8 8 5 3*pi + 8*i*besselj(---,pi)*sin(------)*pi 8 8 3 3*pi 5*pi - 3*i*besselj(---,pi)*cos(------)*sin(------) 8 8 8 3 5*pi 3*pi 5*pi - 3*besselj(---,pi)*sin(------)*sin(------))/(8*sin(------) 8 8 8 8 3*pi *sin(------)*pi) 8 on rounded; *** Domain mode complex changed to complex-rounded r; - 0.119366207318921501826662822529 *hankel1(0.375,3.14159265358979323846264338328) + hankel1( - 0.625,3.14159265358979323846264338328) off complex, rounded; *** Domain mode complex-rounded changed to rounded df(df(besselj(pp,qq)+rr * hankel1(pp*2,qq) * bessely(pp-qq,qq),qq),qq); 2 (besselj(pp - 2,qq)*qq - 2*besselj(pp - 1,qq)*pp*qq 2 + besselj(pp - 1,qq)*qq + besselj(pp,qq)*pp + besselj(pp,qq)*pp 2 + bessely(pp - qq,qq)*hankel1(2*pp - 2,qq)*qq *rr - 4*bessely(pp - qq,qq)*hankel1(2*pp - 1,qq)*pp*qq*rr + bessely(pp - qq,qq)*hankel1(2*pp - 1,qq)*qq*rr 2 + 4*bessely(pp - qq,qq)*hankel1(2*pp,qq)*pp *rr + 2*bessely(pp - qq,qq)*hankel1(2*pp,qq)*pp*rr 2 + df(bessely(pp - qq,qq),qq,2)*hankel1(2*pp,qq)*qq *rr 2 + 2*df(bessely(pp - qq,qq),qq)*hankel1(2*pp - 1,qq)*qq *rr 2 - 4*df(bessely(pp - qq,qq),qq)*hankel1(2*pp,qq)*pp*qq*rr)/qq % Possible values for real args t!*bes!*vr := { 1, pi, -pi, 26 }$ % Possible values for real and imaginary parts of complex args t!*bes!*vc := { 0, 3, -41/2 }$ array s!*bes(4)$ s!*bes(1) := "BesselJ"$ s!*bes(2) := "BesselY"$ s!*bes(3) := "BesselI"$ s!*bes(4) := "BesselK"$ pre := 16; pre := 16 precision pre; 30 preord := 10**pre; preord := 10000000000000000 t!*b!*c := 3; t*b*c := 3 algebraic procedure do!*one!*bessel(s,n,z); (if s = 1 then besselj(n,z) else if s = 2 then bessely(n,z) else if s = 3 then besseli(n,z) else besselk(n,z)); do*one*bessel algebraic procedure pr!*bessel(s,n,z,k); << if t!*b!*c = k then << on rounded; bes1 := do!*one!*bessel(s,n,z); precision(pre+5); bes2 := do!*one!*bessel(s,n,z); if bes1 neq 0 then disc := floor abs(100*(bes2-bes1)*preord/bes1) else disc := 0; precision pre; write s!*bes(s), "(", n, ",", z, ") = ", bes1; if not numberp disc then << precom := !*complex; on complex; disc := disc; if not precom then off complex >>; if disc neq 0 then write " (discrepancy ", disc, "% of one s.f.)"; if numberp disc and disc > 200 then << write "***** WARNING Significant Inaccuracy."; write " Lower precision result:"; write " ", bes1; write " Higher precision result:"; precision(pre+5); write " ", bes2; precision pre >>; off rounded; t!*b!*c := 0 >>; t!*b!*c := t!*b!*c + 1 >>; pr*bessel % About to begin Bessel test. We have a list of possible % values, and we test every Bessel, with every value on the % list as both order and argument. Every Bessel is computed % twice, to different precisions (differing by 3), and any % discrepancy is reported. The value reported is the diff- % erence between the two computed values, expressed as a % percentage of the unit of the least significant displayed % digit. A discrepancy between 100% and 200% means that the % last digit of the displayed value was found to differ at % higher precision; values greater than 200% are cause for % concern. An ideal discrepancy would be between about 1% % and 20%; if the value is found to be zero, no discrepancy % is reported. off msg; for s := 1:4 do << write(" ... Testing ", s!*bes(s), " for real domains ... "); for each n in t!*bes!*vr do for each z in t!*bes!*vr do pr!*bessel(s, n, z, 5) >>; ... Testing BesselJ for real domains ... BesselJ(1, - 3.141592653589793) = - 0.2846153431797528 BesselJ(3.141592653589793,26) = - 0.006989220174690161 (discrepancy 5% of one s.f.) BesselJ(26,1) = 3.660826744416803e-35 ... Testing BesselY for real domains ... BesselY(1,3.141592653589793) = 0.358872916776719 BesselY(3.141592653589793, - 3.141592653589793) = ( 6.283185307179586 0.1545613960392598*( - 1.570796326794897) 3.141592653589793 - 4.829362563540275)/( - 1.570796326794897) BesselY( - 3.141592653589793,26) = - 0.1386083623177915 ... Testing BesselI for real domains ... BesselI(1,1) = 0.565159103992485 BesselI(3.141592653589793,3.141592653589793) = 1.011423335928613 BesselI( - 3.141592653589793, - 3.141592653589793) = -0.8856101155917482 + 0.4221616153281286*i BesselI(26,26) = 68397.86776155122 ... Testing BesselK for real domains ... BesselK(3.141592653589793,1) = 9.025908765806763 BesselK( - 3.141592653589793,3.141592653589793) = 0.1107526602738113 (discrepancy 1% of one s.f.) BesselK(26, - 3.141592653589793) = besselk(26, - 3.141592653589793) on complex; for s := 1:3 do << write (" ... Testing ", s!*bes(s), " for complex domains ... "); for each nr in t!*bes!*vc do for each ni in t!*bes!*vc do for each zr in t!*bes!*vc do for each zi in t!*bes!*vc do pr!*bessel(s, nr+ni*i, zr+zi*i, 9) >>; ... Testing BesselJ for complex domains ... BesselJ(0,-20.5 + 3.0*i) = 1.05389016561334 + 1.410918160335249*i BesselJ(3*i,-20.5 + 3.0*i) = 0.01225787392170983 + 0.01066256817009466*i BesselJ(-20.5*i,-20.5 + 3.0*i) = -6.607837931625446e+38 + 7.203284455482089e+38*i BesselJ(3,-20.5 + 3.0*i) = 1.568613483726435 - 0.7011991107137573*i BesselJ(3 + 3*i,-20.5 + 3.0*i) = 0.007904103001381543 - 0.006566520928092784*i BesselJ(3.0 - 20.5*i,-20.5 + 3.0*i) = -7.069920310202644e+37 - 1.753271554229047e+37*i BesselJ( - 20.5,-20.5 + 3.0*i) = 0.1758742246345278 - 0.332739860634916*i BesselJ(-20.5 + 3.0*i,-20.5 + 3.0*i) = 0.08815299110072903 - 0.1369698556512304*i BesselJ(-20.5 - 20.5*i,-20.5 + 3.0*i) = -5.364748129151297e+46 + 2.608178375230083e+47*i ... Testing BesselY for complex domains ... BesselY(0,-20.5 + 3.0*i) = -1.404746667469566 + 1.060048452645186*i BesselY(3*i,-20.5 + 3.0*i) = 0.4973091982659732 + 0.7985114801567726*i BesselY(-20.5*i,-20.5 + 3.0*i) = -7.203284455482089e+38 - 6.607837931625446e+38*i (discrepancy 10% of one s.f.) BesselY(3,-20.5 + 3.0*i) = 0.6963128100601111 + 1.576222640523309*i BesselY(3 + 3*i,-20.5 + 3.0*i) = -1.117333163968302 + 0.9789575771194796*i BesselY(3.0 - 20.5*i,-20.5 + 3.0*i) = 1.753271554229047e+37 - 7.069920310202644e+37*i BesselY( - 20.5,-20.5 + 3.0*i) = 0.2353954565395826 + 0.144691313932682*i BesselY(-20.5 + 3.0*i,-20.5 + 3.0*i) = -0.1527215881543493 + 0.2371137974094512*i BesselY(-20.5 - 20.5*i,-20.5 + 3.0*i) = -2.608178375230083e+47 - 5.364748129151297e+46*i ... Testing BesselI for complex domains ... BesselI(0,-20.5 + 3.0*i) = -6.891185608459107e+7 - 1.506065792318474e+7*i BesselI(3*i,-20.5 + 3.0*i) = -6879.511500081906 - 1745.250262122384*i BesselI(-20.5*i,-20.5 + 3.0*i) = 4.756052726395977e+40 - 1.836844915663626e+40*i BesselI(3,-20.5 + 3.0*i) = 5.56813934269752e+7 + 1.026348768124686e+7*i BesselI(3 + 3*i,-20.5 + 3.0*i) = 5917.619873410601 - 1360.623977225443*i BesselI(3.0 - 20.5*i,-20.5 + 3.0*i) = 1.125266362182384e+40 - 3.487618645712341e+39*i BesselI( - 20.5,-20.5 + 3.0*i) = -0.0005935206748136158 - 0.001045512507192035*i BesselI(-20.5 + 3.0*i,-20.5 + 3.0*i) = 0.0001367212493519011 - 0.00008899048591346719*i BesselI(-20.5 - 20.5*i,-20.5 + 3.0*i) = -4.705541285798881e+45 + 6.968045629188714e+45*i off complex; on msg; write (" ..."); ... write ("Bessel test complete."); Bessel test complete. end; (TIME: specfn 175193 188193) End of Lisp run after 175.22+14.41 seconds