MODULE quadruple_precision ! Using Lahey's LF90, this MUST be compiled with optimization OFF. ! i.e. using LF90 quad -o0 ! OK with ELF90 (no optimization). ! Use quad_df for the Lahey/Fujitsu LF95 compiler. ! N.B. This is quadruple precision implemented in SOFTWARE, hence it is SLOW. ! This package has NOT been tested with other Fortran 90 compilers. ! The operations +-*/ & sqrt will probably work correctly with other compilers ! for PCs. Functions and routines which initialize quadruple-precision ! constants, particularly the trigonometric and exponential functions will ! only work if the compiler initializes double-precision constants from ! decimal code in exactly the same way as the Lahey compilers. ! The basic algorithms come from: ! Linnainmma, S. Software for doubled-precision floating-point computations, ! ACM Trans, Math. Software, vol. 7, pp.272-283, 1981. ! Programmer : Alan Miller ! e-mail: amiller@bigpond.net.au URL: www.ozemail.com.au/~milleraj ! Latest revision - 18 September 2002 ! 4 Aug 97. Added new algorithm for exponential. Takes half the time of the ! previous Taylor series, but errors 2.5 times larger. The old ! exponential is included here under the name exp_taylor for anyone ! who needs the extra accuracy. To use it instead of the new ! algorithm, change the module procedure name in INTERFACE exp from ! longexp to exp_taylor. ! 5 Aug 97. Found way to reduce cancellation errors in yesterday's algorithm ! for the exponential. Removed exp_taylor. ! 18 Aug 97. Added table of quadruple-precision constants. ! 8 Sept 97. Added str_quad which reads a character string and converts it to ! quadruple-precision form. ! 15 Oct 97. Added quad_str which takes a quadruple-precision value and outputs ! a character string containing its decimal value to 30 significant ! digits. ! Added overlays to the ** operator for quadruple-precision values. ! 15 Jan 98. Added ATAN2. ! 19 Jan 98. Added <, <=, ==, >= and >. ! 27 Dec 98. Added quad/real, quad/integer, integer/quad, real/quad, dp/quad, ! int+quad, real+quad, dp+quad, int-quad, real-quad, dp-quad, ! SUM, DOT_PRODUCT & MATMUL. ! 29 Dec 98. Correction to routine string_quad for strings of < 5 characters. ! 10 May 99. Added EPSILON for quad type. ! 5 Oct 99. log10e corrected. ! 15 Oct 99. Corrected function quad_pow_int. ! 18 Oct 99. Rewrote quad_pow_int to use the binary power method. ! 11 Nov 99. Added overlaying of assignments, e.g. quad = int, etc. ! 21 Jan 00. Added interface for EPSILON. ! 17 Feb 02. Corrected ACOS & ASIN for arguments close to +1 or -1. ! 18 Feb 02. Further improvement to ACOS. ! 18 Sep 02. Added reference to Linnainmaa in comments at beginning. IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 100) INTEGER, PARAMETER :: nbits = DIGITS(1.0_dp) REAL (dp), PARAMETER :: constant = 2._dp**(nbits + 11 - nbits/2) + 1._dp ! ^^ ! N.B. IEEE format for double precision uses 53 bits ! for the mantissa. PC floating-point processors ! use an extra 11 bits - hence that 11. REAL (dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp PRIVATE :: zero, one, constant, nbits TYPE :: quad REAL (dp) :: hi, lo END TYPE quad INTERFACE OPERATOR (*) MODULE PROCEDURE longmul MODULE PROCEDURE mult_quad_int MODULE PROCEDURE mult_int_quad MODULE PROCEDURE mult_quad_dp MODULE PROCEDURE mult_dp_quad MODULE PROCEDURE mult_quad_real MODULE PROCEDURE mult_real_quad END INTERFACE INTERFACE OPERATOR (/) MODULE PROCEDURE longdiv MODULE PROCEDURE div_quad_dp MODULE PROCEDURE div_quad_real MODULE PROCEDURE div_quad_int MODULE PROCEDURE div_int_quad MODULE PROCEDURE div_real_quad MODULE PROCEDURE div_dp_quad END INTERFACE INTERFACE OPERATOR (+) MODULE PROCEDURE longadd MODULE PROCEDURE quad_add_int MODULE PROCEDURE quad_add_real MODULE PROCEDURE quad_add_dp MODULE PROCEDURE int_add_quad MODULE PROCEDURE real_add_quad MODULE PROCEDURE dp_add_quad END INTERFACE INTERFACE OPERATOR (-) MODULE PROCEDURE longsub MODULE PROCEDURE quad_sub_int MODULE PROCEDURE quad_sub_real MODULE PROCEDURE quad_sub_dp MODULE PROCEDURE int_sub_quad MODULE PROCEDURE real_sub_quad MODULE PROCEDURE dp_sub_quad MODULE PROCEDURE negate_quad END INTERFACE INTERFACE ASSIGNMENT (=) MODULE PROCEDURE quad_eq_int MODULE PROCEDURE quad_eq_real MODULE PROCEDURE quad_eq_dp MODULE PROCEDURE int_eq_quad MODULE PROCEDURE real_eq_quad MODULE PROCEDURE dp_eq_quad END INTERFACE INTERFACE OPERATOR (**) MODULE PROCEDURE quad_pow_int MODULE PROCEDURE quad_pow_real MODULE PROCEDURE quad_pow_dp MODULE PROCEDURE quad_pow_quad END INTERFACE INTERFACE OPERATOR (<) MODULE PROCEDURE quad_lt END INTERFACE INTERFACE OPERATOR (<=) MODULE PROCEDURE quad_le END INTERFACE INTERFACE OPERATOR (==) MODULE PROCEDURE quad_eq END INTERFACE INTERFACE OPERATOR (>=) MODULE PROCEDURE quad_ge END INTERFACE INTERFACE OPERATOR (>) MODULE PROCEDURE quad_gt END INTERFACE INTERFACE SCALE MODULE PROCEDURE qscale END INTERFACE INTERFACE ABS MODULE PROCEDURE qabs END INTERFACE INTERFACE SQRT MODULE PROCEDURE longsqrt END INTERFACE INTERFACE LOG MODULE PROCEDURE longlog END INTERFACE INTERFACE EXP MODULE PROCEDURE longexp END INTERFACE INTERFACE SIN MODULE PROCEDURE longsin END INTERFACE INTERFACE COS MODULE PROCEDURE longcos END INTERFACE INTERFACE TAN MODULE PROCEDURE longtan END INTERFACE INTERFACE ASIN MODULE PROCEDURE longasin END INTERFACE INTERFACE ACOS MODULE PROCEDURE longacos END INTERFACE INTERFACE ATAN MODULE PROCEDURE longatan END INTERFACE INTERFACE ATAN2 MODULE PROCEDURE qatan2 END INTERFACE INTERFACE SUM MODULE PROCEDURE quad_sum END INTERFACE INTERFACE DOT_PRODUCT MODULE PROCEDURE quad_dot_product END INTERFACE INTERFACE MATMUL MODULE PROCEDURE q_matmul12 MODULE PROCEDURE q_matmul21 MODULE PROCEDURE q_matmul22 END INTERFACE INTERFACE EPSILON MODULE PROCEDURE q_epsilon END INTERFACE TYPE (quad), PARAMETER :: & pi = quad( 0.3141592653589793D+01, 0.1224646799147353D-15 ), & piby2 = quad( 0.1570796326794897D+01, -.3828568698926950D-15 ), & piby3 = quad( 0.1047197551196598D+01, -.3292527815701405D-15 ), & piby4 = quad( 0.7853981633974484D+00, -.8040613248383182D-16 ), & piby6 = quad( 0.5235987755982990D+00, -.1646263907850702D-15 ), & twopi = quad( 0.6283185307179586D+01, 0.2449293598294707D-15 ), & ln_pi = quad( 0.1144729885849400D+01, 0.2323105560877391D-15 ), & sqrtpi = quad( 0.1772453850905516D+01, -.7666586499825800D-16 ), & fact_pt5 = quad( 0.8862269254527582D+00, -.1493552349616447D-15 ), & sqrt2pi = quad( 0.2506628274631001D+01, -.6273750096546544D-15 ), & lnsqrt2pi = quad( 0.9189385332046728D+00, -.3878294158067242D-16 ), & one_on2pi = quad( 0.1591549430918953D+00, 0.4567181289366658D-16 ), & two_on_rtpi = quad( 0.1128379167095513D+01, -.4287537502368968D-15 ), & deg2rad = quad( 0.1745329251994330D-01, -.3174581724866598D-17 ), & rad2deg = quad( 0.5729577951308232D+02, -.1987849567057628D-14 ), & ln2 = quad( 0.6931471805599454D+00, -.8783183432405266D-16 ), & ln10 = quad( 0.2302585092994046D+01, -.2170756223382249D-15 ), & log2e = quad( 0.1442695040888964D+01, -.6457785410341630D-15 ), & log10e = quad( 0.4342944819032519D+00, -.5552037773430574D-16 ), & log2_10 = quad( 0.3321928094887362D+01, 0.1661617516973592D-15 ), & log10_2 = quad( 0.3010299956639812D+00, -.2803728127785171D-17 ), & euler = quad( 0.5772156649015330D+00, -.1159652176149463D-15 ), & e = quad( 0.2718281828459045D+01, 0.1445646891729250D-15 ), & sqrt2 = quad( 0.1414213562373095D+01, 0.1253716717905022D-15 ), & sqrt3 = quad( 0.1732050807568877D+01, 0.3223954471431004D-15 ), & sqrt10 = quad( 0.3162277660168380D+01, -.6348773795572286D-15 ) CONTAINS FUNCTION exactmul2(a, c) RESULT(ac) ! Procedure exactmul2, translated from Pascal, from: ! Linnainmaa, Seppo (1981). Software for doubled-precision floating-point ! computations. ACM Trans. on Math. Software (TOMS), 7, 272-283. REAL (dp), INTENT(IN) :: a, c TYPE (quad) :: ac ! Local variables REAL (dp) :: a1, a2, c1, c2, t t = constant * a a1 = (a - t) + t ! Lahey's optimization removes the brackets ! and sets a1 = a which defeats the whole point. a2 = a - a1 t = constant * c c1 = (c - t) + t c2 = c - c1 ac%hi = a * c ac%lo = (((a1 * c1 - ac%hi) + a1 * c2) + c1 * a2) + c2 * a2 RETURN END FUNCTION exactmul2 FUNCTION longmul(a, c) RESULT(ac) ! Procedure longmul, translated from Pascal, from: ! Linnainmaa, Seppo (1981). Software for doubled-precision floating-point ! computations. ACM Trans. on Math. Software (TOMS), 7, 272-283. TYPE (quad), INTENT(IN) :: a, c TYPE (quad) :: ac ! Local variables REAL (dp ) :: zz TYPE (quad) :: z z = exactmul2(a%hi, c%hi) zz = ((a%hi + a%lo) * c%lo + a%lo * c%hi) + z%lo ac%hi = z%hi + zz ac%lo = (z%hi - ac%hi) + zz RETURN END FUNCTION longmul FUNCTION mult_quad_int(a, i) RESULT(c) ! Multiply quadruple-precision number (a) by an integer (i). TYPE (quad), INTENT(IN) :: a INTEGER, INTENT(IN) :: i TYPE (quad) :: c IF (i == 0) THEN c = quad(zero, zero) ELSE IF (i == 1) THEN c = a ELSE IF (i == -1) THEN c = -a ELSE c = exactmul2(a%hi, DBLE(i)) + exactmul2(a%lo, DBLE(i)) END IF RETURN END FUNCTION mult_quad_int FUNCTION mult_int_quad(i, a) RESULT(c) ! Multiply quadruple-precision number (a) by an integer (i). INTEGER, INTENT(IN) :: i TYPE (quad), INTENT(IN) :: a TYPE (quad) :: c IF (i == 0) THEN c = quad(zero, zero) ELSE IF (i == 1) THEN c = a ELSE IF (i == -1) THEN c = -a ELSE c = exactmul2(a%hi, DBLE(i)) + exactmul2(a%lo, DBLE(i)) END IF RETURN END FUNCTION mult_int_quad FUNCTION mult_quad_dp(a, b) RESULT(c) ! Multiply a quadruple-precision number (a) by a double-precision number (b). TYPE (quad), INTENT(IN) :: a REAL (dp), INTENT(IN) :: b TYPE (quad) :: c ! Local variables TYPE (quad) :: z REAL (dp) :: zz z = exactmul2(a%hi, b) zz = a%lo * b + z%lo c%hi = z%hi + zz c%lo = (z%hi - c%hi) + zz RETURN END FUNCTION mult_quad_dp FUNCTION mult_quad_real(a, b) RESULT(c) ! Multiply a quadruple-precision number (a) by a real number (b). TYPE (quad), INTENT(IN) :: a REAL, INTENT(IN) :: b TYPE (quad) :: c ! Local variables TYPE (quad) :: z REAL (dp) :: zz z = exactmul2(a%hi, DBLE(b)) zz = a%lo * b + z%lo c%hi = z%hi + zz c%lo = (z%hi - c%hi) + zz RETURN END FUNCTION mult_quad_real FUNCTION mult_dp_quad(b, a) RESULT(c) ! Multiply a quadruple-precision number (a) by a double-precision number (b). TYPE (quad), INTENT(IN) :: a REAL (dp), INTENT(IN) :: b TYPE (quad) :: c ! Local variables TYPE (quad) :: z REAL (dp) :: zz z = exactmul2(a%hi, b) zz = a%lo * b + z%lo c%hi = z%hi + zz c%lo = (z%hi - c%hi) + zz RETURN END FUNCTION mult_dp_quad FUNCTION mult_real_quad(a, b) RESULT(c) ! Multiply a quadruple-precision number (a) by a double-precision number (b). REAL, INTENT(IN) :: a TYPE (quad), INTENT(IN) :: b TYPE (quad) :: c ! Local variables TYPE (quad) :: z REAL (dp) :: zz z = exactmul2(b%hi, DBLE(a)) zz = b%lo * a + z%lo c%hi = z%hi + zz c%lo = (z%hi - c%hi) + zz RETURN END FUNCTION mult_real_quad FUNCTION longdiv(a, c) RESULT(ac) ! Procedure longdiv, translated from Pascal, from: ! Linnainmaa, Seppo (1981). Software for doubled-precision floating-point ! computations. ACM Trans. on Math. Software (TOMS), 7, 272-283. TYPE (quad), INTENT(IN) :: a, c TYPE (quad) :: ac ! Local variables REAL (dp ) :: z, zz TYPE (quad) :: q z = a%hi / c%hi q = exactmul2(c%hi, z) zz = ((((a%hi - q%hi) - q%lo) + a%lo) - z*c%lo) / (c%hi + c%lo) ac%hi = z + zz ac%lo = (z - ac%hi) + zz RETURN END FUNCTION longdiv FUNCTION div_quad_dp(a, b) RESULT(c) ! Divide a quadruple-precision number (a) by a double-precision number (b) TYPE (quad), INTENT(IN) :: a REAL (dp), INTENT(IN) :: b TYPE (quad) :: c ! Local variables REAL (dp ) :: z, zz TYPE (quad) :: q z = a%hi / b q = exactmul2(b, z) zz = (((a%hi - q%hi) - q%lo) + a%lo) / b c%hi = z + zz c%lo = (z - c%hi) + zz RETURN END FUNCTION div_quad_dp FUNCTION div_quad_real(a, b) RESULT(c) ! Divide a quadruple-precision number (a) by a real number (b) TYPE (quad), INTENT(IN) :: a REAL, INTENT(IN) :: b TYPE (quad) :: c ! Local variables REAL (dp ) :: z, zz TYPE (quad) :: q z = a%hi / b q = exactmul2(DBLE(b), z) zz = (((a%hi - q%hi) - q%lo) + a%lo) / b c%hi = z + zz c%lo = (z - c%hi) + zz RETURN END FUNCTION div_quad_real FUNCTION div_quad_int(a, b) RESULT(c) ! Divide a quadruple-precision number (a) by an integer (b) TYPE (quad), INTENT(IN) :: a INTEGER, INTENT(IN) :: b TYPE (quad) :: c ! Local variables REAL (dp ) :: z, zz TYPE (quad) :: q z = a%hi / b q = exactmul2(DBLE(b), z) zz = (((a%hi - q%hi) - q%lo) + a%lo) / b c%hi = z + zz c%lo = (z - c%hi) + zz RETURN END FUNCTION div_quad_int FUNCTION div_int_quad(a, c) RESULT(ac) ! Procedure longdiv, translated from Pascal, from: ! Linnainmaa, Seppo (1981). Software for doubled-precision floating-point ! computations. ACM Trans. on Math. Software (TOMS), 7, 272-283. INTEGER, INTENT(IN) :: a TYPE (quad), INTENT(IN) :: c TYPE (quad) :: ac ! Local variables REAL (dp ) :: z, zz TYPE (quad) :: q z = DBLE(a) / c%hi q = exactmul2(c%hi, z) zz = (((a - q%hi) - q%lo) - z*c%lo) / (c%hi + c%lo) ac%hi = z + zz ac%lo = (z - ac%hi) + zz RETURN END FUNCTION div_int_quad FUNCTION div_real_quad(a, c) RESULT(ac) ! Procedure longdiv, translated from Pascal, from: ! Linnainmaa, Seppo (1981). Software for doubled-precision floating-point ! computations. ACM Trans. on Math. Software (TOMS), 7, 272-283. REAL, INTENT(IN) :: a TYPE (quad), INTENT(IN) :: c TYPE (quad) :: ac ! Local variables REAL (dp ) :: z, zz TYPE (quad) :: q z = DBLE(a) / c%hi q = exactmul2(c%hi, z) zz = (((a - q%hi) - q%lo) - z*c%lo) / (c%hi + c%lo) ac%hi = z + zz ac%lo = (z - ac%hi) + zz RETURN END FUNCTION div_real_quad FUNCTION div_dp_quad(a, c) RESULT(ac) ! Procedure longdiv, translated from Pascal, from: ! Linnainmaa, Seppo (1981). Software for doubled-precision floating-point ! computations. ACM Trans. on Math. Software (TOMS), 7, 272-283. REAL (dp), INTENT(IN) :: a TYPE (quad), INTENT(IN) :: c TYPE (quad) :: ac ! Local variables REAL (dp ) :: z, zz TYPE (quad) :: q z = a / c%hi q = exactmul2(c%hi, z) zz = (((a - q%hi) - q%lo) - z*c%lo) / (c%hi + c%lo) ac%hi = z + zz ac%lo = (z - ac%hi) + zz RETURN END FUNCTION div_dp_quad FUNCTION longadd(a, c) RESULT(ac) ! Procedure longadd, translated from Pascal, from: ! Linnainmaa, Seppo (1981). Software for doubled-precision floating-point ! computations. ACM Trans. on Math. Software (TOMS), 7, 272-283. TYPE (quad), INTENT(IN) :: a, c TYPE (quad) :: ac ! Local variables REAL (dp ) :: z, q, zz z = a%hi + c%hi q = a%hi - z zz = (((q + c%hi) + (a%hi - (q + z))) + a%lo) + c%lo ac%hi = z + zz ac%lo = (z - ac%hi) + zz RETURN END FUNCTION longadd FUNCTION quad_add_int(a, c) RESULT(ac) TYPE (quad), INTENT(IN) :: a INTEGER, INTENT(IN) :: c TYPE (quad) :: ac ! Local variables REAL (dp ) :: z, q, zz z = a%hi + c q = a%hi - z zz = (((q + c) + (a%hi - (q + z))) + a%lo) ac%hi = z + zz ac%lo = (z - ac%hi) + zz RETURN END FUNCTION quad_add_int FUNCTION quad_add_real(a, c) RESULT(ac) TYPE (quad), INTENT(IN) :: a REAL, INTENT(IN) :: c TYPE (quad) :: ac ! Local variables REAL (dp ) :: z, q, zz z = a%hi + c q = a%hi - z zz = (((q + c) + (a%hi - (q + z))) + a%lo) ac%hi = z + zz ac%lo = (z - ac%hi) + zz RETURN END FUNCTION quad_add_real FUNCTION quad_add_dp(a, c) RESULT(ac) TYPE (quad), INTENT(IN) :: a REAL (dp), INTENT(IN) :: c TYPE (quad) :: ac ! Local variables REAL (dp ) :: z, q, zz z = a%hi + c q = a%hi - z zz = (((q + c) + (a%hi - (q + z))) + a%lo) ac%hi = z + zz ac%lo = (z - ac%hi) + zz RETURN END FUNCTION quad_add_dp FUNCTION int_add_quad(c, a) RESULT(ac) INTEGER, INTENT(IN) :: c TYPE (quad), INTENT(IN) :: a TYPE (quad) :: ac ! Local variables REAL (dp ) :: z, q, zz z = a%hi + c q = a%hi - z zz = (((q + c) + (a%hi - (q + z))) + a%lo) ac%hi = z + zz ac%lo = (z - ac%hi) + zz RETURN END FUNCTION int_add_quad FUNCTION real_add_quad(c, a) RESULT(ac) REAL, INTENT(IN) :: c TYPE (quad), INTENT(IN) :: a TYPE (quad) :: ac ! Local variables REAL (dp ) :: z, q, zz z = a%hi + c q = a%hi - z zz = (((q + c) + (a%hi - (q + z))) + a%lo) ac%hi = z + zz ac%lo = (z - ac%hi) + zz RETURN END FUNCTION real_add_quad FUNCTION dp_add_quad(c, a) RESULT(ac) REAL (dp), INTENT(IN) :: c TYPE (quad), INTENT(IN) :: a TYPE (quad) :: ac ! Local variables REAL (dp ) :: z, q, zz z = a%hi + c q = a%hi - z zz = (((q + c) + (a%hi - (q + z))) + a%lo) ac%hi = z + zz ac%lo = (z - ac%hi) + zz RETURN END FUNCTION dp_add_quad FUNCTION longsub(a, c) RESULT(ac) ! Adapted from longadd by changing signs of c. ! Linnainmaa, Seppo (1981). Software for doubled-precision floating-point ! computations. ACM Trans. on Math. Software (TOMS), 7, 272-283. TYPE (quad), INTENT(IN) :: a, c TYPE (quad) :: ac ! Local variables REAL (dp ) :: z, q, zz z = a%hi - c%hi q = a%hi - z zz = (((q - c%hi) + (a%hi - (q + z))) + a%lo) - c%lo ac%hi = z + zz ac%lo = (z - ac%hi) + zz RETURN END FUNCTION longsub FUNCTION quad_sub_int(a, c) RESULT(ac) TYPE (quad), INTENT(IN) :: a INTEGER, INTENT(IN) :: c TYPE (quad) :: ac ! Local variables REAL (dp ) :: z, q, zz z = a%hi - c q = a%hi - z zz = (((q - c) + (a%hi - (q + z))) + a%lo) ac%hi = z + zz ac%lo = (z - ac%hi) + zz RETURN END FUNCTION quad_sub_int FUNCTION quad_sub_real(a, c) RESULT(ac) TYPE (quad), INTENT(IN) :: a REAL, INTENT(IN) :: c TYPE (quad) :: ac ! Local variables REAL (dp ) :: z, q, zz z = a%hi - c q = a%hi - z zz = (((q - c) + (a%hi - (q + z))) + a%lo) ac%hi = z + zz ac%lo = (z - ac%hi) + zz RETURN END FUNCTION quad_sub_real FUNCTION quad_sub_dp(a, c) RESULT(ac) TYPE (quad), INTENT(IN) :: a REAL (dp), INTENT(IN) :: c TYPE (quad) :: ac ! Local variables REAL (dp ) :: z, q, zz z = a%hi - c q = a%hi - z zz = (((q - c) + (a%hi - (q + z))) + a%lo) ac%hi = z + zz ac%lo = (z - ac%hi) + zz RETURN END FUNCTION quad_sub_dp FUNCTION int_sub_quad(a, c) RESULT(ac) INTEGER, INTENT(IN) :: a TYPE (quad), INTENT(IN) :: c TYPE (quad) :: ac ! Local variables REAL (dp ) :: z, q, zz z = DBLE(a) - c%hi q = DBLE(a) - z zz = ((q - c%hi) + (DBLE(a) - (q + z))) - c%lo ac%hi = z + zz ac%lo = (z - ac%hi) + zz RETURN END FUNCTION int_sub_quad FUNCTION real_sub_quad(a, c) RESULT(ac) REAL, INTENT(IN) :: a TYPE (quad), INTENT(IN) :: c TYPE (quad) :: ac ! Local variables REAL (dp ) :: z, q, zz z = a - c%hi q = a - z zz = ((q - c%hi) + (a - (q + z))) - c%lo ac%hi = z + zz ac%lo = (z - ac%hi) + zz RETURN END FUNCTION real_sub_quad FUNCTION dp_sub_quad(a, c) RESULT(ac) REAL (dp), INTENT(IN) :: a TYPE (quad), INTENT(IN) :: c TYPE (quad) :: ac ! Local variables REAL (dp ) :: z, q, zz z = a - c%hi q = a - z zz = ((q - c%hi) + (a - (q + z))) - c%lo ac%hi = z + zz ac%lo = (z - ac%hi) + zz RETURN END FUNCTION dp_sub_quad FUNCTION negate_quad(a) RESULT(b) ! Change the sign of a quadruple-precision number. ! In many cases, a & b will occupy the same locations. TYPE (quad), INTENT(IN) :: a TYPE (quad) :: b b%hi = -a%hi b%lo = -a%lo RETURN END FUNCTION negate_quad SUBROUTINE quad_eq_int(a, i) ! Assignment TYPE (quad), INTENT(OUT) :: a INTEGER, INTENT(IN) :: i a%hi = i a%lo = 0 RETURN END SUBROUTINE quad_eq_int SUBROUTINE quad_eq_real(a, r) ! Assignment TYPE (quad), INTENT(OUT) :: a REAL, INTENT(IN) :: r a%hi = r a%lo = zero RETURN END SUBROUTINE quad_eq_real SUBROUTINE quad_eq_dp(a, d) ! Assignment TYPE (quad), INTENT(OUT) :: a REAL (dp), INTENT(IN) :: d a%hi = d a%lo = zero RETURN END SUBROUTINE quad_eq_dp SUBROUTINE int_eq_quad(i, a) ! Assignment INTEGER, INTENT(OUT) :: i TYPE (quad), INTENT(IN) :: a i = a%hi RETURN END SUBROUTINE int_eq_quad SUBROUTINE real_eq_quad(r, a) ! Assignment REAL, INTENT(OUT) :: r TYPE (quad), INTENT(IN) :: a r = a%hi RETURN END SUBROUTINE real_eq_quad SUBROUTINE dp_eq_quad(d, a) ! Assignment REAL (dp), INTENT(OUT) :: d TYPE (quad), INTENT(IN) :: a d = a%hi RETURN END SUBROUTINE dp_eq_quad FUNCTION quad_lt(x, y) RESULT(is_it) ! Comparison of 2 logical numbers TYPE (quad), INTENT(IN) :: x, y LOGICAL :: is_it ! Local variable TYPE (quad) :: diff diff = x - y is_it = (diff%hi < zero) RETURN END FUNCTION quad_lt FUNCTION quad_le(x, y) RESULT(is_it) ! Comparison of 2 logical numbers TYPE (quad), INTENT(IN) :: x, y LOGICAL :: is_it ! Local variable TYPE (quad) :: diff diff = x - y is_it = .NOT. (diff%hi > zero) RETURN END FUNCTION quad_le FUNCTION quad_eq(x, y) RESULT(is_it) ! Comparison of 2 logical numbers TYPE (quad), INTENT(IN) :: x, y LOGICAL :: is_it ! Local variable TYPE (quad) :: diff diff = x - y is_it = (diff%hi == zero) RETURN END FUNCTION quad_eq FUNCTION quad_ge(x, y) RESULT(is_it) ! Comparison of 2 logical numbers TYPE (quad), INTENT(IN) :: x, y LOGICAL :: is_it ! Local variable TYPE (quad) :: diff diff = x - y is_it = .NOT. (diff%hi < zero) RETURN END FUNCTION quad_ge FUNCTION quad_gt(x, y) RESULT(is_it) ! Comparison of 2 logical numbers TYPE (quad), INTENT(IN) :: x, y LOGICAL :: is_it ! Local variable TYPE (quad) :: diff diff = x - y is_it = (diff%hi > zero) RETURN END FUNCTION quad_gt FUNCTION quad_pow_int(a, i) RESULT(b) ! Raise a quadruple=precision number (a) to a power. TYPE (quad), INTENT(IN) :: a INTEGER, INTENT(IN) :: i TYPE (quad) :: b ! Local variables INTEGER :: ia, j, emax TYPE (quad) :: power LOGICAL :: first SELECT CASE (i) CASE (0) b = quad(one, zero) CASE (-1023:-1, 1:1023) ia = ABS(i) first = .TRUE. power = a emax = MAXEXPONENT(one) DO j = 0, 9 IF (BTEST(ia, j)) THEN IF (first) THEN b = power first = .FALSE. ELSE IF (EXPONENT(power%hi) + EXPONENT(b%hi) > emax) THEN WRITE(*, *) '** Exponential overflow in routine QUAD_POW_INT **' RETURN END IF b = b * power END IF ia = IBCLR(ia, j) IF (ia == 0) EXIT END IF IF (EXPONENT(power%hi) > emax/2) THEN WRITE(*, *) '** Exponential overflow in routine QUAD_POW_INT **' RETURN END IF power = power * power END DO IF (i < 0) b = quad(one, zero) / b CASE DEFAULT IF (a%hi < zero) THEN IF (i * LOG(-a%hi) > LOG( HUGE(one) )) THEN WRITE(*, *) '** Exponential overflow in routine QUAD_POW_INT **' RETURN ELSE IF (i * LOG(-a%hi) < LOG( TINY(one) )) THEN b = quad(zero, zero) ELSE ia = ABS(i) IF (BTEST(ia,0)) THEN b = - EXP( LOG(-a)*i ) ! i is odd (test least sign. bit) ELSE b = EXP( LOG(-a)*i ) ! i is even END IF END IF ELSE IF (a%hi > zero) THEN IF (i * LOG(a%hi) > LOG( HUGE(one) )) THEN WRITE(*, *) '** Exponential overflow in routine QUAD_POW_INT **' RETURN ELSE IF (i * LOG(a%hi) < LOG( TINY(one) )) THEN b = quad(zero, zero) ELSE b = EXP( LOG(a)*i ) END IF ELSE b = quad(zero, zero) END IF END SELECT RETURN END FUNCTION quad_pow_int FUNCTION quad_pow_real(a, r) RESULT(b) ! Raise a quadruple=precision number (a) to a power. TYPE (quad), INTENT(IN) :: a REAL, INTENT(IN) :: r TYPE (quad) :: b IF (a%hi < zero) THEN WRITE(*, *) & ' *** Error: attempt to raise negative quad. prec. number to a real power ***' b = quad(zero, zero) ELSE IF (a%hi > zero) THEN b = EXP( LOG(a)*r ) ELSE b = quad(zero, zero) END IF RETURN END FUNCTION quad_pow_real FUNCTION quad_pow_dp(a, d) RESULT(b) ! Raise a quadruple=precision number (a) to a power. TYPE (quad), INTENT(IN) :: a REAL (dp), INTENT(IN) :: d TYPE (quad) :: b IF (a%hi < zero) THEN WRITE(*, *) & ' *** Error: attempt to raise negative quad. prec. number to a real power ***' b = quad(zero, zero) ELSE IF (a%hi > zero) THEN b = EXP( LOG(a)*d ) ELSE b = quad(zero, zero) END IF RETURN END FUNCTION quad_pow_dp FUNCTION quad_pow_quad(a, q) RESULT(b) ! Raise a quadruple=precision number (a) to a power. TYPE (quad), INTENT(IN) :: a, q TYPE (quad) :: b IF (a%hi < zero) THEN WRITE(*, *) & ' *** Error: attempt to raise negative quad. prec. number to a real power ***' b = quad(zero, zero) ELSE IF (a%hi > zero) THEN b = EXP( LOG(a)*q ) ELSE b = quad(zero, zero) END IF RETURN END FUNCTION quad_pow_quad FUNCTION qscale(a, i) RESULT(b) ! Multiply a by 2^i TYPE (quad), INTENT(IN) :: a INTEGER, INTENT(IN) :: i TYPE (quad) :: b b%hi = SCALE(a%hi, i) b%lo = SCALE(a%lo, i) RETURN END FUNCTION qscale FUNCTION qabs(a) RESULT(b) ! Absolute value of a quadruple-precision number TYPE (quad), INTENT(IN) :: a TYPE (quad) :: b IF (a%hi < zero) THEN b%hi = - a%hi b%lo = - a%lo ELSE b%hi = a%hi b%lo = a%lo END IF RETURN END FUNCTION qabs FUNCTION longsqrt(a) RESULT(b) ! This is modified from procedure sqrt2 of: ! Dekker, T.J. (1971). 'A floating-point technique for extending the ! available precision', Numer. Math., 18, 224-242. TYPE (quad), INTENT(IN) :: a TYPE (quad) :: b ! Local variables REAL (dp) :: t, res TYPE (quad) :: tt ! Check that ahi >= 0. IF (a%hi < 0._dp) THEN WRITE(*, *) ' *** Negative argument for longsqrt ***' RETURN ELSE IF (a%hi == 0.d0) THEN b%hi = 0._dp b%lo = 0._dp RETURN END IF ! First approximation is t = sqrt(a). t = SQRT(a%hi) tt = exactmul2(t, t) res = t + (((a%hi - tt%hi) - tt%lo) + a%lo) * 0.5_dp / t b%lo = (t - res) + (((a%hi - tt%hi) - tt%lo) + a%lo) * 0.5_dp / t b%hi = res RETURN END FUNCTION longsqrt FUNCTION longlog(x) RESULT(y) ! Quadruple-precision logarithm to base e ! Halley's algorithm using double-precision logarithm as starting value. ! Solves: y ! f(y) = e - x = 0 TYPE (quad), INTENT(IN) :: x TYPE (quad) :: y ! Local variables TYPE (quad) :: expy, f y%hi = LOG(x%hi) y%lo = 0._dp expy = EXP(y) f = expy - x f = SCALE(f, 1) y = y - f / (expy + x) RETURN END FUNCTION longlog FUNCTION longexp(x) RESULT(y) ! Calculate a quadruple-precision exponential ! Method: ! x x.log2(e) nint[x.log2(e)] + frac[x.log2(e)] ! e = 2 = 2 ! ! iy fy ! = 2 . 2 ! Then ! fy y.ln(2) ! 2 = e ! ! Now y.ln(2) will be less than 0.3466 in absolute value. ! This is halved and a Pade approximation is used to approximate e^x over ! the region (-0.1733, +0.1733). This approximation is then squared. ! WARNING: No overflow checks! TYPE (quad), INTENT(IN) :: x TYPE (quad) :: y ! Local variables TYPE (quad) :: temp, ysq, sum1, sum2 INTEGER :: iy y = x / ln2 iy = NINT(REAL(y%hi)) y = (y - DBLE(iy)) * ln2 y = SCALE(y, -1) ! The Pade series is: ! p0 + p1.y + p2.y^2 + p3.y^3 + ... + p9.y^9 ! ------------------------------------------ ! p0 - p1.y + p2.y^2 - p3.y^3 + ... - p9.y^9 ! ! sum1 is the sum of the odd powers, sum2 is the sum of the even powers ysq = y * y sum1 = y * ((((ysq + 3960.)*ysq + 2162160._dp)*ysq + 302702400._dp)*ysq + & 8821612800._dp) sum2 = (((90.*ysq + 110880.)*ysq + 30270240._dp)*ysq + 2075673600._dp)*ysq + & 17643225600._dp ! sum2 + sum1 2.sum1 ! Now approximation = ----------- = 1 + ----------- = 1 + 2.temp ! sum2 - sum1 sum2 - sum1 ! ! Then (1 + 2.temp)^2 = 4.temp.(1 + temp) + 1 temp = sum1 / (sum2 - sum1) y = temp * (temp + 1._dp) y = SCALE(y, 2) y = y + 1._dp y = SCALE(y, iy) RETURN END FUNCTION longexp SUBROUTINE longmodr(a, b, n, rem) ! Extended arithmetic calculation of the 'rounded' modulus: ! a = n.b + rem ! where all quantities are in quadruple-precision, except the integer ! number of multiples, n. The absolute value of the remainder (rem) ! is not greater than b/2. ! The result is exact. rem may occupy the same location as either input. ! Programmer: Alan Miller ! Latest revision - 11 September 1986 ! Fortran version - 4 December 1996 TYPE (quad), INTENT(IN) :: a, b INTEGER, INTENT(OUT) :: n TYPE (quad), INTENT(OUT) :: rem ! Local variables TYPE (quad) :: temp ! Check that b%hi .ne. 0 IF (b%hi == 0._dp) THEN WRITE(*, *) ' *** Error in longmodr - 3rd argument zero ***' RETURN END IF ! Calculate n. temp = a / b n = NINT(REAL(temp%hi)) ! Calculate remainder preserving full accuracy. temp = exactmul2(DBLE(n), b%hi) rem%hi = a%hi rem%lo = zero temp = rem - temp rem%hi = a%lo rem%lo = zero temp = rem + temp rem = exactmul2(DBLE(n), b%lo) rem = temp - rem RETURN END SUBROUTINE longmodr ! Extended accuracy arithmetic sine, cosine & tangent (about 31 decimals). ! Calculates b = sin, cos or tan (a), where all quantities are in ! quadruple-precision, using table look-up and a Taylor series expansion. ! The result may occupy the same locations as the input value. ! Much of the code is common to all three functions, and this is in a ! subroutine longcst. FUNCTION longsin(a) RESULT(b) TYPE (quad), INTENT(IN) :: a TYPE (quad) :: b ! Local variables LOGICAL :: sine, cosine, tangent ! Set logical variables for sine function. sine = .true. cosine = .false. tangent = .false. CALL longcst(a, b, sine, cosine, tangent) RETURN END FUNCTION longsin FUNCTION longcos(a) RESULT(b) TYPE (quad), INTENT(IN) :: a TYPE (quad) :: b ! Local variables LOGICAL :: sine, cosine, tangent ! Set logical variables for sine function. sine = .false. cosine = .true. tangent = .false. CALL longcst(a, b, sine, cosine, tangent) RETURN END FUNCTION longcos FUNCTION longtan(a) RESULT(b) TYPE (quad), INTENT(IN) :: a TYPE (quad) :: b ! Local variables LOGICAL :: sine, cosine, tangent ! Set logical variables for sine function. sine = .false. cosine = .false. tangent = .true. CALL longcst(a, b, sine, cosine, tangent) RETURN END FUNCTION longtan SUBROUTINE longcst(a, b, sine, cosine, tangent) TYPE (quad), INTENT(IN) :: a TYPE (quad), INTENT(OUT) :: b LOGICAL, INTENT(IN) :: sine, cosine, tangent ! Local variables LOGICAL :: pos TYPE (quad) :: d, term, temp, angle, piby40, sum1, sum2, sin INTEGER :: npi, ipt, i REAL (dp) :: tol15 = 1.E-15_dp, tol30 = 1.E-30_dp ! sin(i.pi/40), i = 0(1)20 REAL (dp) :: table(2, 0:20) = RESHAPE( (/ & 0.0000000000000000E+00_dp, 0.0000000000000000E+00_dp, & 0.7845909572784494E-01_dp, 0.1464397249532491E-17_dp, & 0.1564344650402309E+00_dp, -.2770509565052586E-16_dp, & 0.2334453638559054E+00_dp, 0.2058612230858154E-16_dp, & 0.3090169943749475E+00_dp, -.8267172724967036E-16_dp, & 0.3826834323650898E+00_dp, -.1005077269646159E-16_dp, & 0.4539904997395468E+00_dp, -.1292033036231312E-16_dp, & 0.5224985647159488E+00_dp, 0.6606794454708078E-16_dp, & 0.5877852522924732E+00_dp, -.1189570533007057E-15_dp, & 0.6494480483301838E+00_dp, -.1134903961116171E-15_dp, & 0.7071067811865476E+00_dp, -.4833646656726458E-16_dp, & 0.7604059656000310E+00_dp, -.1036987135483477E-15_dp, & 0.8090169943749476E+00_dp, -.1381828784809282E-15_dp, & 0.8526401643540922E+00_dp, 0.4331886637554353E-16_dp, & 0.8910065241883680E+00_dp, -.1474714419679880E-15_dp, & 0.9238795325112868E+00_dp, -.9337725537817898E-16_dp, & 0.9510565162951536E+00_dp, -.7008780156242836E-16_dp, & 0.9723699203976766E+00_dp, 0.4478912629332321E-16_dp, & 0.9876883405951378E+00_dp, -.4416018005989794E-16_dp, & 0.9969173337331280E+00_dp, 0.1235153006196267E-16_dp, & 0.1000000000000000E+01_dp, 0.0000000000000000E+00_dp /), & (/ 2, 21 /) ) ! pi/40 piby40%hi = 0.7853981633974484E-01_dp piby40%lo = -.1081617080994607E-16_dp ! Reduce angle to range (-pi/2, +pi/2) by subtracting an integer multiple of pi. CALL longmodr(a, pi, npi, angle) ! Find nearest multiple of pi/40 to angle. CALL longmodr(angle, piby40, ipt, d) ! Sum 1 = 1 - d**2/2! + d**4/4! - d**6/6! + ... ! Sum 2 = d - d**3/3! + d**5/5! - d**7/7! + ... sum1%hi = zero sum1%lo = zero sum2%hi = zero sum2%lo = zero pos = .false. term = d i = 2 20 IF (ABS(term%hi) > tol15) THEN term = term * d ! Use quad. precision IF (i == 2 .OR. i == 4 .OR. i == 8) THEN term%hi = term%hi / i term%lo = term%lo / i ELSE term = term / DBLE(i) END IF IF (pos) THEN sum1 = sum1 + term ELSE sum1 = sum1 - term END IF ELSE term%hi = term%hi * d%hi / DBLE(i) ! Double prec. adequate IF (pos) THEN sum1%lo = sum1%lo + term%hi ELSE sum1%lo = sum1%lo - term%hi END IF END IF ! Repeat for sum2 i = i + 1 IF (ABS(term%hi) > tol15) THEN term = term * d / DBLE(i) ! Use quad. precision IF (pos) THEN sum2 = sum2 + term ELSE sum2 = sum2 - term END IF ELSE term%hi = term%hi * d%hi / DBLE(i) ! Double prec. adequate IF (pos) THEN sum2%lo = sum2%lo + term%hi ELSE sum2%lo = sum2%lo - term%hi END IF END IF i = i + 1 pos = .NOT. pos IF (ABS(term%hi) > tol30) GO TO 20 sum1 = sum1 + 1._dp ! Now add the 1st terms sum2 = sum2 + d ! for max. accuracy ! Construct sine, cosine or tangent. ! Sine first. sin(angle + d) = sin(angle).cos(d) + cos(angle).sin(d) IF (sine .OR. tangent) THEN IF (ipt >= 0) THEN temp%hi = table(1, ipt) temp%lo = table(2, ipt) ELSE temp%hi = - table(1, -ipt) temp%lo = - table(2, -ipt) END IF b = sum1 * temp IF (ipt >= 0) THEN temp%hi = table(1, 20-ipt) temp%lo = table(2, 20-ipt) ELSE temp%hi = table(1, 20+ipt) temp%lo = table(2, 20+ipt) END IF b = b + sum2 * temp IF (npi /= 2*(npi/2)) THEN b = -b END IF IF (tangent) THEN sin = b END IF END IF ! Cosine or tangent. IF (cosine .OR. tangent) THEN IF (ipt >= 0) THEN temp%hi = table(1, ipt) temp%lo = table(2, ipt) ELSE temp%hi = - table(1, -ipt) temp%lo = - table(2, -ipt) END IF b = sum2 * temp IF (ipt >= 0) THEN temp%hi = table(1, 20-ipt) temp%lo = table(2, 20-ipt) ELSE temp%hi = table(1, 20+ipt) temp%lo = table(2, 20+ipt) END IF b = sum1 * temp - b IF (npi /= 2*(npi/2)) THEN b = -b END IF END IF ! Tangent. IF (tangent) THEN ! Check that bhi .ne. 0 IF (b%hi == 0.d0) THEN WRITE(*, *) ' *** Infinite tangent - routine longcst ***' b%hi = HUGE(one) b%lo = 0.d0 RETURN END IF b = sin / b END IF RETURN END SUBROUTINE longcst FUNCTION longasin(a) RESULT(b) ! Quadratic-precision arc sine (about 31 decimals). ! One Newton-Raphson iteration to solve: f(b) = sin(b) - a = 0, ! except when a close to -1 or +1. ! The result (b) may occupy the same location as the input values (a). ! Use ACOS when |a| is close to 1. TYPE (quad), INTENT(IN) :: a TYPE (quad) :: b ! Local variables TYPE (quad) :: y, c ! Check that -1 <= a%hi <= +1. IF (a%hi < -1._dp .OR. a%hi > 1._dp) THEN WRITE(*, *) ' *** Argument outside range for longasin ***' RETURN END IF IF (ABS(a%hi) < 0.866) THEN ! First approximation is y = asin(a). ! Quadruple-precision result is y - [sin(y) - a]/cos(y). y%hi = ASIN(a%hi) y%lo = zero b = y + (a - SIN(y)) / COS(y%hi) ELSE ! Calculate acos(c) where c = sqrt(1 - a^2) c = SQRT(one - a*a) y%hi = ACOS(c%hi) y%lo = zero b = y + (COS(y) - c) / SIN(y%hi) IF (a%hi < zero) b = -b END IF RETURN END FUNCTION longasin FUNCTION longacos(a) RESULT(b) ! Quadratic-precision arc cosine (about 31 decimals). ! Newton-Raphson iteration to solve: f(b) = cos(b) - a = 0. ! The result (b) may occupy the same location as the input values (a). ! When |a| is near 1, use formula from p.175 of ! `Software Manual for the Elementary Functions' by W.J. Cody, Jr. & ! W. Waite, Prentice-Hall, 1980. TYPE (quad), INTENT(IN) :: a TYPE (quad) :: b ! Local variables TYPE (quad) :: y, c ! Check that -1 <= a%hi <= +1. IF (a%hi < -1._dp .OR. a%hi > 1._dp) THEN WRITE(*, *) ' *** Argument outside range for longacos ***' RETURN END IF IF (ABS(a%hi) < 0.866) THEN ! First approximation is y = acos(a). ! Quadruple-precision result is y + [cos(y) - a]/sin(y). y%hi = ACOS(a%hi) y%lo = 0._dp b = y + (COS(y) - a) / SIN(y%hi) ELSE ! Calculate 2.asin(c) where c = sqrt([1 - |a|]/2) c = SQRT((one - ABS(a))/2) y%hi = ASIN(c%hi) y%lo = zero b = (y - (SIN(y) - c) / COS(y%hi))*2 IF (a%hi < zero) b = pi - b END IF RETURN END FUNCTION longacos FUNCTION longatan(a) RESULT(b) ! Quadratic-precision arc tangent (about 31 decimals). ! Newton-Raphson iteration to solve: f(b) = tan(b) - a = 0. ! The result (b) may occupy the same location as the input values (a). TYPE (quad), INTENT(IN) :: a TYPE (quad) :: b ! Local variables TYPE (quad) :: y ! First approximation is y = atan(a). ! Quadruple-precision result is y - [tan(y) - a] * cos(y)**2. y%hi = ATAN(a%hi) y%lo = 0._dp b = y - (TAN(y) - a) * (COS(y%hi))**2 RETURN END FUNCTION longatan FUNCTION qatan2(y, x) RESULT(b) ! Quadratic-precision arc tangent (about 31 decimals). ! As for arc tangent (y/x) except that the result is in the range ! -pi < ATAN2 <= pi. ! The signs of x and y determine the quadrant. TYPE (quad), INTENT(IN) :: y, x TYPE (quad) :: b ! Local variables TYPE (quad) :: z ! First approximation is z = atan2(y, x). ! Quadruple-precision result is z - [tan(z) - (y/x)] * cos(z)**2. z%hi = ATAN2(y%hi, x%hi) z%lo = 0._dp IF (x%hi == zero) THEN b = z ELSE b = z - (TAN(z) - y/x) * (COS(z%hi))**2 END IF RETURN END FUNCTION qatan2 FUNCTION quad_sum(a) RESULT(s) ! Quadruple-precision SUM TYPE (quad), INTENT(IN), DIMENSION(:) :: a TYPE (quad) :: s ! Local variables INTEGER :: i s = quad(zero, zero) DO i = 1, SIZE(a) s = s + a(i) END DO RETURN END FUNCTION quad_sum FUNCTION quad_dot_product(a, b) RESULT(ab) ! Quadruple-precision DOT_PRODUCT TYPE (quad), INTENT(IN), DIMENSION(:) :: a, b TYPE (quad) :: ab ! Local variables INTEGER :: i, n ab = quad(zero, zero) n = SIZE(a) IF (n /= SIZE(b)) THEN WRITE(*, *) ' ** Error invoking DOT_PRODUCT - different argument sizes **' WRITE(*, '(a, i10, a, i10)') ' Size of 1st argument = ', n, & ' Size of 2nd argument = ', SIZE(b) RETURN END IF DO i = 1, n ab = ab + a(i)*b(i) END DO RETURN END FUNCTION quad_dot_product FUNCTION q_matmul12(a, b) RESULT(ab) ! Quadruple-precision MATMUL ! Rank of A = 1, rank of B = 2 TYPE (quad), INTENT(IN), DIMENSION(:) :: a TYPE (quad), INTENT(IN), DIMENSION(:,:) :: b TYPE (quad), DIMENSION( SIZE(b,2) ) :: ab ! Local variables INTEGER :: j, na, nb1, nb2 ! Check dimensions na = SIZE(a) nb1 = SIZE(b, 1) nb2 = SIZE(b, 2) IF (na /= nb1) THEN WRITE(*, *) ' ** Incompatible dimensions for quad-prec. MATMUL' RETURN END IF DO j = 1, nb2 ab(j) = quad_dot_product( a, b(:,j) ) END DO RETURN END FUNCTION q_matmul12 FUNCTION q_matmul21(a, b) RESULT(ab) ! Quadruple-precision MATMUL ! Rank of A = 2, rank of B = 1 TYPE (quad), INTENT(IN), DIMENSION(:,:) :: a TYPE (quad), INTENT(IN), DIMENSION(:) :: b TYPE (quad), DIMENSION( SIZE(a,1) ) :: ab ! Local variables INTEGER :: i, na1, na2, nb ! Check dimensions na1 = SIZE(a, 1) na2 = SIZE(a, 2) nb = SIZE(b) IF (na2 /= nb) THEN WRITE(*, *) ' ** Incompatible dimensions for quad-prec. MATMUL' RETURN END IF DO i = 1, na1 ab(i) = quad_dot_product( a(i,:), b ) END DO RETURN END FUNCTION q_matmul21 FUNCTION q_matmul22(a, b) RESULT(ab) ! Quadruple-precision MATMUL ! Rank of A = 2, rank of B = 2 TYPE (quad), INTENT(IN), DIMENSION(:,:) :: a TYPE (quad), INTENT(IN), DIMENSION(:,:) :: b TYPE (quad), DIMENSION(SIZE(a,1),SIZE(b,2)) :: ab ! Local variables INTEGER :: i, j, na1, na2, nb1, nb2 ! Check dimensions na1 = SIZE(a, 1) na2 = SIZE(a, 2) nb1 = SIZE(b, 1) nb2 = SIZE(b, 2) IF (na2 /= nb1) THEN WRITE(*, *) ' ** Incompatible dimensions for quad-prec. MATMUL' RETURN END IF DO i = 1, na1 DO j = 1, nb2 ab(i,j) = quad_dot_product( a(i,:), b(:,j) ) END DO END DO RETURN END FUNCTION q_matmul22 SUBROUTINE string_quad(string, value, ier) ! Convert a character string to a quadruple-precision quantity. ! Error indicator ier = 0 if value OK ! = 1 if string has > 50 characters and no decimal point ! in about the first 45 ! = 2 if a non-numeric character is read in the mantissa ! part CHARACTER (LEN=*), INTENT(IN) :: string TYPE (quad), INTENT(OUT) :: value INTEGER, INTENT(OUT) :: ier ! Local variables CHARACTER (LEN=50) :: str INTEGER :: length, pos, decpt, status, power10, i, i1, i2 REAL (dp) :: sgn, temp CHARACTER (LEN= 5) :: expnt str = ADJUSTL(string) length = LEN_TRIM( ADJUSTL(string) ) IF (length > 50) THEN ! Truncate string to 50 characters preserving ! the exponent, if present. pos = SCAN(string(length-4:length), 'DdEe') IF (pos > 0) THEN str(45+pos:50) = string(length-5+pos:length) END IF ! Check that the string contains a '.' IF (INDEX(str, '.') == 0) THEN ier = 1 value = quad(zero, zero) RETURN END IF END IF ier = 0 ! Determine sign of mantissa sgn = +1._dp IF (str(1:1) == '-') THEN sgn = -1._dp str = str(2:) ELSE IF (str(1:1) == '+') THEN str = str(2:) END IF ! Separate the exponent, if there is one length = LEN_TRIM(str) IF (length > 4) THEN pos = SCAN(str(length-4:length), 'DdEe') IF (pos > 0) THEN expnt = str(length-5+pos:) str(length-5+pos:) = ' ' ELSE expnt = ' ' END IF ELSE expnt = ' ' END IF ! decpt = position of decimal point decpt = INDEX(str, '.') IF (decpt > 0) str = str(:decpt-1) // str(decpt+1:) ! Read str in blocks of up to 9 digits at a time, as a large integer. i1 = 1 length = LEN_TRIM(str) value = quad(zero, zero) DO i2 = MIN(i1+8, length) READ(str(i1:i2), '(f9.0)', IOSTAT=status) temp IF (status /= 0) THEN ier = 2 RETURN END IF IF (i1 > 1) THEN value = value * (10._dp ** (i2+1-i1)) + temp ELSE value = quad(temp, zero) END IF i1 = i2 + 1 IF (i1 > length) EXIT END DO IF (sgn < zero) value = -value ! Multiply by appropriate power of 10 for position of the decimal point, ! and the exponent. IF (expnt == ' ') THEN power10 = 0 ELSE i = LEN_TRIM(expnt) READ(expnt(2:i), '(i4)') power10 END IF IF (decpt > 0) power10 = power10 - (length + 1 - decpt) ! As 1.D+15 is represented exactly in double-precision IEEE arithmetic, ! multiply by multiples of 1.D+15 or 1.D-15. IF (power10 == 0) RETURN IF (power10 < 0) THEN DO i = 1, -power10/15 value = value / 1.D+15 END DO i = MOD(-power10, 15) IF (i /= 0) value = value / (10._dp ** i) ELSE DO i = 1, power10/15 value = value * 1.D+15 END DO i = MOD(power10, 15) IF (i /= 0) value = value * (10._dp ** i) END IF RETURN END SUBROUTINE string_quad SUBROUTINE quad_string(value, string, ier) ! Convert a quadruple-precision quantity to a decimal character string. ! Error indicator ier = 0 if conversion OK ! = 1 if the length of the string < 36 characters. TYPE (quad), INTENT(IN) :: value CHARACTER (LEN=*), INTENT(OUT) :: string INTEGER, INTENT(OUT) :: ier ! Local variables CHARACTER (LEN= 1) :: sgn CHARACTER (LEN=17) :: str1, str2 TYPE (quad) :: val INTEGER :: dec_expnt, i REAL (dp) :: tmp IF (LEN(string) < 36) THEN ier = 1 string = '***' RETURN END IF ier = 0 ! Check if value = zero. IF (value%hi == zero) THEN string = ' 0.00' RETURN END IF IF (value%hi < zero) THEN sgn = '-' val = - value ELSE sgn = ' ' val = value END IF ! Use LOG10 to set the exponent. dec_expnt = FLOOR( LOG10(val%hi) ) ! Get the first 15 decimal digits IF (dec_expnt /= 14) THEN val = val * EXP( LOG(quad(10._dp, zero)) * REAL(14 - dec_expnt) ) END IF WRITE(str1, '(f16.0)') val%hi ! Calculate the remainder READ(str1, '(f16.0)') tmp val = val - tmp ! If val is -ve, subtract 1 from the last digit of str1, and add 1 to val. IF (val%hi < -0.5D-16) THEN tmp = tmp - one WRITE(str1, '(f16.0)') tmp val = val + one END IF val = val * 1.D15 ! write the second 15 digits WRITE(str2, '(f16.0)') val%hi ! If str2 consists of asterisks, add 1 in the last digit to str1. ! Set str2 to zeroes. IF (str2(2:2) == '*') THEN tmp = tmp + one WRITE(str1, '(f17.0)') tmp IF (str1(1:1) /= ' ') THEN dec_expnt = dec_expnt + 1 ELSE str1 = str1(2:17) END IF str2 = '000000000000000.' END IF ! Replace leading blanks with zeroes DO i = 1, 15 IF (str2(i:i) /= ' ') EXIT str2(i:i) = '0' END DO ! Combine str1 & str2, removing decimal points & adding exponent. i = INDEX(str1, '.') str1(i:i) = ' ' str2(16:16) = ' ' string = '.' // TRIM(ADJUSTL(str1)) // TRIM(ADJUSTL(str2)) // 'E' WRITE(str1, '(i4.2)') dec_expnt+1 string = TRIM(string) // ADJUSTL(str1) ! Restore the sign. IF (sgn == '-') THEN string = '-' // ADJUSTL(string) ELSE string = ADJUSTL(string) END IF RETURN END SUBROUTINE quad_string FUNCTION q_epsilon(x) RESULT(eps) ! Returns the machine accuracy. ! eps is the smallest value such that x.(1 + eps) > x for all x. ! This value is machine dependent. It has been checked only for the ! the following PC compilers: ! Lahey's ELF90 ! Lahey/Fujitsi LF95 ! Compaq (Digital) DVF5 TYPE (quad), INTENT(IN) :: x TYPE (quad) :: eps eps = quad(0.6163e-32_dp, 0.0_dp) RETURN END FUNCTION q_epsilon END MODULE quadruple_precision