SUBROUTINE BITOUT(
C THIS ROUTINE WRITES OUT THE BITS IN IA
CHARACTER*16 CTEMP
CHARACTER*4 COUT(16)
INTEGER IIN(8)
DATA
COUT/'0000','0001','0010','0011','0100','0101','0110',
$
'0111','1000','1001','1010','1011','1100','1101','1110','1111'/
WRITE(CTEMP,'(Z8)')IA
READ(CTEMP,'(8Z1)')IIN
WRITE(LUN,*)(COUT(IIN(I)+1),I=1,8)
RETURN
END
SUBROUTINE DplMP(DX,IXE)
PARAMETER (IT=63)
C ADDS DOUBLE PRECISION DX TO MULTIPLE PRECISION IXE
DIMENSION IXE(*),IZ(2)
REAL*8 DX,TEMP
EQUIVALENCE (TEMP,IZ(1))
CWC DATA IAA/Z80000000/
C First find the bit equivalent of DX
IF(DX.EQ.0D0)RETURN
TEMP=DX
ITZ1=IZ(1)
IZ1L=32
ITZ2=IZ(2)
IZ2L=32
IF(ITZ2.LT.0)THEN
ISIGN=-1
ITZ2=IBCLR(ITZ2,31)
C ITZ2=ITZ2+IAA
ELSE
ISIGN=1
ENDIF
IZ2L=31
C EXAMINE DOUBLE PRECISION EXPONENT
ITT=ISHFT(ITZ2,-20)
IP2=-1023+ITT
IPB=IP2/14+1
IF(IP2.LT.0)IPB=IPB-1
C *** NORMALLY IPB IS LESS THAN IXE(1)
IOFF=IXE(1)-IPB
IF(IOFF.GT.IT-1)RETURN
C IX IXEIN
IXEOUT
C -1 1 1
C 0 B2 B2 - UNCHANGED
C 0 B3 B3
C A2 B4 A2+B4
C A3 B5 A3+B5
C A4 B6 A4+B6
C 0 B7 B7 - UNCHANGED
C ...
C 0 BIT+1 BIT+1 - UNCHANGED
C When This is not true, will
need to
C *** discard some of IXE, which requires renormalization first
IF(IOFF.LT.0)THEN
CALL MPNORM(IXE,0)
IDIFF=IPB-IXE(1)
IF(IDIFF.GT.0)THEN
IXE(1)=IPB
C *** IPB > IXE(1)
C IX IXEIN
IXEOUT
C 1 -1 1
C A2 0 A2
C A3 0 A3
C A4 B2 A4+B2
C A5 B3 A5+B3
C 0 B4 B4
C ...
C 0 BIT+1-2 BIT+1-2
DO I=IT+1,IDIFF+1,-1
IXE(I)=IXE(I-IDIFF)
ENDDO
DO
I=2,MIN(IDIFF+1,IT+1)
IXE(I)=0
ENDDO
IOFF=0
ELSE
IOFF=IXE(1)-IPB
IF(IOFF.GT.IT-1)RETURN
ENDIF
ENDIF
JFACT=IP2-14*(IPB-1)
C LEN IS THE LENGTH OF IX(3) IN IZ(2)
LEN=JFACT+1
c jfact is between 0 and 13
C ITT=ITT+2**20
C ! !20-jfact -- LEN IS JFACT+1 (CAN BE
1 TO 14)
C !31 !20 0!31 !IBB !0 !IB
C
0011111111110001100110011001100110011001100110011001100110011010--------
C REMOVE THE FIRST 11 SPACES FROM ITZ2
ITZ2=ISHFT(ITZ2,11)
IZ2L=32-11
C Double precision re-uses what is now position 31 with an assumed 1
ITZ2=IBSET(ITZ2,31)
C IF(ITZ2.GE.0)ITZ2=ITZ2+IAA
C !<--jfact+1--->!<---14---->!
C !20 0!31 !IBB !0 !IB
C 10001100110011001100110011001100110011001100110011010--------
C sign is reversed for a right shift, note that there are still 32
positions
IRS=LEN-32
IT1=ISHFT(ITZ2,IRS)
IXE(2+IOFF)=IXE(2+IOFF)+ISIGN*IT1
C NOW REMOVE LEN FROM ITZ2
ITZ2=ISHFT(ITZ2,LEN)
IZ2L=IZ2L-LEN
C !<---14---->!
C 0!31 !IBB !0 !IB
C 01100110011001100110011001100110011010--------
DO I=3,IT+1
C *** ARE THERE 14 SPACES LEFT IN IZ(2)
C *** AND ANY INTEGERS LEFT IN IXE
IF((IZ1L.GT.0).AND.(I+IOFF.LE.IT+1))THEN
IF(IZ2L.LT.14)THEN
C NEED TO RIGHT SHIFT
ITZ1 BY IZ2L
IRS=IZ2L
IT1=ISHFT(ITZ1,-IRS)
ITZ2=ITZ2+IT1
IZ2L=32
C REMOVE THE USED BITS FROM ITIZ1
ITZ1=ISHFT(ITZ1,32-IRS)
IZ1L=IZ1L-IRS
ENDIF
IRS=-18
IT1=ISHFT(ITZ2,IRS)
IXE(I+IOFF)=IXE(I+IOFF)+ISIGN*IT1
ITZ2=ISHFT(ITZ2,14)
IZ2L=IZ2L-14
ELSE
RETURN
ENDIF
IF(I+IOFF.GT.IT)RETURN
ENDDO
END
FUNCTION DPMP(IX)
C INPUT IS MULTIPLE PRECISION IX, RETURN IS A DOUBLE PRECISION DP
C CALL AS MPMDB(IX(2),DP) T USE BRENDT'S NUMBERS, There is a need for
C At most 5 terms after normalization.
DIMENSION IX(*)
REAL*8 BDP,DPMP
DATA IB/16384/
CALL MPNORM(IX,0)
IF(IX(1).LT.-66)THEN
DPMP=0D0
RETURN
ENDIF
IF(IX(1).GT.72)THEN
DPMP=1D307
RETURN
ENDIF
BDP=(1D0*IB)**(IX(1)-6)
DPMP=0
DO I=7,3,-1
DPMP=DPMP+IX(I)*BDP
BDP=BDP*IB
ENDDO
DPMP=DPMP+IX(2)*BDP
RETURN
END
SUBROUTINE DxMP(DX,IXE,IY)
PARAMETER (IT=63)
C Multiplies DOUBLE PRECISION DX times MULTIPLE PRECISION IXE
DIMENSION
IXE(*),IXT(2*IT+1),IZ(2),IDP(IT+1),IY(*)
REAL*8 DX,TEMP
EQUIVALENCE (TEMP,IZ(1))
IF(DX.EQ.0D0)THEN
CALL MPINIT(IY)
RETURN
ENDIF
TEMP=DX
ITZ1=IZ(1)
IZ1L=32
ITZ2=IZ(2)
IZ2L=32
IF(ITZ2.LT.0)THEN
ISIGN=-1
ITZ2=IBCLR(ITZ2,31)
ELSE
ISIGN=1
ENDIF
IZ2L=31
C EXAMINE DOUBLE PRECISION EXPONENT
ITT=ISHFT(ITZ2,-20)
IP2=-1023+ITT
IPB=IP2/14+1
IF(IP2.LT.0)IPB=IPB-1
JFACT=IP2-14*(IPB-1)
IDP(1)=IPB
C LEN IS THE LENGTH OF IX(3) IN IZ(2)
LEN=JFACT+1
c jfact is between 0 and 13
C ITT=ITT+2**20
C REMOVE THE FIRST 11 SPACES FROM ITZ2
ITZ2=ISHFT(ITZ2,11)
IZ2L=32-11
C Double precision re-uses what is now position 31 with an assumed 1
ITZ2=IBSET(ITZ2,31)
IRS=LEN-32
IT1=ISHFT(ITZ2,IRS)
IDP(2)=IT1
C NOW REMOVE LEN FROM ITZ2
ITZ2=ISHFT(ITZ2,LEN)
IZ2L=IZ2L-LEN
C !<---14---->!
C 0!31 !IBB !0 !IB
C 01100110011001100110011001100110011010--------
DO I=3,IT+1
C *** ARE THERE 14 SPACES LEFT IN IZ(2)
IF(IZ1L.GT.0)THEN
IF(IZ2L.LT.14)THEN
C NEED TO RIGHT SHIFT
ITZ1 BY IZ2L
IRS=IZ2L
IT1=ISHFT(ITZ1,-IRS)
ITZ2=ITZ2+IT1
IZ2L=32
C REMOVE THE USED BITS FROM ITIZ1
ITZ1=ISHFT(ITZ1,32-IRS)
IZ1L=IZ1L-IRS
ENDIF
IRS=-18
IT1=ISHFT(ITZ2,IRS)
IDP(I)=IT1
ITZ2=ISHFT(ITZ2,14)
IZ2L=IZ2L-14
ELSE
MDP=I
GOTO 15
ENDIF
ENDDO
MDP=IT+1
15 CONTINUE
DO J=2,MDP-1
IDP(J)=ISIGN*IDP(J)
ENDDO
DO J=MDP,IT+1
IDP(J)=0
ENDDO
CALL MPNORM(IDP,0)
CALL MPNORM(IDP,0)
CALL MPNORM(IXE,0)
IXT(1)=IXE(1)+IDP(1)-1
DO I=2,2*IT+1
IXT(I)=0
ENDDO
DO I=1,MDP-1
DO J=1,IT
K=I+J
IXT(K)=IXT(K)+IDP(I+1)*IXE(J+1)
ENDDO
ENDDO
CALL MPNORM(IXT,1)
DO I=1,IT+1
IY(I)=IXT(I)
ENDDO
RETURN
END
SUBROUTINE MPCLEAN(IX1)
C *** produces an ix1 with 16 decimal digits and the rest zeroes
PARAMETER (IT=63)
REAL*8 R1,DPMP
INTEGER
IX1(*),IX2(IT+1),IXT(IT+1)
CALL MPNORM(IX1,0)
R1=DPMP(IX1)
DO I=1,IT+1
IXT(I)=IX1(I)
ENDDO
C THE LOG IS OF X1, BUT THIS IS LIKELY TO BE OUTSIDE DP RANGE
C 4.21... IS LOG10(IB)
IF(R1.EQ.0D0)THEN
PRINT*,' R1 = 0, NO
CLEANING NEEDED '
RETURN
ENDIF
AL10=LOG10(ABS(R1))
IF(AL10.GT.0)THEN
JE=AL10+1
ELSE
JE=AL10
ENDIF
C MAKE R1 0.XXXXXX
IF(JE.GT.0)THEN
DO I=1,JE
CALL MPDIVI(IXT,IXT,10)
ENDDO
CALL MPNORM(IXT,0)
ELSEIF(JE.LT.0)THEN
DO I=1,-JE
CALL MPXI(IXT,IXT,10)
IF(I.EQ.4*(I/4))CALL
MPNORM(IXT,0)
ENDDO
ENDIF
DO I=1,4
CALL MPXI(IXT,IXT,10000)
R1=DPMP(IXT)
IC=R1
IX2(1)=1
IX2(2)=-IC
DO J=3,IT+1
IX2(J)=0
ENDDO
CALL MPPLMP(IX2,IXT,IXT)
CALL MPNORM(IXT,0)
ENDDO
DO I=1,4
CALL
MPDIVI(IXT,IXT,10000)
ENDDO
R1=DPMP(IXT)
IF(JE.LT.0)THEN
DO I=1,-JE
CALL MPDIVI(IXT,IXT,10)
ENDDO
CALL MPNORM(IXT,0)
ELSEIF(JE.GT.0)THEN
DO I=1,JE
CALL MPXI(IXT,IXT,10)
IF(I.EQ.4*(I/4))CALL
MPNORM(IXT,0)
ENDDO
ENDIF
DO I=2,IT+1
IXT(I)=-IXT(I)
ENDDO
CALL MPPLMP(IX1,IXT,IX1)
RETURN
END
SUBROUTINE
MPdivI(IX,IY,IDIV)
PARAMETER (IT=63,IB=16384)
DIMENSION
IX(*),IXT(IT+5),IY(*)
IF(ABS(IDIV).GT.32768)THEN
CALL MPINIT(IXT)
CALL MPPLI(IXT,IXT,IDIV)
CALL MPINV(IXT,IXT)
CALL MPXMP(IX,IXT,IY)
RETURN
ENDIF
IST=1
IF(IDIV.LT.0)IST=-1
CALL MPNORM(IX,0)
IF(IX(2).LT.0)THEN
IST=-IST
DO I=2,IT+1
IXT(I-1)=-IX(I)
ENDDO
ELSE
DO I=2,IT+1
IXT(I-1)=IX(I)
ENDDO
ENDIF
DO I=IT+1,IT+5
IXT(I)=0
ENDDO
DO I=1,IT+4
IT1=ABS(IXT(I))/ABS(IDIV)
IREM=ABS(IXT(I))-IT1*ABS(IDIV)
IXT(I)=SIGN(IT1,IXT(I)*IST)
IXT(I+1)=IXT(I+1)+SIGN(IREM*IB,IXT(I)*IST)
ENDDO
ISMALL=IT+6
DO I=1,IT+5
IF(IXT(I).NE.0)THEN
ISMALL=I
GOTO 15
ENDIF
ENDDO
15 CONTINUE
IF(ISMALL.EQ.IT+6)THEN
IY(1)=1
DO I=2,IT+1
IY(I)=0
ENDDO
RETURN
ENDIF
IY(1)=IX(1)+1-ISMALL
ITO=MIN(IT+1-ISMALL,IT-1)
DO I=2,ITO+2
IY(I)=IXT(I+ISMALL-2)
ENDDO
DO I=ITO+3,IT+1
IY(I)=0
ENDDO
RETURN
END
SUBROUTINE
MPEXP(MPXP,MPEXPP)
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (IT=63)
DIMENSION
MPEXPP(*),MPXP(*)
DIMENSION
MPTEMP(IT+1),MPXM(IT+1)
DIMENSION MPLN2(64)
C MPLN2 WAS CALCLUATED WITH IT=63 USING CODE IN MPLOG63.ZIP - LOG2.OUT
DATA MPLN2/1, 1,-5027,-7809, 8007, 3962,-5390, 7739,
3680,
2 1011, -596,-3057, 3341,-6541,-6610,-7463,
3 -6092, 5980,-7445,-4189,-4193,-1930, 2075,
7867,-3487, 5525,
4 5311,-2896,-6036, 4333, 2988,-7332,
5 1249,-7868, 2518,-3149,-7098,-5704,
2384,-5779,-6055,-7605,
6 -1487, 4525,-4248, 4401, 7950, 2028,
7 -6230,
7108,-5042,-6652,-4970,-4191,-1250,-6785, 7470,-3377,
8 -7700,
6614, 5411, 2806,-2070,-5164/
XP=DPMP(MPXP)
IF(ABS(XP).GT.5.209827825D9)THEN
MPEXPP(1)=536870912
IF(XP.GT.0)THEN
MPEXPP(2)=1
ELSE
MPEXPP(2)=-1
ENDIF
DO I=3,IT+1
MPEXPP(I)=0
ENDDO
RETURN
ENDIF
KLN2=ABS(XP)/LOG(2D0)+.5D0
IF(XP.LT.0)KLN2=-KLN2
CALL
MPXI(MPLN2,MPTEMP,KLN2)
CALL
MPmiMP(MPXP,MPTEMP,MPXM)
CALL MPEXPS(MPXM,MPEXPP)
KMOD=KLN2/14
KMULT=KLN2-KMOD*14
MMULT=2**ABS(KMULT)
MPEXPP(1)=MPEXPP(1)+KMOD
IF(KMULT.GT.0)THEN
CALL
MPXI(MPEXPP,MPEXPP,MMULT)
ELSEIF(KMULT.LT.0)THEN
CALL MPDIVI(MPEXPP,MPEXPP,MMULT)
ENDIF
CALL MPNORM(MPEXPP,0)
RETURN
END
SUBROUTINE MPEXPS(MPXP,MPEXP)
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (IT=63)
DIMENSION
MPEXP(*),MPXP(*),MPADD(IT+1),MPX0(IT+1),
2
MP2C(IT+1),MPCSQ(IT+1),MPC(IT+1)
LOGICAL LTEST
AN=SQRT((IT-1)*14D0)
N=AN+1
XP=DPMP(MPXP)
IF(ABS(XP).NE.0)THEN
AK=(1+(IT-1)*14)/AN+(LOG(ABS(XP))-LOG(AN))/LOG(2D0)
K=AK+1
K=MAX(1,K)
ELSE
K=0
LTEST=.TRUE.
CALL MPINIT(MPEXP)
MPEXP(1)=1
MPEXP(2)=1
DO I=3,IT+1
MPEXP(I)=0
ENDDO
DO I=1,IT+1
MPX0(I)=MPXP(I)
IF(MPX0(I).NE.0)LTEST=.FALSE.
ENDDO
IF(LTEST)RETURN
ENDIF
KMOD=K/14
KDIV=K-KMOD*14
MDIV=2**KDIV
CALL
MPDIVI(MPXP,MPX0,MDIV)
MPX0(1)=MPX0(1)-KMOD
CALL MPINIT(MPADD)
CALL
MPPLMP(MPX0,MPADD,MPADD)
CALL MPINIT(MPEXP)
CALL
MPPLMP(MPADD,MPEXP,MPEXP)
CALL MPNORM(MPEXP,0)
DO I=2,N
CALL
MPDIVI(MPADD,MPADD,I)
CALL
MPXMP(MPX0,MPADD,MPADD)
CALL MPNORM(MPADD,0)
CALL
MPPLMP(MPEXP,MPADD,MPEXP)
CALL MPNORM(MPEXP,0)
ENDDO
C THIS LOOP IS OF (1+C)**2 --> 1+2C+C**2
DO I=1,IT+1
MPC(I)=MPEXP(I)
ENDDO
DO I=1,K
CALL
MPXMP(MPC,MPC,MPCSQ)
CALL MPXI(MPC,MP2C,2)
CALL
MPPLMP(MP2C,MPCSQ,MPC)
CALL MPNORM(MPC,0)
ENDDO
CALL MPPLI(MPC,MPEXP,1)
CALL MPNORM(MPEXP,0)
RETURN
END
SUBROUTINE MPINIT(IX)
PARAMETER (IT=63)
DIMENSION IX(*)
IX(1)=-536870912
DO I=2,IT+1
IX(I)=0
ENDDO
RETURN
END
SUBROUTINE MPINV(IXP,IXIP)
C F(IXI)=1-IXI*IX=ERR
C F'(IXI)=-IX
C IXI=IXI*(2-IXI*IX) -- SEE
MPINV.DOC
PARAMETER (IT=63)
DIMENSION
IXP(*),IXIP(*),IX(IT+1),IXT(IT+1),IXI(IT+1)
DIMENSION
I2(IT+1),IXIX(IT+1),I2MXIX(IT+1)
REAL*8 DPMP,XI,X
IX(1)=0
DO I=2,IT+1
IX(I)=IXP(I)
ENDDO
X=DPMP(IX)
IF(X.EQ.0D0)THEN
PRINT*,' NO INVERSE
DELIBERATE CRASH'
IXI(1)=536870912
IXI(2)=1
DO I=3,IT+1
IXI(I)=0
ENDDO
XI=1/X
RETURN
ENDIF
XI=1/X
CALL MPINIT(IXT)
CALL DPLMP(XI,IXT)
CALL MPINIT(I2)
CALL DPLMP(2D0,I2)
DO I=1,5
CALL MPXMP(IXT,IX,IXIX)
CALL
MPmiMP(I2,IXIX,I2MXIX)
CALL
MPXMP(IXT,I2MXIX,IXT)
ENDDO
IXIP(1)=IXT(1)-IXP(1)
DO I=2,IT+1
IXIP(I)=IXT(I)
ENDDO
RETURN
END
SUBROUTINE MPLOG(IXP,ILOGX)
IMPLICIT REAL*8 (A-H,O-Z)
C description in mpkroot2.doc and ln.doc
C AM=2d0**K <-- AM is probably beyond the range of the integers
C (IC+1)**M=IX+1
C F'(IC=-(M)*(IC+1)**(M-1)
C IC=IC-(1+IC)(1/(M+1))(((1+IKR)**2)(K TIMES)-IX-1)
PARAMETER (IT=63,IB=16384)
DIMENSION IXP(*),ILOGX(*)
DIMENSION
IXM1(IT+1),ILOGT(IT+1)
DIMENSION
IC(IT+1),IX(IT+1),ILOGB(64)
C THE LOG OF IB IS CALCLUATED WITH IT=63 IN MPLOG63.ZIP
DATA ILOGB/ 1,
10,-4849, 5369,-2587, 6311, 6467,-6339,
2 2369,-2231, 8037, 6357,-2384, 6724,
5758,-6183,
3 -3363, 1794,-5930, 6886,
6832, 5750,-3711,-4553, 339,-4565,
4
-7568,-7781,-2580,-4871,-7326,-4343,
5 1095, 4538, 2481, 5060,-1073, 2066, 603, 1009,-2856,-8167,
6 -4430,-2190,
6068,-3915,-3386,-4381,
7 -5294, 1204,-5058,
5172,-4048, 6861,-1122, 3320, 6273, 1867,
8 6894,-5703,-6164, 6514,
3785,-2536/
Mest=10*SQRT(1D0*IT)
MTEMP=IXP(1)
IX(1)=0
DO I=2,IT+1
IX(I)=IXP(I)
ENDDO
CALL MPPLI(IX,IXM1,-1)
C estimate the number of operations
C1=DPMP(IXM1)
IF(C1.GT.2.75D303)THEN
XSUM=0
AMULT=1/(1D0*IB)
AMULTP=AMULT
DO I=1,IT
XSUM=XSUM+AMULTP*IXM1(I+1)
AMULTP=AMULT*AMULT
ENDDO
ALABC1=IXM1(1)*LOG(1D0*IB)+LOG(ABS(XSUM))
ELSEIF(ABS(C1).EQ.0)THEN
CALL MPINIT(ILOGX)
RETURN
ELSE
ALABC1=LOG(ABS(C1))
ENDIF
IF(ABS(C1).LT.1D-1)THEN
CALL MPLNS(IXM1,ILOGX)
RETURN
ELSEIF(ABS(C1).LT.1D6)THEN
K=14
ELSE
K=LOG(ALABC1*100)/LOG(2d0)
ENDIF
CALL MPKROOT(IC,IXM1,K)
CALL MPLNS(IC,ILOGX)
KMOD=K/14
KDIV=K-KMOD*14
MDIV=2**KDIV
ILOGX(1)=ILOGX(1)+KMOD
CALL MPXI(ILOGX,ILOGX,MDIV)
IX(1)=MTEMP
CALL
MPXI(ILOGB,ILOGT,MTEMP)
CALL
MPPLMP(ILOGX,ILOGT,ILOGX)
RETURN
END
SUBROUTINE MPKROOT(IC,IX,K)
IMPLICIT REAL*8 (A-H,O-Z)
C description in mpkroot2.doc
C SOLVES FOR IC SUCH THAT (1+IC)**2**K = 1+IX
C AM=2d0**K <-- AM is probably beyond the range of the integers
C (IKR+1)**M=IX+1
C F'(IKR)=-(M)*(IKR+1)**(M-1)
C IKR=IKR-(1+IKR)(1/(M+1))(((1+IKR)**2)(K TIMES)-IX-1)
PARAMETER (IT=63,IB=16384)
DIMENSION IX(*),IC(*)
DIMENSION
ICT(IT+1),ICT2(IT+1),IFC(IT+1),IXP1I(IT+1),ICP1(IT+1)
CALL MPPLI(IX,IXP1I,1)
CALL MPINV(IXP1I,IXP1I)
AM=2D0**K
c M=2**K
X=DPMP(IX)
IF(X+1.GT.1D-15)THEN
ALNM=LOG(X+1)
ELSE
XSUM=0
AMULT=1/(1D0*IB)
AMULTP=AMULT
DO I=1,IT
XSUM=XSUM+AMULTP*IX(I+1)
AMULTP=AMULT*AMULT
ENDDO
ALNM=IX(1)*LOG(1D0*IB)+LOG(ABS(XSUM))
ENDIF
RAT=ALNM/AM
IF(ABS(X).LT.1D-8)THEN
C1=X/AM+X**2/AM
ELSEIF(ABS(RAT).LT.1D-8)THEN
C1=RAT*(1-.5D0*(RAT*(1-RAT/3)))
ELSEIF(1+X.LT.1)THEN
C1=(1+X)**(1/AM)-1
ELSE
C1=EXP(ALNM/AM)-1
ENDIF
EPS=ABS(C1*(1D0*IB)**(-IT))
OPC=1+C1
KMOD=K/14
KDIV=K-KMOD*14
MDIV=2**KDIV
CALL MPINIT(IC)
CALL DPLMP(C1,IC)
DO I=1,5
C *** FORM (1+C)**2**k -1 =
C c --> 2c+c**2 first squaring 1+c is now (1+c)**2 -- first 2
to the 1
C c --> 2c+c**2 second squaring 1+c is now (1+c)**2**2 2 to the 2
c ...
c c -->
C OPC=DPMP(IC)+1
CALL MPPLI(IC,ICP1,1)
C AMULT=1/AM
DO J=1,IT+1
ICT(J)=IC(J)
ENDDO
DO J=1,K
CALL
MPXMP(ICT,ICT,ICT2)
CALL MPXI(ICT,ICT,2)
CALL
MPPLMP(ICT,ICT2,ICT)
ENDDO
CALL MPMIMP(ICT,IX,IFC)
CALL MPXMP(IXP1I,IFC,IFC)
CALL MPXMP(ICP1,IFC,IFC)
IFC(1)=IFC(1)-KMOD
IF(MDIV.NE.1)CALL
MPDIVI(IFC,IFC,MDIV)
c FC=DPMP(IFC)
c PRINT*,' from mpkroot2
i,fc',i,fc
CALL MPMIMP(IC,IFC,IC)
ENDDO
RETURN
END
SUBROUTINE MPLNS(IC,ILOG)
IMPLICIT REAL*8 (A-H,O-Z)
C description in mpkroot2.doc
C code finds ln(1+c) by expanding as c-c*c/2 etc.
C *** THE GAIN PER STEP IS ~ 16 DIGITS, THUS isteps=4.2*it/16
PARAMETER (IT=63,IB=16384)
DIMENSION IC(*),ILOG(*)
DIMENSION
ICT(IT+1),ICT2(IT+1),ILOGT(IT+1)
C1=DPMP(IC)
M1=(-IT*LOG(1D0*IB)+LOG(2D0))/LOG(ABS(C1))
M=M1-LOG(1D0*(M1+1))/LOG(ABS(C1))
C1I=ABS(1/C1)
DO I=1,IT+1
ILOGT(I)=IC(I)
ICT(I)=IC(I)
ENDDO
DO I=2,M,2
CALL MPXMP(IC,ICT,ICT)
CALL MPDIVI(ICT,ICT2,I)
CALL
MPMIMP(ILOGT,ICT2,ILOGT)
CALL MPXMP(IC,ICT,ICT)
CALL MPDIVI(ICT,ICT2,I+1)
CALL
MPPLMP(ILOGT,ICT2,ILOGT)
ENDDO
DO I=1,IT+1
ILOG(I)=ILOGT(I)
ENDDO
RETURN
END
SUBROUTINE
MPNLOGS(MPX,MPLOG)
C see nrlog.doc
PARAMETER (IT=63)
DIMENSION MPX(*),MPLOG(*),MPEXP1(IT+1),MPLOGN(IT+1)
DIMENSION MPTEMP(IT+1)
REAL*8 DPMP
CALL MPINIT(MPLOG)
CALL
DPLMP(LOG(DPMP(MPX)),MPLOG)
DO I=1,5
MPLOGN(1)=MPLOG(1)
DO J=2,IT+1
MPLOGN(J)=-MPLOG(J)
ENDDO
CALL
MPEXPS(MPLOGN,MPEXP1)
CALL
MPXMP(MPEXP1,MPX,MPTEMP)
CALL
MPPLI(MPTEMP,MPTEMP,-1)
PRINT*,' I, ERR
',I,DPMP(MPTEMP)
CALL
MPPLMP(MPLOG,MPTEMP,MPLOG)
ENDDO
RETURN
END
SUBROUTINE MPNLOG(MPX,MPLOG)
C see nrlog.doc
PARAMETER (IT=63)
DIMENSION
MPX(*),MPLOG(*),MPEXP1(IT+1),MPLOGN(IT+1)
DIMENSION MPTEMP(IT+1)
REAL*8 DPMP
CALL MPINIT(MPLOG)
CALL
DPLMP(LOG(DPMP(MPX)),MPLOG)
DO I=1,5
MPLOGN(1)=MPLOG(1)
DO J=2,IT+1
MPLOGN(J)=-MPLOG(J)
ENDDO
CALL MPEXP(MPLOGN,MPEXP1)
CALL
MPXMP(MPEXP1,MPX,MPTEMP)
CALL
MPPLI(MPTEMP,MPTEMP,-1)
PRINT*,' I, ERR
',I,DPMP(MPTEMP)
CALL MPPLMP(MPLOG,MPTEMP,MPLOG)
ENDDO
RETURN
END
SUBROUTINE MPNORM(IX,IDIR)
PARAMETER
(IT=63,IB=16384,IBS2=8192)
C *** IX MUST BE DIMENSIONED AT LEAST 6 MORE THAN IT
DIMENSION IX(*),IXT(2*IT+6)
ILIM=IT
IF(IDIR.EQ.1)ILIM=2*IT
DO I=2,ILIM+1
IXT(I+5)=IX(I)
ENDDO
DO I=1,6
IXT(I)=0
ENDDO
DO I=ILIM+6,2,-1
IF(IXT(I).GE.IBS2)THEN
IT1=IXT(I)/IB
IXT(I)=IXT(I)-IT1*IB
IXT(I-1)=IXT(I-1)+IT1
IF(IXT(I).GE.IBS2)THEN
IXT(I)=IXT(I)-IB
IXT(I-1)=IXT(I-1)+1
ENDIF
ELSEIF(IXT(I).LT.-IBS2)THEN
IT1=ABS(IXT(I))/IB
IXT(I)=IXT(I)+IT1*IB
IXT(I-1)=IXT(I-1)-IT1
IF(IXT(I).LT.-IBS2)THEN
IXT(I)=IXT(I)+IB
IXT(I-1)=IXT(I-1)-1
ENDIF
ENDIF
ENDDO
ISMALL=ILIM+6
DO I=1,ILIM+5
IF(IXT(I).NE.0)THEN
ISMALL=I
GOTO 15
ENDIF
ENDDO
15 CONTINUE
IF(ISMALL.EQ.ILIM+6)THEN
IX(1)=-536870912
DO I=2,IT+1
IX(I)=0
ENDDO
RETURN
ENDIF
IX(1)=IX(1)+7-ISMALL
ILIM=MIN(IT,IT+7-ISMALL)
DO I=2,ILIM+1
IX(I)=IXT(I+ISMALL-2)
ENDDO
DO I=ILIM+2,IT+1
IX(I)=0
ENDDO
RETURN
END
SUBROUTINE
MPOUT(MPVAL,CRET,NCHAR)
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (IT=63)
DIMENSION
MPLVAL(IT+1),MPFVAL(IT+1),MPL10(64),MPL10I(64)
DIMENSION
MPTEMP(IT+1),MPVAL(*)
CHARACTER*(*) CRET
DATA MPL10/ 1,
2, 4958,-7305, 6827,-7504, 5866, 5558,
2 -7496, 2892, 2601,-1884,
-304,-6295, 7673, 929,
3 -926,-1528,
496, 2909,-5789,-1013, 6066,-5628, 1485,-3645,
4 3008, 2840,-3881,-8014,-7644, 1552,
5 4286, 2924,-1137, 7713, 46, 5211, 5752, 2878, 2601,-5582,
6 4382, 7575,-5524, 3779, 6194, -260,
7 -2451,-1343,
2740,-3649,-6890, 1773, 2514, 1650,-1091,-5729,
8 -3237, 2264, 1288, 5690,
4919,-2926/
DATA MPL10I/ 0, 7115, 7877, 4718, 5177,-3417,-5411,
5471,
2 5786,-1146, 4599, 1676, 291,-1740, 4114,-2946,
3 5794, 7922,-2774,-1813,-7897,
8003,-8168,-6171, 2322,-5454,
4 5682, 3220, 2700, 6599, 7888,-1422,
5 -4514, 3572,
5575,-6766,-6444,-2456, 6838,-1485,-7064, 3583,
6 2560, 3648,
710, 953,-6445, 2839,
7 -3851, 2137, 619, 5281, 6583, -602,-1606,-3157,-5608,
1061,
8 5940, 6154,-5582, 4570,-7923, 2542/
CALL MPNORM(MPVAL,0)
DO I=2,IT+1
IF(MPVAL(I).NE.0)GOTO 7
ENDDO
CRET='0.0'
RETURN
7 CONTINUE
ISIGN=1
IF(MPVAL(2).LT.0)THEN
ISIGN=-1
MPTEMP(1)=MPVAL(1)
DO I=2,IT+1
MPTEMP(I)=-MPVAL(I)
ENDDO
ELSE
DO I=1,IT+1
MPTEMP(I)=MPVAL(I)
ENDDO
ENDIF
CALL MPLOG(MPTEMP,MPLVAL)
CALL
MPXMP(MPLVAL,MPL10I,MPLVAL)
VAL=DPMP(MPLVAL)
IF(VAL.GT.0)THEN
JE=VAL
ELSE
JE=-ABS(VAL-1)
ENDIF
CALL
MPPLI(MPLVAL,MPLVAL,-JE)
CALL MPNORM(MPLVAL,0)
CALL
MPXMP(MPLVAL,MPL10,MPLVAL)
CALL MPEXP(MPLVAL,MPFVAL)
CALL MPXI(MPFVAL,MPFVAL,100)
CALL MPNORM(MPFVAL,0)
MSUB=MPFVAL(2)
MDIG=NCHAR/3-1
DO J=3,MDIG+1
IF(J-MPFVAL(1).GT.IT+1)GOTO 15
IF(MPFVAL(J).GT.0)GOTO 15
IF(MPFVAL(J).LT.0)THEN
MSUB=MSUB-1
GOTO 15
ENDIF
ENDDO
15 CONTINUE
AOUT=MSUB
AOUT=AOUT/100
AOUT=AOUT*ISIGN
WRITE(CRET,'(F5.2)')AOUT
CALL
MPPLI(MPFVAL,MPFVAL,-MSUB)
DO I=1,MDIG
CALL
MPXI(MPFVAL,MPFVAL,1000)
CALL MPNORM(MPFVAL,0)
IF(MPFVAL(1).EQ.1)THEN
MSUB=MPFVAL(2)
ELSE
MSUB=0
GOTO 25
ENDIF
DO J=3,MDIG
IF(J-MPFVAL(1).GT.MDIG+1)GOTO 25
IF(MPFVAL(J).GT.0)GOTO
25
IF(MPFVAL(J).LT.0)THEN
MSUB=MSUB-1
GOTO 25
ENDIF
ENDDO
25 CONTINUE
NW=3*I+3
WRITE(CRET(NW:),'(I3)')MSUB
CALL
MPPLI(MPFVAL,MPFVAL,-MSUB)
CALL MPNORM(MPFVAL,0)
ENDDO
DO I=6,6+NCHAR
IF(CRET(I:I).EQ.'
')CRET(I:I)='0'
ENDDO
CRET(NCHAR:)=' '
CRET(5+NCHAR+1:5+NCHAR+3)='
E'
IF(ABS(JE).LT.1000)THEN
WRITE(CRET(5+NCHAR+4:),'(I5)')JE
ELSEIF(ABS(JE).LT.10000)THEN
WRITE(CRET(5+NCHAR+4:),'(I6)')JE
ELSE
WRITE(CRET(5+NCHAR+4:),*)JE
ENDIF
RETURN
END
SUBROUTINE MPOUT0(IXT,CRET)
PARAMETER (IT=63)
REAL*8 R1,DPMP
CHARACTER*300 CRET
CHARACTER*300 CTEMP
CHARACTER*15 CHEAD
INTEGER
IC(291),IX1(IT+1),IX2(IT+1),IXT(IT+1)
DO I=1,IT+1
IX1(I)=IXT(I)
ENDDO
call mpnorm(ix1,0)
R1=DPMP(IX1)
CRET=' '
ISIGN=1
IF(R1.LT.0)THEN
ISIGN=-1
DO I=2,IT+1
IX1(I)=-IX1(I)
ENDDO
ENDIF
C THE LOG IS OF X1, BUT THIS IS LIKELY TO BE OUTSIDE DP RANGE
C 4.21... IS LOG10(IB)
IF(R1.EQ.0D0)THEN
CTEMP='0.0'
RETURN
ENDIF
AL10=LOG10(ABS(R1))
IF(AL10.GT.0)THEN
JE=AL10+1
ELSE
JE=AL10
ENDIF
C MAKE R1 0.XXXXXX
IF(JE.GT.0)THEN
DO I=1,JE
CALL MPDIVI(IX1,IX1,10)
ENDDO
CALL MPNORM(IX1,0)
ELSEIF(JE.LT.0)THEN
DO I=1,-JE
CALL MPXI(IX1,IX1,10)
IF(I.EQ.4*(I/4))CALL
MPNORM(IX1,0)
ENDDO
ENDIF
DO I=1,290
CALL MPXI(IX1,IX1,10)
R1=DPMP(IX1)
IC(I)=R1
IX2(1)=1
IX2(2)=-IC(I)
DO J=3,IT+1
IX2(J)=0
ENDDO
CALL MPPLMP(IX2,IX1,IX1)
CALL MPNORM(IX1,0)
ENDDO
IF(IC(290).GE.5)IC(289)=IC(289)+1
DO I=289,2,-1
IF(IC(I).GE.10)THEN
IC(I)=IC(I)-10
IC(I-1)=IC(I-1)+1
ENDIF
ENDDO
IF(IC(1).EQ.10)THEN
IC(1)=1
JE=JE+1
ENDIF
WRITE(CTEMP,'(289I1)')(IC(I),I=1,289)
IF(JE.GE.0)THEN
CHEAD(1:2)=' E'
ELSE
CHEAD(1:2)='-E'
ENDIF
WRITE(CHEAD(3:),'(I9)')ABS(JE)
CHEAD(13:14)=' '
IF(ISIGN.EQ.-1)CHEAD(15:15)='-'
CRET=CHEAD//CTEMP(1:284)
RETURN
END
SUBROUTINE
MPPLI(IX,IY,IADD)
PARAMETER (IT=63)
DIMENSION IX(*),IY(*),IYT(IT+1)
REAL*8 DP
IF(ABS(IADD).GE.2147475456)THEN
DP=IADD
CALL DPLMP(DP,IY)
ELSE
IYT(1)=1
IYT(2)=IADD
DO I=3,IT+1
IYT(I)=0
ENDDO
CALL MPPLMP(IX,IYT,IY)
ENDIF
RETURN
END
SUBROUTINE
MPmiMP(IX1,IX2,IX3)
c subtracts ix2 from ix1
PARAMETER (IT=63)
DIMENSION
IX1(*),IX2(*),IX3(*),IXT(IT+1)
C IX1 IX2
IX3
C 3 1 3
C A2 -
A2
C A3 -
A3
C A4 B2
A4+B2
C ...
C A IT1+1 B IT1-2
A IT1+1 + B IT1 - 2
C B IT1-1 B IT1 -1
C ... ...
C B(IT2) B(IT2)
C 0
C This routine uses signed integers
IF(IX1(1).GE.IX2(1))THEN
IXT(1)=IX1(1)
I2OFF=IX1(1)-IX2(1)
ILIM=MIN(I2OFF+1,IT+1)
DO I=2,ILIM
IXT(I)=IX1(I)
ENDDO
ILIM2=MIN(IT+1,IT+I2OFF+1)
DO I=ILIM+1,ILIM2
IXT(I)=IX1(I)-IX2(I-I2OFF)
ENDDO
ILIM3=MIN(IT+I2OFF+1,IT)
DO I=ILIM2+1,ILIM3
IXT(I)=-IX2(I-I2OFF)
ENDDO
DO I=ILIM2+1,IT+1
IXT(I)=0
ENDDO
ELSE
IXT(1)=IX2(1)
I1OFF=IX2(1)-IX1(1)
ILIM=MIN(I1OFF+1,IT+1)
DO I=2,ILIM
IXT(I)=-IX2(I)
ENDDO
ILIM2=MIN(IT+1,IT+I1OFF+1)
DO I=ILIM+1,ILIM2
IXT(I)=IX1(I-I1OFF)-IX2(I)
ENDDO
DO I=ILIM2+1,IT+1
IXT(I)=0
ENDDO
ENDIF
DO I=1,IT+1
IX3(I)=IXT(I)
ENDDO
RETURN
END
SUBROUTINE
MPplMP(IX1,IX2,IX3)
PARAMETER (IT=63)
DIMENSION
IX1(*),IX2(*),IX3(*),IXT(IT+1)
C IX1 IX2
IX3
C 3 1 3
C A2 -
A2
C A3 -
A3
C A4 B2
A4+B2
C ...
C A IT1+1 B IT1-2
A IT1+1 + B IT1 - 2
C B IT1-1 B IT1 -1
C ...
...
C B(IT2) B(IT2)
C 0
C This routine uses signed integers
IF(IX1(1).GE.IX2(1))THEN
IXT(1)=IX1(1)
I2OFF=IX1(1)-IX2(1)
C INTEGER OVERFLOW CAN CAUSE I2OFF TO BE < 0
IF(I2OFF.LT.0)THEN
DO I=1,IT+1
IX3(I)=IX1(I)
ENDDO
RETURN
ENDIF
ILIM=MIN(I2OFF+1,IT+1)
DO I=2,ILIM
IXT(I)=IX1(I)
ENDDO
ILIM2=MIN(IT+1,IT+I2OFF+1)
DO I=ILIM+1,ILIM2
IXT(I)=IX1(I)+IX2(I-I2OFF)
ENDDO
ILIM3=MIN(IT+I2OFF+1,IT)
DO I=ILIM2+1,ILIM3
IXT(I)=IX2(I-I2OFF)
ENDDO
DO I=ILIM2+1,IT+1
IXT(I)=0
ENDDO
ELSE
IXT(1)=IX2(1)
I1OFF=IX2(1)-IX1(1)
C INTEGER OVERFLOW CAN CAUSE I1OFF TO BE < 0
IF(I1OFF.LT.0)THEN
DO I=1,IT+1
IX3(I)=IX2(I)
ENDDO
RETURN
ENDIF
ILIM=MIN(I1OFF+1,IT+1)
DO I=2,ILIM
IXT(I)=IX2(I)
ENDDO
ILIM2=MIN(IT+1,IT+I1OFF+1)
DO I=ILIM+1,ILIM2
IXT(I)=IX1(I-I1OFF)+IX2(I)
ENDDO
DO I=ILIM2+1,IT+1
IXT(I)=0
ENDDO
ENDIF
DO I=1,IT+1
IX3(I)=IXT(I)
ENDDO
RETURN
END
SUBROUTINE
MPSQRT(IX,ISQR,IFL)
C F(ISQR)=X-ISQR*ISQR=ERR
C F'(ISQR)=-2*ISQR
C ISQR=ISQR+(1/2)(X-ISQR*ISQR)/ISQR
C ISQR=ISQR(1+(1/2)XINV(X-ISQR*ISQR)
C ***
C ISQR=ISQR*(3-IXI*ISQR*ISQR)/2
-- SEE MPINV.DOC
PARAMETER (IT=63)
DIMENSION IX(*),ISQR(*)
DIMENSION
I3(IT+1),IXI(IT+1),IXT(IT+1),IX1(IT+1)
REAL*8 DSQR,DPMP,X
ITEMP=IX(1)/2
IX1(1)=IX(1)-2*ITEMP
DO I=2,IT+1
IX1(I)=IX(I)
ENDDO
X=DPMP(IX1)
IFL=1
IF(X.LE.0D0)THEN
PRINT*,' NO SQRT '
IFL=-1
RETURN
ENDIF
DSQR=SQRT(X)
CALL MPINIT(ISQR)
CALL DPLMP(DSQR,ISQR)
CALL MPINV(IX1,IXI)
CALL MPINIT(I3)
CALL DPLMP(3D0,I3)
DO I=1,5
C *** FORM C**2
CALL MPXMP(ISQR,ISQR,IXT)
C *** MULTIPLY BY 1/X
CALL MPXMP(IXI,IXT,IXT)
C *** SUBTRACT FROM 3
CALL MPMIMP(I3,IXT,IXT)
C *** MULTIPY BY THE CURRENT ESTIMATOR
CALL MPXMP(ISQR,IXT,IXT)
C *** DIVIDE BY 2
CALL MPDIVI(IXT,ISQR,2)
ENDDO
ISQR(1)=ISQR(1)+ITEMP
RETURN
END
SUBROUTINE MPxI(IX,IY,IMUL)
PARAMETER (IT=63)
DIMENSION IX(*),IY(*),IYT(IT+1)
CALL MPNORM(IX,0)
C NOTE THAT THE RETURNED IY IS NOT NORMALIZED.
IF(ABS(IMUL).GT.131072)THEN
c following is from "Mult by large integer.doc"
K1=IMUL/16384
K2=IMUL-K1*16384
IYT(1)=IX(1)+1
DO I=2,IT+1
IYT(I)=K1*IX(I)
ENDDO
CALL MPNORM(IYT,0)
IY(1)=IX(1)
DO I=2,IT+1
IY(I)=K2*IX(I)
ENDDO
CALL MPNORM(IY,0)
CALL MPPLMP(IY,IYT,IY)
RETURN
ENDIF
IY(1)=IX(1)
DO I=2,IT+1
IY(I)=IMUL*IX(I)
ENDDO
RETURN
END
SUBROUTINE MPxMP(IX,IY,IZ)
PARAMETER (IT=63)
DIMENSION
IX(*),IY(*),IZ(*),IZT(2*IT+1)
CALL MPNORM(IX,0)
CALL MPNORM(IY,0)
IF(IX(1).GE.536870912.AND.IY(1).GT.0)THEN
C *** INFINITY = IB**536870912 --> INFINITY CAN ONLY MOVE DOWN
DO I=1,IT+1
IZ(I)=IX(I)
ENDDO
RETURN
ENDIF
IF(IY(1).GE.536870912.AND.IX(1).GT.0)THEN
C *** INFINITY = IB**536870912 --> INFINITY CAN ONLY MOVE DOWN
DO I=1,IT+1
IZ(I)=IY(I)
ENDDO
RETURN
ENDIF
IZT(1)=IX(1)+IY(1)-1
DO I=2,2*IT+1
IZT(I)=0
ENDDO
DO I=1,IT
DO J=1,IT
K=I+J
IZT(K)=IZT(K)+IX(I+1)*IY(J+1)
ENDDO
ENDDO
CALL MPNORM(IZT,1)
DO I=1,IT+1
IZ(I)=IZT(I)
ENDDO
RETURN
END