Extended Precision Arithmetic
This section describes a set of routines for mixed precision arithmetic. The routines are designed to allow the computation and use of the full quadruple precision result from the multiplication of two double precision numbers. An array called the accumulator stores the result of this multiplication. The result of the multiplication is added to the current contents of the accumulator. It is also possible to add a double precision number to the accumulator or to store a double precision approximation in the accumulator.
The mixed double precision arithmetic routines are described below. The accumulator array, QACC, is a double precision array of length 2. Double precision variables are denoted by DA and DB. Available operations are:
Initialize a real accumulator, QACC ← DA.
CALL DQINI (DA, QACC)
Store a real accumulator, DA ← QACC.
CALL DQSTO (QACC, DA)
Add to a real accumulator, QACC ← QACC + DA.
CALL DQADD (DA, QACC)
Add a product to a real accumulator, QACC ← QACC + DA*DB.
CALL DQMUL (DA, DB, QACC)
There are also mixed double complex arithmetic versions of the above routines. The accumulator, ZACC, is a double precision array of length 4. Double complex variables are denoted by ZA and ZB. Available operations are:
Initialize a complex accumulator, ZACC ← ZA.
CALL ZQINI (ZA, ZACC)
Store a complex accumulator, ZA ← ZACC.
CALL ZQSTO (ZACC, ZA)
Add to a complex accumulator, ZACC ← ZACC + ZA.
CALL ZQADD (ZA, ZACC)
Add a product to a complex accumulator, ZACC ← ZACC + ZA * ZB.
CALL ZQMUL (ZA, ZB, ZACC)
Example
In this example, the value of 1.0D0/3.0D0 is computed in quadruple precision using Newton’s method. Four iterations of
with a = 3 are taken. The error ax ‑ 1 is then computed. The results are accurate to approximately twice the usual double precision accuracy, as given by the IMSL routine DMACH(4), in the Reference Material section of this manual;. Since DMACH is machine dependent, the actual accuracy obtained is also machine dependent.
USE IMSL_LIBRARIES
IMPLICIT NONE
INTEGER I, NOUT
DOUBLE PRECISION A, DACC(2), DMACH, ERROR, SACC(2), X(2), X1, X2, EPSQ
!
CALL UMACH (2, NOUT)
A = 3.0D0
CALL DQINI (1.0001D0/A, X)
! Compute X(K+1) = X(K) - A*X(K)*X(K)
! + X(K)
DO 10 I=1, 4
X1 = X(1)
X2 = X(2)
! Compute X + X
CALL DQADD (X1, X)
CALL DQADD (X2, X)
! Compute X*X
CALL DQINI (0.0D0, DACC)
CALL DQMUL (X1, X1, DACC)
CALL DQMUL (X1, X2, DACC)
CALL DQMUL (X1, X2, DACC)
CALL DQMUL (X2, X2, DACC)
! Compute -A*(X*X)
CALL DQINI (0.0D0, SACC)
CALL DQMUL (-A, DACC(1), SACC)
CALL DQMUL (-A, DACC(2), SACC)
! Compute -A*(X*X) + (X + X)
CALL DQADD (SACC(1), X)
CALL DQADD (SACC(2), X)
10 CONTINUE
! Compute A*X - 1
CALL DQINI (0.0D0, SACC)
CALL DQMUL (A, X(1), SACC)
CALL DQMUL (A, X(2), SACC)
CALL DQADD (-1.0D0, SACC)
CALL DQSTO (SACC, ERROR)
! ERROR should be less than MACHEPS**2
EPSQ = AMACH(4)
EPSQ = EPSQ * EPSQ
WRITE (NOUT,99999) ERROR, ERROR/EPSQ
!
99999 FORMAT (' A*X - 1 = ', D15.7, ' = ', F10.5, '*MACHEPS**2')
END
Output
A*X - 1 = 0.6162976D-32 = 0.12500*MACHEPS**2
