Changing default RPM built to GCL.
[maxima.git] / tests / rtest_elliptic.mac
blobe85792fba6daaef58899f49cdc5a9cb5607042ac
1 /* Tests of Jacobi elliptic functions and elliptic integrals */
3 kill(all);
4 done$
6 /* derivatives */
7 diff(jacobi_sn(u,m),u);
8 jacobi_cn(u,m)*jacobi_dn(u,m);
10 diff(jacobi_sn(u,m),m);
11 jacobi_cn(u,m)*jacobi_dn(u,m)*(u-elliptic_e(asin(jacobi_sn(u,m)),m)/(1-m))
12   /(2*m)
13  +jacobi_cn(u,m)^2*jacobi_sn(u,m)/(2*(1-m));
15 diff(jacobi_cn(u,m),u);
16 -jacobi_dn(u,m)*jacobi_sn(u,m);
18 diff(jacobi_cn(u,m),m);
19 -(jacobi_dn(u,m)*jacobi_sn(u,m)*(u-elliptic_e(asin(jacobi_sn(u,m)),m)/(1-m))/(2*m))
20   -(jacobi_cn(u,m)*jacobi_sn(u,m)^2/(2*(1-m)));
22 diff(jacobi_dn(u,m),u);
23 -m*jacobi_cn(u,m)*jacobi_sn(u,m);
25 diff(jacobi_dn(u,m),m);
26 -(jacobi_cn(u,m)*jacobi_sn(u,m)*(u-elliptic_e(asin(jacobi_sn(u,m)),m)/(1-m))/2)
27   -(jacobi_dn(u,m)*'jacobi_sn(u,m)^2/(2*(1-m)));
29 diff(inverse_jacobi_sn(u,m),u);
30 1/(sqrt(1-u^2)*sqrt(1-m*u^2));
32 diff(inverse_jacobi_sn(u,m),m);
33 ((elliptic_e(asin(u),m)-(1-m)*elliptic_f(asin(u),m))/m-(u*sqrt(1-u^2)/sqrt(1-m*u^2)))/(1-m);
35 diff(inverse_jacobi_cn(u,m),u);
36 -(1/(sqrt(1-u^2)*sqrt(m*u^2-m+1)));
38 diff(inverse_jacobi_cn(u,m),m);
39 ((elliptic_e(asin(sqrt(1-u^2)),m)-(1-m)*elliptic_f(asin(sqrt(1-u^2)),m))/m
40    -(sqrt(1-u^2)*abs(u)/sqrt(1-m*(1-u^2))))/(1-m);
42 diff(inverse_jacobi_dn(u,m),u);
43 1/(sqrt(1-u^2)*sqrt(u^2+m-1));
45 diff(inverse_jacobi_dn(u,m),m);
46 ((elliptic_e(asin(sqrt(1-u^2)/sqrt(m)),m)-(1-m)*elliptic_f(asin(sqrt(1-u^2)/sqrt(m)),m))/m
47   -sqrt(1-(1-u^2)/m)*sqrt(1-u^2)/(sqrt(m)*abs(u)))/(1-m)
48 -(sqrt(1-u^2)/(2*m^(3/2)*sqrt(1-(1-u^2)/m)*abs(u)));
50 diff(elliptic_e(phi,m),phi);
51 sqrt(1-m*sin(phi)^2);
53 diff(elliptic_e(phi,m),m);
54 (elliptic_e(phi,m)-elliptic_f(phi,m))/(2*m);
56 diff(elliptic_f(phi,m),phi);
57 1/sqrt(1-m*sin(phi)^2);
59 diff(elliptic_f(phi,m),m);
60 ((elliptic_e(phi,m)-(1-m)*elliptic_f(phi,m))/m
61    -(cos(phi)*sin(phi)/sqrt(1-m*sin(phi)^2)))
62  /(2*(1-m));
64 diff(elliptic_pi(n,phi,m),n);
65 ((-(n*sqrt(1-m*sin(phi)^2)*sin(2*phi))/(2*(1-n*sin(phi)^2)))
66  +((m-n)*elliptic_f(phi,m))/n+elliptic_e(phi,m)
67  +((n^2-m)*elliptic_pi(n,phi,m))/n)
68  /(2*(m-n)*(n-1));
70 diff(elliptic_pi(n,phi,m),phi);
71 1/(sqrt(1-m*sin(phi)^2)*(1-n*sin(phi)^2));
73 diff(elliptic_pi(n,phi,m),m);
74 ((-(m*sin(2*phi))/(2*(m-1)*sqrt(1-m*sin(phi)^2)))
75  +elliptic_e(phi,m)/(m-1)+elliptic_pi(n,phi,m))
76  /(2*(n-m));
78 diff(elliptic_kc(m),m);
79 (elliptic_ec(m)-(1-m)*elliptic_kc(m))/(2*(1-m)*m);
81 diff(elliptic_ec(m),m);
82 (elliptic_ec(m)-elliptic_kc(m))/(2*m);
84 /* Integrals */
86 integrate(jacobi_sn(u,m),u); /* A&S 16.24.1 */
87 log(jacobi_dn(u,m)-sqrt(m)*jacobi_cn(u,m))/sqrt(m);
89 integrate(jacobi_cn(u,m),u); /* A&S 16.24.2 */
90 acos(jacobi_dn(u,m))/sqrt(m);
92 integrate(jacobi_dn(u,m),u); /* A&S 16.24.3 */
93 asin(jacobi_sn(u,m));
95 integrate(jacobi_cd(u,m),u); /* A&S 16.24.4 */
96 log(sqrt(m)*jacobi_sd(u,m)+jacobi_nd(u,m))/sqrt(m);
98 /* Use functions.wolfram.com 09.35.21.001.01, not A&S 16.24.5 */
99 integrate(jacobi_sd(u,m),u);
100 -(sqrt(1-m*jacobi_cd(u,m)^2)*jacobi_dn(u,m)*asin(sqrt(m)*jacobi_cd(u,m))
101   /((1-m)*sqrt(m)));
103 /* Use functions.wolfram.com 09.32.21.0001.01, not A&S 16.24.6 */
104 integrate(jacobi_nd(u,m),u); 
105 sqrt(1-jacobi_cd(u,m)^2)*acos(jacobi_cd(u,m))/((1-m)*jacobi_sd(u,m));
107 integrate(jacobi_dc(u,m),u); /* A&S 16.24.7 */
108 log(jacobi_sc(u,m)+jacobi_nc(u,m));
110 integrate(jacobi_nc(u,m),u); /* A&S 16.24.8 */
111 log(sqrt(1-m)*jacobi_sc(u,m)+jacobi_dc(u,m))/sqrt(1-m);
113 integrate(jacobi_sc(u,m),u); /* A&S 16.24.9 */
114 log(sqrt(1-m)*jacobi_nc(u,m)+jacobi_dc(u,m))/sqrt(1-m);
116 integrate(jacobi_ns(u,m),u); /* A&S 16.24.10 */
117 log(jacobi_ds(u,m)-jacobi_cs(u,m));
119 /* Use functions.wolfram.com 09.30.21.0001.01, not A&S 16.24.11 */
120 integrate(jacobi_ds(u,m),u);
121 log((1-jacobi_cn(u,m))/jacobi_sn(u,m));
123 integrate(jacobi_cs(u,m),u); /* A&S 16.24.12 */
124 log(jacobi_ns(u,m)-jacobi_ds(u,m));
126 /* Check the integrals and derivatives by confirming
128        f(x,m)-diff(integral(x,m),x),x) = constant
130   Look at Taylor expansion about zero, rather than messing about with 
131   elliptic function.  This was sufficient to find several errors  */
132 (te(f,n):=ratsimp(
133            taylor(f(x,m)
134            -diff(integrate(f(x,m),x),x),x,0,n)),
135 ti(f,n):=ratsimp(
136            taylor(integrate(f(x,m),x),x,0,n)
137             -integrate(taylor(f(x,m),x,0,n-1),x)),
138 td(f,n):=ratsimp(
139            taylor(diff(f(x,m),x),x,0,n)
140             -diff(taylor(f(x,m),x,0,n+1),x)),
141 /* Compare analytic and numerical integral */
142 ni(f,x1,x2,m):= block(
143   [x,n,I,In,Ia,key:1,eps:1.0e-14],
144   I:integrate(f(x,n),x),
145   In:quad_qag(f(x,m),x,x1,x2,key),
146   Ia:subst([x=float(x2),n=float(m)],I)-subst([x=float(x1),n=float(m)],I),
147   if ( abs(Ia-In[1]) < eps ) then
148     true
149   else
150     [Ia,In[1]]
152 done);
153 done;
155 te(jacobi_sn,4);
158 te(jacobi_cn,4);
161 te(jacobi_dn,4);
164 te(jacobi_cd,4);
167 assume(m>0,m<1); /* also ok for m<0 and m>0 */
168 [m > 0, m < 1];
169 te(jacobi_sd,4);
171 forget(m>0,m<1);
172 [m > 0, m < 1];
174 te(jacobi_nd,4);
177 te(jacobi_dc,4);
180 te(jacobi_nc,4);
183 te(jacobi_sc,4);
186 /* jacobi_ns, jacobi_ds and jacobi_cs are singular at x=0 
187    Compare numerical and analytic integrals for a single case.
190 ni(jacobi_ns,1,2,1/2);
191 true;
193 ni(jacobi_ds,1,2,1/2);
194 true;
196 ni(jacobi_cs,1,2,1/2);
197 true;
199 kill(te,ti,td,ni);
200 done;
202 /* Slightly modified version of test_table taken from rtest_expintegral.mac */
204 (test_table(func,table,eps) :=
205 block([badpoints : [],
206        abserr    : 0,
207        maxerr    : -1],
208   for entry in table do
209     block([z : entry[1],
210            result, answer],
211       z : expand(rectform(bfloat(entry[1]))),
212       result : rectform(apply(func, z)),
213       answer : expand(rectform(bfloat(entry[2]))),
214       abserr : abs(result-answer),
215       maxerr : max(maxerr,abserr),
216       if abserr > eps then
217         badpoints : cons ([z,result,answer,abserr],badpoints)
218     ),
219   if badpoints # [] then
220     cons(maxerr,badpoints)
221   else
222     badpoints
223 ),done);
224 done;
226 /* These test values come from http://getnet.net/~cherry/testrf.mac */
228 (mrf:[[[1,2,0],1.3110287771461b0],
229      [[%i,-%i,0],1.8540746773014b0],
230      [[0.5,1,0],1.8540746773014b0],
231      [[%i-1,%i,0],0.79612586584234b0-%i*(1.2138566698365b0)],
232      [[2,3,4],0.58408284167715b0],
233      [[%i,-%i,2],1.0441445654064b0],
234      [[%i-1,%i,1-%i],0.93912050218619b0-%i*(0.53296252018635b0)]],
235 done);
236 done;
238 test_table('carlson_rf, mrf, 1.5b-13);
241 (mrc:[[[0,1/4],bfloat(%pi)],
242       [[9/4,2],log(2b0)],
243       [[0,%i],(1-%i)*(1.1107207345396b0)],
244       [[-%i,%i],1.2260849569072b0-%i*(0.34471136988768b0)],
245       [[1/4,-2],log(2b0)/3],
246       [[%i,-1],0.77778596920447b0+%i*(0.19832484993429b0)],
247       [[0,1/4],%pi],
248       [[9/4,2],log(2)],
249       [[2,1],-log(sqrt(2)-1)],
250       [[-%i,%i],-log(sqrt(2)-1)/2+ %pi/4-%i*(log(sqrt(2)-1)/2+%pi/4)],
251       [[1/4,-2],log(2)/3],
252       [[%i,-1],sqrt(sqrt(2)/4-1/4)*atan(sqrt(sqrt(2)-1))-
253                sqrt(sqrt(2)/16+1/16)*log(-sqrt(2*sqrt(2)+2)+sqrt(2)+1)+
254                %i*(sqrt(sqrt(2)/4+1/4)*atan(sqrt(sqrt(2)-1))+sqrt(sqrt(2)/16-
255                1/16)*log(-sqrt(2*sqrt(2)+2)+sqrt(2)+1))],
256       [[0,1],%pi/2],
257       [[%i,%i+1],%pi/4+%i*log(sqrt(2)-1)/2]],
258 done);
259 done;
261 test_table('carlson_rc, mrc, 2.5b-14);
264 (mrj:[[[0,1,2,3],0.77688623778582b0],
265       [[2,3,4,5],0.14297579667157b0],
266       [[2,3,4,-1+%i],0.13613945827771b0-%i*(0.38207561624427b0)],
267       [[%i,-%i,0,2],1.6490011662711b0],
268       [[-1+%i,-1-%i,1,2],0.9414835884122b0],
269       [[%i,-%i,0,1-%i],1.8260115229009b0+%i*(1.2290661908643b0)],
270       [[-1+%i,-1-%i,1,-3+%i],-0.61127970812028b0-%i*(1.0684038390007b0)],
271       [[-1+%i,-2-%i,-%i,-1+%i],1.8249027393704b0-%i*(1.2218475784827b0)],
272       [[2,3,4,-0.5],0.24723819703052b0],
273       [[2,3,4,-5],-0.12711230042964b0]],
274 done);
275 done;
277 test_table('carlson_rj, mrj, 1b-13);
280 (mrd:[[[0,2,1],1.7972103521034b0],
281       [[2,3,4],0.16510527294261b0],
282       [[%i,-%i,2],0.6593385415422b0],
283       [[0,%i,-%i],1.270819627191b0+%i*(2.7811120159521b0)],
284       [[0,%i-1,%i],-1.8577235439239b0-%i*(0.96193450888839b0)],
285       [[-2-%i,-%i,-1+%i],1.8249027393704b0-%i*(1.2218475784827b0)]],
286 done);
287 done;
289 test_table('carlson_rd, mrd, 6e-14);
292 /* Some tests of the Jacobian elliptic functions.  
293  * Just some tests at random points, to verify that we are accurate.
294  * The reference values were obtained from Mathematica, but we could 
295  * also just compute the values using a much larger fpprec.
296  */
298 test_table('jacobi_sn,
299            [[[1b0+%i*1b0, .7b0], 1.134045971912365274954b0 + 0.352252346922494477621b0*%i],
300             [[1b0+%i*1b0, 2b0], 0.98613318109123804740b0 + 0.09521910819727230780b0*%i],
301             [[1b0+%i*1b0, 2b0+3b0*%i], 0.94467077879445294981b0 - 0.19586410083100945528b0* %i],
302             [[1.785063352082689082581887b0 *%i + 8.890858759528509578661528b-1,
303               9.434463451695984398149033b-1 - 1.476052973708684178844821b-1 * %i],
304              1.345233534539675700312281b0 - 7.599023926615176284214056b-2 * %i]],
305            1b-15);
308 test_table('jacobi_cn,
309            [[[1b0+%i*1b0, .7b0], 0.571496591371764254029b0 - 0.698989917271916772991b0*%i ],
310             [[1b0+%i*1b0, 2b0], 0.33759463268642431412b0 - 0.27814044708010806377b0*%i],
311             [[1b0+%i*1b0, 2b0+3b0*%i], -0.52142079485827170824b0 - 0.35485177134179601850b0*%i],
312             [[100b0, .7b0], 0.93004753815774770476196b0]], 
313            6b-14);
316 test_table('jacobi_dn,
317            [[[1b0+%i*1b0, .7b0], 0.62297154372331777630564880787568b0 - 0.448863598031509643267241389621738b0 *%i ],
318             [[1b0+%i*1b0, 2b0], 0.1913322443206041462495606602242b0 - 0.9815253294150083432282549919753b0 * %i],
319             [[1b0+%i*1b0, 2b0+3b0*%i], 0.6147387452173944656984656771134b0 - 1.4819401302071697495918834416787b0 * %i],
320             [[100b0, .7b0], 0.95157337933724324055428565654872978b0],
321             [[1.785063352082689082581887b0 *%i + 8.890858759528509578661528b-1,
322               9.434463451695984398149033b-1 - 1.476052973708684178844821b-1 * %i],
323              -8.617730683333292717095686b-1 *%i - 2.663978258141280808361839b-1]], 
324            2b-15);
327 /* These routines for cn and dn work well for small (<= 1?) values of
328  * u and m.  They have known issues for large real values of u.
329  */
330 (ascending_transform(u,m) :=
331   block([root_m : expand(rectform(sqrt(m))), mu, root_mu1, v],
332     mu : expand(rectform(4*root_m/(1+root_m)^2)),
333     root_mu1 : expand(rectform((1-root_m)/(1+root_m))),
334     v : expand(rectform(u/(1+root_mu1))),
335     [v, mu, root_mu1]),
336  elliptic_dn_ascending(u,m) :=
337   if is(abs(m-1) < 4*10^(-fpprec)) then sech(u)
338   else
339     block([v, mu, root_mu1, new],
340       [v, mu, root_mu1] : ascending_transform(u,m),
341       new : elliptic_dn_ascending(v, mu),
342       expand(rectform((1-root_mu1)/mu*(new^2 + root_mu1)/new))),
343  elliptic_cn_ascending(u,m) :=
344   if is(abs(m-1) < 4*10^(-fpprec)) then sech(u)
345   else
346     block([v, mu, root_mu1, new],
347       [v, mu, root_mu1] : ascending_transform(u,m),
348       new : elliptic_dn_ascending(v, mu),
349       expand(rectform((1+root_mu1)/mu*(new^2-root_mu1)/new))),
350  done);
351 done;
353 /* Test with random values for the argument and parameter. */
354 (test_random(n, eps, testf, truef) :=
355   block([badpoints : [], maxerr : -1],
356     for k : 1 thru n do
357       block([z : bfloat(1-2*random(1d0) + %i * (1-2*random(1d0))),
358              m : bfloat(1-2*random(1d0) + %i * (1-2*random(1d0))),
359              ans, expected, abserr],
360         ans : testf(z, m),
361         expected : truef(z, m),
362         abserr : abs(ans-expected),
363         maxerr : max(maxerr, abserr),
364         if abserr > eps then
365           badpoints : cons([[z,m], ans, expected, abserr], badpoints)),
366     if badpoints # [] then
367       cons(maxerr, badpoints)
368     else
369       badpoints),
370  done);
371 done;
373 test_random(50, 2b-15, 'jacobi_dn, 'elliptic_dn_ascending);
376 test_random(50, 2b-15, 'jacobi_cn, 'elliptic_cn_ascending);
379 /* Test elliptic_f by using the fact that
381 inverse_jacobi_sn(x,m) = elliptic_f(asin(x), m)
384 (test_ef(x,m) := jacobi_sn(elliptic_f(asin(x),m), m), id(x,m):=x, done);
385 done;
387 test_random(100, 6b-13, 'test_ef, 'id);
390 /* Test elliptic_kc These values are from 
392  * http://functions.wolfram.com/EllipticIntegrals/EllipticK/03/01/
393  */
395 block([oldfpprec : fpprec, fpprec:100],
396   test_table('elliptic_kc,
397              [[[1/2], 8*%pi^(3/2)/gamma(-1/4)^2],
398               [[17-12*sqrt(2)], 2*(2+sqrt(2))*%pi^(3/2)/gamma(-1/4)^2],
399               [[-1], gamma(1/4)^2/4/sqrt(2*%pi)]],
400              2b-100));
403 /* Some tests for specific values */
404 inverse_jacobi_sn(1,m);
405 elliptic_kc(m);
407 inverse_jacobi_sn(x,0);
408 asin(x);
410 inverse_jacobi_sn(x,1);
411 log(tan(asin(x)/2 + %pi/4));
413 inverse_jacobi_sd(1/sqrt(1-m),m);
414 elliptic_kc(m);
416 inverse_jacobi_ds(sqrt(1-m),m);
417 elliptic_kc(m);
419 /* elliptic_kc(1) is undefined */
420 errcatch(elliptic_kc(1));
422 errcatch(elliptic_kc(1.0));
424 errcatch(elliptic_kc(1.0b0));
427 /* Test noun/verb results from elliptic functions */
428 diff(inverse_jacobi_sn(x,0),x);
429 1/sqrt(1-x^2);
431 diff(elliptic_pi(4/3,phi,0),phi);
432 sec(phi)^2/(1-tan(phi)^2/3);
434 diff(elliptic_pi(3/4,phi,0),phi);
435 sec(phi)^2/(1+tan(phi)^2/4);
437 diff(elliptic_pi(1,phi,0),phi);
438 sec(phi)^2;
440 /* signbfloat:false provokes "Is 0 positive, negative, or zero?" */
441 (signbfloat:false, 
442  diff(elliptic_pi(1,phi,0),phi));
443 sec(phi)^2;
445 (reset (signbfloat), 0);
448 /* Special values */
449 elliptic_ec(1);
452 elliptic_ec(0);
453 %pi/2;
455 elliptic_ec(1/2);
456 gamma(3/4)^2/(2*sqrt(%pi))+%pi^(3/2)/(4*gamma(3/4)^2);
458 elliptic_ec(-1);
459 sqrt(2)*elliptic_ec(1/2);
461 /* Test periodicity of elliptic_e */
462 elliptic_e(x, 1);
463 2*round(x/%pi)+sin(x);
465 elliptic_e(3,1/3);
466 elliptic_e(3-%pi,1/3)+2*elliptic_ec(1/3);
468 /* Bug #2629: elliptic_kc(3.0) not accurate */
469 test_table('elliptic_kc, [[[3.0], elliptic_kc(3b0)]], 1b-15);
472 /* Bug #2630: inverse_jacobi_cn(-2.0, 3.0) generates an error */
473 test_table('inverse_jacobi_cn, [[[-2.0, 3.0], 2.002154760912212-3.202503914656527*%i]],
474   1b-15);
477 /* Bug #2615: Numeric evaluation of inverse Jacobi elliptic functions is wrong for some inputs 
479 is(abs(jacobi_dn(inverse_jacobi_dn(-2.0,3.0), 3.0) + 2) < 1d-14);
480 true;
482 /* verify that elliptical functions handle non-rectangular complex numbers
483  * mailing list 2016-07-14 "Jacobi elliptic functions, maxima 5.38.1"
484  */
486 /* list of functions here is everything that has a SIMP-[$%](ELLIPTIC|JACOBI)_FOO in sr/ellipt.lisp */
488 elliptic_fcns :
489   [ elliptic_e, elliptic_ec, elliptic_eu, elliptic_f, elliptic_kc, elliptic_pi,
490     inverse_jacobi_cd, inverse_jacobi_cn, inverse_jacobi_cs, inverse_jacobi_dc,
491     inverse_jacobi_dn, inverse_jacobi_ds, inverse_jacobi_nc, inverse_jacobi_nd,
492     inverse_jacobi_ns, inverse_jacobi_sc, inverse_jacobi_sd, inverse_jacobi_sn,
493     jacobi_am,
494     jacobi_cd, jacobi_cn, jacobi_cs, jacobi_dc, jacobi_dn, jacobi_ds,
495     jacobi_nc, jacobi_nd, jacobi_ns, jacobi_sc, jacobi_sd, jacobi_sn ];
496  */
498 /* code to generate test cases below
500 elliptic_fcns_1 : sublist (elliptic_fcns, lambda ([f], errcatch(f('a)) # []));
501 elliptic_fcns_2 : sublist (elliptic_fcns, lambda ([f], errcatch(f('a, 'b)) # []));
502 elliptic_fcns_3 : sublist (elliptic_fcns, lambda ([f], errcatch(f('a, 'b, 'c)) # []));
504 kill (u1, k1, n1);
506 with_stdout ("/tmp/foo.mac",
508   for f in elliptic_fcns_1
509     do block ([a : f(u1), b : f('rectform(u1))],
510         print ('if 'complex_floatp(a) and 'complex_floatp(b) and abs(a - b) < 1e-12 then true else a # b, ";"),
511         print (true, ";"),
512         print ("")),
513   
514   for f in elliptic_fcns_2
515     do block ([a : f(u1, k1), b : f('rectform(u1), k1)],
516         print ('if 'complex_floatp(a) and 'complex_floatp(b) and abs(a - b) < 1e-12 then true else a # b, ";"),
517         print (true, ";"),
518         print ("")),
519   
520   for f in elliptic_fcns_3
521     do block ([a : f(n1, u1, k1), b : f(n1, 'rectform(u1), k1)],
522         print ('if 'complex_floatp(a) and 'complex_floatp(b) and abs(a - b) < 1e-12 then true else a # b, ";"),
523         print (true, ";"),
524         print ("")));
525  */
527 (u1 : 0.5*(1.0 - 0.25*%i)^2,
528  k1 : 0.775,
529  n1 : -1.375,
530  complex_floatp(e) := floatnump(realpart(e)) and floatnump(imagpart(e)),
531  0);
534 if complex_floatp(elliptic_ec(u1))
535  and complex_floatp(elliptic_ec(rectform(u1)))
536  and (abs(elliptic_ec(rectform(u1)) - elliptic_ec(u1)) < 1.0E-12) then true
537  else elliptic_ec(u1) # elliptic_ec(rectform(u1)) ; 
538 true ; 
540 if complex_floatp(elliptic_kc(u1))
541  and complex_floatp(elliptic_kc(rectform(u1)))
542  and (abs(elliptic_kc(rectform(u1)) - elliptic_kc(u1)) < 1.0E-12) then true
543  else elliptic_kc(u1) # elliptic_kc(rectform(u1)) ; 
544 true ; 
546 if complex_floatp(elliptic_e(u1, k1))
547  and complex_floatp(elliptic_e(rectform(u1), k1))
548  and (abs(elliptic_e(rectform(u1), k1) - elliptic_e(u1, k1)) < 1.0E-12)
549  then true else elliptic_e(u1, k1) # elliptic_e(rectform(u1), k1) ; 
550 true ; 
552 if complex_floatp(elliptic_eu(u1, k1))
553  and complex_floatp(elliptic_eu(rectform(u1), k1))
554  and (abs(elliptic_eu(rectform(u1), k1) - elliptic_eu(u1, k1)) < 1.0E-12)
555  then true else elliptic_eu(u1, k1) # elliptic_eu(rectform(u1), k1) ; 
556 true ; 
558 if complex_floatp(elliptic_f(u1, k1))
559  and complex_floatp(elliptic_f(rectform(u1), k1))
560  and (abs(elliptic_f(rectform(u1), k1) - elliptic_f(u1, k1)) < 1.0E-12)
561  then true else elliptic_f(u1, k1) # elliptic_f(rectform(u1), k1) ; 
562 true ; 
564 if complex_floatp(inverse_jacobi_cd(u1, k1))
565  and complex_floatp(inverse_jacobi_cd(rectform(u1), k1))
566  and (abs(inverse_jacobi_cd(rectform(u1), k1) - inverse_jacobi_cd(u1, k1)) < 1.0E-12)
567  then true else inverse_jacobi_cd(u1, k1) # inverse_jacobi_cd(rectform(u1), k1)
568  ; 
569 true ; 
571 if complex_floatp(inverse_jacobi_cn(u1, k1))
572  and complex_floatp(inverse_jacobi_cn(rectform(u1), k1))
573  and (abs(inverse_jacobi_cn(rectform(u1), k1) - inverse_jacobi_cn(u1, k1)) < 1.0E-12)
574  then true else inverse_jacobi_cn(u1, k1) # inverse_jacobi_cn(rectform(u1), k1)
575  ; 
576 true ; 
578 if complex_floatp(inverse_jacobi_cs(u1, k1))
579  and complex_floatp(inverse_jacobi_cs(rectform(u1), k1))
580  and (abs(inverse_jacobi_cs(rectform(u1), k1) - inverse_jacobi_cs(u1, k1)) < 1.0E-12)
581  then true else inverse_jacobi_cs(u1, k1) # inverse_jacobi_cs(rectform(u1), k1)
582  ; 
583 true ; 
585 if complex_floatp(inverse_jacobi_dc(u1, k1))
586  and complex_floatp(inverse_jacobi_dc(rectform(u1), k1))
587  and (abs(inverse_jacobi_dc(rectform(u1), k1) - inverse_jacobi_dc(u1, k1)) < 1.0E-12)
588  then true else inverse_jacobi_dc(u1, k1) # inverse_jacobi_dc(rectform(u1), k1)
589  ; 
590 true ; 
592 if complex_floatp(inverse_jacobi_dn(u1, k1))
593  and complex_floatp(inverse_jacobi_dn(rectform(u1), k1))
594  and (abs(inverse_jacobi_dn(rectform(u1), k1) - inverse_jacobi_dn(u1, k1)) < 1.0E-12)
595  then true else inverse_jacobi_dn(u1, k1) # inverse_jacobi_dn(rectform(u1), k1)
596  ; 
597 true ; 
599 if complex_floatp(inverse_jacobi_ds(u1, k1))
600  and complex_floatp(inverse_jacobi_ds(rectform(u1), k1))
601  and (abs(inverse_jacobi_ds(rectform(u1), k1) - inverse_jacobi_ds(u1, k1)) < 1.0E-12)
602  then true else inverse_jacobi_ds(u1, k1) # inverse_jacobi_ds(rectform(u1), k1)
603  ; 
604 true ; 
606 if complex_floatp(inverse_jacobi_nc(u1, k1))
607  and complex_floatp(inverse_jacobi_nc(rectform(u1), k1))
608  and (abs(inverse_jacobi_nc(rectform(u1), k1) - inverse_jacobi_nc(u1, k1)) < 1.0E-12)
609  then true else inverse_jacobi_nc(u1, k1) # inverse_jacobi_nc(rectform(u1), k1)
610  ; 
611 true ; 
613 if complex_floatp(inverse_jacobi_nd(u1, k1))
614  and complex_floatp(inverse_jacobi_nd(rectform(u1), k1))
615  and (abs(inverse_jacobi_nd(rectform(u1), k1) - inverse_jacobi_nd(u1, k1)) < 1.0E-12)
616  then true else inverse_jacobi_nd(u1, k1) # inverse_jacobi_nd(rectform(u1), k1)
617  ; 
618 true ; 
620 if complex_floatp(inverse_jacobi_ns(u1, k1))
621  and complex_floatp(inverse_jacobi_ns(rectform(u1), k1))
622  and (abs(inverse_jacobi_ns(rectform(u1), k1) - inverse_jacobi_ns(u1, k1)) < 1.0E-12)
623  then true else inverse_jacobi_ns(u1, k1) # inverse_jacobi_ns(rectform(u1), k1)
624  ; 
625 true ; 
627 if complex_floatp(inverse_jacobi_sc(u1, k1))
628  and complex_floatp(inverse_jacobi_sc(rectform(u1), k1))
629  and (abs(inverse_jacobi_sc(rectform(u1), k1) - inverse_jacobi_sc(u1, k1)) < 1.0E-12)
630  then true else inverse_jacobi_sc(u1, k1) # inverse_jacobi_sc(rectform(u1), k1)
631  ; 
632 true ; 
634 if complex_floatp(inverse_jacobi_sd(u1, k1))
635  and complex_floatp(inverse_jacobi_sd(rectform(u1), k1))
636  and (abs(inverse_jacobi_sd(rectform(u1), k1) - inverse_jacobi_sd(u1, k1)) < 1.0E-12)
637  then true else inverse_jacobi_sd(u1, k1) # inverse_jacobi_sd(rectform(u1), k1)
638  ; 
639 true ; 
641 if complex_floatp(inverse_jacobi_sn(u1, k1))
642  and complex_floatp(inverse_jacobi_sn(rectform(u1), k1))
643  and (abs(inverse_jacobi_sn(rectform(u1), k1) - inverse_jacobi_sn(u1, k1)) < 1.0E-12)
644  then true else inverse_jacobi_sn(u1, k1) # inverse_jacobi_sn(rectform(u1), k1)
645  ; 
646 true ; 
648 if complex_floatp(jacobi_am(u1, k1))
649  and complex_floatp(jacobi_am(rectform(u1), k1))
650  and (abs(jacobi_am(rectform(u1), k1) - jacobi_am(u1, k1)) < 1.0E-12) then true
651  else jacobi_am(u1, k1) # jacobi_am(rectform(u1), k1) ; 
652 true ; 
654 if complex_floatp(jacobi_cd(u1, k1))
655  and complex_floatp(jacobi_cd(rectform(u1), k1))
656  and (abs(jacobi_cd(rectform(u1), k1) - jacobi_cd(u1, k1)) < 1.0E-12) then true
657  else jacobi_cd(u1, k1) # jacobi_cd(rectform(u1), k1) ; 
658 true ; 
660 if complex_floatp(jacobi_cn(u1, k1))
661  and complex_floatp(jacobi_cn(rectform(u1), k1))
662  and (abs(jacobi_cn(rectform(u1), k1) - jacobi_cn(u1, k1)) < 1.0E-12) then true
663  else jacobi_cn(u1, k1) # jacobi_cn(rectform(u1), k1) ; 
664 true ; 
666 if complex_floatp(jacobi_cs(u1, k1))
667  and complex_floatp(jacobi_cs(rectform(u1), k1))
668  and (abs(jacobi_cs(rectform(u1), k1) - jacobi_cs(u1, k1)) < 1.0E-12) then true
669  else jacobi_cs(u1, k1) # jacobi_cs(rectform(u1), k1) ; 
670 true ; 
672 if complex_floatp(jacobi_dc(u1, k1))
673  and complex_floatp(jacobi_dc(rectform(u1), k1))
674  and (abs(jacobi_dc(rectform(u1), k1) - jacobi_dc(u1, k1)) < 1.0E-12) then true
675  else jacobi_dc(u1, k1) # jacobi_dc(rectform(u1), k1) ; 
676 true ; 
678 if complex_floatp(jacobi_dn(u1, k1))
679  and complex_floatp(jacobi_dn(rectform(u1), k1))
680  and (abs(jacobi_dn(rectform(u1), k1) - jacobi_dn(u1, k1)) < 1.0E-12) then true
681  else jacobi_dn(u1, k1) # jacobi_dn(rectform(u1), k1) ; 
682 true ; 
684 if complex_floatp(jacobi_ds(u1, k1))
685  and complex_floatp(jacobi_ds(rectform(u1), k1))
686  and (abs(jacobi_ds(rectform(u1), k1) - jacobi_ds(u1, k1)) < 1.0E-12) then true
687  else jacobi_ds(u1, k1) # jacobi_ds(rectform(u1), k1) ; 
688 true ; 
690 if complex_floatp(jacobi_nc(u1, k1))
691  and complex_floatp(jacobi_nc(rectform(u1), k1))
692  and (abs(jacobi_nc(rectform(u1), k1) - jacobi_nc(u1, k1)) < 1.0E-12) then true
693  else jacobi_nc(u1, k1) # jacobi_nc(rectform(u1), k1) ; 
694 true ; 
696 if complex_floatp(jacobi_nd(u1, k1))
697  and complex_floatp(jacobi_nd(rectform(u1), k1))
698  and (abs(jacobi_nd(rectform(u1), k1) - jacobi_nd(u1, k1)) < 1.0E-12) then true
699  else jacobi_nd(u1, k1) # jacobi_nd(rectform(u1), k1) ; 
700 true ; 
702 if complex_floatp(jacobi_ns(u1, k1))
703  and complex_floatp(jacobi_ns(rectform(u1), k1))
704  and (abs(jacobi_ns(rectform(u1), k1) - jacobi_ns(u1, k1)) < 1.0E-12) then true
705  else jacobi_ns(u1, k1) # jacobi_ns(rectform(u1), k1) ; 
706 true ; 
708 if complex_floatp(jacobi_sc(u1, k1))
709  and complex_floatp(jacobi_sc(rectform(u1), k1))
710  and (abs(jacobi_sc(rectform(u1), k1) - jacobi_sc(u1, k1)) < 1.0E-12) then true
711  else jacobi_sc(u1, k1) # jacobi_sc(rectform(u1), k1) ; 
712 true ; 
714 if complex_floatp(jacobi_sd(u1, k1))
715  and complex_floatp(jacobi_sd(rectform(u1), k1))
716  and (abs(jacobi_sd(rectform(u1), k1) - jacobi_sd(u1, k1)) < 1.0E-12) then true
717  else jacobi_sd(u1, k1) # jacobi_sd(rectform(u1), k1) ; 
718 true ; 
720 if complex_floatp(jacobi_sn(u1, k1))
721  and complex_floatp(jacobi_sn(rectform(u1), k1))
722  and (abs(jacobi_sn(rectform(u1), k1) - jacobi_sn(u1, k1)) < 1.0E-12) then true
723  else jacobi_sn(u1, k1) # jacobi_sn(rectform(u1), k1) ; 
724 true ; 
726 if complex_floatp(elliptic_pi(n1, u1, k1))
727  and complex_floatp(elliptic_pi(n1, rectform(u1), k1))
728  and (abs(elliptic_pi(n1, rectform(u1), k1) - elliptic_pi(n1, u1, k1)) < 1.0E-12)
729  then true else elliptic_pi(n1, u1, k1) # elliptic_pi(n1, rectform(u1), k1) ; 
730 true ; 
733 with_stdout ("/tmp/bar.mac",
735   for f in elliptic_fcns_1
736     do block ([a : f(u1), b : f('rectform(u1))],
737         print ('if 'complex_bfloatp(a) and 'complex_bfloatp(b) and abs(a - b) < 1b-12 then true else a # b, ";"),
738         print (true, ";"),
739         print ("")),
740   
741   for f in elliptic_fcns_2
742     do block ([a : f(u1, k1), b : f('rectform(u1), k1)],
743         print ('if 'complex_bfloatp(a) and 'complex_bfloatp(b) and abs(a - b) < 1b-12 then true else a # b, ";"),
744         print (true, ";"),
745         print ("")),
746   
747   for f in elliptic_fcns_3
748     do block ([a : f(n1, u1, k1), b : f(n1, 'rectform(u1), k1)],
749         print ('if 'complex_bfloatp(a) and 'complex_bfloatp(b) and abs(a - b) < 1b-12 then true else a # b, ";"),
750         print (true, ";"),
751         print ("")));
752  */
754 (u1 : 0.5b0*(1.0b0 - 0.25b0*%i)^2,
755  k1 : 0.775b0,
756  n1 : -1.375b0,
757  complex_bfloatp(e) := bfloatp(realpart(e)) and bfloatp(imagpart(e)),
758  0);
761 if complex_bfloatp(elliptic_ec(u1))
762  and complex_bfloatp(elliptic_ec(rectform(u1)))
763  and (abs(elliptic_ec(rectform(u1)) - elliptic_ec(u1)) < 1.0b-12) then true
764  else elliptic_ec(u1) # elliptic_ec(rectform(u1)) ; 
765 true ; 
767 if complex_bfloatp(elliptic_kc(u1))
768  and complex_bfloatp(elliptic_kc(rectform(u1)))
769  and (abs(elliptic_kc(rectform(u1)) - elliptic_kc(u1)) < 1.0b-12) then true
770  else elliptic_kc(u1) # elliptic_kc(rectform(u1)) ; 
771 true ; 
773 if complex_bfloatp(elliptic_e(u1, k1))
774  and complex_bfloatp(elliptic_e(rectform(u1), k1))
775  and (abs(elliptic_e(rectform(u1), k1) - elliptic_e(u1, k1)) < 1.0b-12)
776  then true else elliptic_e(u1, k1) # elliptic_e(rectform(u1), k1) ; 
777 true ; 
779 if complex_bfloatp(elliptic_eu(u1, k1))
780  and complex_bfloatp(elliptic_eu(rectform(u1), k1))
781  and (abs(elliptic_eu(rectform(u1), k1) - elliptic_eu(u1, k1)) < 1.0b-12)
782  then true else elliptic_eu(u1, k1) # elliptic_eu(rectform(u1), k1) ; 
783 true ; 
785 if complex_bfloatp(elliptic_f(u1, k1))
786  and complex_bfloatp(elliptic_f(rectform(u1), k1))
787  and (abs(elliptic_f(rectform(u1), k1) - elliptic_f(u1, k1)) < 1.0b-12)
788  then true else elliptic_f(u1, k1) # elliptic_f(rectform(u1), k1) ; 
789 true ; 
791 if complex_bfloatp(inverse_jacobi_cd(u1, k1))
792  and complex_bfloatp(inverse_jacobi_cd(rectform(u1), k1))
793  and (abs(inverse_jacobi_cd(rectform(u1), k1) - inverse_jacobi_cd(u1, k1)) < 1.0b-12)
794  then true else inverse_jacobi_cd(u1, k1) # inverse_jacobi_cd(rectform(u1), k1)
795  ; 
796 true ; 
798 if complex_bfloatp(inverse_jacobi_cn(u1, k1))
799  and complex_bfloatp(inverse_jacobi_cn(rectform(u1), k1))
800  and (abs(inverse_jacobi_cn(rectform(u1), k1) - inverse_jacobi_cn(u1, k1)) < 1.0b-12)
801  then true else inverse_jacobi_cn(u1, k1) # inverse_jacobi_cn(rectform(u1), k1)
802  ; 
803 true ; 
805 if complex_bfloatp(inverse_jacobi_cs(u1, k1))
806  and complex_bfloatp(inverse_jacobi_cs(rectform(u1), k1))
807  and (abs(inverse_jacobi_cs(rectform(u1), k1) - inverse_jacobi_cs(u1, k1)) < 1.0b-12)
808  then true else inverse_jacobi_cs(u1, k1) # inverse_jacobi_cs(rectform(u1), k1)
809  ; 
810 true ; 
812 if complex_bfloatp(inverse_jacobi_dc(u1, k1))
813  and complex_bfloatp(inverse_jacobi_dc(rectform(u1), k1))
814  and (abs(inverse_jacobi_dc(rectform(u1), k1) - inverse_jacobi_dc(u1, k1)) < 1.0b-12)
815  then true else inverse_jacobi_dc(u1, k1) # inverse_jacobi_dc(rectform(u1), k1)
816  ; 
817 true ; 
819 if complex_bfloatp(inverse_jacobi_dn(u1, k1))
820  and complex_bfloatp(inverse_jacobi_dn(rectform(u1), k1))
821  and (abs(inverse_jacobi_dn(rectform(u1), k1) - inverse_jacobi_dn(u1, k1)) < 1.0b-12)
822  then true else inverse_jacobi_dn(u1, k1) # inverse_jacobi_dn(rectform(u1), k1)
823  ; 
824 true ; 
826 if complex_bfloatp(inverse_jacobi_ds(u1, k1))
827  and complex_bfloatp(inverse_jacobi_ds(rectform(u1), k1))
828  and (abs(inverse_jacobi_ds(rectform(u1), k1) - inverse_jacobi_ds(u1, k1)) < 1.0b-12)
829  then true else inverse_jacobi_ds(u1, k1) # inverse_jacobi_ds(rectform(u1), k1)
830  ; 
831 true ; 
833 if complex_bfloatp(inverse_jacobi_nc(u1, k1))
834  and complex_bfloatp(inverse_jacobi_nc(rectform(u1), k1))
835  and (abs(inverse_jacobi_nc(rectform(u1), k1) - inverse_jacobi_nc(u1, k1)) < 1.0b-12)
836  then true else inverse_jacobi_nc(u1, k1) # inverse_jacobi_nc(rectform(u1), k1)
837  ; 
838 true ; 
840 if complex_bfloatp(inverse_jacobi_nd(u1, k1))
841  and complex_bfloatp(inverse_jacobi_nd(rectform(u1), k1))
842  and (abs(inverse_jacobi_nd(rectform(u1), k1) - inverse_jacobi_nd(u1, k1)) < 1.0b-12)
843  then true else inverse_jacobi_nd(u1, k1) # inverse_jacobi_nd(rectform(u1), k1)
844  ; 
845 true ; 
847 if complex_bfloatp(inverse_jacobi_ns(u1, k1))
848  and complex_bfloatp(inverse_jacobi_ns(rectform(u1), k1))
849  and (abs(inverse_jacobi_ns(rectform(u1), k1) - inverse_jacobi_ns(u1, k1)) < 1.0b-12)
850  then true else inverse_jacobi_ns(u1, k1) # inverse_jacobi_ns(rectform(u1), k1)
851  ; 
852 true ; 
854 if complex_bfloatp(inverse_jacobi_sc(u1, k1))
855  and complex_bfloatp(inverse_jacobi_sc(rectform(u1), k1))
856  and (abs(inverse_jacobi_sc(rectform(u1), k1) - inverse_jacobi_sc(u1, k1)) < 1.0b-12)
857  then true else inverse_jacobi_sc(u1, k1) # inverse_jacobi_sc(rectform(u1), k1)
858  ; 
859 true ; 
861 if complex_bfloatp(inverse_jacobi_sd(u1, k1))
862  and complex_bfloatp(inverse_jacobi_sd(rectform(u1), k1))
863  and (abs(inverse_jacobi_sd(rectform(u1), k1) - inverse_jacobi_sd(u1, k1)) < 1.0b-12)
864  then true else inverse_jacobi_sd(u1, k1) # inverse_jacobi_sd(rectform(u1), k1)
865  ; 
866 true ; 
868 if complex_bfloatp(inverse_jacobi_sn(u1, k1))
869  and complex_bfloatp(inverse_jacobi_sn(rectform(u1), k1))
870  and (abs(inverse_jacobi_sn(rectform(u1), k1) - inverse_jacobi_sn(u1, k1)) < 1.0b-12)
871  then true else inverse_jacobi_sn(u1, k1) # inverse_jacobi_sn(rectform(u1), k1)
872  ; 
873 true ; 
875 if complex_bfloatp(jacobi_am(u1, k1))
876  and complex_bfloatp(jacobi_am(rectform(u1), k1))
877  and (abs(jacobi_am(rectform(u1), k1) - jacobi_am(u1, k1)) < 1.0b-12) then true
878  else jacobi_am(u1, k1) # jacobi_am(rectform(u1), k1) ; 
879 true ; 
881 if complex_bfloatp(jacobi_cd(u1, k1))
882  and complex_bfloatp(jacobi_cd(rectform(u1), k1))
883  and (abs(jacobi_cd(rectform(u1), k1) - jacobi_cd(u1, k1)) < 1.0b-12) then true
884  else jacobi_cd(u1, k1) # jacobi_cd(rectform(u1), k1) ; 
885 true ; 
887 if complex_bfloatp(jacobi_cn(u1, k1))
888  and complex_bfloatp(jacobi_cn(rectform(u1), k1))
889  and (abs(jacobi_cn(rectform(u1), k1) - jacobi_cn(u1, k1)) < 1.0b-12) then true
890  else jacobi_cn(u1, k1) # jacobi_cn(rectform(u1), k1) ; 
891 true ; 
893 if complex_bfloatp(jacobi_cs(u1, k1))
894  and complex_bfloatp(jacobi_cs(rectform(u1), k1))
895  and (abs(jacobi_cs(rectform(u1), k1) - jacobi_cs(u1, k1)) < 1.0b-12) then true
896  else jacobi_cs(u1, k1) # jacobi_cs(rectform(u1), k1) ; 
897 true ; 
899 if complex_bfloatp(jacobi_dc(u1, k1))
900  and complex_bfloatp(jacobi_dc(rectform(u1), k1))
901  and (abs(jacobi_dc(rectform(u1), k1) - jacobi_dc(u1, k1)) < 1.0b-12) then true
902  else jacobi_dc(u1, k1) # jacobi_dc(rectform(u1), k1) ; 
903 true ; 
905 if complex_bfloatp(jacobi_dn(u1, k1))
906  and complex_bfloatp(jacobi_dn(rectform(u1), k1))
907  and (abs(jacobi_dn(rectform(u1), k1) - jacobi_dn(u1, k1)) < 1.0b-12) then true
908  else jacobi_dn(u1, k1) # jacobi_dn(rectform(u1), k1) ; 
909 true ; 
911 if complex_bfloatp(jacobi_ds(u1, k1))
912  and complex_bfloatp(jacobi_ds(rectform(u1), k1))
913  and (abs(jacobi_ds(rectform(u1), k1) - jacobi_ds(u1, k1)) < 1.0b-12) then true
914  else jacobi_ds(u1, k1) # jacobi_ds(rectform(u1), k1) ; 
915 true ; 
917 if complex_bfloatp(jacobi_nc(u1, k1))
918  and complex_bfloatp(jacobi_nc(rectform(u1), k1))
919  and (abs(jacobi_nc(rectform(u1), k1) - jacobi_nc(u1, k1)) < 1.0b-12) then true
920  else jacobi_nc(u1, k1) # jacobi_nc(rectform(u1), k1) ; 
921 true ; 
923 if complex_bfloatp(jacobi_nd(u1, k1))
924  and complex_bfloatp(jacobi_nd(rectform(u1), k1))
925  and (abs(jacobi_nd(rectform(u1), k1) - jacobi_nd(u1, k1)) < 1.0b-12) then true
926  else jacobi_nd(u1, k1) # jacobi_nd(rectform(u1), k1) ; 
927 true ; 
929 if complex_bfloatp(jacobi_ns(u1, k1))
930  and complex_bfloatp(jacobi_ns(rectform(u1), k1))
931  and (abs(jacobi_ns(rectform(u1), k1) - jacobi_ns(u1, k1)) < 1.0b-12) then true
932  else jacobi_ns(u1, k1) # jacobi_ns(rectform(u1), k1) ; 
933 true ; 
935 if complex_bfloatp(jacobi_sc(u1, k1))
936  and complex_bfloatp(jacobi_sc(rectform(u1), k1))
937  and (abs(jacobi_sc(rectform(u1), k1) - jacobi_sc(u1, k1)) < 1.0b-12) then true
938  else jacobi_sc(u1, k1) # jacobi_sc(rectform(u1), k1) ; 
939 true ; 
941 if complex_bfloatp(jacobi_sd(u1, k1))
942  and complex_bfloatp(jacobi_sd(rectform(u1), k1))
943  and (abs(jacobi_sd(rectform(u1), k1) - jacobi_sd(u1, k1)) < 1.0b-12) then true
944  else jacobi_sd(u1, k1) # jacobi_sd(rectform(u1), k1) ; 
945 true ; 
947 if complex_bfloatp(jacobi_sn(u1, k1))
948  and complex_bfloatp(jacobi_sn(rectform(u1), k1))
949  and (abs(jacobi_sn(rectform(u1), k1) - jacobi_sn(u1, k1)) < 1.0b-12) then true
950  else jacobi_sn(u1, k1) # jacobi_sn(rectform(u1), k1) ; 
951 true ; 
953 if complex_bfloatp(elliptic_pi(n1, u1, k1))
954  and complex_bfloatp(elliptic_pi(n1, rectform(u1), k1))
955  and (abs(elliptic_pi(n1, rectform(u1), k1) - elliptic_pi(n1, u1, k1)) < 1.0b-12)
956  then true else elliptic_pi(n1, u1, k1) # elliptic_pi(n1, rectform(u1), k1) ; 
957 true ; 
959 /* Test derivatives at a few random points using numerical differentiation
960    
961    func  - a function
962    df - deriviative of func wrt p-th argument
963    p - derivative wrt p-th argument 
964    vars - variables in expression deriv
965    table - table of points to evaluate
966    eps - report errors > eps
967    delta - offset for numerical derivative
969 (test_deriv(func,df,vars,p,table,eps,delta) :=
970 block([badpoints : [],
971        abserr    : 0,
972        maxerr    : -1],
973   for %z in table do
974     block([%z0, %z1, deriv, nderiv],
975       %z1:makelist(if i=p then %z[i]+delta else %z[i],i,1,length(%z)),
976       %z0:makelist(if i=p then %z[i]-delta else %z[i],i,1,length(%z)),
977       nderiv: float((apply(func,%z1)-apply(func,%z0))/(2*delta)),
978       deriv : float(subst(maplist("=",vars,%z),df)),
979       abserr : abs(nderiv-deriv),
980       maxerr : max(maxerr,abserr),
981       if abserr > eps then
982         badpoints : cons ([%z,deriv,nderiv,abserr],badpoints)
983     ),
984   if badpoints # [] then
985     cons(maxerr,badpoints)
986   else
987     badpoints
988     ),done);
989 done;
991 /* test points for two-arg functions */
992 (l2: [[0.2, 0.3], [0.2, 0.5], [1.5, 0.7], [0.5, 0.99], [1.5, 0.8],
993     [1.57, 0.8], [1.5707, 0.8], [-1.0, 0.8], [2.0, 0.8]],done);
994 done;
996 /* test points for three-arg functions */
997 /* bug #3221 elliptic_pi() wrong for 2nd arg > float(%pi/2) */
998 (l3: [[0.1,0.2,0.3],[0.1,0.2,0.5],[0.2,1.5,0.7],[0.001,0.5,0.99],
999     [0.2,1.5,0.8],[0.2,1.57,0.8],[0.2,1.5707,0.8]],done);
1000 done;
1002 test_deriv('elliptic_f,diff(elliptic_f(z,m),z),[z,m],1,l2,2.0e-7,1.0e-8);
1005 test_deriv('elliptic_f,diff(elliptic_f(z,m),m),[z,m],2,l2,2.0e-7,1.0e-8);
1008 test_deriv('elliptic_e,diff(elliptic_e(z,m),z),[z,m],1,l2,2.0e-7,1.0e-8);
1011 test_deriv('elliptic_e,diff(elliptic_e(z,m),m),[z,m],2,l2,2.0e-7,1.0e-8);
1014 test_deriv('elliptic_pi,diff(elliptic_pi(n,z,m),n),[n,z,m],1,l3,2.0e-7,1.0e-8);
1017 test_deriv('elliptic_pi,diff(elliptic_pi(n,z,m),z),[n,z,m],2,l3,2.0e-7,1.0e-8);
1020 test_deriv('elliptic_pi,diff(elliptic_pi(n,z,m),m),[n,z,m],3,l3,2.0e-7,1.0e-8);
1023 closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
1024 closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
1027  * Some additional tests for complex values. Use the fact that
1028  * elliptic_pi(0, phi, m) = elliptic_f(phi,m);
1029  */
1030 closeto(elliptic_pi(0, 0.5, 0.25)- elliptic_f(0.5, 0.25), 1e-15);
1031 true;
1032 closeto(elliptic_pi(0, 0.5*%i, 0.25)- elliptic_f(0.5*%i, 0.25), 1e-15);
1033 true;
1035 closeto(elliptic_pi(0, 0.5b0, 0.25b0)- elliptic_f(0.5b0, 0.25b0), 1e-15);
1036 true;
1037 closeto(elliptic_pi(0, 0.5b0*%i, 0.25b0)- elliptic_f(0.5b0*%i, 0.25b0), 1e-15);
1038 true;
1040 /* Test for Bug 3221 */
1042 closeto(elliptic_pi(0.2,1.57,0.8) - 2.5770799919605668, 1e-14);
1043 true;
1044 closeto(elliptic_pi(0.2,1.58,0.8) - 2.60502920656026151, 1e-14);
1045 true;
1046 closeto(elliptic_pi(0.20,2.00,0.80) - 3.6543835090829025, 7.47e-11);
1047 true;
1049 (oldfpprec:fpprec,fpprec:32);
1051 closeto(elliptic_pi(0.2b0,1.57b0,0.8b0) - 2.5770799919605668058849721013196b0, 3.09b-32);
1052 true;
1053 closeto(elliptic_pi(0.2b0,1.58b0,0.8b0) - 2.6050292065602615145481924917132b0, 1.24b-32);
1054 true;
1055 closeto(elliptic_pi(0.2b0,2.0b0,0.8b0) - 3.6543835090082902501670624830938b0, 1b-32);
1056 true;
1058 (fpprec:oldfpprec,kill(l2,l3,test_deriv,oldfpprec),done);
1059 done;