*/////////////////////////////////////////////////////////////////////// * 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 ++++++++++++++++++++++++++