PROGRAM SAMPLE C C David M. Smith 6-17-96 C C This is a test program for FMLIB 1.1, a multiple-precision real C arithmetic package. A few example FM calculations are carried C out using 60 significant digit precision. C C The output is saved in file FMSAMPLE.LOG. A comparison file, C FMSAMPLE.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 four 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 various array sizes. C C INTEGER NDIGMX,NBITS,LPACK,LUNPCK,LMWA,LJSUMS,LMBUFF 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 ) 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 FM variables. All are in C unpacked format. C DOUBLE PRECISION MA(0:LUNPCK),MB(0:LUNPCK),MC(0:LUNPCK), * MD(0:LUNPCK) C C Character strings used for input and output. C CHARACTER *80 ST1 INTEGER ITER,J,K,KLOG,NERROR LOGICAL FMCOMP,IMCOMP C C Set precision to give at least 60 significant digits C and initialize the FMLIB package. C C Note that any program using the FM package MUST call C FMSET before using the package. C CALL FMSET(60) C NERROR = 0 C C Write output to the standard FM output (unit KW, defined C in subroutine FMSET), and also to the file FMSAMPLE.LOG. C KLOG = 18 OPEN (KLOG,FILE='FMSAMPLE.LOG') C C 1. Find a root of the equation C f(x) = x**5 - 3x**4 + x**3 - 4x**2 + x - 6 = 0. C C Use Newton's method with initial guess x = 3.12. C This version is not tuned for speed. See the FMSQRT 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 FMST2M('3.12',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 FMFORM('F65.60',MA,ST1) WRITE (KW,120) 0,ST1(1:65) WRITE (KLOG,120) 0,ST1(1:65) 120 FORMAT(/I10,4X,A) C DO 130 ITER = 1, 10 C C MC is f(MA). C CALL FMEQ(MA,MC) CALL FMADDI(MC,-3) CALL FMMPY(MC,MA,MC) CALL FMADDI(MC,1) CALL FMMPY(MC,MA,MC) CALL FMADDI(MC,-4) CALL FMMPY(MC,MA,MC) CALL FMADDI(MC,1) CALL FMMPY(MC,MA,MC) CALL FMADDI(MC,-6) C C MD is f'(MA). C CALL FMMPYI(MA,5,MD) CALL FMADDI(MD,-12) CALL FMMPY(MD,MA,MD) CALL FMADDI(MD,3) CALL FMMPY(MD,MA,MD) CALL FMADDI(MD,-8) CALL FMMPY(MD,MA,MD) CALL FMADDI(MD,1) C CALL FMDIV(MC,MD,MB) CALL FMSUB(MA,MB,MB) C C Print each iteration. C CALL FMFORM('F65.60',MB,ST1) WRITE (KW,120) ITER,ST1(1:65) WRITE (KLOG,120) ITER,ST1(1:65) C C Stop iterating if MA and MB agree to over C 60 places. C CALL FMSUB(MA,MB,MD) CALL FMABS(MD,MD) CALL FMST2M('1.0E-61',MC) IF (FMCOMP(MD,'LT',MC)) GO TO 140 C C Set MA = MB for the next iteration. C CALL FMEQ(MB,MA) 130 CONTINUE C C Check the answer. C 140 ST1 = '3.120656215326726500470956013523797484654623935599066014'// * '9888284358' CALL FMST2M(ST1,MC) CALL FMSUB(MC,MB,MD) CALL FMABS(MD,MD) CALL FMST2M('1.0E-61',MC) IF (FMCOMP(MD,'GT',MC)) THEN NERROR = NERROR + 1 WRITE (KW,150) WRITE (KLOG,150) 150 FORMAT(/' Error in sample case number 1.'/) ENDIF C C 2. Compute the Riemann Zeta function for s=3. C C Use Gosper's formula Zeta(3) = C (5/4)*Sum[ (-1)**k * (k!)**2 / ((k+1)**2 * (2k+1)!) ] C while k = 0, 1, .... C C MA is the current partial sum. C MB is the current term. C MC is k! C MD is (2k+1)! C CALL FMI2M(1,MA) CALL FMEQ(MA,MC) CALL FMEQ(MA,MD) DO 170 K = 1, 200 CALL FMMPYI(MC,K,MC) J = 2*K*(2*K+1) CALL FMMPYI(MD,J,MD) CALL FMSQR(MC,MB) J = (K+1)*(K+1) CALL FMDIVI(MB,J,MB) CALL FMDIV(MB,MD,MB) IF (MOD(K,2).EQ.0) THEN CALL FMADD(MA,MB,MA) ELSE CALL FMSUB(MA,MB,MA) ENDIF 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'/) GO TO 180 ENDIF 170 CONTINUE C C Print the result. C 180 CALL FMMPYI(MA,5,MA) CALL FMDIVI(MA,4,MA) CALL FMFORM('F65.60',MA,ST1) WRITE (KW,190) ST1(1:65) WRITE (KLOG,190) ST1(1:65) 190 FORMAT(' Zeta(3) = ',A) C C Check the answer. C ST1 = '1.20205690315959428539973816151144999076498629234049888'// * '1792271555' CALL FMST2M(ST1,MC) CALL FMSUB(MA,MC,MD) CALL FMABS(MD,MD) CALL FMST2M('1.0E-61',MC) IF (FMCOMP(MD,'GT',MC)) THEN NERROR = NERROR + 1 WRITE (KW,200) WRITE (KLOG,200) 200 FORMAT(/' Error in sample case number 2.'/) ENDIF C C 3. Integer multiple precision calculations. C C Fermat's theorem says x**(p-1) mod p = 1 C when p is prime and x is not a multiple of p. C If x**(p-1) mod p gives 1 for some p with C several different x's, then it is very likely C that p is prime (but it is not certain until C further tests are done). C C Find a 70-digit number p that is "probably" prime. C C MA is the value p being tested. C CALL IMI2M(10,MA) CALL IMI2M(69,MB) CALL IMPWR(MA,MB,MA) C C To speed up the search, test only values that are C not multiples of 2, 3, 5, 7, 11, 13. C K = 2*3*5*7*11*13 CALL IMDIVI(MA,K,MA) CALL IMMPYI(MA,K,MA) CALL IMI2M(K,MB) CALL IMADD(MA,MB,MA) CALL IMI2M(1,MD) CALL IMADD(MA,MD,MA) CALL IMI2M(3,MC) C DO 220 J = 1, 100 C C Compute 3**(p-1) mod p C CALL IMSUB(MA,MD,MB) CALL IMPMOD(MC,MB,MA,MC) IF (IMCOMP(MC,'EQ',MD)) THEN C C Check that 7**(p-1) mod p is also 1. C CALL IMI2M(7,MC) CALL IMPMOD(MC,MB,MA,MC) IF (IMCOMP(MC,'EQ',MD)) THEN WRITE (KW,210) J WRITE (KLOG,210) J 210 FORMAT(///' Sample 3.',8X,I5,' values were tested'/) GO TO 230 ENDIF ENDIF C CALL IMI2M(3,MC) CALL IMI2M(K,MB) CALL IMADD(MA,MB,MA) 220 CONTINUE C C Print the result. C 230 CALL IMFORM('I72',MA,ST1) WRITE (KW,240) ST1(1:72) WRITE (KLOG,240) ST1(1:72) 240 FORMAT(' p = ',A) C C Check the answer. C ST1 = '1000000000000000000000000000000000000000000000000000'// * '000000000000659661' CALL IMST2M(ST1,MC) IF (IMCOMP(MA,'NE',MC)) THEN NERROR = NERROR + 1 WRITE (KW,250) WRITE (KLOG,250) 250 FORMAT(/' Error in sample case number 3.'/) ENDIF C IF (NERROR.EQ.0) THEN WRITE (KW,260) ' All results were ok.' WRITE (KLOG,260) ' All results were ok.' 260 FORMAT(//A/) ENDIF STOP END