tests/rtest_elliptic.mac

   1 /* Tests of Jacobi elliptic functions and elliptic integrals */
   2
   3 kill(all);
   4 done$
   5
   6 /* derivatives */
   7 diff(jacobi_sn(u,m),u);
   8 jacobi_cn(u,m)*jacobi_dn(u,m);
   9
  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));
  14
  15 diff(jacobi_cn(u,m),u);
  16 -jacobi_dn(u,m)*jacobi_sn(u,m);
  17
  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)));
  21
  22 diff(jacobi_dn(u,m),u);
  23 -m*jacobi_cn(u,m)*jacobi_sn(u,m);
  24
  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)));
  28
  29 diff(inverse_jacobi_sn(u,m),u);
  30 1/(sqrt(1-u^2)*sqrt(1-m*u^2));
  31
  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);
  34
  35 diff(inverse_jacobi_cn(u,m),u);
  36 -(1/(sqrt(1-u^2)*sqrt(m*u^2-m+1)));
  37
  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);
  41
  42 diff(inverse_jacobi_dn(u,m),u);
  43 1/(sqrt(1-u^2)*sqrt(u^2+m-1));
  44
  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)));
  49
  50 diff(elliptic_e(phi,m),phi);
  51 sqrt(1-m*sin(phi)^2);
  52
  53 diff(elliptic_e(phi,m),m);
  54 (elliptic_e(phi,m)-elliptic_f(phi,m))/(2*m);
  55
  56 diff(elliptic_f(phi,m),phi);
  57 1/sqrt(1-m*sin(phi)^2);
  58
  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));
  63
  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));
  69
  70 diff(elliptic_pi(n,phi,m),phi);
  71 1/(sqrt(1-m*sin(phi)^2)*(1-n*sin(phi)^2));
  72
  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));
  77
  78 diff(elliptic_kc(m),m);
  79 (elliptic_ec(m)-(1-m)*elliptic_kc(m))/(2*(1-m)*m);
  80
  81 diff(elliptic_ec(m),m);
  82 (elliptic_ec(m)-elliptic_kc(m))/(2*m);
  83
  84 /* Integrals */
  85
  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);
  88
  89 integrate(jacobi_cn(u,m),u); /* A&S 16.24.2 */
  90 acos(jacobi_dn(u,m))/sqrt(m);
  91
  92 integrate(jacobi_dn(u,m),u); /* A&S 16.24.3 */
  93 asin(jacobi_sn(u,m));
  94
  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);
  97
  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)));
 102
 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));
 106
 107 integrate(jacobi_dc(u,m),u); /* A&S 16.24.7 */
 108 log(jacobi_sc(u,m)+jacobi_nc(u,m));
 109
 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);
 112
 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);
 115
 116 integrate(jacobi_ns(u,m),u); /* A&S 16.24.10 */
 117 log(jacobi_ds(u,m)-jacobi_cs(u,m));
 118
 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));
 122
 123 integrate(jacobi_cs(u,m),u); /* A&S 16.24.12 */
 124 log(jacobi_ns(u,m)-jacobi_ds(u,m));
 125
 126 /* Check the integrals and derivatives by confirming
 127
 128        f(x,m)-diff(integral(x,m),x),x) = constant
 129
 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]]
 151 ),
 152 done);
 153 done;
 154
 155 te(jacobi_sn,4);
 156 0;
 157
 158 te(jacobi_cn,4);
 159 0;
 160
 161 te(jacobi_dn,4);
 162 0;
 163
 164 te(jacobi_cd,4);
 165 0;
 166
 167 assume(m>0,m<1); /* also ok for m<0 and m>0 */
 168 [m > 0, m < 1];
 169 te(jacobi_sd,4);
 170 0;
 171 forget(m>0,m<1);
 172 [m > 0, m < 1];
 173
 174 te(jacobi_nd,4);
 175 0;
 176
 177 te(jacobi_dc,4);
 178 0;
 179
 180 te(jacobi_nc,4);
 181 0;
 182
 183 te(jacobi_sc,4);
 184 0;
 185
 186 /* jacobi_ns, jacobi_ds and jacobi_cs are singular at x=0
 187    Compare numerical and analytic integrals for a single case.
 188 */
 189
 190 ni(jacobi_ns,1,2,1/2);
 191 true;
 192
 193 ni(jacobi_ds,1,2,1/2);
 194 true;
 195
 196 ni(jacobi_cs,1,2,1/2);
 197 true;
 198
 199 kill(te,ti,td,ni);
 200 done;
 201
 202 /* Slightly modified version of test_table taken from rtest_expintegral.mac */
 203
 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,reverse(badpoints))
 221   else
 222     badpoints
 223 ),done);
 224 done;
 225
 226 /* These test values come from http://getnet.net/~cherry/testrf.mac */
 227
 228 /*
 229  * rf(1,2,0) = (gamma(1/4)^2/4/sqrt(2*%pi)
 230  * rf(%i,-%i,0) = 1/2*integrate(1/sqrt(t^2+1)/sqrt(t),t,0,inf) = beta(1/4,1/4)/4;
 231  * rf(1/2,1,0) = 1/sqrt(1/2)*rf(1,2,0) (See https://dlmf.nist.gov/19.20.E1)
 232  */
 233 (mrf:[[[1,2,0], gamma(1/4)^2/4/sqrt(2*%pi)],
 234      [[0.5,1,0], gamma(1/4)^2/(4*sqrt(%pi))]],
 235 done);
 236 done;
 237
 238 (mrf2:[[[%i,-%i,0], beta(1/4,1/4)/4],
 239      [[%i-1,%i,0],0.79612586584234b0-%i*(1.2138566698365b0)],
 240      [[%i,-%i,2],1.0441445654064b0],
 241      [[2,3,4],0.58408284167715b0],
 242      [[%i-1,%i,1-%i],0.93912050218619b0-%i*(0.53296252018635b0)]],
 243 done);
 244 done;
 245
 246 test_table('carlson_rf, mrf, 8.3267b-17);
 247 [];
 248
 249 test_table('carlson_rf, mrf2, 1.2102b-8);
 250 [];
 251
 252 /*
 253  * rc(0,1/4) = 1/2*integrate(1/sqrt(t)/(t+1/4), t, 0, inf)
 254  *           = %pi
 255  *
 256  * rc(9/4,2) = 1/2*integrate(1/sqrt(t+9/4)/(t+2), t, 0, inf)
 257  *           = log(2)
 258  * After doing a logcontract.
 259  *
 260  * rc(0,%i) = 1/2*integrate(1/sqrt(t)/(t+%i), t, 0, inf);
 261  *          = (1-%i)*%pi/2^(3/2)
 262  *
 263  * rc(2,1) = 1/2*integrate(1/sqrt(t+2)/(t+1), t, 0, inf)
 264  *         = (log(sqrt(2)+1)-log(sqrt(2)-1))/2
 265  *         = -log(sqrt(2)-1)
 266  *
 267  * After doing a logcontract, ratsimp/algebraic, and logcontract with
 268  *  logconcoeffp, and a sqrtdenest
 269  *
 270  * rc(0,1) = 1/2*integrate(1/sqrt(t)/(t+1), t, 0, inf)
 271  *         = %pi/2
 272  * rc(%i, %i+1) = 1/2*integrate(1/sqrt(t+%i)/(t+%i+1), t, 0, inf)
 273  *     = (%pi-2*atan((-1)^(1/4)))/2
 274  *     = %pi/4+%i/2*log(sqrt(2)-1)
 275  *
 276  * After applying rectform, ratsimp, logcontract, then another
 277  * logcontract with logconcoeffp set to featurep(m) or ratnump(m).
 278  *
 279  */
 280 (mrc:[[[0b0,1/4],bfloat(%pi)],
 281       [[9/4,2b0],log(2b0)],
 282       [[0b0,%i],(1-%i)*%pi/2^(3/2)],
 283       [[-%i,%i],1.2260849569072b0-%i*(0.34471136988768b0)],
 284       [[1/4,-2],log(2b0)/3],
 285       [[%i,-1],0.77778596920447b0+%i*(0.19832484993429b0)],
 286       [[0,1/4],%pi],
 287       [[9/4,2],log(2)],
 288       [[2,1],-log(sqrt(2)-1)],
 289       [[-%i,%i],-log(sqrt(2)-1)/2+ %pi/4-%i*(log(sqrt(2)-1)/2+%pi/4)],
 290       [[1/4,-2],log(2)/3],
 291       [[%i,-1],sqrt(sqrt(2)/4-1/4)*atan(sqrt(sqrt(2)-1))-
 292                sqrt(sqrt(2)/16+1/16)*log(-sqrt(2*sqrt(2)+2)+sqrt(2)+1)+
 293                %i*(sqrt(sqrt(2)/4+1/4)*atan(sqrt(sqrt(2)-1))+sqrt(sqrt(2)/16-
 294                1/16)*log(-sqrt(2*sqrt(2)+2)+sqrt(2)+1))],
 295       [[0,1],%pi/2],
 296       [[%i,%i+1],%pi/4+%i*log(sqrt(2)-1)/2]],
 297 done);
 298 done;
 299
 300 test_table('carlson_rc, mrc, 2.5b-14);
 301 [];
 302
 303 (mrj:[[[0,1,2,3],0.77688623778582b0],
 304       [[2,3,4,5],0.14297579667157b0],
 305       [[2,3,4,-1+%i],0.13613945827771b0-%i*(0.38207561624427b0)],
 306       [[-1+%i,-1-%i,1,2],0.9414835884122b0],
 307       [[-1+%i,-1-%i,1,-3+%i],-0.61127970812028b0-%i*(1.0684038390007b0)],
 308       [[-1+%i,-2-%i,-%i,-1+%i],1.8249027393704b0-%i*(1.2218475784827b0)],
 309       [[2,3,4,-0.5],0.24723819703052b0],
 310       [[2,3,4,-5],-0.12711230042964b0]],
 311 done);
 312 done;
 313
 314 (mrj2:[[[%i,-%i,0,2],1.6490011662711b0],
 315       [[%i,-%i,0,1-%i],1.8260115229009b0+%i*(1.2290661908643b0)]],
 316 done);
 317 done;
 318
 319 test_table('carlson_rj, mrj, 1b-13);
 320 [];
 321
 322 test_table('carlson_rj, mrj2, 2.3452b-8);
 323 [];
 324
 325 /*
 326  * rd(0,2,1) = (3*gamma(3/4)^2)/(sqrt(2)*sqrt(%pi))
 327 */
 328 (mrd:[[[0,2,1],(3*gamma(3/4)^2)/(sqrt(2)*sqrt(%pi))],
 329       [[2,3,4],0.16510527294261b0],
 330       [[-2-%i,-%i,-1+%i],1.8249027393704b0-%i*(1.2218475784827b0)]],
 331 done);
 332 done;
 333
 334 (mrd2:[[[%i,-%i,2],0.6593385415422b0],
 335       [[0,%i,-%i],1.270819627191b0+%i*(2.7811120159521b0)],
 336       [[0,%i-1,%i],-1.8577235439239b0-%i*(0.96193450888839b0)]],
 337 done);
 338 done;
 339
 340 test_table('carlson_rd, mrd, 6e-14);
 341 [];
 342
 343 test_table('carlson_rd, mrd2, 4.7818b-8);
 344 [];
 345
 346 /* Some tests of the Jacobian elliptic functions.
 347  * Just some tests at random points, to verify that we are accurate.
 348  * The reference values were obtained from Mathematica, but we could
 349  * also just compute the values using a much larger fpprec.
 350  */
 351
 352 test_table('jacobi_sn,
 353            [[[1b0+%i*1b0, .7b0], 1.134045971912365274954b0 + 0.352252346922494477621b0*%i],
 354             [[1b0+%i*1b0, 2b0], 0.98613318109123804740b0 + 0.09521910819727230780b0*%i],
 355             [[1b0+%i*1b0, 2b0+3b0*%i], 0.94467077879445294981b0 - 0.19586410083100945528b0* %i],
 356             [[1.785063352082689082581887b0 *%i + 8.890858759528509578661528b-1,
 357               9.434463451695984398149033b-1 - 1.476052973708684178844821b-1 * %i],
 358              1.345233534539675700312281b0 - 7.599023926615176284214056b-2 * %i]],
 359            1b-15);
 360 [];
 361
 362 test_table('jacobi_cn,
 363            [[[1b0+%i*1b0, .7b0], 0.571496591371764254029b0 - 0.698989917271916772991b0*%i ],
 364             [[1b0+%i*1b0, 2b0], 0.33759463268642431412b0 - 0.27814044708010806377b0*%i],
 365             [[1b0+%i*1b0, 2b0+3b0*%i], -0.52142079485827170824b0 - 0.35485177134179601850b0*%i],
 366             [[100b0, .7b0], 0.93004753815774770476196b0]],
 367            6b-14);
 368 [];
 369
 370 test_table('jacobi_dn,
 371            [[[1b0+%i*1b0, .7b0], 0.62297154372331777630564880787568b0 - 0.448863598031509643267241389621738b0 *%i ],
 372             [[1b0+%i*1b0, 2b0], 0.1913322443206041462495606602242b0 - 0.9815253294150083432282549919753b0 * %i],
 373             [[1b0+%i*1b0, 2b0+3b0*%i], 0.6147387452173944656984656771134b0 - 1.4819401302071697495918834416787b0 * %i],
 374             [[100b0, .7b0], 0.95157337933724324055428565654872978b0],
 375             [[1.785063352082689082581887b0 *%i + 8.890858759528509578661528b-1,
 376               9.434463451695984398149033b-1 - 1.476052973708684178844821b-1 * %i],
 377              -8.617730683333292717095686b-1 *%i - 2.663978258141280808361839b-1]],
 378            2b-15);
 379 [];
 380
 381 /* These routines for cn and dn work well for small (<= 1?) values of
 382  * u and m.  They have known issues for large real values of u.
 383  */
 384 (ascending_transform(u,m) :=
 385   block([root_m : expand(rectform(sqrt(m))), mu, root_mu1, v],
 386     mu : expand(rectform(4*root_m/(1+root_m)^2)),
 387     root_mu1 : expand(rectform((1-root_m)/(1+root_m))),
 388     v : expand(rectform(u/(1+root_mu1))),
 389     [v, mu, root_mu1]),
 390  elliptic_dn_ascending(u,m) :=
 391   if is(abs(m-1) < 4*10^(-fpprec)) then sech(u)
 392   else
 393     block([v, mu, root_mu1, new],
 394       [v, mu, root_mu1] : ascending_transform(u,m),
 395       new : elliptic_dn_ascending(v, mu),
 396       expand(rectform((1-root_mu1)/mu*(new^2 + root_mu1)/new))),
 397  elliptic_cn_ascending(u,m) :=
 398   if is(abs(m-1) < 4*10^(-fpprec)) then sech(u)
 399   else
 400     block([v, mu, root_mu1, new],
 401       [v, mu, root_mu1] : ascending_transform(u,m),
 402       new : elliptic_dn_ascending(v, mu),
 403       expand(rectform((1+root_mu1)/mu*(new^2-root_mu1)/new))),
 404  done);
 405 done;
 406
 407 /* Test with random values for the argument and parameter. */
 408 (test_random(n, eps, testf, truef) :=
 409   block([badpoints : [], maxerr : -1],
 410     for k : 1 thru n do
 411       block([z : bfloat(1-2*random(1d0) + %i * (1-2*random(1d0))),
 412              m : bfloat(1-2*random(1d0) + %i * (1-2*random(1d0))),
 413              ans, expected, abserr],
 414         ans : testf(z, m),
 415         expected : truef(z, m),
 416         abserr : abs(ans-expected),
 417         maxerr : max(maxerr, abserr),
 418         if abserr > eps then
 419           badpoints : cons([[z,m], ans, expected, abserr], badpoints)),
 420     if badpoints # [] then
 421       cons(maxerr, badpoints)
 422     else
 423       badpoints),
 424  done);
 425 done;
 426
 427 test_random(50, 2b-15, 'jacobi_dn, 'elliptic_dn_ascending);
 428 [];
 429
 430 test_random(50, 2b-15, 'jacobi_cn, 'elliptic_cn_ascending);
 431 [];
 432
 433 /* Test elliptic_f by using the fact that
 434
 435 inverse_jacobi_sn(x,m) = elliptic_f(asin(x), m)
 436
 437 */
 438 (test_ef(x,m) := jacobi_sn(elliptic_f(asin(x),m), m), id(x,m):=x, done);
 439 done;
 440
 441 test_random(100, 6b-13, 'test_ef, 'id);
 442 [];
 443
 444 /* Test elliptic_kc These values are from
 445  *
 446  * http://functions.wolfram.com/EllipticIntegrals/EllipticK/03/01/
 447  */
 448
 449 block([oldfpprec : fpprec, fpprec:100],
 450   test_table('elliptic_kc,
 451              [[[1/2], 8*%pi^(3/2)/gamma(-1/4)^2],
 452               [[17-12*sqrt(2)], 2*(2+sqrt(2))*%pi^(3/2)/gamma(-1/4)^2],
 453               [[-1], gamma(1/4)^2/4/sqrt(2*%pi)]],
 454              2b-100));
 455 [];
 456
 457 /* Some tests for specific values */
 458 inverse_jacobi_sn(1,m);
 459 elliptic_kc(m);
 460
 461 inverse_jacobi_sn(x,0);
 462 asin(x);
 463
 464 inverse_jacobi_sn(x,1);
 465 log(tan(asin(x)/2 + %pi/4));
 466
 467 inverse_jacobi_sd(1/sqrt(1-m),m);
 468 elliptic_kc(m);
 469
 470 inverse_jacobi_ds(sqrt(1-m),m);
 471 elliptic_kc(m);
 472
 473 /* elliptic_kc(1) is undefined */
 474 errcatch(elliptic_kc(1));
 475 [];
 476 errcatch(elliptic_kc(1.0));
 477 [];
 478 errcatch(elliptic_kc(1.0b0));
 479 [];
 480
 481 /* Test noun/verb results from elliptic functions */
 482 diff(inverse_jacobi_sn(x,0),x);
 483 1/sqrt(1-x^2);
 484
 485 diff(elliptic_pi(4/3,phi,0),phi);
 486 sec(phi)^2/(1-tan(phi)^2/3);
 487
 488 diff(elliptic_pi(3/4,phi,0),phi);
 489 sec(phi)^2/(1+tan(phi)^2/4);
 490
 491 diff(elliptic_pi(1,phi,0),phi);
 492 sec(phi)^2;
 493
 494 /* signbfloat:false provokes "Is 0 positive, negative, or zero?" */
 495 (signbfloat:false,
 496  diff(elliptic_pi(1,phi,0),phi));
 497 sec(phi)^2;
 498
 499 (reset (signbfloat), 0);
 500 0;
 501
 502 /* Special values */
 503 elliptic_ec(1);
 504 1;
 505
 506 elliptic_ec(0);
 507 %pi/2;
 508
 509 elliptic_ec(1/2);
 510 gamma(3/4)^2/(2*sqrt(%pi))+%pi^(3/2)/(4*gamma(3/4)^2);
 511
 512 elliptic_ec(-1);
 513 sqrt(2)*elliptic_ec(1/2);
 514
 515 /* Test periodicity of elliptic_e */
 516 elliptic_e(x, 1);
 517 2*round(x/%pi)+sin(x-%pi*round(x/%pi));
 518
 519 elliptic_e(3,1/3);
 520 elliptic_e(3-%pi,1/3)+2*elliptic_ec(1/3);
 521
 522 /* Bug #2629: elliptic_kc(3.0) not accurate */
 523 test_table('elliptic_kc, [[[3.0], elliptic_kc(3b0)]], 1b-15);
 524 [];
 525
 526 /* Bug #2630: inverse_jacobi_cn(-2.0, 3.0) generates an error */
 527 test_table('inverse_jacobi_cn, [[[-2.0, 3.0], 2.002154760912212-3.202503914656527*%i]],
 528   1b-15);
 529 [];
 530
 531 /* Bug #2615: Numeric evaluation of inverse Jacobi elliptic functions is wrong for some inputs
 532 */
 533 is(abs(jacobi_dn(inverse_jacobi_dn(-2.0,3.0), 3.0) + 2) < 1d-14);
 534 true;
 535
 536 /* elliptical functions handling of non-rectangular complex numbers
 537  * mailing list 2016-07-14 "Jacobi elliptic functions, maxima 5.38.1"
 538  * mailing list 2020-05-28 "Unexpected result from ev"
 539  */
 540
 541 /* list of functions here is everything that has a SIMP-[$%](ELLIPTIC|JACOBI)_FOO in sr/ellipt.lisp */
 542 /*
 543 elliptic_fcns :
 544   [ elliptic_e, elliptic_ec, elliptic_eu, elliptic_f, elliptic_kc, elliptic_pi,
 545     inverse_jacobi_cd, inverse_jacobi_cn, inverse_jacobi_cs, inverse_jacobi_dc,
 546     inverse_jacobi_dn, inverse_jacobi_ds, inverse_jacobi_nc, inverse_jacobi_nd,
 547     inverse_jacobi_ns, inverse_jacobi_sc, inverse_jacobi_sd, inverse_jacobi_sn,
 548     jacobi_am,
 549     jacobi_cd, jacobi_cn, jacobi_cs, jacobi_dc, jacobi_dn, jacobi_ds,
 550     jacobi_nc, jacobi_nd, jacobi_ns, jacobi_sc, jacobi_sd, jacobi_sn ];
 551  */
 552
 553 /* code to generate test cases below
 554
 555 elliptic_fcns_1 : sublist (elliptic_fcns, lambda ([f], errcatch(f('a)) # []));
 556 elliptic_fcns_2 : sublist (elliptic_fcns, lambda ([f], errcatch(f('a, 'b)) # []));
 557 elliptic_fcns_3 : sublist (elliptic_fcns, lambda ([f], errcatch(f('a, 'b, 'c)) # []));
 558
 559 kill (u1, k1, n1);
 560
 561 with_stdout ("/tmp/foo.mac",
 562
 563   for f in elliptic_fcns_1
 564     do block ([a : f('rectform(u1))],
 565         print ('complex_floatp(a), ";"),
 566         print (true, ";"),
 567         print ("")),
 568
 569   for f in elliptic_fcns_2
 570     do block ([a : f('rectform(u1), k1)],
 571         print ('complex_floatp(a), ";"),
 572         print (true, ";"),
 573         print ("")),
 574
 575   for f in elliptic_fcns_3
 576     do block ([a : f(n1, 'rectform(u1), k1)],
 577         print ('complex_floatp(a), ";"),
 578         print (true, ";"),
 579         print ("")));
 580  */
 581
 582 (u1 : 0.5*(1.0 - 0.25*%i)^2,
 583  k1 : 0.775,
 584  n1 : -1.375,
 585  complex_floatp(e) := floatnump(realpart(e)) and floatnump(imagpart(e)),
 586  0);
 587 0;
 588
 589 complex_floatp(elliptic_ec(rectform(u1))) ;
 590 true ;
 591
 592 complex_floatp(elliptic_kc(rectform(u1))) ;
 593 true ;
 594
 595 complex_floatp(elliptic_e(rectform(u1), k1)) ;
 596 true ;
 597
 598 complex_floatp(elliptic_eu(rectform(u1), k1)) ;
 599 true ;
 600
 601 complex_floatp(elliptic_f(rectform(u1), k1)) ;
 602 true ;
 603
 604 complex_floatp(inverse_jacobi_cd(rectform(u1), k1)) ;
 605 true ;
 606
 607 complex_floatp(inverse_jacobi_cn(rectform(u1), k1)) ;
 608 true ;
 609
 610 complex_floatp(inverse_jacobi_cs(rectform(u1), k1)) ;
 611 true ;
 612
 613 complex_floatp(inverse_jacobi_dc(rectform(u1), k1)) ;
 614 true ;
 615
 616 complex_floatp(inverse_jacobi_dn(rectform(u1), k1)) ;
 617 true ;
 618
 619 complex_floatp(inverse_jacobi_ds(rectform(u1), k1)) ;
 620 true ;
 621
 622 complex_floatp(inverse_jacobi_nc(rectform(u1), k1)) ;
 623 true ;
 624
 625 complex_floatp(inverse_jacobi_nd(rectform(u1), k1)) ;
 626 true ;
 627
 628 complex_floatp(inverse_jacobi_ns(rectform(u1), k1)) ;
 629 true ;
 630
 631 complex_floatp(inverse_jacobi_sc(rectform(u1), k1)) ;
 632 true ;
 633
 634 complex_floatp(inverse_jacobi_sd(rectform(u1), k1)) ;
 635 true ;
 636
 637 complex_floatp(inverse_jacobi_sn(rectform(u1), k1)) ;
 638 true ;
 639
 640 complex_floatp(jacobi_am(rectform(u1), k1)) ;
 641 true ;
 642
 643 complex_floatp(jacobi_cd(rectform(u1), k1)) ;
 644 true ;
 645
 646 complex_floatp(jacobi_cn(rectform(u1), k1)) ;
 647 true ;
 648
 649 complex_floatp(jacobi_cs(rectform(u1), k1)) ;
 650 true ;
 651
 652 complex_floatp(jacobi_dc(rectform(u1), k1)) ;
 653 true ;
 654
 655 complex_floatp(jacobi_dn(rectform(u1), k1)) ;
 656 true ;
 657
 658 complex_floatp(jacobi_ds(rectform(u1), k1)) ;
 659 true ;
 660
 661 complex_floatp(jacobi_nc(rectform(u1), k1)) ;
 662 true ;
 663
 664 complex_floatp(jacobi_nd(rectform(u1), k1)) ;
 665 true ;
 666
 667 complex_floatp(jacobi_ns(rectform(u1), k1)) ;
 668 true ;
 669
 670 complex_floatp(jacobi_sc(rectform(u1), k1)) ;
 671 true ;
 672
 673 complex_floatp(jacobi_sd(rectform(u1), k1)) ;
 674 true ;
 675
 676 complex_floatp(jacobi_sn(rectform(u1), k1)) ;
 677 true ;
 678
 679 complex_floatp(elliptic_pi(n1, rectform(u1), k1)) ;
 680 true ;
 681
 682 /*
 683 with_stdout ("/tmp/bar.mac",
 684
 685   for f in elliptic_fcns_1
 686     do block ([a : f('rectform(u1))],
 687         print ('complex_bfloatp(a), ";"),
 688         print (true, ";"),
 689         print ("")),
 690
 691   for f in elliptic_fcns_2
 692     do block ([a : f('rectform(u1), k1)],
 693         print ('complex_bfloatp(a), ";"),
 694         print (true, ";"),
 695         print ("")),
 696
 697   for f in elliptic_fcns_3
 698     do block ([a : f(n1, 'rectform(u1), k1)],
 699         print ('complex_bfloatp(a), ";"),
 700         print (true, ";"),
 701         print ("")));
 702  */
 703
 704 (u1 : 0.5b0*(1.0b0 - 0.25b0*%i)^2,
 705  k1 : 0.775b0,
 706  n1 : -1.375b0,
 707  complex_bfloatp(e) := bfloatp(realpart(e)) and bfloatp(imagpart(e)),
 708  0);
 709 0;
 710
 711 complex_bfloatp(elliptic_ec(rectform(u1))) ;
 712 true ;
 713
 714 complex_bfloatp(elliptic_kc(rectform(u1))) ;
 715 true ;
 716
 717 complex_bfloatp(elliptic_e(rectform(u1), k1)) ;
 718 true ;
 719
 720 complex_bfloatp(elliptic_eu(rectform(u1), k1)) ;
 721 true ;
 722
 723 complex_bfloatp(elliptic_f(rectform(u1), k1)) ;
 724 true ;
 725
 726 complex_bfloatp(inverse_jacobi_cd(rectform(u1), k1)) ;
 727 true ;
 728
 729 complex_bfloatp(inverse_jacobi_cn(rectform(u1), k1)) ;
 730 true ;
 731
 732 complex_bfloatp(inverse_jacobi_cs(rectform(u1), k1)) ;
 733 true ;
 734
 735 complex_bfloatp(inverse_jacobi_dc(rectform(u1), k1)) ;
 736 true ;
 737
 738 complex_bfloatp(inverse_jacobi_dn(rectform(u1), k1)) ;
 739 true ;
 740
 741 complex_bfloatp(inverse_jacobi_ds(rectform(u1), k1)) ;
 742 true ;
 743
 744 complex_bfloatp(inverse_jacobi_nc(rectform(u1), k1)) ;
 745 true ;
 746
 747 complex_bfloatp(inverse_jacobi_nd(rectform(u1), k1)) ;
 748 true ;
 749
 750 complex_bfloatp(inverse_jacobi_ns(rectform(u1), k1)) ;
 751 true ;
 752
 753 complex_bfloatp(inverse_jacobi_sc(rectform(u1), k1)) ;
 754 true ;
 755
 756 complex_bfloatp(inverse_jacobi_sd(rectform(u1), k1)) ;
 757 true ;
 758
 759 complex_bfloatp(inverse_jacobi_sn(rectform(u1), k1)) ;
 760 true ;
 761
 762 complex_bfloatp(jacobi_am(rectform(u1), k1)) ;
 763 true ;
 764
 765 complex_bfloatp(jacobi_cd(rectform(u1), k1)) ;
 766 true ;
 767
 768 complex_bfloatp(jacobi_cn(rectform(u1), k1)) ;
 769 true ;
 770
 771 complex_bfloatp(jacobi_cs(rectform(u1), k1)) ;
 772 true ;
 773
 774 complex_bfloatp(jacobi_dc(rectform(u1), k1)) ;
 775 true ;
 776
 777 complex_bfloatp(jacobi_dn(rectform(u1), k1)) ;
 778 true ;
 779
 780 complex_bfloatp(jacobi_ds(rectform(u1), k1)) ;
 781 true ;
 782
 783 complex_bfloatp(jacobi_nc(rectform(u1), k1)) ;
 784 true ;
 785
 786 complex_bfloatp(jacobi_nd(rectform(u1), k1)) ;
 787 true ;
 788
 789 complex_bfloatp(jacobi_ns(rectform(u1), k1)) ;
 790 true ;
 791
 792 complex_bfloatp(jacobi_sc(rectform(u1), k1)) ;
 793 true ;
 794
 795 complex_bfloatp(jacobi_sd(rectform(u1), k1)) ;
 796 true ;
 797
 798 complex_bfloatp(jacobi_sn(rectform(u1), k1)) ;
 799 true ;
 800
 801 complex_bfloatp(elliptic_pi(n1, rectform(u1), k1)) ;
 802 true ;
 803
 804 /* Test derivatives at a few random points using numerical differentiation
 805
 806    func  - a function
 807    df - deriviative of func wrt p-th argument
 808    p - derivative wrt p-th argument
 809    vars - variables in expression deriv
 810    table - table of points to evaluate
 811    eps - report errors > eps
 812    delta - offset for numerical derivative
 813 */
 814 (test_deriv(func,df,vars,p,table,eps,delta) :=
 815 block([badpoints : [],
 816        abserr    : 0,
 817        maxerr    : -1],
 818   for %z in table do
 819     block([%z0, %z1, deriv, nderiv],
 820       %z1:makelist(if i=p then %z[i]+delta else %z[i],i,1,length(%z)),
 821       %z0:makelist(if i=p then %z[i]-delta else %z[i],i,1,length(%z)),
 822       nderiv: float((apply(func,%z1)-apply(func,%z0))/(2*delta)),
 823       deriv : float(subst(maplist("=",vars,%z),df)),
 824       abserr : abs(nderiv-deriv),
 825       maxerr : max(maxerr,abserr),
 826       if abserr > eps then
 827         badpoints : cons ([%z,deriv,nderiv,abserr],badpoints)
 828     ),
 829   if badpoints # [] then
 830     cons(maxerr,badpoints)
 831   else
 832     badpoints
 833     ),done);
 834 done;
 835
 836 /* test points for two-arg functions */
 837 (l2: [[0.2, 0.3], [0.2, 0.5], [1.5, 0.7], [0.5, 0.99], [1.5, 0.8],
 838     [1.57, 0.8], [1.5707, 0.8], [-1.0, 0.8], [2.0, 0.8]],done);
 839 done;
 840
 841 /* test points for three-arg functions */
 842 /* bug #3221 elliptic_pi() wrong for 2nd arg > float(%pi/2) */
 843 (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],
 844     [0.2,1.5,0.8],[0.2,1.57,0.8],[0.2,1.5707,0.8]],done);
 845 done;
 846
 847 test_deriv('elliptic_f,diff(elliptic_f(z,m),z),[z,m],1,l2,2.0e-7,1.0e-8);
 848 [];
 849
 850 test_deriv('elliptic_f,diff(elliptic_f(z,m),m),[z,m],2,l2,2.0e-7,1.0e-8);
 851 [];
 852
 853 test_deriv('elliptic_e,diff(elliptic_e(z,m),z),[z,m],1,l2,2.0e-7,1.0e-8);
 854 [];
 855
 856 test_deriv('elliptic_e,diff(elliptic_e(z,m),m),[z,m],2,l2,2.0e-7,1.0e-8);
 857 [];
 858
 859 test_deriv('elliptic_pi,diff(elliptic_pi(n,z,m),n),[n,z,m],1,l3,2.0e-7,1.0e-8);
 860 [];
 861
 862 test_deriv('elliptic_pi,diff(elliptic_pi(n,z,m),z),[n,z,m],2,l3,2.0e-7,1.0e-8);
 863 [];
 864
 865 test_deriv('elliptic_pi,diff(elliptic_pi(n,z,m),m),[n,z,m],3,l3,2.0e-7,1.0e-8);
 866 [];
 867
 868 closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
 869 closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
 870
 871 /*
 872  * Some additional tests for complex values. Use the fact that
 873  * elliptic_pi(0, phi, m) = elliptic_f(phi,m);
 874  */
 875 closeto(elliptic_pi(0, 0.5, 0.25)- elliptic_f(0.5, 0.25), 1e-15);
 876 true;
 877 closeto(elliptic_pi(0, 0.5*%i, 0.25)- elliptic_f(0.5*%i, 0.25), 1e-15);
 878 true;
 879
 880 closeto(elliptic_pi(0, 0.5b0, 0.25b0)- elliptic_f(0.5b0, 0.25b0), 1e-15);
 881 true;
 882 closeto(elliptic_pi(0, 0.5b0*%i, 0.25b0)- elliptic_f(0.5b0*%i, 0.25b0), 1e-15);
 883 true;
 884
 885 /* Test for Bug 3221 */
 886
 887 closeto(elliptic_pi(0.2,1.57,0.8) - 2.5770799919605668, 1e-14);
 888 true;
 889 closeto(elliptic_pi(0.2,1.58,0.8) - 2.60502920656026151, 1e-14);
 890 true;
 891 closeto(elliptic_pi(0.20,2.00,0.80) - 3.6543835090829025, 7.47e-11);
 892 true;
 893
 894 (oldfpprec:fpprec,fpprec:32);
 895 32;
 896 closeto(elliptic_pi(0.2b0,1.57b0,0.8b0) - 2.5770799919605668058849721013196b0, 3.09b-32);
 897 true;
 898 closeto(elliptic_pi(0.2b0,1.58b0,0.8b0) - 2.6050292065602615145481924917132b0, 1.24b-32);
 899 true;
 900 closeto(elliptic_pi(0.2b0,2.0b0,0.8b0) - 3.6543835090082902501670624830938b0, 1b-32);
 901 true;
 902
 903 /* Test for #3733 $gamma vs %gamma confusion */
 904 elliptic_ec(1/2) - gamma(3/4)^2/(2*sqrt(%pi))- %pi^(3/2)/(4*gamma(3/4)^2);
 905 0$
 906
 907 /* #3745 Quoting either elliptic_f [or elliptic_e] inhibits simplification */
 908 'elliptic_f(5,0);
 909 5$
 910
 911 'elliptic_e(5,0);
 912 5$
 913
 914 /* #3967 elliptic_e(5*%pi/4,1) inconsistent with numerical evaluation */
 915 elliptic_e(5*%pi/4, 1);
 916 1/sqrt(2) + 2;
 917
 918 closeto(elliptic_e(5*%pi/4,1) - elliptic_e(float(5*%pi/4),1), 1e-16);
 919 true;
 920
 921 /* #3746 derivative of inverse_jacobi_sn is noun/verb confused.
 922    The same is true for inverse_jacobi_cn and inverse_jacobi_dn, but the
 923    remaining nine inverse Jacobi functions ns, nc, nd, sc, cs, sd, ds, cd,
 924    and dc don't have defined m derivatives.*/
 925 subst(u=0, diff(inverse_jacobi_sn(u,m),m));
 926 0$
 927
 928 subst(u=1, diff(inverse_jacobi_cn(u,m),m));
 929 0$
 930
 931 subst(u=1, diff(inverse_jacobi_dn(u,m),m));
 932 0$
 933
 934 /* Carlson's integrals */
 935
 936 (assume(xp > 0), done);
 937 done$
 938
 939 carlson_rc(xp,xp);
 940 1/sqrt(xp)$
 941
 942 carlson_rc(0,1);
 943 %pi/2$
 944
 945 carlson_rc(0, 1/4);
 946 %pi$
 947
 948 carlson_rc(2, 1);
 949 log(sqrt(2)+1)$
 950
 951 carlson_rc(%i, %i+1);
 952 %pi/4 + %i*log(sqrt(2)-1)/2$
 953
 954 carlson_rc(0, %i);
 955 ((1-%i)*%pi)/2^(3/2)$
 956
 957 carlson_rc(((1+xp)/2)^2, xp);
 958 log(xp)/(xp-1)$
 959
 960 /* log(x) = (x-1)*carlson_rc(((1+x)/2)^2,x) */
 961 closeto((2-1)*carlson_rc(((1+2)/2)^2,2.0) - log(2.0), 1e-16);
 962 true$
 963
 964 /* asin(x) = x*carlson_rc(1-x^2,1) */
 965 closeto(1/2*carlson_rc(1-1/4,1.0) - asin(0.5), 2.2205e-16);
 966 true$
 967
 968 /* acos(x) = sqrt(1-x^2)*carlson_rc(x^2,1) */
 969 closeto(sqrt(1.0-1/4)*carlson_rc(1/4,1.0) - acos(0.5), 2.2205e-16);
 970 true$
 971
 972 /* atan(x) = x*carlson_rc(1,1+x^2) */
 973 closeto(1*carlson_rc(1.0,1+1^2) - atan(1.0), 1.1103e-16);
 974 true$
 975
 976 /* asinh(x) = x*carlson_rc(1+x^2,1) */
 977 closeto(1*carlson_rc(1+1^2,1.0) - asinh(1.0), 2.2205e-16);
 978 true$
 979
 980 /* acosh(x) = sqrt(x^2-1)*carlson_rc(x^2,1) */
 981 closeto(sqrt(2^2-1.0)*carlson_rc(2^2,1.0) - acosh(2.0), 2.2205e-16);
 982 true$
 983
 984 /* atanh(x) = x*carlson_rc(1,1-x^2) */
 985 closeto(1/2*carlson_rc(1.0,1-(1/2)^2) - atanh(0.5), 1.1103e-16);
 986 true$
 987
 988 carlson_rd(x, x, x);
 989 1/x^(3/2)$
 990
 991 carlson_rd(0,2,1);
 992 3*sqrt(%pi)*gamma(3/4)/gamma(1/4)$
 993
 994 carlson_rd(2,0,1);
 995 3*sqrt(%pi)*gamma(3/4)/gamma(1/4)$
 996
 997 closeto(carlson_rd(0, 2, 1.0) - float(3*sqrt(%pi)*gamma(3/4)/gamma(1/4)), 6.66134e-16);
 998 true$
 999
1000 carlson_rf(x,x,x);
1001 1/sqrt(x)$
1002
1003 /*
1004  * Rf(1,2,0) has an analytic solution.  Check some of the permutations too.
1005  */
1006 carlson_rf(1,2,0);
1007 gamma(1/4)^2/(4*sqrt(2*%pi))$
1008
1009 carlson_rf(2,1,0);
1010 gamma(1/4)^2/(4*sqrt(2*%pi))$
1011
1012 carlson_rf(0, 1,2);
1013 gamma(1/4)^2/(4*sqrt(2*%pi))$
1014
1015 /* Check permutations of Rf(%i, -%i, 0) too */
1016 carlson_rf(%i,-%i,0);
1017 gamma(1/4)^2/(4*sqrt(%pi))$
1018
1019 carlson_rf(-%i,%i,0);
1020 gamma(1/4)^2/(4*sqrt(%pi))$
1021
1022 carlson_rf(0, %i,-%i);
1023 gamma(1/4)^2/(4*sqrt(%pi))$
1024
1025
1026 carlson_rj(x,x,x,x);
1027 1/x^(3/2)$
1028
1029 carlson_rj(x,y,z,z);
1030 carlson_rd(x,y,z)$
1031
1032 (assume(pp>0,yp>0), done);
1033 done$
1034
1035 carlson_rj(xp,xp,xp,pp);
1036 3*(1/sqrt(xp)-carlson_rc(xp,pp))/(pp-xp)$
1037
1038 carlson_rj(0,yp,yp,pp);
1039 3*%pi/(2*(yp*sqrt(pp)+pp*sqrt(yp)))$
1040
1041 carlson_rj(x,yp,yp,pp);
1042 3/(pp-yp)*(carlson_rc(x,yp) - carlson_rc(x,pp))$
1043
1044 carlson_rj(x,yp,yp,yp);
1045 (3*(carlson_rc(x,yp)-sqrt(x)/yp))/(2*(yp-x))$
1046
1047 (fpprec:oldfpprec,kill(l2,l3,test_deriv,oldfpprec),done);
1048 done;