Extended Precision Arithmetic

This section describes a set of functions for mixed precision arithmetic. The functions 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 functions 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