\ elip     Complete Elliptic Integral         ACM Algorithm #149

\ Forth Scientific Library Algorithm #2

\ Evaluates the Complete Elliptic Integral,
\     Elip[a, b] = int_0^{\pi/2} 1/Sqrt{a^2 cos^2(t) + b^2 sin^2(t)} dt

\ This function can be used to evaluate the complete elliptic integral
\ of the first kind, by using the relation K[m] = a Elip[a,b],  m = 1 - b^2/a^2
 
\ This code conforms with ANS requiring:
\      1. The Floating-Point word set
\      2. The FCONSTANT PI (3.1415926536...)
 
\ Both a recursive form and an iterative form are given, but because of the
\ large stack consumption the recursive form is probably not of much
\ practical use.

\ Caution: this code can potentially go into an endless loop
\          for certain values of the parameters.

\ Collected Algorithms from ACM, Volume 1 Algorithms 1-220,
\ 1980; Association for Computing Machinery Inc., New York,
\ ISBN 0-89791-017-6

\ (c) Copyright 1994 Everett F. Carter.  Permission is granted by the
\ author to use this software for any application provided the
\ copyright notice is preserved.

\ ( ELIP     V1.2                  21 September 1994   EFC )
\ include lib/ansfloat.4th

[UNDEFINED] felip [IF]
[UNDEFINED] float [IF] [ABORT] [THEN]
[UNDEFINED] pi    [IF] include lib/fpconst.4th [THEN]
: felip ( --, f: a b -- elip[a,b] )     \ nonrecursive version

    BEGIN
      FOVER FOVER F+ 2 S>F F/
      FROT  FROT  F* FSQRT

      FOVER FOVER FOVER F- FABS
      FSWAP S" 1.0e-8" S>FLOAT F*
   F< UNTIL

      FDROP

      pi 2 S>F F/ FSWAP F/
;

\ test driver,  calculates the complete elliptic integral of the first
\ kind (K(m)) using the relation: K[m] = a Elip[a,b], m = 1 - b^2/a^2
\ compare with Abramowitz & Stegun, Handbook of Mathematical Functions,
\ Table 17.1
false [IF]
: felipr ( --, f: a b -- elip[a,b] )     \ recursive form

     FOVER FOVER FOVER F- FABS
     FSWAP S" 1.0e-8" S>FLOAT F*
     F< IF
          FDROP
          pi 2 S>F F/ FSWAP F/
        ELSE
          FOVER FOVER F+ 2 S>F F/
          FROT  FROT  F* FSQRT
          RECURSE
        THEN
;

: elip_test ( -- )
        fclear
        10 set-precision

        CR
        ."  m     K(m) exact    a Elip1[a,b]    a Elip2[a,b] " CR

       ." 0.0    1.57079633   "
       S" 1000.0" S>FLOAT S" 1000.0" S>FLOAT felipr S" 1000.0" S>FLOAT  F* F. ."     "
       S" 1000.0" S>FLOAT S" 1000.0" S>FLOAT felip S" 1000.0" S>FLOAT  F* F. CR

      ." 0.44   1.80632756   "
      S" 400.0" S>FLOAT S" 299.33259" S>FLOAT felipr  S" 400.0" S>FLOAT  F* F. ."     "
      S" 400.0" S>FLOAT S" 299.33259" S>FLOAT felip  S" 400.0" S>FLOAT  F* F. CR

      ." 0.75   2.15651565   "
      S" 1000.0" S>FLOAT S" 500.0" S>FLOAT felipr S" 1000.0" S>FLOAT  F* F. ."     "
      S" 1000.0" S>FLOAT S" 500.0" S>FLOAT felip S" 1000.0" S>FLOAT  F* F. CR

      ." 0.96   3.01611249   "
      S" 500.0" S>FLOAT S" 100.0" S>FLOAT felipr S" 500.0" S>FLOAT  F* F.  ."     "
      S" 500.0" S>FLOAT S" 100.0" S>FLOAT felip S" 500.0" S>FLOAT  F* F. CR
;
elip_test
[THEN]
[THEN]