module QUADRUPLE_PRECISION
!
! Quadruple precision package for the F90 subset - F
! N.B. This has only been tested using the Windows 95 version of F.
!      On other computers you will probably need to change the value of
!      CONSTANT below.
!
! N.B. This is quadruple precision implemented in SOFTWARE, hence it is SLOW.

! N.B. Overlaying of some operators and intrinsic functions is supported.
!      For instance if I is an integer, and Q is quad.-precision, then
!      Q * I is allowed in a program, but at the moment Q / I is not.
!      If I, R, D, Q represent respectively integer, real, double precision
!      and quadruple, then there are 7 combinations for multiplication:
!      Q*I, Q*R, Q*D, Q*Q, I*Q, R*Q, D*Q.   When you consider the number of
!      operators and intrinsic functions, this requires a very large number
!      of functions to support them all.   Check with the interfaces below
!      to see if the combination you want is supported.   If it is not, then
!      you can:
!      1. Write your own function to do it and add it as a MODULE PROCEDURE,
!      2. Reverse the order, e.g. Q+I instead of I+Q, or
!      3. Use Q / quad( real(I,DP), 0.0_dp ) instead of Q/I, i.e. convert the
!         integer to quad.-precision.

! 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 >.
!  2 Feb 98. Converted to F-compatibility
! 16 Feb 98. Added complex operations
! 5 Oct 99.  log10e corrected.
! 15 Oct 99. Corrected 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.
! 18 Sep 02. Added reference to Linnainmaa in comments at beginning.

implicit none

public :: LONGMUL, MULT_QUAD_INT, MULT_INT_QUAD, MULT_QUAD_DP, MULT_DP_QUAD,  &
          MULT_QUAD_REAL, MULT_REAL_QUAD, LONGDIV, DIV_QUAD_DP, DIV_QUAD_REAL, &
          DIV_QUAD_INT, DIV_INT_QUAD, DIV_REAL_QUAD,  DIV_DP_QUAD,  &
          LONGADD, QUAD_ADD_INT, QUAD_ADD_REAL, QUAD_ADD_DP,  &
          INT_ADD_QUAD, REAL_ADD_QUAD, DP_ADD_QUAD, LONGSUB, QUAD_SUB_INT,  &
          QUAD_SUB_REAL, QUAD_SUB_DP, INT_SUB_QUAD, REAL_SUB_QUAD,  &
          DP_SUB_QUAD, NEGATE_QUAD, QUAD_EQ_INT, QUAD_EQ_REAL, QUAD_EQ_DP,  &
          INT_EQ_QUAD, REAL_EQ_QUAD, DP_EQ_QUAD, QUAD_POW_INT, QUAD_POW_REAL, &
          QUAD_POW_DP, QUAD_POW_QUAD, QUAD_LT, QUAD_LE, QUAD_EQ, QUAD_GE,     &
          QUAD_GT, QSCALE, QABS, LONGSQRT, LONGLOG, LONGEXP, LONGSIN, LONGCOS, &
          LONGTAN, LONGASIN, LONGACOS, LONGATAN, QATAN2, EXACTMUL2, LONGMODR, &
          QUAD_SUM, QUAD_DOT_PRODUCT, Q_MATMUL12, Q_MATMUL21, Q_MATMUL22,  &
          STRING_QUAD, LONGCST, QUAD_STRING, Q_EPSILON
private :: QC_ADD, QC_SUB, NEGATE_QC, QC_MULT, QC_DIV, QC_ABS, QC_CMPLX, &
           QC_CONJG, QC_AIMAG, QC_SQRT, QC_EXP, QC_LOG

public :: scale, abs, sqrt, log, exp, sin, cos, tan, asin, acos, atan, atan2, &
          sum, dot_product, matmul, epsilon
public :: cmplx, conjg, aimag

public :: operator (*), operator (/), operator (+), operator (-),    &
          operator (**), operator (<), operator (<=), operator (==), &
          operator (>=), operator (>), assignment (=)

intrinsic scale, abs, sqrt, log, exp, sin, cos, tan, asin, acos, atan, atan2, &
          sum, dot_product, matmul, epsilon
intrinsic cmplx, conjg, aimag

integer, parameter, public    :: DP = selected_real_kind(12, 100)
integer, parameter, private   :: NBITS = digits(1.0_DP)
real (kind=DP), parameter, private :: CONSTANT = 2.0_DP**(NBITS+11-NBITS/2) + 1.0_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.
! For other computers, try:
! REAL (kind=DP), parameter, private :: CONSTANT = 2._DP**(NBITS-NBITS/2) + 1._dp


real (kind=DP), parameter, private :: ZERO = 0.0_DP, ONE = 1.0_DP

type, public :: QUAD
  real (kind=DP) :: HI, LO
end type QUAD

type, public  :: QC
  type (QUAD) :: QR, QI                ! quad_real, quad_imaginary
end type QC

type, public :: FUNC_RESULT           ! Structure for the output from LONGMODR
  type (QUAD)  :: remainder
  integer      :: nmult
end type FUNC_RESULT

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
  module procedure QC_MULT
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
  module procedure QC_DIV
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
  module procedure QC_ADD
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
  module procedure QC_SUB
  module procedure NEGATE_QC
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
  module procedure QC_ABS
end interface

interface cmplx
  module procedure QC_CMPLX
end interface

interface conjg
  module procedure QC_CONJG
end interface

interface aimag
  module procedure QC_AIMAG
end interface

interface sqrt
  module procedure LONGSQRT
  module procedure QC_SQRT
end interface

interface log
  module procedure LONGLOG
  module procedure QC_LOG
end interface

interface exp
  module procedure LONGEXP
  module procedure QC_EXP
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, public ::  &
          PI = QUAD( 0.3141592653589793e+01_DP,  0.1224646799147353e-15_DP ), &
       PIBY2 = QUAD( 0.1570796326794897e+01_DP, -0.3828568698926950e-15_DP ), &
       PIBY3 = QUAD( 0.1047197551196598e+01_DP, -0.3292527815701405e-15_DP ), &
       PIBY4 = QUAD( 0.7853981633974484e+00_DP, -0.8040613248383182e-16_DP ), &
       PIBY6 = QUAD( 0.5235987755982990e+00_DP, -0.1646263907850702e-15_DP ), &
       TWOPI = QUAD( 0.6283185307179586e+01_DP,  0.2449293598294707e-15_DP ), &
       LN_PI = QUAD( 0.1144729885849400e+01_DP,  0.2323105560877391e-15_DP ), &
      SQRTPI = QUAD( 0.1772453850905516e+01_DP, -0.7666586499825800e-16_DP ), &
    FACT_PT5 = QUAD( 0.8862269254527582e+00_DP, -0.1493552349616447e-15_DP ), &
     SQRT2PI = QUAD( 0.2506628274631001e+01_DP, -0.6273750096546544e-15_DP ), &
   LNSQRT2PI = QUAD( 0.9189385332046728e+00_DP, -0.3878294158067242e-16_DP ), &
   ONE_ON2PI = QUAD( 0.1591549430918953e+00_DP,  0.4567181289366658e-16_DP ), &
 TWO_ON_RTPI = QUAD( 0.1128379167095513e+01_DP, -0.4287537502368968e-15_DP ), &
     DEG2RAD = QUAD( 0.1745329251994330e-01_DP, -0.3174581724866598e-17_DP ), &
     RAD2DEG = QUAD( 0.5729577951308232e+02_DP, -0.1987849567057628e-14_DP ), &
         LN2 = QUAD( 0.6931471805599454e+00_DP, -0.8783183432405266e-16_DP ), &
        LN10 = QUAD( 0.2302585092994046e+01_DP, -0.2170756223382249e-15_DP ), &
       LOG2E = QUAD( 0.1442695040888964e+01_DP, -0.6457785410341630e-15_DP ), &
      LOG10E = QUAD( 0.4342944819032519e+00_DP, -0.5552037773430574e-16_DP ), &
     LOG2_10 = QUAD( 0.3321928094887362e+01_DP,  0.1661617516973592e-15_DP ), &
     LOG10_2 = QUAD( 0.3010299956639812e+00_DP, -0.2803728127785171e-17_DP ), &
       EULER = QUAD( 0.5772156649015330e+00_DP, -0.1159652176149463e-15_DP ), &
           E = QUAD( 0.2718281828459045e+01_DP,  0.1445646891729250e-15_DP ), &
       SQRT2 = QUAD( 0.1414213562373095e+01_DP,  0.1253716717905022e-15_DP ), &
       SQRT3 = QUAD( 0.1732050807568877e+01_DP,  0.3223954471431004e-15_DP ), &
      SQRT10 = QUAD( 0.3162277660168380e+01_DP, -0.6348773795572286e-15_DP )

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 (kind=DP), intent(in) :: A, C
type (QUAD)                :: AC

!     Local variables
real (kind=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 (kind=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, real(I,DP)) + EXACTMUL2(A%LO, real(I,DP))
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, real(I,DP)) + EXACTMUL2(A%LO, real(I,DP))
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 (kind=DP), intent(in) :: B
type (QUAD)                :: C

!     Local variables
type (QUAD)             :: Z
real (kind=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 (kind=DP)          :: ZZ

Z = EXACTMUL2(A%HI, real(B,DP))
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 (kind=DP), intent(in) :: B
type (QUAD)                :: C

!     Local variables
type (QUAD)             :: Z
real (kind=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 (kind=DP)          :: ZZ

Z = EXACTMUL2(B%HI, real(A,DP))
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 (kind=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 (kind=DP), intent(in) :: B
type (QUAD)                :: C

!     Local variables
real (kind=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 (kind=DP) :: Z, ZZ
type (QUAD)    :: Q

Z = A%HI / B
Q = EXACTMUL2(real(B,kind=DP), 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 number (b)

type (QUAD), intent(in)  :: A
integer, intent(in)      :: B
type (QUAD)              :: C

!     Local variables
real (kind=DP) :: Z, ZZ
type (QUAD)    :: Q

Z = A%HI / B
Q = EXACTMUL2(real(B,kind=DP), 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, B) result(C)
!     Divide a quadruple-precision number (a) by an integer number (b)

integer, intent(in)      :: A
type (QUAD), intent(in)  :: B
type (QUAD)              :: C

!     Local variables
real (kind=DP) :: Z, ZZ
type (QUAD)    :: Q

Z = real(A,kind=DP) / B%HI
Q = EXACTMUL2(B%HI, Z)
ZZ = (((A - Q%HI) - Q%LO) - Z*B%LO) / (B%HI + B%LO)
C%HI = Z + ZZ
C%LO = (Z - C%HI) + ZZ

return
end function DIV_INT_QUAD



function DIV_DP_QUAD(A, B) result(C)
!     Divide a quadruple-precision number (a) by an integer number (b)

real (kind=DP), intent(in)  :: A
type (QUAD), intent(in)     :: B
type (QUAD)                 :: C

!     Local variables
real (kind=DP) :: Z, ZZ
type (QUAD)    :: Q

Z = A / B%HI
Q = EXACTMUL2(B%HI, Z)
ZZ = (((A - Q%HI) - Q%LO) - Z*B%LO) / (B%HI + B%LO)
C%HI = Z + ZZ
C%LO = (Z - C%HI) + ZZ

return
end function DIV_DP_QUAD



function DIV_REAL_QUAD(A, B) result(C)
!     Divide a quadruple-precision number (a) by an integer number (b)

real, intent(in)         :: A
type (QUAD), intent(in)  :: B
type (QUAD)              :: C

!     Local variables
real (kind=DP) :: Z, ZZ
type (QUAD)    :: Q

Z = real(A,kind=DP) / B%HI
Q = EXACTMUL2(B%HI, Z)
ZZ = (((A - Q%HI) - Q%LO) - Z*B%LO) / (B%HI + B%LO)
C%HI = Z + ZZ
C%LO = (Z - C%HI) + ZZ

return
end function DIV_REAL_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 (kind=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 (kind=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 (kind=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 (kind=DP), intent(in) :: C
type (QUAD)                :: AC

!     Local variables
real (kind=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(A, C) result(AC)

integer, intent(in)      :: A
type (QUAD), intent(in)  :: C
type (QUAD)              :: AC

!     Local variables
real (kind=DP) :: Z, Q, ZZ

Z = C%HI + A
Q = C%HI - Z
ZZ = (((Q + A) + (C%HI - (Q + Z))) + C%LO)
AC%HI = Z + ZZ
AC%LO = (Z - AC%HI) + ZZ

return
end function INT_ADD_QUAD



function REAL_ADD_QUAD(A, C) result(AC)

real, intent(in)         :: A
type (QUAD), intent(in)  :: C
type (QUAD)              :: AC

!     Local variables
real (kind=DP) :: Z, Q, ZZ

Z = C%HI + A
Q = C%HI - Z
ZZ = (((Q + A) + (C%HI - (Q + Z))) + C%LO)
AC%HI = Z + ZZ
AC%LO = (Z - AC%HI) + ZZ

return
end function REAL_ADD_QUAD



function DP_ADD_QUAD(A, C) result(AC)

real (kind=DP), intent(in)  :: A
type (QUAD), intent(in)     :: C
type (QUAD)                 :: AC

!     Local variables
real (kind=DP) :: Z, Q, ZZ

Z = C%HI + A
Q = C%HI - Z
ZZ = (((Q + A) + (C%HI - (Q + Z))) + C%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 (kind=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 (kind=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 (kind=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 (kind=DP), intent(in) :: C
type (QUAD)                :: AC

!     Local variables
real (kind=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 (kind=DP) :: Z, Q, ZZ

Z = real(A,kind=DP) - C%HI
Q = real(A,kind=DP) - Z
ZZ = (((Q - C%HI) + (real(A,kind=DP) - (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 (kind=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 (kind=DP), intent(in)  :: A
type (QUAD), intent(in)     :: C
type (QUAD)                 :: AC

!     Local variables
real (kind=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 (kind=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 (kind=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(unit=*, fmt=*) "** Exponential overflow in routine QUAD_POW_INT **"
            return
          end if
          B = B * POWER
        end if
        IA = ibclr(IA, J)
        if (IA == 0) then
          exit
        end if
      end if
      if (exponent(POWER%HI) > EMAX/2) then
        write(unit=*, fmt=*) "** Exponential overflow in routine QUAD_POW_INT **"
        return
      end if
      POWER = POWER * POWER
    end do
    if (I < 0) then
      B = QUAD(ONE, ZERO) / B
    end if
  case default
    if (A%HI < ZERO) then
      if (I * log(-A%HI) > log( huge(ONE) )) then
        write(unit=*, fmt=*) "** 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(unit=*, fmt=*) "** 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(unit=*, fmt=*)   &
  " *** 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 (kind=DP), intent(in) :: D
type (QUAD)                :: B

if (A%HI < ZERO) then
  write(unit=*, fmt=*)   &
  " *** 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(unit=*, fmt=*)   &
  " *** 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 (kind=DP) :: T, RES
type (QUAD)    :: TT

! Check that ahi >= 0.

if (A%HI < 0.0_DP) then
  write(unit=*, fmt=*) " *** Negative argument for longsqrt ***"
  return
else if (A%HI == ZERO) then
  B%HI = ZERO
  B%LO = ZERO
  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 = ZERO
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 - real(IY,DP)) * 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.0)*YSQ + 2162160.0_DP)*YSQ + 302702400.0_DP)*YSQ +   &
               8821612800.0_DP)
SUM2 = (((90.0*YSQ + 110880.0)*YSQ + 30270240.0_DP)*YSQ + 2075673600.0_DP)*YSQ +  &
       17643225600.0_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 + ONE)
Y = scale(Y, 2)
Y = Y + ONE
Y = scale(Y, IY)

return
end function LONGEXP



function LONGMODR(A, B) result(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
type (FUNC_RESULT)       :: REM

! Local variables

type (QUAD) :: TEMP

! Check that b%hi .ne. 0

if (B%HI == ZERO) then
  write(unit=*, fmt=*) " *** Error in longmodr - 3rd argument zero ***"
  return
end if

! Calculate NMULT.

TEMP = A / B
REM % NMULT = nint(real(TEMP % HI))

! Calculate remainder preserving full accuracy.

TEMP = EXACTMUL2(real(REM % NMULT,DP), B % HI)
REM % REMAINDER % HI = A % HI
REM % REMAINDER % LO = ZERO
TEMP = REM % REMAINDER - TEMP
REM % REMAINDER % HI = A % LO
REM % REMAINDER % LO = ZERO
TEMP = REM % REMAINDER + TEMP
REM % REMAINDER = EXACTMUL2(real(REM % NMULT,DP), B % LO)
REM % REMAINDER = TEMP - REM % REMAINDER

return
end function 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.
B = LONGCST(A, 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.
B = LONGCST(A, 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.
B = LONGCST(A, SINE, COSINE, TANGENT)

return
end function LONGTAN



function LONGCST(A, SINE, COSINE, TANGENT) result(B)

type (QUAD), intent(in)  :: A
logical, intent(in)      :: SINE, COSINE, TANGENT
type (QUAD)              :: B

! Local variables

logical                   :: POS
type (QUAD)               :: D, TERM, TEMP, ANGLE, PIBY40, SUM1, SUM2, SSIN
integer                   :: NPI, IPT, I
type (FUNC_RESULT)        :: REM
real (kind=DP), parameter :: TOL15 = 1.0e-15_DP, TOL30 = 1.0e-30_DP

! sin(i.pi/40), i = 0(1)20
real (kind=DP), parameter, dimension(2, 0:20) :: TABLE = reshape( (/   &
               0.0000000000000000e+00_DP,  0.0000000000000000e+00_DP,  &
               0.7845909572784494e-01_DP,  0.1464397249532491e-17_DP,  &
               0.1564344650402309e+00_DP, -0.2770509565052586e-16_DP,  &
               0.2334453638559054e+00_DP,  0.2058612230858154e-16_DP,  &
               0.3090169943749475e+00_DP, -0.8267172724967036e-16_DP,  &
               0.3826834323650898e+00_DP, -0.1005077269646159e-16_DP,  &
               0.4539904997395468e+00_DP, -0.1292033036231312e-16_DP,  &
               0.5224985647159488e+00_DP,  0.6606794454708078e-16_DP,  &
               0.5877852522924732e+00_DP, -0.1189570533007057e-15_DP,  &
               0.6494480483301838e+00_DP, -0.1134903961116171e-15_DP,  &
               0.7071067811865476e+00_DP, -0.4833646656726458e-16_DP,  &
               0.7604059656000310e+00_DP, -0.1036987135483477e-15_DP,  &
               0.8090169943749476e+00_DP, -0.1381828784809282e-15_DP,  &
               0.8526401643540922e+00_DP,  0.4331886637554353e-16_DP,  &
               0.8910065241883680e+00_DP, -0.1474714419679880e-15_DP,  &
               0.9238795325112868e+00_DP, -0.9337725537817898e-16_DP,  &
               0.9510565162951536e+00_DP, -0.7008780156242836e-16_DP,  &
               0.9723699203976766e+00_DP,  0.4478912629332321e-16_DP,  &
               0.9876883405951378e+00_DP, -0.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 = -0.1081617080994607e-16_DP

! Reduce angle to range (-pi/2, +pi/2) by subtracting an integer multiple of pi.

REM = LONGMODR(A, PI)
ANGLE = REM % REMAINDER
NPI = REM % NMULT

! Find nearest multiple of pi/40 to angle.

REM = LONGMODR(ANGLE, PIBY40)
D = REM % REMAINDER
IPT = REM % NMULT

! 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
do
  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 / real(I,DP)
    end if
    if (POS) then
      SUM1 = SUM1 + TERM
    else
      SUM1 = SUM1 - TERM
    end if
  else
    TERM%HI = TERM%HI * D%HI / real(I,DP)             ! 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 / real(I,DP)                      ! Use quad. precision
    if (POS) then
      SUM2 = SUM2 + TERM
    else
      SUM2 = SUM2 - TERM
    end if
  else
    TERM%HI = TERM%HI * D%HI / real(I,DP)             ! 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) then
    exit
  end if
end do

SUM1 = SUM1 + 1.0_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
    SSIN = 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 == ZERO) then
    write(unit=*, fmt=*) " *** Infinite tangent - routine longcst ***"
    B%HI = huge(ONE)
    B%LO = ZERO
    return
  end if
  B = SSIN / B
end if

return
end function LONGCST



function LONGASIN(A) result(B)

! Quadratic-precision arc sine (about 31 decimals).
! Newton-Raphson iteration to solve:  f(b) = sin(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

! Check that -1 <= a%hi <= +1.

if (A%HI < -1.0_DP .or. A%HI > ONE) then
  write(unit=*, fmt=*) " *** Argument outside range for longasin ***"
  return
end if

! 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)

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).

type (QUAD), intent(in) :: A
type (QUAD)             :: B

! Local variables
type (QUAD)  :: Y

! Check that -1 <= a%hi <= +1.

if (A%HI < -1.0_DP .or. A%HI > ONE) then
  write(unit=*, fmt=*) " *** Argument outside range for longacos ***"
  return
end if

! First approximation is  y = acos(a).
! Quadruple-precision result is  y + [cos(y) - a]/sin(y).

Y%HI = acos(A%HI)
Y%LO = ZERO
B = Y + (cos(Y) - A) / sin(Y%HI)

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 = ZERO
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 = ZERO
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(unit=*, fmt=*) " ** Error invoking DOT_PRODUCT - different argument sizes **"
  write(unit=*, fmt="(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(unit=*, fmt=*) " ** 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(unit=*, fmt=*) " ** 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(unit=*, fmt=*) " ** 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 (kind=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 = ONE
if (STR(1:1) == "-") then
  SGN = -1.0_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)
POS = scan(STR(LENGTH-4:LENGTH), "DdEe")
if (POS > 0) then
  EXPNT = STR(LENGTH-5+POS:)
  STR(LENGTH-5+POS:) = "     "
else
  EXPNT = " "
end if

! decpt = position of decimal point

DECPT = index(STR, ".")
if (DECPT > 0) then
  STR = STR(:DECPT-1) // STR(DECPT+1:)
end if

! 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(unit=STR(I1:I2), fmt="(f9.0)", iostat=STATUS) TEMP
  if (STATUS /= 0) then
    IER = 2
    return
  end if
  if (I1 > 1) then
    VALUE = VALUE * (10.0_DP ** (I2+1-I1)) + TEMP
  else
    VALUE = QUAD(TEMP, ZERO)
  end if
  I1 = I2 + 1
  if (I1 > LENGTH) then
    exit
  end if
end do

if (SGN < ZERO) then
  VALUE = -VALUE
end if

! 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(unit=EXPNT(2:I), fmt="(i4)") POWER10
end if
if (DECPT > 0) then
  POWER10 = POWER10 - (LENGTH + 1 - DECPT)
end if

! 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) then
  return
end if
if (POWER10 < 0) then
  do I = 1, -POWER10/15
    VALUE = VALUE / 1.0e+15_DP
  end do
  I = modulo(-POWER10, 15)
  if (I /= 0) then
    VALUE = VALUE / (10.0_DP ** I)
  end if
else
  do I = 1, POWER10/15
    VALUE = VALUE * 1.0e+15_DP
  end do
  I = modulo(POWER10, 15)
  if (I /= 0) then
    VALUE = VALUE * (10.0_DP ** I)
  end if
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 (kind=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.0_DP, ZERO)) * real(14 - DEC_EXPNT) )
end if
write(unit=STR1, fmt="(f16.0)") VAL%HI

! Calculate the remainder
read(unit=STR1, fmt="(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.5e-16) then
  TMP = TMP - ONE
  write(unit=STR1, fmt="(f16.0)") TMP
  VAL = VAL + ONE
end if
VAL = VAL * 1.0e15_DP

! write the second 15 digits
write(unit=STR2, fmt="(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(unit=STR1, fmt="(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) /= " ") then
    exit
  end if
  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(unit=STR1, fmt="(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.

type (QUAD), intent(in)  :: X
type (QUAD)              :: EPS

EPS = QUAD(0.6163e-32_DP, 0.0_DP)

return
end function Q_EPSILON




! Quad. prec. complex functions & routines

function QC_ADD(X, Y) result(Z)
! Quadruple-precision complex addition

type (QC), intent(in) :: X, Y
type (QC)             :: Z

Z % QR = X % QR + Y % QR
Z % QI = X % QI + Y % QI

return
end function QC_ADD



function QC_SUB(X, Y) result(Z)
! Quadruple-precision complex subtraction

type (QC), intent(in) :: X, Y
type (QC)             :: Z

Z%QR = X%QR - Y%QR
Z%QI = X%QI - Y%QI

return
end function QC_SUB



function NEGATE_QC(X) result(Z)
! Quadruple-precision complex negate sign

type (QC), intent(in) :: X
type (QC)             :: Z

Z%QR = - X%QR
Z%QI = - X%QI

return
end function NEGATE_QC



function QC_MULT(X, Y) result(Z)
! Quadruple-precision complex multiplication

type (QC), intent(in) :: X, Y
type (QC)             :: Z

Z%QR = X%QR * Y%QR - X%QI * Y%QI
Z%QI = X%QI * Y%QR + X%QR * Y%QI

return
end function QC_MULT



function QC_DIV(X, Y) result(Z)
! Quadruple-precision complex division

type (QC), intent(in) :: X, Y
type (QC)             :: Z

type (QUAD) :: T, D

if ( abs(Y%QR%HI) > abs(Y%QI%HI) ) then
  T = Y%QI / Y%QR
  D = Y%QR + T * Y%QI
  Z%QR = (X%QR + T * X%QI) / D
  Z%QI = (X%QI - T * X%QR) / D
else
  T = Y%QR / Y%QI
  D = Y%QI + T * Y%QR
  Z%QR = (T * X%QR + X%QI) / D
  Z%QI = (T * X%QI - X%QR) / D
end if

return
end function QC_DIV



function QC_CMPLX(XR, XI) result(Z)
! Convert pair of quad.precision numbers to qc.

type (QUAD), intent(in) :: XR, XI
type (QC)               :: Z

Z%QR = XR
Z%QI = XI

return
end function QC_CMPLX



function QC_AIMAG(X) result(Z)
! Convert pair of quad.precision numbers to qc.

type (QC), intent(in)   :: X
type (QUAD)             :: Z

Z = X%QI

return
end function QC_AIMAG



function QC_CONJG(X) result(Z)
! Quadruple-precision complex square root

type (QC), intent(in) :: X
type (QC)             :: Z

Z%QR = X%QR
Z%QI = - X%QI

return
end function QC_CONJG



function QC_SQRT(X) result(Z)
! Quadruple-precision complex square root

type (QC), intent(in) :: X
type (QC)             :: Z

! Local variables
type (QUAD)                :: R
real (kind=DP), parameter  :: VSMALL = tiny(1.0_DP)

R = sqrt(X%QR **2 + X%QI **2)
if (R%HI < VSMALL) then
  Z%QR = QUAD( 0.0_DP, 0.0_DP )
  Z%QI = QUAD( 0.0_DP, 0.0_DP )
  return
end if

Z%QR = sqrt( scale( R + abs(X%QR), -1 ) )
Z%QI = abs( scale( X%QI, -1 ) ) / Z%QR
if (X%QR%HI < 0.0_DP) then
  R = Z%QR
  Z%QR = Z%QI
  Z%QI = R
end if
if (X%QI%HI < 0.0_DP) then
  Z%QI = - Z%QI
end if

return
end function QC_SQRT



function QC_ABS(X) result(Z)
! Quadruple-precision complex absolute (modulus)
! SQRT( x%qr **2 + x%qi **2 ) avoiding overflow
! Using either:   X^2 + Y^2 = X^2.(1 + (Y/X)^2)
!           or:             = Y^2.(1 + (X/Y)^2)

type (QC), intent(in) :: X
type (QUAD)           :: Z

type (QUAD) :: A

if ( abs(X%QR%HI) > abs(X%QI%HI) ) then
  A = X%QI / X%QR
  Z = abs( X%QR ) * sqrt( QUAD(1.0_DP, 0.0_DP) + A**2 )
else if ( X%QI%HI == 0.0_DP ) then
  Z = QUAD( 0.0_DP, 0.0_DP )
else
  A = X%QR / X%QI
  Z = abs( X%QI ) * sqrt( QUAD(1.0_DP, 0.0_DP) + A**2 )
end if

return
end function QC_ABS



function QC_EXP(X) result(Z)
! Quadruple-precision complex exponential

type (QC), intent(in) :: X
type (QC)             :: Z

type (QUAD)           :: EXPX

EXPX = exp(X%QR)
Z%QR = EXPX * cos(X%QI)
Z%QI = EXPX * sin(X%QI)

return
end function QC_EXP



function QC_LOG(X) result(Z)
! Quadruple-precision complex logarithm

type (QC), intent(in) :: X
type (QC)             :: Z

Z%QR = log( abs(X) )
Z%QI = atan2( X%QI, X%QR )

return
end function QC_LOG

end module QUADRUPLE_PRECISION