/*
Implementation of methods given in
B. C. Carlson
Computation of elliptic integrals
Numerical Algorithms 10, 13-26 (1995)

Copyright (c) Charles Karney (2009-2013) <charles@karney.com> and
licensed under the MIT/X11 License.  For more information, see
https://geographiclib.sourceforge.io/
*/

/* fpprec:120$  Should be set outside */
etol:0.1b0^fpprec$ /* For Carlson */
tolRF : (3 * etol)^(1/8b0)$
tolRD : (etol/5)^(1/8b0)$
tolRG0 : 2.7b0 * sqrt(etol)$
tolJAC:sqrt(etol)$  /* For Bulirsch */
pi:bfloat(%pi)$
atan2x(y,x) := -atan2(-y, x)$  /* fix atan2(0b0,-1b0) = -pi */
rf(x,y,z) := block([a0:(x+y+z)/3, q,x0:x,y0:y,z0:z,
  an,ln,xx,yy,zz,e2,e3,mul:1b0],
  q:max(abs(a0-x),abs(a0-y),abs(a0-z))/tolRF,
  an:a0,
  while q >= mul * abs(an) do (
    ln:sqrt(x0)*sqrt(y0)+sqrt(y0)*sqrt(z0)+sqrt(z0)*sqrt(x0),
    an:(an+ln)/4,
    x0:(x0+ln)/4,
    y0:(y0+ln)/4,
    z0:(z0+ln)/4,
    mul:mul*4),
  xx:(a0-x)/(mul*an),
  yy:(a0-y)/(mul*an),
  zz:-xx-yy,
  e2:xx*yy-zz^2,
  e3:xx*yy*zz,
  (e3 * (6930 * e3 + e2 * (15015 * e2 - 16380) + 17160) +
    e2 * ((10010 - 5775 * e2) * e2 - 24024) + 240240) /
  (240240 * sqrt(an)))$
rd(x,y,z) := block([a0:(x+y+3*z)/5, q,x0:x,y0:y,z0:z,
  an,ln,xx,yy,zz,e2,e3,e4,e5,s:0b0,mul:1b0],
  q:max(abs(a0-x),abs(a0-y),abs(a0-z))/tolRD,
  an:a0,
  while q >= mul*abs(an) do (
    ln:sqrt(x0)*sqrt(y0)+sqrt(y0)*sqrt(z0)+sqrt(z0)*sqrt(x0),
    s:s+1/(mul*sqrt(z0)*(z0+ln)),
    an:(an+ln)/4,
    x0:(x0+ln)/4,
    y0:(y0+ln)/4,
    z0:(z0+ln)/4,
    mul:4*mul),
  xx:(a0-x)/(mul*an),
  yy:(a0-y)/(mul*an),
  zz:-(xx+yy)/3,
  e2:xx*yy-6*zz^2,
  e3:(3*xx*yy-8*zz^2)*zz,
  e4:3*(xx*yy-zz^2)*zz^2,
  e5:xx*yy*zz^3,
  ((471240 - 540540 * e2) * e5 +
    (612612 * e2 - 540540 * e3 - 556920) * e4 +
    e3 * (306306 * e3 + e2 * (675675 * e2 - 706860) + 680680) +
    e2 * ((417690 - 255255 * e2) * e2 - 875160) + 4084080) /
  (4084080 * mul * an * sqrt(an)) + 3 * s)$

/* R_G(x,y,0) */
rg0(x,y) := block([x0:sqrt(max(x,y)),y0:sqrt(min(x,y)),xn,yn,
  t,s:0b0,mul:0.25b0],
  xn:x0,
  yn:y0,
  while abs(xn-yn) > tolRG0 * xn do (
    t:(xn+yn)/2,
    yn:sqrt(xn*yn),
    xn:t,
    mul : 2*mul,
    s:s+mul*(xn-yn)^2),
  ((x0+y0)^2/4 - s)*pi/(2*(xn+yn)) )$

rg(x, y, z) := (
  if z = 0b0 then (z:y, y:0b0),
  (z * rf(x, y, z) - (x-z) * (y-z) * rd(x, y, z) / 3
    + sqrt(x * y / z)) / 2)$

rj(x, y, z, p) := block([A0:(x + y + z + 2*p)/5,An,delta:(p-x) * (p-y) * (p-z),
  Q,x0:x,y0:y,z0:z,p0:p,mul:1b0,mul3:1b0,s:0b0,X,Y,Z,P,E2,E3,E4,E5],
  An : A0,
  Q : max(abs(A0-x), abs(A0-y), abs(A0-z), abs(A0-p)) / tolRD,
  while Q >= mul * abs(An) do block([lam,d0,e0],
    lam : sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0),
    d0 : (sqrt(p0)+sqrt(x0)) * (sqrt(p0)+sqrt(y0)) * (sqrt(p0)+sqrt(z0)),
    e0 : delta/(mul3 * d0^2),
    s : s + rc(1b0, 1b0 + e0)/(mul * d0),
    An : (An + lam)/4,
    x0 : (x0 + lam)/4,
    y0 : (y0 + lam)/4,
    z0 : (z0 + lam)/4,
    p0 : (p0 + lam)/4,
    mul : 4*mul,
    mul3 : 64*mul3),
  X : (A0 - x) / (mul * An),
  Y : (A0 - y) / (mul * An),
  Z : (A0 - z) / (mul * An),
  P : -(X + Y + Z) / 2,
  E2 : X*Y + X*Z + Y*Z - 3*P*P,
  E3 : X*Y*Z + 2*P * (E2 + 2*P*P),
  E4 : (2*X*Y*Z + P * (E2 + 3*P*P)) * P,
  E5 : X*Y*Z*P*P,
  ((471240 - 540540 * E2) * E5 +
    (612612 * E2 - 540540 * E3 - 556920) * E4 +
    E3 * (306306 * E3 + E2 * (675675 * E2 - 706860) + 680680) +
    E2 * ((417690 - 255255 * E2) * E2 - 875160) + 4084080) /
  (4084080 * mul * An * sqrt(An)) + 6 * s)$

rd(x, y, z) := block([A0:(x + y + 3*z)/5,An,Q,x0:x,y0:y,z0:z,
  mul:1b0,s:0b0,X,Y,Z,E2,E3,E4,E5],
  An : A0,
  Q : max(abs(A0-x), abs(A0-y), abs(A0-z)) / tolRD,
  while Q >= mul * abs(An) do block([lam],
    lam : sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0),
    s : s + 1/(mul * sqrt(z0) * (z0 + lam)),
    An : (An + lam)/4,
    x0 : (x0 + lam)/4,
    y0 : (y0 + lam)/4,
    z0 : (z0 + lam)/4,
    mul : 4*mul),
  X : (A0 - x) / (mul * An),
  Y : (A0 - y) / (mul * An),
  Z : -(X + Y) / 3,
  E2 : X*Y - 6*Z*Z,
  E3 : (3*X*Y - 8*Z*Z)*Z,
  E4 : 3 * (X*Y - Z*Z) * Z*Z,
  E5 : X*Y*Z*Z*Z,
  ((471240 - 540540 * E2) * E5 +
    (612612 * E2 - 540540 * E3 - 556920) * E4 +
    E3 * (306306 * E3 + E2 * (675675 * E2 - 706860) + 680680) +
    E2 * ((417690 - 255255 * E2) * E2 - 875160) + 4084080) /
  (4084080 * mul * An * sqrt(An)) + 3 * s)$

rc(x, y):=if x < y then atan(sqrt((y - x) / x)) / sqrt(y - x)
else ( if x = y and y > 0b0 then 1 / sqrt(y) else
  atanh( if y > 0b0 then sqrt((x - y) / x)
    else sqrt(x / (x - y)) ) / sqrt(x - y) )$

/* k^2 = m */
ec(k2):=2*rg0(1b0-k2,1b0)$
kc(k2):=rf(0b0,1b0-k2,1b0)$
dc(k2):=rd(0b0,1b0-k2,1b0)/3$
pic(k2,alpha2):=if alpha2 # 0b0 then
kc(k2)+alpha2*rj(0b0,1b0-k2,1b0,1b0-alpha2)/3
else kc(k2)$
gc(k2,alpha2):=if alpha2 # 0b0 then (if k2 # 1b0 then
kc(k2) + (alpha2-k2)*rj(0b0,1b0-k2,1b0,1b0-alpha2)/3
else rc(1b0,1b0-alpha2))
else ec(k2)$
hc(k2,alpha2):=if alpha2 # 0b0 then (if k2 # 1b0 then
kc(k2) - (1b0-alpha2)*rj(0b0,1b0-k2,1b0,1b0-alpha2)/3
else rc(1b0,1b0-alpha2))
else kc(k2) - dc(k2)$

/* Implementation of methods given in

Roland Bulirsch
Numerical Calculation of Elliptic Integrals and Elliptic Functions
Numericshe Mathematik 7, 78-90 (1965)

*/

sncndn(x,mc):=block([bo, a, b, c, d, l, sn, cn, dn, m, n],
  local(m, n),
  if mc # 0 then (
    bo:is(mc < 0b0),
    if bo then (
      d:1-mc,
      mc:-mc/d,
      d:sqrt(d),
      x:d*x),
    dn:a:1,
    for i:0 thru 12 do (
      l:i,
      m[i]:a,
      n[i]:mc:sqrt(mc),
      c:(a+mc)/2,
      if abs(a-mc)<=tolJAC*a then return(false),
      mc:a*mc,
      a:c
      ),
    x:c*x,
    sn:sin(x),
    cn:cos(x),
    if sn#0b0 then (
      a:cn/sn,
      c:a*c,
      for i:l step -1 thru 0 do (
        b:m[i],
        a:c*a,
        c:dn*c,
        dn:(n[i]+a)/(b+a),
        a:c/b
        ),
      a:1/sqrt(c*c+1b0),
      sn:if sn<0b0 then -a else a,
      cn:c*sn
      ),
    if bo then (
      a:dn,
      dn:cn,
      cn:a,
      sn:sn/d
      )
    ) else /* mc = 0 */ (
    sn:tanh(x),
    dn:cn:sech(x)
/*    d:exp(x), a:1/d, b:a+d, cn:dn:2/b,
    if x < 0.3b0 then (
      d:x*x*x*x,
      d:(d*(d*(d*(d+93024b0)+3047466240b0)+24135932620800b0)+
        20274183401472000b0)/60822550204416000b0,
      sn:cn*(x*x*x*d+sin(x))
      ) else
    sn:(d-a)/b */
    ),
  [sn,cn,dn]
  )$

Delta(sn, k2):=sqrt(1 - k2 * sn*sn)$

/* Versions of incomplete functions in terms of Jacobi elliptic function
with u = am(phi) real and in [0,K(k2)] */

eirx(sn,cn,dn,k2,ec):=block([t,cn2:cn*cn,dn2:dn*dn,sn2:sn*sn,ei,kp2:1b0-k2],
  ei : (if  k2 <= 0b0 then
    rf(cn2, dn2, 1b0) - k2 * sn2 * rd(cn2, dn2, 1b0) / 3 else
    (if kp2 >= 0b0 then
      kp2 * rf(cn2, dn2, 1b0) +
      k2 * kp2 * sn2 * rd(cn2, 1b0, dn2) / 3 +
      k2 * abs(cn) / dn else
      - kp2 * sn2 * rd(dn2, 1b0, cn2) / 3 + dn / abs(cn) ) ),
  ei : ei * abs(sn),
  if cn < 0b0 then ei : 2 * ec - ei,
  if sn < 0 then ei : -ei,
  ei)$

dirx(sn, cn, dn, k2, dc):=block(
  [di:abs(sn) * sn*sn * rd(cn*cn, dn*dn, 1b0) / 3],
  if cn < 0b0 then di : 2 * dc - di,
  if sn < 0b0 then di : -di,
  di)$

hirx(sn, cn, dn, k2, alpha2, hc):=block(
  [cn2 : cn*cn, dn2 : dn*dn, sn2 : sn*sn, hi, alphap2:1b0-alpha2],
  hi : abs(sn) * (rf(cn2, dn2, 1b0) -
    alphap2 * sn2 *
    rj(cn2, dn2, 1b0, cn2 + alphap2 * sn2) / 3),
  if cn < 0 then hi : 2 * hc - hi,
  if sn < 0 then hi : -hi,
  hi)$

deltae(sn, cn, dn, k2, ec):=(
  if cn < 0 then (cn : -cn, sn : -sn),
  eirx(sn, cn, dn, k2, ec) * (pi/2) / ec - atan2x(sn, cn))$
deltad(sn, cn, dn, k2, dc):=(
  if cn < 0 then (cn : -cn, sn : -sn),
  dirx(sn, cn, dn, k2, ec) * (pi/2) / dc - atan2x(sn, cn))$
deltah(sn, cn, dn, k2, alpha2, hc):=(
  if cn < 0 then (cn : -cn, sn : -sn),
  hirx(sn, cn, dn, k2, alpha2, hc) * (pi/2) / hc - atan2x(sn, cn))$

einv(x, k2, ec):=block(
  [n : floor(x / (2 * ec) + 0.5b0), phi, eps : k2/(sqrt(1b0-k2) + 1b0)^2,
  err:1b0],
  x : x - 2 * ec * n,
  phi : pi * x / (2 * ec),
  phi : phi - eps * sin(2 * phi) / 2,
  while abs(err) > tolJAC do block([sn:sin(phi), cn:cos(phi), dn],
    dn : Delta(sn, k2),
    err : (eirx(sn, cn, dn, k2, ec) - x)/dn,
    phi : phi - err),
  n * pi + phi)$

deltaeinv(stau, ctau, k2, ec) := block([tau],
  if ctau < 0 then (ctau : -ctau, stau : -stau),
  tau : atan2x(stau, ctau),
  einv(tau * ec / (pi/2), k2, ec) - tau)$