| 1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798 |
- *DECK MPMUL
- SUBROUTINE MPMUL (X, Y, Z)
- C***BEGIN PROLOGUE MPMUL
- C***SUBSIDIARY
- C***PURPOSE Subsidiary to DQDOTA and DQDOTI
- C***LIBRARY SLATEC
- C***TYPE ALL (MPMUL-A)
- C***AUTHOR (UNKNOWN)
- C***DESCRIPTION
- C
- C Multiplies X and Y, returning result in Z, for 'mp' X, Y and Z.
- C The simple o(t**2) algorithm is used, with four guard digits and
- C R*-rounding. Advantage is taken of zero digits in X, but not in Y.
- C Asymptotically faster algorithms are known (see Knuth, VOL. 2),
- C but are difficult to implement in FORTRAN in an efficient and
- C machine-independent manner. In comments to other 'mp' routines,
- C M(t) is the time to perform t-digit 'mp' multiplication. Thus
- C M(t) = o(t**2) with the present version of MPMUL, but
- C M(t) = o(t.log(t).log(log(t))) is theoretically possible.
- C
- C The arguments X(*), Y(*), and Z(*), and the variable R in COMMON are
- C all INTEGER arrays of size 30. See the comments in the routine
- C MPBLAS for the reason for this choice.
- C
- C***SEE ALSO DQDOTA, DQDOTI, MPBLAS
- C***ROUTINES CALLED MPCHK, MPERR, MPMLP, MPNZR
- C***COMMON BLOCKS MPCOM
- C***REVISION HISTORY (YYMMDD)
- C 791001 DATE WRITTEN
- C ?????? Modified for use with BLAS. Blank COMMON changed to named
- C COMMON. R given dimension 12.
- C 890531 Changed all specific intrinsics to generic. (WRB)
- C 891214 Prologue converted to Version 4.0 format. (BAB)
- C 900402 Added TYPE section. (WRB)
- C 930124 Increased Array size in MPCON for SUN -r8. (RWC)
- C***END PROLOGUE MPMUL
- COMMON /MPCOM/ B, T, M, LUN, MXR, R(30)
- INTEGER B, T, R, X(*), Y(*), Z(*), RS, RE, XI, C, RI
- C***FIRST EXECUTABLE STATEMENT MPMUL
- CALL MPCHK (1, 4)
- I2 = T + 4
- I2P = I2 + 1
- C FORM SIGN OF PRODUCT
- RS = X(1)*Y(1)
- IF (RS.NE.0) GO TO 10
- C SET RESULT TO ZERO
- Z(1) = 0
- RETURN
- C FORM EXPONENT OF PRODUCT
- 10 RE = X(2) + Y(2)
- C CLEAR ACCUMULATOR
- DO 20 I = 1, I2
- 20 R(I) = 0
- C PERFORM MULTIPLICATION
- C = 8
- DO 40 I = 1, T
- XI = X(I+2)
- C FOR SPEED, PUT THE NUMBER WITH MANY ZEROS FIRST
- IF (XI.EQ.0) GO TO 40
- CALL MPMLP (R(I+1), Y(3), XI, MIN (T, I2 - I))
- C = C - 1
- IF (C.GT.0) GO TO 40
- C CHECK FOR LEGAL BASE B DIGIT
- IF ((XI.LT.0).OR.(XI.GE.B)) GO TO 90
- C PROPAGATE CARRIES AT END AND EVERY EIGHTH TIME,
- C FASTER THAN DOING IT EVERY TIME.
- DO 30 J = 1, I2
- J1 = I2P - J
- RI = R(J1) + C
- IF (RI.LT.0) GO TO 70
- C = RI/B
- 30 R(J1) = RI - B*C
- IF (C.NE.0) GO TO 90
- C = 8
- 40 CONTINUE
- IF (C.EQ.8) GO TO 60
- IF ((XI.LT.0).OR.(XI.GE.B)) GO TO 90
- C = 0
- DO 50 J = 1, I2
- J1 = I2P - J
- RI = R(J1) + C
- IF (RI.LT.0) GO TO 70
- C = RI/B
- 50 R(J1) = RI - B*C
- IF (C.NE.0) GO TO 90
- C NORMALIZE AND ROUND RESULT
- 60 CALL MPNZR (RS, RE, Z, 0)
- RETURN
- 70 WRITE (LUN, 80)
- 80 FORMAT (' *** INTEGER OVERFLOW IN MPMUL, B TOO LARGE ***')
- GO TO 110
- 90 WRITE (LUN, 100)
- 100 FORMAT (' *** ILLEGAL BASE B DIGIT IN CALL TO MPMUL,',
- 1 ' POSSIBLE OVERWRITING PROBLEM ***')
- 110 CALL MPERR
- Z(1) = 0
- RETURN
- END
|
