PROGRAM SAMPLE C C David M. Smith 9-17-96 C C This is a test program for ZMLIB 1.1, a multiple-precision real C arithmetic package. A few example ZM calculations are carried C out using 30 significant digit precision. C C This program uses both ZMLIB.f and FMLIB.f. C C The output is saved in file ZMSAMPLE.LOG. A comparison file, C ZMSAMPLE.CHK, is provided showing the expected output from 32-bit C (IEEE arithmetic) machines. When run on other computers, all the C numerical results should still be the same, but the number of terms C needed for some of the results might be slightly different. The C program checks all the results and the last line of the log file C should be "All results were ok." C C----------------------------------------------------------------------- C C These five common blocks contain information that must be saved C between calls, so they should be declared in the main program. C The parameter statement defines array sizes and pointers, and C contains the FMLIB parameters, followed by ZMLIB parameters. C C INTEGER NDIGMX,NBITS,LPACK,LUNPCK,LMWA,LJSUMS,LMBUFF,LPACKZ, * LUNPKZ,KPTIMP,KPTIMU,LMBUFZ C PARAMETER ( NDIGMX=256 , NBITS=64 , * LPACK = (NDIGMX+1)/2+1 , LUNPCK = (6*NDIGMX)/5+20, * LMWA = 2*LUNPCK , LJSUMS = 8*(LUNPCK+2), * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000+6, * LPACKZ = 2*LPACK+1 , LUNPKZ = 2*LUNPCK+1, * KPTIMP = LPACK+1 , KPTIMU = LUNPCK+1, * LMBUFZ = 2*LMBUFF+10 ) C INTEGER JFORMZ,JPRNTZ C COMMON /ZMUSER/ JFORMZ,JPRNTZ C INTEGER NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC,KFLAG, * KWARN,KROUND,KSWIDE,KESWCH,KDEBUG DOUBLE PRECISION MBASE C COMMON /FMUSER/ MBASE,NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC, * KFLAG,KWARN,KROUND,KSWIDE,KESWCH,KDEBUG C DOUBLE PRECISION DPMAX INTEGER NCALL,KACCSW,IUNKNO,NDG2MX,INTMAX,KSUB REAL RUNKNO,SPMAX DOUBLE PRECISION MWA,MXEXP,MXEXP2,MEXPUN,MEXPOV,MUNKNO,MXBASE, * MAXINT C COMMON /FM/ MWA(LMWA),NCALL,KACCSW,MXEXP,MXEXP2,MEXPUN,MEXPOV, * MUNKNO,IUNKNO,RUNKNO,MXBASE,NDG2MX,SPMAX,DPMAX, * MAXINT,INTMAX,KSUB C INTEGER LHASH1,LHASH2 PARAMETER (LHASH1=0 , LHASH2=256) DOUBLE PRECISION DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB INTEGER NDIGPI,NDIGE,NDIGLB,NDIGLI,KHASHT,KHASHV, * NGRD21,NGRD52,NGRD22 DOUBLE PRECISION MBSPI,MBSE,MBSLB,MBSLI,MPISAV,MESAV,MLBSAV, * MLN1,MLN2,MLN3,MLN4,MBLOGS,MEXPAB REAL ALOGMB,ALOGMT,ALOGM2,ALOGMX C COMMON /FMSAVE/ NDIGPI,NDIGE,NDIGLB,NDIGLI,MBSPI,MBSE,MBSLB,MBSLI, * MPISAV(0:LUNPCK),MESAV(0:LUNPCK),MLBSAV(0:LUNPCK), * MLN1(0:LUNPCK),MLN2(0:LUNPCK),MLN3(0:LUNPCK), * MLN4(0:LUNPCK),MBLOGS,MEXPAB,ALOGMB,ALOGM2,ALOGMX, * ALOGMT,DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB,KHASHT(LHASH1:LHASH2), * KHASHV(LHASH1:LHASH2),NGRD21,NGRD52,NGRD22 C CHARACTER CMBUFF,CMCHAR CHARACTER *6 NAMEST C COMMON /FMBUFF/ CMBUFF(LMBUFF),NAMEST(0:50),CMCHAR C C----------------------------------------------------------------------- C C Declare arrays for ZM variables. All are in C unpacked format. C DOUBLE PRECISION MA(0:LUNPKZ),MB(0:LUNPKZ),MC(0:LUNPKZ), * MD(0:LUNPKZ),MAFM(0:LUNPCK),MBFM(0:LUNPCK) C C Character strings used for input and output. C CHARACTER *80 ST1 INTEGER ITER,K,KLOG,NERROR LOGICAL FMCOMP C C Set precision to give at least 30 significant digits C and initialize both the ZMLIB and FMLIB packages. C C Note that any program using the ZM package MUST call C ZMSET before using the package. C CALL ZMSET(30) C NERROR = 0 C C Write output to the standard FM output (unit KW, defined C in subroutine FMSET), and also to the file ZMSAMPLE.LOG. C KLOG = 18 OPEN (KLOG,FILE='ZMSAMPLE.LOG') C C 1. Find a complex root of the equation C f(x) = x**5 - 3x**4 + x**3 - 4x**2 + x - 6 = 0. C C Newton's method with initial guess x = .56 + 1.06 i. C This version is not tuned for speed. See the ZMSQRT C routine for possible ways to increase speed. C Horner's rule is used to evaluate the function: C f(x) = ((((x-3)*x+1)*x-4)*x+1)*x-6. C C MA is the previous iterate. C MB is the current iterate. C CALL ZMST2M('.56 + 1.06 i',MA) C C Print the first iteration. C WRITE (KW,110) WRITE (KLOG,110) 110 FORMAT(//' Sample 1. Find a root of f(x) = x**5 - 3x**4 + ', * 'x**3 - 4x**2 + x - 6 = 0.'/// * ' Iteration Newton Approximation') CALL ZMFORM('F32.30','F32.30',MA,ST1) WRITE (KW,120) 0,ST1(1:69) WRITE (KLOG,120) 0,ST1(1:69) 120 FORMAT(/I6,4X,A) C DO 130 ITER = 1, 10 C C MC is f(MA). C CALL ZMEQ(MA,MC) CALL ZMADDI(MC,-3) CALL ZMMPY(MC,MA,MC) CALL ZMADDI(MC,1) CALL ZMMPY(MC,MA,MC) CALL ZMADDI(MC,-4) CALL ZMMPY(MC,MA,MC) CALL ZMADDI(MC,1) CALL ZMMPY(MC,MA,MC) CALL ZMADDI(MC,-6) C C MD is f'(MA). C CALL ZMMPYI(MA,5,MD) CALL ZMADDI(MD,-12) CALL ZMMPY(MD,MA,MD) CALL ZMADDI(MD,3) CALL ZMMPY(MD,MA,MD) CALL ZMADDI(MD,-8) CALL ZMMPY(MD,MA,MD) CALL ZMADDI(MD,1) C CALL ZMDIV(MC,MD,MB) CALL ZMSUB(MA,MB,MB) C C Print each iteration. C CALL ZMFORM('F32.30','F32.30',MB,ST1) WRITE (KW,120) ITER,ST1(1:69) WRITE (KLOG,120) ITER,ST1(1:69) C C Stop iterating if MA and MB agree to over C 30 places. C CALL ZMSUB(MA,MB,MD) CALL ZMABS(MD,MAFM) C C The ABS result is real -- do a real (FM) compare. C CALL FMST2M('1.0E-31',MBFM) IF (FMCOMP(MAFM,'LT',MBFM)) GO TO 140 C C Set MA = MB for the next iteration. C CALL ZMEQ(MB,MA) 130 CONTINUE C C Check the answer. C 140 ST1 = '0.561958308335403235498111195347453 +'// * '1.061134679604332556983391239058885 i' CALL ZMST2M(ST1,MC) CALL ZMSUB(MC,MB,MD) CALL ZMABS(MD,MAFM) CALL FMST2M('1.0E-31',MBFM) IF (FMCOMP(MAFM,'GT',MBFM)) THEN NERROR = NERROR + 1 WRITE (KW,150) WRITE (KLOG,150) 150 FORMAT(/' Error in sample case number 1.'/) ENDIF C C 2. Compute Exp(1.23-2.34i). C C Use the direct Taylor series. See the ZMEXP routine C for a faster way to get Exp(x). C C MA is x. C MB is the current term, x**n/n!. C MC is the current partial sum. C CALL ZMST2M('1.23-2.34i',MA) CALL ZMI2M(1,MB) CALL ZMEQ(MB,MC) DO 170 K = 1, 100 CALL ZMMPY(MB,MA,MB) CALL ZMDIVI(MB,K,MB) CALL ZMADD(MC,MB,MC) C C Test for convergence. KFLAG will be 1 if the result C of the last add or subtract is the same as one of the C input arguments. C IF (KFLAG.EQ.1) THEN WRITE (KW,160) K WRITE (KLOG,160) K 160 FORMAT(///' Sample 2.',8X,I5,' terms were added to get ', * 'Exp(1.23-2.34i)'/) GO TO 180 ENDIF 170 CONTINUE C C Print the result. C 180 CALL ZMFORM('F33.30','F32.30',MC,ST1) WRITE (KW,190) ST1(1:70) WRITE (KLOG,190) ST1(1:70) 190 FORMAT(' Result= ',A) C C Check the answer. C ST1 = '-2.379681796854777515745457977696745 -'// * '2.458032970832342652397461908326042 i' CALL ZMST2M(ST1,MD) CALL ZMSUB(MD,MC,MD) CALL ZMABS(MD,MAFM) CALL FMST2M('1.0E-31',MBFM) IF (FMCOMP(MAFM,'GT',MBFM)) THEN NERROR = NERROR + 1 WRITE (KW,200) WRITE (KLOG,200) 200 FORMAT(/' Error in sample case number 2.'/) ENDIF C IF (NERROR.EQ.0) THEN WRITE (KW,210) ' All results were ok.' WRITE (KLOG,210) ' All results were ok.' 210 FORMAT(//A/) ENDIF STOP END