*///////////////////////////////////////////////////////////////////////
* thrift.f - cheap imitations of intrinsic functions
*   by Andy Allinger, 2013, released to the public domain
*   Be sure to read the comments regarding the range and precision.
*   Also be warned that these functions may be SLOWER than the intrinsics
*    on some compilers.
*///////////////////////////////////////////////////////////////////////
*  Quartic rational algebraic functions are used to replace function calls.
*  Given:  y = c0 + (c1 @ x) + (c2 @ x^2) + (c3 @ x^3) + (c4 @ x^4)
*  Let:  P = c4^(1/4)
*        Q = (c3 - P^3) / (4 @ P^3)
*        S = 3 @ Q^2 + 8 @ Q^3 + (c1 @ P - 2 @ c2 @ Q) / P^2
*        R = c2 / P^2 - (2 @ Q) - (6 @ Q^2) - S
*        T = c0 - Q^4 - (Q^2) @ (R + S) - (R @ S)
*  Thus:  y = ((P@x + Q)^2 + P@x + R)((P@x + Q)^2 + S) + T
*
*    attributed to Donald Knuth, published in:
*///////////////////////////////////////////////////////////////////////
* William H. Press, Brian P. Flannery, Saul A. Teukolsky & Willam T. Vetterling
* Numerical Recipes in C:  the Art of Scientific Computing
* Cambridge University Press, 1988
*///////////////////////////////////////////////////////////////////////

*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
*     quartic rational algebraic function approximation of e^x
*        relative error < .000014 for x < 0
*        not for positive numbers
*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      FUNCTION SNAP (t)      ! swift Napier      
       IMPLICIT NONE
       REAL t, x, SNAP
       REAL C(10), ax, axb, axb2, numer, denom, f
       SAVE C
       DATA C / 0.1037642555, 0.6703349337, 1.075752351, 0.2436912224, 
     &  -0.5696177849E-01, 0.2477126235, -0.7225424107, 0.9719112436, 
     &  0.2440869124, -0.1446184861 /

       x = t
       f = 1.0
  50   IF (x < -7.) THEN           ! reduce argument to preferred domain
         x = x + 7.
         f = f * .00091188196555
         GO TO 50
       END IF
       
*          fast quartic
       ax = C(1) * x
       axb = ax + C(2)
       axb2 = axb * axb
       numer = (axb2 + ax + C(3)) * (axb2 + C(4)) + C(5)
 
       ax = C(6) * x
       axb = ax + C(7)
       axb2 = axb * axb
       denom = (axb2 + ax + C(8)) * (axb2 + C(9)) + C(10) 

       SNAP = f * numer / denom
       RETURN
      END  ! of SNAP


*#######################################################################
*     quartic rational algebraic function approximation of ln(x)
*        absolute error < .00025 for x < 1
*        not for x > 1
*#######################################################################
      FUNCTION SLOG (t)         ! swift logarithm
       IMPLICIT NONE
       REAL C(10), t, f, x, SLOG
       SAVE C

       DATA C / 30064.81789, -9296.414198, -19466.37761, -1294.772254,
     &  -7.253826679, 11954.83362, 31605.20508, 8339.845481, 
     &  294.5602191,  1.000000000 /

*           validate
       IF (t .LE. 0.) THEN
         SLOG = -999.
         RETURN
       ELSE IF (t .GE. 1.) THEN
         SLOG = 0.
         RETURN
       END IF

*         reduce argument to preferred domain
       x = t
       f = 0.
  50   IF (x < .00673794699908) THEN
         x = x * 148.41315910257
         f = f - 5.
         GO TO 50
       END IF

*                synthetic division
         SLOG = f + ((((C(1) * x + C(2)) * x + C(3)) 
     &      * x + C(4)) * x + C(5))
     &   / ((((C(6) * x + C(7)) * x + C(8)) 
     &      * x + C(9)) * x + C(10))
       RETURN
      END  ! of SLOG


*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
*  quartic rational algebraic function approximation of SIN(x)
*   valid in the range 0...2 @ Pi, accuracy is a couple parts per million
*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
      FUNCTION SSIN (t)         ! swift sine
       IMPLICIT NONE
       REAL t, x, x2, SSIN
       REAL C(5), hpi, pi, piah, twopi     
       SAVE C, hpi, pi, piah, twopi
       DATA C / -0.1038742270,   0.9999915683, 
     &             0.002196769622, 0.06274478232, 1.0 /       
       DATA hpi, pi, piah, twopi / 1.5707963268, 3.1415926536,
     &  4.7123889804, 6.2831853072 /
               
*          choose what quadrant
       IF (t < hpi) THEN
           x = t
       ELSE IF (t < piah) THEN
           x = pi - t
       ELSE
           x = t - twopi
       END IF
      
*          apply rational algebraic function
       x2 = x **2
       SSIN = ((C(1) * x2 + C(2)) * x) 
     &    / ((C(3) * x2 + C(4)) * x2 + C(5))
       RETURN
      END  ! of SSIN


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC       
C  quartic rational algebraic function approximation of COS(x)
C   valid in the range 0...2 @ Pi
C   accurate better than a part per million !
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC       
      FUNCTION SCOS (t)         ! swift cosine
       IMPLICIT NONE
       REAL t, x, x2, SCOS
       REAL C(6), hpi, pi, piah, twopi
       SAVE C, hpi, pi, piah, twopi
       DATA C / 0.2025147905E-01, -0.4552531975, 0.9999999316, 
     &  0.9623103039E-03, 0.4474536302E-01, 1.000000000 /
       DATA hpi, pi, piah, twopi / 1.5707963268, 3.1415926536,
     &  4.7123889804, 6.2831853072 /
       
*          choose what quadrant
       IF (t < hpi) THEN
         x = t
       ELSE IF (t < piah) THEN
         x = pi - t
       ELSE
         x = t - twopi
       END IF
      
*          apply rational algebraic function
       x2 = x **2
       SCOS = ((C(1) * x2 + C(2)) * x2 + C(3)) 
     &    / ((C(4) * x2 + C(5)) * x2 + C(6))

*              fix sign
       IF (ABS(t - pi) < hpi) SCOS = -SCOS
       RETURN
      END  ! of SCOS


*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
* a rational function approximation to the arctangent, for all x and y
* accuracy is a few parts per million
*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
      FUNCTION SARC (y, x)  ! swift arctangent
       IMPLICIT NONE
       REAL y, x, SARC
       
       REAL o, p, q, q2, C(5), hpi, pi, piah, twopi
       SAVE C, hpi, pi, piah, twopi
       DATA C / 0.4224727843, 0.9999528872, 0.5615783790E-01, 
     &  0.7549406108, 1.000000000 /
       DATA hpi, pi, piah, twopi / 1.5707963268, 3.1415926536,
     &  4.7123889804, 6.2831853072 /
     
*           Java programmers, this is for you       
       IF (x) 10, 20, 30
  10   IF (y) 11, 12, 13
  11   IF (x < y) GO TO 50    
       GO TO 60
  12   SARC = pi 
       RETURN
  13   IF (-x < y) GO TO 40
       GO TO 50
  20   IF (y) 21, 22, 23
  21   SARC = piah
       RETURN
  22   SARC = 0.
       RETURN
  23   SARC = hpi
       RETURN
  30   IF (y) 31, 32, 33
  31   IF (x < -y) GO TO 60
       o = +1.
       p = twopi
       q = y / x
       GO TO 70
  32   SARC = 0.
       RETURN
  33   IF (x < y) GO TO 40
       o = +1.
       p = 0.
       q = y / x
       GO TO 70
  40   o = -1.
       p = hpi
       q = x / y
       GO TO 70
  50   o = +1.
       p = pi
       q = y / x
       GO TO 70
  60   o = -1.
       p = piah
       q = x / y
       GO TO 70

  70   CONTINUE 
       q2 = q **2
       q = ((C(1) * q2 + C(2)) * q) 
     &   / ((C(3) * q2 + C(4)) * q2 + C(5))
       SARC = SIGN(q, o) + p
       RETURN 
      END ! of SARC
      
      
*-----------------------------------------------------------------------
* RANG - Random Approximate Normal Generator
* 
* refer to:
*  Hand-Book on Statistical Distributions for experimentalists
*  Christian Walck, Particle Physics Group, University of Stockholm
*  10 September 2007
*
* origin:  G.E.P. Box and M.E. Muller, 1950.
*
* The random number generator is assumed to already be initialized.
*
*                            I/O LIST
*__Name_____Type_______I/O_______Description____________________________
*  z1       Real       Out       Standard normal random variable
*  z2       Real       Out       Standard normal random variable
*-----------------------------------------------------------------------
      FUNCTION RANG ()
       IMPLICIT NONE
       REAL RANG
       REAL URAND
       
       REAL R, S,               ! temp scalars
     &      Z1, Z2              ! pair of standard random variables
       REAL SLOG, SSIN, SCOS    ! function approximations
       LOGICAL ODD              ! toggle of how many times f'n has been called
       REAL TWOPI
       PARAMETER (TWOPI = 6.2831853071796)
       SAVE Z1, Z2, ODD
       DATA ODD / .FALSE. /

*          begin
       ODD = .NOT. ODD
       IF (ODD) THEN
  10     CONTINUE
         R = SQRT(-2. * SLOG(URAND()))
         IF (R .NE. R) GO TO 10         ! in case SLOG is inaccurate
         S = TWOPI * URAND()
         Z1 = R * SSIN(S)
         Z2 = R * SCOS(S)
         RANG = Z1
       ELSE
         RANG = Z2
       END IF
       RETURN
      END ! of RANG
*+++++++++++++++++++++++ end of file thrift.f ++++++++++++++++++++++++++
