C ALGORITHM 786, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 24,NO. 4, December, 1998, P. 359--367. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Fortran90/ # Fortran90/Doc/ # Fortran90/Doc/Makefile # Fortran90/Doc/readme # Fortran90/Drivers/ # Fortran90/Drivers/Dp/ # Fortran90/Drivers/Dp/RES1 # Fortran90/Drivers/Dp/RES2 # Fortran90/Drivers/Dp/RES3 # Fortran90/Drivers/Dp/RES4 # Fortran90/Drivers/Dp/RES5 # Fortran90/Drivers/Dp/RES6 # Fortran90/Drivers/Dp/driver1.f90 # Fortran90/Drivers/Dp/driver2.f90 # Fortran90/Drivers/Dp/driver3.f90 # Fortran90/Drivers/Dp/driver4.f90 # Fortran90/Drivers/Dp/driver5.f90 # Fortran90/Drivers/Dp/driver6.f90 # Fortran90/Src/ # Fortran90/Src/Dp/ # Fortran90/Src/Dp/fmlib.f90 # Fortran90/Src/Dp/fmzm90.f90 # Fortran90/Src/Dp/fmzmcomm.f90 # Fortran90/Src/Dp/zmlib.f90 # This archive created: Thu Mar 25 10:55:09 1999 export PATH; PATH=/bin:$PATH if test ! -d 'Fortran90' then mkdir 'Fortran90' fi cd 'Fortran90' if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'Makefile' then echo shar: will not over-write existing file "'Makefile'" else cat << SHAR_EOF > 'Makefile' # Define EPC F90 compiler and flags #FC = epcf90 #FFLAGS = -C -d1 -g -temp=/tmp -u #FFLAGS = -temp=/tmp -O # Define Nag f90 compiler and flags FC = f90 FFLAGS = -g # Define rule for .f to .o and .f90 to .o .SUFFIXES : .f .f90 .o .f.o: $(FC) $(FFLAGS) -c $< .f90.o: $(FC) $(FFLAGS) -c $< all: res1 res2 res3 res4 res5 res6 res1: zmlib.o fmlib.o driver1.o $(FC) $(FFLAGS) zmlib.o fmlib.o driver1.o -o driver1 driver1 > res1 res2: zmlib.o fmlib.o driver2.o $(FC) $(FFLAGS) zmlib.o fmlib.o driver2.o -o driver2 driver2 > res2 res3: fmlib.o driver3.o $(FC) $(FFLAGS) fmlib.o driver3.o -o driver3 driver3 > res3 res4: driver4.o fmlib.o $(FC) $(FFLAGS) driver4.o fmlib.o -o driver4 driver4 > res4 res5: fmzmcomm.o fmzm90.o fmlib.o zmlib.o driver5.o $(FC) $(FFLAGS) fmzmcomm.o fmzm90.o fmlib.o zmlib.o driver5.o -o driver5 driver5 > res5 res6: fmzmcomm.o fmzm90.o fmlib.o zmlib.o driver6.o $(FC) $(FFLAGS) fmzmcomm.o fmzm90.o fmlib.o zmlib.o driver6.o -o driver6 driver6 > res6 clean: rm -rf driver4 driver6 driver3 driver5 driver1 driver2 rm -rf *.o *.LOG res* SHAR_EOF fi # end of overwriting check if test -f 'readme' then echo shar: will not over-write existing file "'readme'" else cat << SHAR_EOF > 'readme' This is a list of the files for version 1.1 of FMLIB and ZMLIB. 1. zmlib.f90 Subroutine library for complex operations 2. fmlib.f90 Subroutine library for real operations 3. testzm.f90 Test program for most of the ZM routines 4. zmsample.f90 Small sample program using ZM 5. zmsample.chk Expected output file from zmsample.f90 6. testfm.f90 Test program for most of the FM routines 7. fmsample.f90 Small sample program using FM 8. fmsample.chk Expected output file from fmsample.f90 9. fmzm90.f90 Fortran-90 interface module 10. fmzmcomm.f90 Fortran-90 module for common blocks 11. Test90.f90 Test program for fmzm90 12. Sample90.f90 Small sample program using fmzm90 13. SAMPLE90.CHK Expected output file from sample90.f A makefile detailing the building of the ZM and FM libraries and the running of all the test examples is provided. My web site contains copies of other related papers and files. In 1998 it was located at "http://cse.eng.lmu.edu/~dsmith/FMLIB.html". If that location changes in the future, try searching for the keyword "dsmithfmlibrary" to find the site. =========================================================================== =========================================================================== USER'S GUIDE FOR THE FM PACKAGE The various lists of available multiple precision operations and routines have been collected here, along with some general advice on using the package. See the programs fmsample.f90, zmsample.f90, and sample90.f90 for some examples of initializing and using the package. This version of the package uses code with the names of routines, variables, and files in lower case, but in this file as well as in comment lines in the code such names are emphasized by writing them in upper case. INITIALIZATION: Before ANY part of the FM package can be used, the base and precision to be used must be defined, along with several other saved parameters. If any complex arithmetic is to be used, put CALL ZMSET(N) in the main program before any multiple precision operations are done, with N replaced by the number of decimal digits of accuracy to be used. This will initialize both FMLIB and ZMLIB packages. If only real arithmetic is to be used, put CALL FMSET(N) in the main program before any multiple precision operations are done, with N replaced by the number of decimal digits of accuracy to be used. This will initialize the FMLIB package. One of these calls must be present whether the FM/ZM routines are to be called directly by the user, or the Fortran-90 interface routines are to be used. For compatibility when the interface module is used, the derived type routine names FM_SET or ZM_SET may be used in place of FMSET or ZMSET. MODULE/COMMON: Some common blocks used for saved parameters must be declared in the main program. If the Fortran-90 interface is used, put USE FMZM at the beginning of the main program and also in each routine that uses type FM, IM, or ZM variables. If the Fortran-90 interface is not used, put the common blocks given in zmsample.f90 at the top of the main program if complex arithmetic is used, or put the common blocks given in fmsample.f90 at the top of the main program if only real arithmetic is used. ROUTINE NAMES: For each multiple precision operation there are several routines with related names that perform variations of that operation. For example, the addition operation has these forms: Using the Fortran-90 interface module, to perform real (floating-point) multiple precision addition, declare the variables with TYPE ( FM ) A,B,C and then add using C = A + B Normally, using the interface module avoids the need to know the name of the FM routine being called. For some operations, usually those that are not Fortran-90 functions (such as formatting a number), a direct call may be needed. The addition above can be done as CALL FM_ADD(A,B,C) If fmlib.f90 is used without the interface module, then the multiple precision numbers are declared as arrays DOUBLE PRECISION A(0:LUNPCK),B(0:LUNPCK),C(0:LUNPCK) where LUNPCK is defined in the PARAMETER statement included with the FM common blocks. The numbers are then added by calling the FMLIB routine where the arguments are assumed to be arrays, not TYPE (FM) derived types: CALL FMADD(A,B,C) For each of the routines listed below (like FMADD), there is a version that assumes the arguments have the appropriate derived type. These have the same names, except "_" has been inserted after the first two letters of the name (like FM_ADD). If direct calls are done instead of using the interface module, there is another form for these routine names that is used when the arrays have been declared in a packed format that takes roughly half as much space: DOUBLE PRECISION A(0:LPACK),B(0:LPACK),C(0:LPACK) The routines that work with packed arrays have names where the second letter has been changed from M to P: CALL FPADD(A,B,C) The packed versions are slower. For multiple precision integer or complex operations there are similar Fortran-90 derived types and the various routines: USE FMZM ... TYPE ( IM ) A,B,C TYPE ( ZM ) X,Y,Z ... C = A + B ... Z = X + Y with explicit calls of the form CALL IM_ADD(A,B,C) CALL ZM_ADD(X,Y,Z) Without using the interface module: DOUBLE PRECISION A(0:LUNPCK),B(0:LUNPCK),C(0:LUNPCK) DOUBLE PRECISION X(0:LUNPKZ),Y(0:LUNPKZ),Z(0:LUNPKZ) ... CALL IMADD(A,B,C) ... CALL ZMADD(X,Y,Z) Packed format without the interface module: DOUBLE PRECISION A(0:LPACK),B(0:LPACK),C(0:LPACK) DOUBLE PRECISION X(0:LPACKZ),Y(0:LPACKZ),Z(0:LPACKZ) ... CALL IPADD(A,B,C) ... CALL ZPADD(X,Y,Z) ------------------------------------------------------------------------ ------------------- Fortran-90 Interface Notes --------------------- There are three multiple precision data types: FM (multiple precision real) IM (multiple precision integer) ZM (multiple precision complex) Some the the interface routines assume that the precision chosen in the calling program (using FM_SET or ZM_SET) represents more significant digits than does the machine's double precision. All the functions defined in this module are standard Fortran-90 functions, except for several direct conversion functions: TO_FM is a function for converting other types of numbers to type FM. Note that TO_FM(3.12) converts the REAL constant to FM, but it is accurate only to single precision. TO_FM(3.12D0) agrees with 3.12 to double precision accuracy, and TO_FM('3.12') or TO_FM(312)/TO_FM(100) agrees to full FM accuracy. TO_IM converts to type IM, and TO_ZM converts to type ZM. Functions are also supplied for converting the three multiple precision types to the other numeric data types: TO_INT converts to machine precision integer TO_SP converts to single precision TO_DP converts to double precision TO_SPZ converts to single precision complex TO_DPZ converts to double precision complex WARNING: When multiple precision type declarations are inserted in an existing program, take care in converting functions like DBLE(X), where X has been declared as a multiple precision type. If X was single precision in the original program, then replacing the DBLE(X) by TO_DP(X) in the new version could lose accuracy. For this reason, the Fortran type-conversion functions defined in this module assume that results should be multiple precision whenever inputs are. Examples: DBLE(TO_FM('1.23E+123456')) is type FM REAL(TO_FM('1.23E+123456')) is type FM REAL(TO_ZM('3.12+4.56i')) is type FM = TO_FM('3.12') INT(TO_FM('1.23')) is type IM = TO_IM(1) INT(TO_IM('1E+23')) is type IM CMPLX(TO_FM('1.23'),TO_FM('4.56')) is type ZM Programs using this module may sometimes need to call FM, IM, or ZM routines directly. This is normally the case when routines are needed that are not Fortran-90 intrinsics, such as the formatting subroutine FMFORM. In a program using this module, suppose MAFM has been declared with TYPE ( FM ) MAFM. To use the routine FMFORM, which expects the second argument to be an array and not a derived type, the call would have to be CALL FMFORM('F65.60',MAFM%MFM,ST1) so that the array contained in MAFM is passed. As an alternative so the user can refer directly to the FM-, IM-, and ZM-type variables and avoid the cumbersome "%MFM" suffixes, this module contains a collection of interface routines to supply any needed argument conversions. For each FM, IM, and ZM routine that is designed to be called by the user, there is also a version that assumes any multiple-precision arguments are derived types instead of arrays. Each interface routine has the same name as the original with an underscore after the first two letters of the routine name. To convert the number to a character string with F65.60 format, use CALL FM_FORM('F65.60',MAFM,ST1) if MAFM is of TYPE ( FM ), or use CALL FMFORM('F65.60',MA,ST1) if MA is declared as an array. All the routines shown below may be used this way. For each of the operations =, +, -, *, /, **, .EQ., .NE., .GT., .GE., .LT., and .LE., the interface module defines all mixed mode variations involving one of the three multiple precision derived types and another argument having one of the types: { integer, real, double, complex, complex double, FM, IM, ZM }. So mixed mode expressions such as MAFM = 12 MAFM = MAFM + 1 IF (ABS(MAFM).LT.1.0D-23) THEN are handled correctly. Not all the named functions are defined for all three multiple precision derived types, so the list below shows which can be used. The labels "real", "integer", and "complex" refer to types FM, IM, and ZM respectively, "string" means the function accepts character strings (e.g., TO_FM('3.45')), and "other" means the function can accept any of the machine precision data types integer, real, double, complex, or complex double. For functions that accept two or more arguments, like ATAN2 or MAX, all the arguments must be of the same type. AVAILABLE OPERATIONS: = + - * / ** .EQ. .NE. .GT. .GE. .LT. .LE. ABS real integer complex ACOS real complex AIMAG complex AINT real complex ANINT real complex ASIN real complex ATAN real complex ATAN2 real BTEST integer CEILING real complex CMPLX real integer CONJ complex COS real complex COSH real complex DBLE real integer complex DIGITS real integer complex DIM real integer DINT real complex DOTPRODUCT real integer complex EPSILON real EXP real complex EXPONENT real FLOOR real integer complex FRACTION real complex HUGE real integer complex INT real integer complex LOG real complex LOG10 real complex MATMUL real integer complex MAX real integer MAXEXPONENT real MIN real integer MINEXPONENT real MOD real integer MODULO real integer NEAREST real NINT real integer complex PRECISION real complex RADIX real integer complex RANGE real integer complex REAL real integer complex RRSPACING real SCALE real complex SETEXPONENT real SIGN real integer SIN real complex SINH real complex SPACING real SQRT real complex TAN real complex TANH real complex TINY real integer complex TO_FM real integer complex string other TO_IM real integer complex string other TO_ZM real integer complex string other TO_INT real integer complex TO_SP real integer complex TO_DP real integer complex TO_SPZ real integer complex TO_DPZ real integer complex ------------------------------------------------------------------------ ----------- Routines for Real Floating-Point Operations ------------ These are the FM routines that are designed to be called by the user. All are subroutines except logical function FMCOMP. MA, MB, MC refer to FM format numbers. In each case it is permissible to use the same array more than once in the calling sequence. The statement MA = MA*MA can be written CALL FMMPY(MA,MA,MA). For each of these routines there is also a version available for which the argument list is the same but all FM numbers are in packed format. The routines using packed numbers have the same names except 'FM' is replaced by 'FP' at the start of each name. FMABS(MA,MB) MB = ABS(MA) FMACOS(MA,MB) MB = ACOS(MA) FMADD(MA,MB,MC) MC = MA + MB FMADDI(MA,IVAL) MA = MA + IVAL Increment an FM number by a one word integer. Note this call does not have an "MB" result like FMDIVI and FMMPYI. FMASIN(MA,MB) MB = ASIN(MA) FMATAN(MA,MB) MB = ATAN(MA) FMATN2(MA,MB,MC) MC = ATAN2(MA,MB) FMBIG(MA) MA = Biggest FM number less than overflow. FMCHSH(MA,MB,MC) MB = COSH(MA), MC = SINH(MA). Faster than making two separate calls. FMCOMP(MA,LREL,MB) Logical comparison of MA and MB. LREL is a CHARACTER*2 value identifying which comparison is made. Example: IF (FMCOMP(MA,'GE',MB)) ... FMCONS Set several saved constants that depend on MBASE, the base being used. FMCONS should be called immediately after changing MBASE. FMCOS(MA,MB) MB = COS(MA) FMCOSH(MA,MB) MB = COSH(MA) FMCSSN(MA,MB,MC) MB = COS(MA), MC = SIN(MA). Faster than making two separate calls. FMDIG(NSTACK,KST) Find a set of precisions to use during Newton iteration for finding a simple root starting with about double precision accuracy. FMDIM(MA,MB,MC) MC = DIM(MA,MB) FMDIV(MA,MB,MC) MC = MA/MB FMDIVI(MA,IVAL,MB) MB = MA/IVAL IVAL is a one word integer. FMDP2M(X,MA) MA = X Convert from double precision to FM. FMDPM(X,MA) MA = X Convert from double precision to FM. Much faster than FMDP2M, but MA agrees with X only to D.P. accuracy. See the comments in the two routines. FMEQ(MA,MB) MB = MA Both have precision NDIG. This is the version to use for standard B = A statements. FMEQU(MA,MB,NA,NB) MB = MA Version for changing precision. MA has NA digits (i.e., MA was computed using NDIG = NA), and MB will be defined having NB digits. MB is zero-padded if NB.GT.NA MB is rounded if NB.LT.NA FMEXP(MA,MB) MB = EXP(MA) FMFORM(FORM,MA,STRING) MA is converted to a character string using format FORM and returned in STRING. FORM can represent I, F, E, or 1PE formats. Example: CALL FMFORM('F60.40',MA,STRING) FMFPRT(FORM,MA) Print MA on unit KW using FORM format. FMI2M(IVAL,MA) MA = IVAL Convert from one word integer to FM. FMINP(LINE,MA,LA,LB) MA = LINE Input conversion. Convert LINE(LA) through LINE(LB) from characters to FM. FMINT(MA,MB) MB = INT(MA) Integer part of MA. FMIPWR(MA,IVAL,MB) MB = MA**IVAL Raise an FM number to a one word integer power. FMLG10(MA,MB) MB = LOG10(MA) FMLN(MA,MB) MB = LOG(MA) FMLNI(IVAL,MA) MA = LOG(IVAL) Natural log of a one word integer. FMM2DP(MA,X) X = MA Convert from FM to double precision. FMM2I(MA,IVAL) IVAL = MA Convert from FM to integer. FMM2SP(MA,X) X = MA Convert from FM to single precision. FMMAX(MA,MB,MC) MC = MAX(MA,MB) FMMIN(MA,MB,MC) MC = MIN(MA,MB) FMMOD(MA,MB,MC) MC = MA mod MB FMMPY(MA,MB,MC) MC = MA*MB FMMPYI(MA,IVAL,MB) MB = MA*IVAL Multiply by a one word integer. FMNINT(MA,MB) MB = NINT(MA) Nearest FM integer. FMOUT(MA,LINE,LB) LINE = MA Convert from FM to character. LINE is a character array of length LB. FMPI(MA) MA = pi FMPRNT(MA) Print MA on unit KW using current format. FMPWR(MA,MB,MC) MC = MA**MB FMREAD(KREAD,MA) MA is returned after reading one (possibly multi-line) FM number on unit KREAD. This routine reads numbers written by FMWRIT. FMRPWR(MA,K,J,MB) MB = MA**(K/J) Rational power. Faster than FMPWR for functions like the cube root. FMSET(NPREC) Set default values and machine-dependent variables to give at least NPREC base 10 digits plus three base 10 guard digits. Must be called to initialize FM package. FMSIGN(MA,MB,MC) MC = SIGN(MA,MB) Sign transfer. FMSIN(MA,MB) MB = SIN(MA) FMSINH(MA,MB) MB = SINH(MA) FMSP2M(X,MA) MA = X Convert from single precision to FM. FMSQR(MA,MB) MB = MA*MA Faster than FMMPY. FMSQRT(MA,MB) MB = SQRT(MA) FMST2M(STRING,MA) MA = STRING Convert from character string to FM. Often more convenient than FMINP, which converts an array of CHARACTER*1 values. Example: CALL FMST2M('123.4',MA). FMSUB(MA,MB,MC) MC = MA - MB FMTAN(MA,MB) MB = TAN(MA) FMTANH(MA,MB) MB = TANH(MA) FMULP(MA,MB) MB = One Unit in the Last Place of MA. FMWRIT(KWRITE,MA) Write MA on unit KWRITE. Multi-line numbers will have '&' as the last nonblank character on all but the last line. These numbers can then be read easily using FMREAD. ------------------------------------------------------------------------ ----------------- Routines for Integer Operations ------------------ These are the integer routines that are designed to be called by the user. All are subroutines except logical function IMCOMP. MA, MB, MC refer to IM format numbers. In each case the version of the routine to handle packed IM numbers has the same name, with 'IM' replaced by 'IP'. IMABS(MA,MB) MB = ABS(MA) IMADD(MA,MB,MC) MC = MA + MB IMBIG(MA) MA = Biggest IM number less than overflow. IMCOMP(MA,LREL,MB) Logical comparison of MA and MB. LREL is a CHARACTER*2 value identifying which comparison is made. Example: IF (IMCOMP(MA,'GE',MB)) ... IMDIM(MA,MB,MC) MC = DIM(MA,MB) IMDIV(MA,MB,MC) MC = int(MA/MB) Use IMDIVR if the remainder is also needed. IMDIVI(MA,IVAL,MB) MB = int(MA/IVAL) IVAL is a one word integer. Use IMDVIR to get the remainder also. IMDIVR(MA,MB,MC,MD) MC = int(MA/MB), MD = MA mod MB When both the quotient and remainder are needed, this routine is twice as fast as calling both IMDIV and IMMOD. IMDVIR(MA,IVAL,MB,IREM) MB = int(MA/IVAL), IREM = MA mod IVAL IVAL and IREM are one word integers. IMEQ(MA,MB) MB = MA IMFM2I(MAFM,MB) MB = MAFM Convert from real (FM) format to integer (IM) format. IMFORM(FORM,MA,STRING) MA is converted to a character string using format FORM and returned in STRING. FORM can represent I, F, E, or 1PE formats. Example: CALL IMFORM('I70',MA,STRING) IMFPRT(FORM,MA) Print MA on unit KW using FORM format. IMGCD(MA,MB,MC) MC = greatest common divisor of MA and MB. IMI2FM(MA,MBFM) MBFM = MA Convert from integer (IM) format to real (FM) format. IMI2M(IVAL,MA) MA = IVAL Convert from one word integer to IM. IMINP(LINE,MA,LA,LB) MA = LINE Input conversion. Convert LINE(LA) through LINE(LB) from characters to IM. IMM2DP(MA,X) X = MA Convert from IM to double precision. IMM2I(MA,IVAL) IVAL = MA Convert from IM to one word integer. IMMAX(MA,MB,MC) MC = MAX(MA,MB) IMMIN(MA,MB,MC) MC = MIN(MA,MB) IMMOD(MA,MB,MC) MC = MA mod MB IMMPY(MA,MB,MC) MC = MA*MB IMMPYI(MA,IVAL,MB) MB = MA*IVAL Multiply by a one word integer. IMMPYM(MA,MB,MC,MD) MD = MA*MB mod MC Slightly faster than calling IMMPY and IMMOD separately, and it works for cases where IMMPY would return OVERFLOW. IMOUT(MA,LINE,LB) LINE = MA Convert from IM to character. LINE is a character array of length LB. IMPMOD(MA,MB,MC,MD) MD = MA**MB mod MC IMPRNT(MA) Print MA on unit KW. IMPWR(MA,MB,MC) MC = MA**MB IMREAD(KREAD,MA) MA is returned after reading one (possibly multi-line) IM number on unit KREAD. This routine reads numbers written by IMWRIT. IMSIGN(MA,MB,MC) MC = SIGN(MA,MB) Sign transfer. IMSQR(MA,MB) MB = MA*MA Faster than IMMPY. IMST2M(STRING,MA) MA = STRING Convert from character string to IM. Often more convenient than IMINP, which converts an array of CHARACTER*1 values. Example: CALL IMST2M('12345678901',MA). IMSUB(MA,MB,MC) MC = MA - MB IMWRIT(KWRITE,MA) Write MA on unit KWRITE. Multi-line numbers will have '&' as the last nonblank character on all but the last line. These numbers can then be read easily using IMREAD. Many of the IM routines call FM routines, but none of the FM routines call IM routines, so the IM routines can be omitted if none are called explicitly from a program. ------------------------------------------------------------------------ ---------- Routines for Complex Floating-Point Operations ---------- These are the routines in ZMLIB that are designed to be called by the user. All are subroutines, and in each case the version of the routine to handle packed ZM numbers has the same name, with 'ZM' replaced by 'ZP'. MA, MB, MC refer to ZM format complex numbers. MAFM, MBFM, MCFM refer to FM format real numbers. INTEG is a Fortran INTEGER variable. ZVAL is a Fortran COMPLEX variable. In each case it is permissible to use the same array more than once in the calling sequence. The statement MA = MA*MA may be written CALL ZMMPY(MA,MA,MA). ZMABS(MA,MBFM) MBFM = ABS(MA) Result is real. ZMACOS(MA,MB) MB = ACOS(MA) ZMADD(MA,MB,MC) MC = MA + MB ZMADDI(MA,INTEG) MA = MA + INTEG Increment an ZM number by a one word integer. Note this call does not have an "MB" result like ZMDIVI and ZMMPYI. ZMARG(MA,MBFM) MBFM = Argument(MA) Result is real. ZMASIN(MA,MB) MB = ASIN(MA) ZMATAN(MA,MB) MB = ATAN(MA) ZMCHSH(MA,MB,MC) MB = COSH(MA), MC = SINH(MA). Faster than 2 calls. ZMCMPX(MAFM,MBFM,MC) MC = CMPLX(MAFM,MBFM) ZMCONJ(MA,MB) MB = CONJG(MA) ZMCOS(MA,MB) MB = COS(MA) ZMCOSH(MA,MB) MB = COSH(MA) ZMCSSN(MA,MB,MC) MB = COS(MA), MC = SIN(MA). Faster than 2 calls. ZMDIV(MA,MB,MC) MC = MA / MB ZMDIVI(MA,INTEG,MB) MB = MA / INTEG ZMEQ(MA,MB) MB = MA ZMEQU(MA,MB,NDA,NDB) MB = MA Version for changing precision. (NDA and NDB are as in FMEQU) ZMEXP(MA,MB) MB = EXP(MA) ZMFORM(FORM1,FORM2,MA,STRING) STRING = MA MA is converted to a character string using format FORM1 for the real part and FORM2 for the imaginary part. The result is returned in STRING. FORM1 and FORM2 can represent I, F, E, or 1PE formats. Example: CALL ZMFORM('F20.10','F15.10',MA,STRING) ZMFPRT(FORM1,FORM2,MA) Print MA on unit KW using formats FORM1 and FORM2. ZMI2M(INTEG,MA) MA = CMPLX(INTEG,0) ZM2I2M(INTEG1,INTEG2,MA) MA = CMPLX(INTEG1,INTEG2) ZMIMAG(MA,MBFM) MBFM = IMAG(MA) Imaginary part. ZMINP(LINE,MA,LA,LB) MA = LINE Input conversion. Convert LINE(LA) through LINE(LB) from characters to ZM. LINE is a character array of length at least LB. ZMINT(MA,MB) MB = INT(MA) Integer part of both Real and Imaginary parts of MA. ZMIPWR(MA,INTEG,MB) MB = MA ** INTEG Integer power function. ZMLG10(MA,MB) MB = LOG10(MA) ZMLN(MA,MB) MB = LOG(MA) ZMM2I(MA,INTEG) INTEG = INT(REAL(MA)) ZMM2Z(MA,ZVAL) ZVAL = MA ZMMPY(MA,MB,MC) MC = MA * MB ZMMPYI(MA,INTEG,MB) MB = MA * INTEG ZMNINT(MA,MB) MB = NINT(MA) Nearest integer of both Real and Imaginary. ZMOUT(MA,LINE,LB,LAST1,LAST2) LINE = MA Convert from FM to character. LINE is the returned character array. LB is the dimensioned size of LINE. LAST1 is returned as the position in LINE of the last character of REAL(MA). LAST2 is returned as the position in LINE of the last character of AIMAG(MA). ZMPRNT(MA) Print MA on unit KW using current format. ZMPWR(MA,MB,MC) MC = MA ** MB ZMREAD(KREAD,MA) MA is returned after reading one (possibly multi-line) ZM number on unit KREAD. This routine reads numbers written by ZMWRIT. ZMREAL(MA,MBFM) MBFM = REAL(MA) Real part. ZMRPWR(MA,IVAL,JVAL,MB) MB = MA ** (IVAL/JVAL) ZMSET(NPREC) Initialize ZM package. Set precision to the equivalent of at least NPREC base 10 digits. ZMSIN(MA,MB) MB = SIN(MA) ZMSINH(MA,MB) MB = SINH(MA) ZMSQR(MA,MB) MB = MA*MA Faster than ZMMPY. ZMSQRT(MA,MB) MB = SQRT(MA) ZMST2M(STRING,MA) MA = STRING Convert from character string to ZM. Often more convenient than ZMINP, which converts an array of CHARACTER*1 values. Example: CALL ZMST2M('123.4+5.67i',MA). ZMSUB(MA,MB,MC) MC = MA - MB ZMTAN(MA,MB) MB = TAN(MA) ZMTANH(MA,MB) MB = TANH(MA) ZMWRIT(KWRITE,MA) Write MA on unit KWRITE. Multi-line numbers are formatted for automatic reading with ZMREAD. ZMZ2M(ZVAL,MA) MA = ZVAL ------------------------------------------------------------------------ -------------------------- fmlib.f90 Notes --------------------------- The FM routines in this package perform floating-point multiple-precision arithmetic, and the IM routines perform integer multiple-precision arithmetic. 1. INITIALIZING THE PACKAGE Before calling any routine in the package, several variables in the common blocks /FMUSER/, /FM/, /FMBUFF/, and /FMSAVE/ must be initialized. These four common blocks contain information that is saved between calls, so they should be declared in the main program. Subroutine FMSET initializes these variables to default values and defines all machine-dependent values in the package. After calling FMSET once at the start of a program, the user may sometimes want to reset some of the variables in these common blocks. These variables are described below. 2. REPRESENTATION OF FM NUMBERS MBASE is the base in which the arithmetic is done. MBASE must be bigger than one, and less than or equal to the square root of the largest representable integer. For best efficiency MBASE should be large, but no more than about 1/4 of the square root of the largest representable integer. Input and output conversions are much faster when MBASE is a power of ten. NDIG is the number of base MBASE digits that are carried in the multiple precision numbers. NDIG must be at least two. The upper limit for NDIG is defined in the PARAMETER statement at the top of each routine and is restricted only by the amount of memory available. Sometimes it is useful to dynamically vary NDIG during the program. Use FMEQU to round numbers to lower precision or zero-pad them to higher precision when changing NDIG. It is rare to need to change MBASE during a program. Use FMCONS to reset some saved constants that depend on MBASE. FMCONS should be called immediately after changing MBASE. There are two representations for a floating multiple precision number. The unpacked representation used by the routines while doing the computations is base MBASE and is stored in NDIG+2 words. A packed representation is available to store the numbers in the user's program in compressed form. In this format, the NDIG (base MBASE) digits of the mantissa are packed two per word to conserve storage. Thus the external, packed form of a number requires (NDIG+1)/2+2 words. This version uses double precision arrays to hold the numbers. Version 1.0 of FM used integer arrays, which are faster on some machines. The package can easily be changed to use integer arrays -- see section 11 on EFFICIENCY below. The unpacked format of a floating multiple precision number is as follows. A number MA is kept in an array with MA(1) containing the exponent and MA(2) through MA(NDIG+1) containing one digit of the mantissa, expressed in base MBASE. The array is dimensioned to start at MA(0), with the approximate number of bits of precision stored in MA(0). This precision value is intended to be used by FM functions that need to monitor cancellation error in addition and subtraction. The cancellation monitor code is usually disabled for user calls, and FM functions only check for cancellation when they must. Tracking cancellation causes most routines to run slower, with addition and subtraction being affected the most. The exponent is a power of MBASE and the implied radix point is immediately before the first digit of the mantissa. Every nonzero number is normalized so that the second array element (the first digit of the mantissa) is nonzero. In both representations the sign of the number is carried on the second array element only. Elements 3,4,... are always nonnegative. The exponent is a signed integer and may be as large in magnitude as MXEXP (defined in FMSET). For MBASE = 10,000 and NDIG = 4, the number -pi would have these representations: Word 1 2 3 4 5 Unpacked: 1 -3 1415 9265 3590 Packed: 1 -31415 92653590 Word 0 would be 42 in both formats, indicating that the mantissa has about 42 bits of precision. Because of normalization in a large base, the equivalent number of base 10 significant digits for an FM number may be as small as LOG10(MBASE)*(NDIG-1) + 1. The integer routines use the FMLIB format to represent numbers, without the number of digits (NDIG) being fixed. Integers in IM format are essentially variable precision, using the minimum number of words to represent each value. For programs using both FM and IM numbers, FM routines should not be called with IM numbers, and IM routines should not be called with FM numbers, since the implied value of NDIG used for an IM number may not match the explicit NDIG expected by an FM routine. Use the conversion routines IMFM2I and IMI2FM to change between the FM and IM formats. 3. INPUT/OUTPUT ROUTINES All versions of the input routines perform free-format conversion from characters to FM numbers. a. Conversion to or from a character array FMINP converts from a character*1 array to an FM number. FMOUT converts an FM number to base 10 and formats it for output as an array of type character*1. The output is left justified in the array, and the format is defined by two variables in common, so that a separate format definition does not have to be provided for each output call. The user sets JFORM1 and JFORM2 to determine the output format. JFORM1 = 0 E format ( .314159M+6 ) = 1 1PE format ( 3.14159M+5 ) = 2 F format ( 314159.000 ) JFORM2 is the number of significant digits to display (if JFORM1 = 0 or 1). If JFORM2.EQ.0 then a default number of digits is chosen. The default is roughly the full precision of the number. JFORM2 is the number of digits after the decimal point (if JFORM1 = 2). See the FMOUT documentation for more details. b. Conversion to or from a character string FMST2M converts from a character string to an FM number. FMFORM converts an FM number to a character string according to a format provided in each call. The format description is more like that of a Fortran FORMAT statement, and integer or fixed-point output is right justified. c. Direct read or write FMPRNT uses FMOUT to print one FM number. FMFPRT uses FMFORM to print one FM number. FMWRIT writes FM numbers for later input using FMREAD. FMREAD reads FM numbers written by FMWRIT. The values given to JFORM1 and JFORM2 can be used to define a default output format when FMOUT or FMPRNT are called. The explicit format used in a call to FMFORM or FMFPRT overrides the settings of JFORM1 and JFORM2. KW is the unit number to be used for standard output from the package, including error and warning messages, and trace output. For multiple precision integers, the corresponding routines IMINP, IMOUT, IMST2M, IMFORM, IMPRNT, IMFPRT, IMWRIT, and IMREAD provide similar input and output conversions. For output of IM numbers, JFORM1 and JFORM2 are ignored and integer format (JFORM1=2, JFORM2=0) is used. For further description of these routines, see sections 9 and 10 below. 4. ARITHMETIC TRACING NTRACE and LVLTRC control trace printout from the package. NTRACE = 0 No printout except warnings and errors. = 1 The result of each call to one of the routines is printed in base 10, using FMOUT. = -1 The result of each call to one of the routines is printed in internal base MBASE format. = 2 The input arguments and result of each call to one of the routines is printed in base 10, using FMOUT. = -2 The input arguments and result of each call to one of the routines is printed in base MBASE format. LVLTRC defines the call level to which the trace is done. LVLTRC = 1 means only FM routines called directly by the user are traced, LVLTRC = 2 also prints traces for FM routines called by other FM routines called directly by the user, etc. In the above description, internal MBASE format means the number is printed as it appears in the array --- an exponent followed by NDIG base MBASE digits. 5. ERROR CONDITIONS KFLAG is a condition parameter returned by the package after each call to one of the routines. Negative values indicate conditions for which a warning message will be printed unless KWARN = 0. Positive values indicate conditions that may be of interest but are not errors. No warning message is printed if KFLAG is nonnegative. KFLAG = 0 Normal operation. = 1 One of the operands in FMADD or FMSUB was insignificant with respect to the other, so that the result was equal to the argument of larger magnitude. = 2 In converting an FM number to a one word integer in FMM2I, the FM number was not exactly an integer. The next integer toward zero was returned. = -1 NDIG was less than 2 or more than NDIGMX. = -2 MBASE was less than 2 or more than MXBASE. = -3 An exponent was out of range. = -4 Invalid input argument(s) to an FM routine. UNKNOWN was returned. = -5 + or - OVERFLOW was generated as a result from an FM routine. = -6 + or - UNDERFLOW was generated as a result from an FM routine. = -7 The input string (array) to FMINP was not legal. = -8 The character array was not large enough in an input or output routine. = -9 Precision could not be raised enough to provide all requested guard digits. Increasing NDIGMX in all the PARAMETER statements may fix this. UNKNOWN was returned. = -10 An FM input argument was too small in magnitude to convert to the machine's single or double precision in FMM2SP or FMM2DP. Check that the definitions of SPMAX and DPMAX in FMSET are correct for the current machine. Zero was returned. When a negative KFLAG condition is encountered, the value of KWARN determines the action to be taken. KWARN = 0 Execution continues and no message is printed. = 1 A warning message is printed and execution continues. = 2 A warning message is printed and execution stops. The default setting is KWARN = 1. When an overflow or underflow is generated for an operation in which an input argument was already an overflow or underflow, no additional message is printed. When an unknown result is generated and an input argument was already unknown, no additional message is printed. In these cases the negative KFLAG value is still returned. IM routines handle exceptions like OVERFLOW or UNKNOWN in the same way as FM routines. When using IMMPY, the product of two large positive integers will return +OVERFLOW. The routine IMMPYM can be used to obtain a modular result without overflow. The largest representable IM integer is MBASE**NDIGMX - 1. For example, if MBASE is 10**7 and NDIGMX is set to 256, integers less than 10**1792 can be used. 6. OTHER PARAMETERS KRAD = 0 All angles in the trigonometric functions and inverse functions are measured in degrees. = 1 All angles are measured in radians. (Default) KROUND = 0 All final results are chopped (rounded toward zero). Intermediate results are rounded. = 1 All results are rounded to the nearest FM number, or to the value with an even last digit if the result is halfway between two FM numbers. (Default) KSWIDE defines the maximum screen width to be used for all unit KW output. Default is 80. KESWCH controls the action taken in FMINP and other input routines for strings like 'E7' that have no digits before the exponent field. Default is for 'E7' to translate like '1.0E+7'. CMCHAR defines the exponent letter to be used for FM variable output. Default is 'M', as in 1.2345M+678. KDEBUG = 0 Error checking is not done for valid input arguments and parameters like NDIG and MBASE upon entry to each routine. (Default) = 1 Some error checking is done. (Slower speed) See FMSET for additional description of these and other variables defining various FM conditions. 7. ARRAY DIMENSIONS The dimensions of the arrays in the FM package are defined using a PARAMETER statement at the top of each routine. The size of these arrays depends on the values of parameters NDIGMX and NBITS. NDIGMX is the maximum value the user may set for NDIG. NBITS is the number of bits used to represent integers for a given machine. See the EFFICIENCY discussion below. The standard version of FMLIB sets NDIGMX = 256, so on a 32-bit machine using MBASE = 10**7 the maximum precision is about 7*255+1 = 1786 significant digits. To change dimensions so that 10,000 significant digit calculation can be done, NDIGMX needs to be at least 10**4/7 + 5 = 1434. This allows for a few user guard digits to be defined when the package is initialized using CALL FMSET(10000). Changing 'NDIGMX = 256' to 'NDIGMX = 1434' everywhere in the package and the user's calling program will define all the new array sizes. If NDIG much greater than 256 is to be used and elementary functions will be needed, they will be faster if array MJSUMS is larger. The parameter defining the size of MJSUMS is set in the standard version by LJSUMS = 8*(LUNPCK+2). The 8 means that up to eight concurrent sums can be used by the elementary functions. The approximate number needed for best speed is given by the formula 0.051*Log(MBASE)*NDIG**(1/3) + 1.85 For example, with MBASE=10**7 and NDIG=1434 this gives 11. Changing 'LJSUMS = 8*(LUNPCK+2)' to 'LJSUMS =11*(LUNPCK+2)' everywhere in the package and the user's calling program will give slightly better speed. FM numbers in packed format have dimension 0:LPACK, and those in unpacked format have dimension 0:LUNPCK. 8. PORTABILITY In FMSET there is some machine-dependent code that attempts to approximate the largest representable integer value. The current code works on all machines tested, but if an FM run fails, check the MAXINT and INTMAX loops in FMSET. Values for SPMAX and DPMAX are also defined in FMSET that should be set to values near overflow for single precision and double precision. Setting KDEBUG = 1 may also identify some errors if a run fails. Some compilers object to a function like FMCOMP with side effects such as changing KFLAG or other common variables. Blocks of code in FMCOMP and IMCOMP that modify common are identified so they may be removed or commented out to produce a function without side effects. This disables trace printing in FMCOMP and IMCOMP, and error codes are not returned in KFLAG. See FMCOMP and IMCOMP for further details. All variables are explicitly declared in each routine. There is a commented IMPLICIT NONE statement in each routine that can be enabled to get more compiler diagnostic information in some testing or debugging situations. 9. NEW FOR VERSION 1.1 Version 1.0 used integer arrays and integer arithmetic internally to perform the multiple precision operations. Version 1.1 uses double precision arithmetic and arrays internally. This is usually faster at higher precisions, and on many machines it is also faster at lower precisions. Version 1.1 is written so that the arithmetic used can easily be changed from double precision to integer, or any other available arithmetic type. This permits the user to make the best use of a given machine's arithmetic hardware. See the EFFICIENCY discussion below. Several routines have undergone minor modification, but only a few changes should affect programs that used FM 1.0. Many of the routines are faster in version 1.1, because code has been added to take advantage of special cases for individual functions instead of using general formulas that are more compact. For example, there are separate routines using series for SINH and COSH instead of just calling EXP. FMEQU was the only routine that required the user to give the value of the current precision. This was to allow automatic rounding or zero-padding when changing precision. Since few user calls change precision, a new routine has been added for this case. FMEQ now handles this case and has a simple argument list that does not include the value of NDIG. FMEQU is used for changing precision. See the list of FM routines above for details. All variable names beginning with M in the package are now declared as double precision, so FM common blocks in the user's program need D.P. declarations, and FM variables (arrays) used in the calling program need to be D.P. /FMUSER/ is a common block holding parameters that define the arithmetic to be used and other user options. Several new variables have been added, including screen width to be used for output. See above for further description. /FMSAVE/ is a common block for saving constants to avoid re-computing them. Several new variables have been added. /FMBUFF/ is a common block containing a character array used to format FM numbers for output. Two new items have been added. New routines: All the IM routines are new for version 1.1. FMADDI increments an FM number by a small integer. It runs in O(1) time, on the average. FMCHSH returns both SINH(MA) and COSH(MA). When both are needed, this is almost twice as fast as making separate calls to FMCOSH and FMSINH. FMCSSN returns both SIN(MA) and COS(MA). When both are needed, this is almost twice as fast as making separate calls to FMCOS and FMSIN. FMFORM uses a format string to convert an FM number to a character string. FMFPRT prints an FM number using a format string. FMREAD reads an FM number written using FMWRIT. FMRPWR computes an FM number raised to a rational power. For cube roots and similar rational powers it is usually much faster than FMPWR. FMSQR squares an FM number. It is faster than using FMMPY. FMST2M converts character strings to FM format. Since FMINP converts character arrays, this routine can be more convenient for easily defining an FM number. For example, CALL FMST2M('123.4',MA). FMWRIT writes an FM number using a format for multi-line numbers with '&' at the end of all but the last line of a multi-line number. This allows automatic reading of FM numbers without needing to know the base, precision or format under which they were written. One extra word has been added to the dimensions of all FM numbers. Word zero in each array contains a value used to monitor cancellation error arising from addition or subtraction. This value approximates the number of bits of precision for an FM value. It allows higher level FM functions to detect cases where too much cancellation has occurred. KACCSW is a switch variable in COMMON /FM/ used internally to enable cancellation error monitoring. 10. EFFICIENCY To take advantage of hardware architecture on different machines, the package has been designed so that the arithmetic used to perform the multiple precision operations can easily be changed. All variables that must be changed to get a different arithmetic have names beginning with 'M' and are declared using REAL (KIND(0.0D0)) :: m.... For example, to change the package to use integer arithmetic internally, make these two changes everywhere in the package: change 'REAL (KIND(0.0D0)) :: m' to 'INTEGER m', change 'DINT(' to 'INT('. On some systems, changing 'DINT(' to '(' may give better speed. When changing to a different type of arithmetic, all FM common blocks and arrays in the user's program must be changed to agree. In a few places in FM, where a DINT function is not supposed to be changed, it is spelled 'DINT (' so the global change will not find it. This version restricts the base used to be also representable in integer variables, so using precision above double usually does not save much time unless integers can also be declared at a higher precision. Using IEEE Extended would allow a base of around 10**9 to be chosen, but the delayed digit-normalization method used for multiplication and division means that a slightly smaller base like 10**8 would usually run faster. This would usually not be much faster than using 10**7 with double precision. The value of NBITS defined as a parameter in most FM routines refers to the number of bits used to represent integers in an M-variable word. Typical values for NBITS are: 24 for IEEE single precision, 32 for integer, 53 for IEEE double precision. NBITS controls only array size, so setting it too high is ok, but then the program will use more memory than necessary. For cases where special compiler directives or minor re-writing of the code may improve speed, several of the most important loops in FM are identified by comments containing the string '(Inner Loop)'. ------------------------------------------------------------------------ -------------------------- zmlib.f90 Notes --------------------------- The ZM routines perform complex floating-point multiple-precision arithmetic. These routines use a Fortran 90 version of the FMLIB package (version 1.1) for real floating-point multiple-precision arithmetic. FMLIB is Algorithm 693, ACM Transactions on Mathematical Software, Vol. 17, No. 2, June 1991, pages 273-283. This package and FMLIB 1.1 use double precision arithmetic and arrays internally. This is usually faster at higher precision, and on many machines it is also faster at lower precision. Both packages are written so that the arithmetic used can easily be changed from double precision to integer, or another available arithmetic type. See the EFFICIENCY discussion in the fmlib.f90 Notes for details. 1. INITIALIZING THE PACKAGE Before calling any routine in the package, several variables in the common blocks /FMUSER/, /FM/, /FMSAVE/, /FMBUFF/, and /ZMUSER/ must be initialized. These common blocks contain information that is saved between calls, so they should be declared in the main program. Subroutine ZMSET initializes these variables to default values and defines all machine-dependent values in the package. After calling ZMSET once at the start of a program, the user may sometimes want to reset some of the variables in common blocks /FMUSER/ or /ZMUSER/. 2. REPRESENTATION OF ZM NUMBERS The format for complex FM numbers (called ZM numbers below) is very similar to that for real FM numbers in FMLIB. Each ZM array holds two FM numbers to represent the real and imaginary parts of a complex number. Each ZM array is twice as long as a corresponding FM array, with the imaginary part starting at the midpoint of the array. As with FM, there are packed and unpacked formats for the numbers. 3. INPUT/OUTPUT ROUTINES All versions of the input routines perform free-format conversion from characters to ZM numbers. a. Conversion to or from a character array ZMINP converts from a character*1 array to an ZM number. ZMOUT converts an ZM number to base 10 and formats it for output as an array of type character*1. The output is left justified in the array, and the format is defined by variables in common, so that a separate format definition does not have to be provided for each output call. For the output format of ZM numbers, JFORM1 and JFORM2 determine the format for the individual parts of a complex number as described in the FMLIB documentation. JFORMZ (in /ZMUSER/) determines the combined output format of the real and imaginary parts. JFORMZ = 1 normal setting : 1.23 - 4.56 i = 2 use capital I : 1.23 - 4.56 I = 3 parenthesis format ( 1.23 , -4.56 ) JPRNTZ (in /ZMUSER/) controls whether to print real and imaginary parts on one line whenever possible. JPRNTZ = 1 print both parts as a single string : 1.23456789M+321 - 9.87654321M-123 i = 2 print on separate lines without the 'i' : 1.23456789M+321 -9.87654321M-123 b. Conversion to or from a character string ZMST2M converts from a character string to an ZM number. ZMFORM converts an ZM number to a character string according to a format provided in each call. The format descriptions are more like that of a Fortran FORMAT statement, and integer or fixed-point output is right justified. c. Direct read or write ZMPRNT uses ZMOUT to print one ZM number. ZMFPRT uses ZMFORM to print one ZM number. ZMWRIT writes ZM numbers for later input using ZMREAD. ZMREAD reads ZM numbers written by ZMWRIT. For further description of these routines, see the list of ZM routines above. 4. ARRAY DIMENSIONS The parameters LPACKZ and LUNPKZ define the size of the packed and unpacked ZM arrays. The real part starts at the beginning of the array, and the imaginary part starts at word KPTIMP for packed format or at word KPTIMU for unpacked format. =========================================================================== =========================================================================== SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'RES1' then echo shar: will not over-write existing file "'RES1'" else cat << SHAR_EOF > 'RES1' 53 cases tested. No errors were found. SHAR_EOF fi # end of overwriting check if test -f 'RES2' then echo shar: will not over-write existing file "'RES2'" else cat << SHAR_EOF > 'RES2' Sample 1. Find a root of f(x) = x**5 - 3x**4 + x**3 - 4x**2 + x - 6 = 0. Iteration Newton Approximation 0 .560000000000000000000000000000 + 1.060000000000000000000000000000 i 1 .561964780980333719745880263787 + 1.061135231152741154895778904059 i 2 .561958308372772219534516409947 + 1.061134679566247415769456345141 i 3 .561958308335403235495113920123 + 1.061134679604332556981397796290 i 4 .561958308335403235498111195347 + 1.061134679604332556983391239059 i 5 .561958308335403235498111195347 + 1.061134679604332556983391239059 i Sample 2. 44 terms were added to get Exp(1.23-2.34i) Result= -2.379681796854777515745457977697 - 2.458032970832342652397461908326 i All results were ok. SHAR_EOF fi # end of overwriting check if test -f 'RES3' then echo shar: will not over-write existing file "'RES3'" else cat << SHAR_EOF > 'RES3' 108 cases tested. No errors were found. SHAR_EOF fi # end of overwriting check if test -f 'RES4' then echo shar: will not over-write existing file "'RES4'" else cat << SHAR_EOF > 'RES4' Sample 1. Find a root of f(x) = x**5 - 3x**4 + x**3 - 4x**2 + x - 6 = 0. Iteration Newton Approximation 0 3.120000000000000000000000000000000000000000000000000000000000 1 3.120656718532108533919391265947916793506741449899073468862023 2 3.120656215327022122238354686569835883519704471397219749798884 3 3.120656215326726500470956115551705969611230193197937042123082 4 3.120656215326726500470956013523797484654623935599078168006617 5 3.120656215326726500470956013523797484654623935599066014988828 6 3.120656215326726500470956013523797484654623935599066014988828 Sample 2. 109 terms were added Zeta(3) = 1.202056903159594285399738161511449990764986292340498881792272 Sample 3. 22 values were tested p = 1000000000000000000000000000000000000000000000000000000000000000659661 All results were ok. SHAR_EOF fi # end of overwriting check if test -f 'RES5' then echo shar: will not over-write existing file "'RES5'" else cat << SHAR_EOF > 'RES5' 603 cases tested. No errors were found. SHAR_EOF fi # end of overwriting check if test -f 'RES6' then echo shar: will not over-write existing file "'RES6'" else cat << SHAR_EOF > 'RES6' Sample 1. Real root of f(x) = x**5 - 3x**4 + x**3 - 4x**2 + x - 6 = 0. Iteration Newton Approximation 0 3.120000000000000000000000000000000000000000000000000000000000 1 3.120656718532108533919391265947916793506741449899073468862023 2 3.120656215327022122238354686569835883519704471397219749798884 3 3.120656215326726500470956115551705969611230193197937042123082 4 3.120656215326726500470956013523797484654623935599078168006617 5 3.120656215326726500470956013523797484654623935599066014988828 6 3.120656215326726500470956013523797484654623935599066014988828 Sample 2. 109 terms were added Zeta(3) = 1.202056903159594285399738161511449990764986292340498881792272 Sample 3. 22 values were tested p = 1000000000000000000000000000000000000000000000000000000000000000659661 Sample 4. Complex root of f(x) = x**5 - 3x**4 + x**3 - 4x**2 + x - 6 = 0. Iteration Newton Approximation 0 .560000000000000000000000000000 + 1.060000000000000000000000000000 i 1 .561964780980333719745880263787 + 1.061135231152741154895778904059 i 2 .561958308372772219534516409947 + 1.061134679566247415769456345141 i 3 .561958308335403235495113920123 + 1.061134679604332556981397796290 i 4 .561958308335403235498111195347 + 1.061134679604332556983391239059 i 5 .561958308335403235498111195347 + 1.061134679604332556983391239059 i Sample 5. 44 terms were added to get Exp(1.23-2.34i) Result= -2.379681796854777515745457977697 - 2.458032970832342652397461908326 i All results were ok. SHAR_EOF fi # end of overwriting check if test -f 'driver1.f90' then echo shar: will not over-write existing file "'driver1.f90'" else cat << SHAR_EOF > 'driver1.f90' PROGRAM test ! David M. Smith 6-14-96 ! This is a test program for ZMLIB 1.1, a multiple-precision complex ! arithmetic package. Most of the ZM routines are tested, and the ! results are checked to 50 significant digits. ! This program uses both ZMLIB.f90 and FMLIB.f90. ! These five common blocks contain information that must be saved ! between calls, so they should be declared in the main program. ! The parameter statement defines array sizes and pointers, and ! contains the FMLIB parameters, followed by ZMLIB parameters. ! .. Parameters .. INTEGER, PARAMETER :: lhash1 = 0, lhash2 = 256, nbits = 64, ndigmx = 256 INTEGER, PARAMETER :: lpack = (ndigmx+1)/2 + 1 INTEGER, PARAMETER :: lpackz = 2*lpack + 1 INTEGER, PARAMETER :: lunpck = (6*ndigmx)/5 + 20 INTEGER, PARAMETER :: lunpkz = 2*lunpck + 1 INTEGER, PARAMETER :: kptimp = lpack + 1 INTEGER, PARAMETER :: kptimu = lunpck + 1 INTEGER, PARAMETER :: ljsums = 8*(lunpck+2) INTEGER, PARAMETER :: lmbuff = ((lunpck+3)*(nbits-1)*301)/2000 + 6 INTEGER, PARAMETER :: lmbufz = 2*lmbuff + 10 INTEGER, PARAMETER :: lmwa = 2*lunpck ! .. ! .. Local Scalars .. INTEGER :: klog, ncase, nerror ! Character strings used for input and output. CHARACTER (160) :: st1, st2 ! Declare arrays for ZM complex variables (MA, MB, MC, MD) ! and for FM real variables (MAFM, MBFM). All are in ! unpacked format. ! .. ! .. Local Arrays .. REAL (KIND(0.0D0)) :: ma(0:lunpkz), mafm(0:lunpck), mb(0:lunpkz), mbfm(0:lunpck), & mc(0:lunpkz), md(0:lunpkz) ! .. ! .. External Subroutines .. EXTERNAL test1, test2, test3, test4, test5, test6, test7, test8, zmset ! .. ! .. Scalars in Common .. REAL :: alogm2, alogmb, alogmt, alogmx, runkno, spmax REAL (KIND(0.0D0)) :: dlogeb, dlogmb, dlogpi, dlogtn, dlogtp, & dlogtw, dpeps, dpmax, dppi REAL (KIND(0.0D0)) :: maxint, mbase, mblogs, mbse, mbslb, mbsli, & mbspi, mexpab, mexpov, mexpun, munkno, mxbase, mxexp, mxexp2 INTEGER :: intmax, iunkno, jform1, jform2, jformz, jprntz, kaccsw, & kdebug, keswch, kflag, krad, kround, ksub, kswide, kw, kwarn, lvltrc, & ncall, ndg2mx, ndig, ndige, ndiglb, ndigli, ndigpi, ngrd21, ngrd22, & ngrd52, ntrace CHARACTER (1) :: cmchar ! .. ! .. Arrays in Common .. REAL (KIND(0.0D0)) :: mesav(0:lunpck), mlbsav(0:lunpck), mln1(0:lunpck), & mln2(0:lunpck), mln3(0:lunpck), mln4(0:lunpck), mpisav(0:lunpck), & mwa(lmwa) INTEGER :: khasht(lhash1:lhash2), khashv(lhash1:lhash2) CHARACTER (1) :: cmbuff(lmbuff) CHARACTER (6) :: namest(0:50) ! .. ! .. Common Blocks .. COMMON /fm/mwa, ncall, kaccsw, mxexp, mxexp2, mexpun, mexpov, munkno, & iunkno, runkno, mxbase, ndg2mx, spmax, dpmax, maxint, intmax, ksub COMMON /fmbuff/cmbuff, namest, cmchar COMMON /fmsave/ndigpi, ndige, ndiglb, ndigli, mbspi, mbse, mbslb, mbsli, & mpisav, mesav, mlbsav, mln1, mln2, mln3, mln4, mblogs, mexpab, alogmb, & alogm2, alogmx, alogmt, dlogmb, dlogtn, dlogtw, dlogtp, dlogpi, dppi, & dpeps, dlogeb, khasht, khashv, ngrd21, ngrd52, ngrd22 COMMON /fmuser/mbase, ndig, jform1, jform2, krad, kw, ntrace, lvltrc, & kflag, kwarn, kround, kswide, keswch, kdebug COMMON /zmuser/jformz, jprntz ! .. ! Set precision to give at least 50 significant digits ! and initialize the FMLIB package. CALL zmset(50) ! Write output to the standard FM output (unit KW, defined ! in subroutine FMSET), and also to the file TESTZM.LOG. klog = 18 OPEN (klog,file='TESTZM.LOG') ! NERROR is the number of errors found. ! NCASE is the number of cases tested. nerror = 0 ! Test input and output conversion. CALL test1(ma,mb,mc,md,mafm,mbfm,st1,st2,ncase,nerror,klog) ! Test add and subtract. CALL test2(ma,mb,mc,md,mafm,mbfm,st1,st2,ncase,nerror,klog) ! Test multiply, divide and square root. CALL test3(ma,mb,mc,md,mafm,mbfm,st1,st2,ncase,nerror,klog) ! Test exponentials. CALL test4(ma,mb,mc,md,mafm,mbfm,st1,st2,ncase,nerror,klog) ! Test logarithms. CALL test5(ma,mb,mc,md,mafm,mbfm,st1,st2,ncase,nerror,klog) ! Test trigonometric functions. CALL test6(ma,mb,mc,md,mafm,mbfm,st1,st2,ncase,nerror,klog) ! Test inverse trigonometric functions. CALL test7(ma,mb,mc,md,mafm,mbfm,st1,st2,ncase,nerror,klog) ! Test hyperbolic functions. CALL test8(ma,mb,mc,md,mafm,mbfm,st1,st2,ncase,nerror,klog) ! End of tests. IF (nerror==0) THEN WRITE (kw,90000) ncase WRITE (klog,90000) ncase ELSE ! Write some of the initialized values in common. WRITE (klog,*) ' NDIG,MBASE,JFORM1,JFORM2,KRAD = ' WRITE (klog,*) ndig, mbase, jform1, jform2, krad WRITE (klog,*) ' KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND = ' WRITE (klog,*) kw, ntrace, lvltrc, kflag, kwarn, kround WRITE (klog,*) ' NCALL,MXEXP,MXEXP2,KACCSW,MEXPUN,MEXPOV' WRITE (klog,*) ncall, mxexp, mxexp2, kaccsw, mexpun, mexpov WRITE (klog,*) ' MUNKNO,IUNKNO,RUNKNO,MXBASE,NDG2MX = ' WRITE (klog,*) munkno, iunkno, runkno, mxbase, ndg2mx WRITE (klog,*) ' MAXINT,INTMAX,SPMAX,DPMAX = ' WRITE (klog,*) maxint, intmax, spmax, dpmax WRITE (klog,*) ' ALOGMB,ALOGM2,ALOGMX,ALOGMT,DLOGMB,DLOGTN =' WRITE (klog,*) alogmb, alogm2, alogmx, alogmt, dlogmb, dlogtn WRITE (klog,*) ' DLOGTW,DLOGTP,DLOGPI,DPPI =' WRITE (klog,*) dlogtw, dlogtp, dlogpi, dppi WRITE (klog,*) ' DPEPS,DLOGEB =' WRITE (klog,*) dpeps, dlogeb WRITE (kw,90010) ncase, nerror WRITE (klog,90010) ncase, nerror END IF WRITE (kw,*) ' End of run.' STOP 90000 FORMAT (///1X,I5,' cases tested. No errors were found.'/) 90010 FORMAT (///1X,I5,' cases tested.',I4,' error(s) found.'/) END PROGRAM test SUBROUTINE test1(ma,mb,mc,md,mafm,mbfm,st1,st2,ncase,nerror,klog) ! Input and output testing. ! Logical function for comparing FM numbers. ! .. Parameters .. INTEGER, PARAMETER :: nbits = 64, ndigmx = 256 INTEGER, PARAMETER :: lpack = (ndigmx+1)/2 + 1 INTEGER, PARAMETER :: lpackz = 2*lpack + 1 INTEGER, PARAMETER :: lunpck = (6*ndigmx)/5 + 20 INTEGER, PARAMETER :: lunpkz = 2*lunpck + 1 INTEGER, PARAMETER :: kptimp = lpack + 1 INTEGER, PARAMETER :: kptimu = lunpck + 1 INTEGER, PARAMETER :: ljsums = 8*(lunpck+2) INTEGER, PARAMETER :: lmbuff = ((lunpck+3)*(nbits-1)*301)/2000 + 6 INTEGER, PARAMETER :: lmbufz = 2*lmbuff + 10 INTEGER, PARAMETER :: lmwa = 2*lunpck ! .. ! .. Scalar Arguments .. INTEGER :: klog, ncase, nerror CHARACTER (160) :: st1, st2 ! .. ! .. Array Arguments .. REAL (KIND(0.0D0)) :: ma(0:lunpkz), mafm(0:lunpck), mb(0:lunpkz), mbfm(0:lunpck), & mc(0:lunpkz), md(0:lunpkz) ! .. ! .. External Functions .. LOGICAL, EXTERNAL :: fmcomp ! .. ! .. External Subroutines .. EXTERNAL errprt, fmi2m, fmipwr, zm2i2m, zmabs, zmdivi, zmform, zmmpyi, & zmst2m, zmsub ! .. ! .. Scalars in Common .. REAL (KIND(0.0D0)) :: mbase INTEGER :: jform1, jform2, kdebug, keswch, kflag, krad, kround, kswide, & kw, kwarn, lvltrc, ndig, ntrace ! .. ! .. Common Blocks .. COMMON /fmuser/mbase, ndig, jform1, jform2, krad, kw, ntrace, lvltrc, & kflag, kwarn, kround, kswide, keswch, kdebug ! .. WRITE (kw,90000) ncase = 1 CALL zmst2m('123 + 456 i',ma) CALL zm2i2m(123,456,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-48,mbfm) ! Use the .NOT. because FMCOMP returns FALSE for special ! cases like MD = UNKNOWN, and these should be treated ! as errors for these tests. IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMST2M',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 2 st1 = '0.3505154639175257731958762886597938144329896907216495 + ' // & '0.7319587628865979381443298969072164948453608247422680 i' CALL zmst2m(st1,ma) CALL zm2i2m(34,71,mc) CALL zmdivi(mc,97,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMST2M',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 3 st1 = '0.3505154639175257731958762886597938144329896907216495E-5 ' // & '+ 0.7319587628865979381443298969072164948453608247422680D-5 i' CALL zmst2m(st1,ma) CALL zm2i2m(34,71,mc) CALL zmdivi(mc,9700000,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-55,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMST2M',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 4 st1 = '7.699115044247787610619469026548672566371681415929204e 03 ' // & '- 5.221238938053097345132743362831858407079646017699115M 03 I' CALL zmst2m(st1,ma) CALL zm2i2m(87,-59,mc) CALL zmdivi(mc,113,mc) CALL zmmpyi(mc,10000,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-47,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMST2M',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 5 st1 = '7.699115044247787610619469026548672566371681415929204e+3 ' // & '- 5.221238938053097345132743362831858407079646017699115M+3 i' CALL zmst2m(st1,ma) CALL zmform('F53.33','F50.30',ma,st2) CALL zmst2m(st2,ma) st1 = '7699.115044247787610619469026548673 ' // & '-5221.238938053097345132743362831858 i' CALL zmst2m(st1,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-30,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMFORM',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 6 st1 = '7.699115044247787610619469026548672566371681415929204e+3 ' // & '- 5.221238938053097345132743362831858407079646017699115M+3 i' CALL zmst2m(st1,ma) CALL zmform('I9','I7',ma,st2) CALL zmst2m(st2,ma) st1 = '7699 -5221 i' CALL zmst2m(st1,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(0,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMFORM',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 7 st1 = '7.699115044247787610619469026548672566371681415929204e+3 ' // & '- 5.221238938053097345132743362831858407079646017699115M+3 i' CALL zmst2m(st1,ma) CALL zmform('E59.50','E58.49',ma,st2) CALL zmst2m(st2,ma) st1 = '7.6991150442477876106194690265486725663716814159292E3' // & '- 5.221238938053097345132743362831858407079646017699E3 i' CALL zmst2m(st1,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-48,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMFORM',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 8 st1 = '7.699115044247787610619469026548672566371681415929204e+3 ' // & '- 5.221238938053097345132743362831858407079646017699115M+3 i' CALL zmst2m(st1,ma) CALL zmform('1PE59.50','1PE58.49',ma,st2) CALL zmst2m(st2,ma) CALL zmst2m(st1,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-44,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMFORM',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF RETURN 90000 FORMAT (/' Testing input and output routines.') END SUBROUTINE test1 SUBROUTINE test2(ma,mb,mc,md,mafm,mbfm,st1,st2,ncase,nerror,klog) ! Test add and subtract. ! .. Parameters .. INTEGER, PARAMETER :: nbits = 64, ndigmx = 256 INTEGER, PARAMETER :: lpack = (ndigmx+1)/2 + 1 INTEGER, PARAMETER :: lpackz = 2*lpack + 1 INTEGER, PARAMETER :: lunpck = (6*ndigmx)/5 + 20 INTEGER, PARAMETER :: lunpkz = 2*lunpck + 1 INTEGER, PARAMETER :: kptimp = lpack + 1 INTEGER, PARAMETER :: kptimu = lunpck + 1 INTEGER, PARAMETER :: ljsums = 8*(lunpck+2) INTEGER, PARAMETER :: lmbuff = ((lunpck+3)*(nbits-1)*301)/2000 + 6 INTEGER, PARAMETER :: lmbufz = 2*lmbuff + 10 INTEGER, PARAMETER :: lmwa = 2*lunpck ! .. ! .. Scalar Arguments .. INTEGER :: klog, ncase, nerror CHARACTER (160) :: st1, st2 ! .. ! .. Array Arguments .. REAL (KIND(0.0D0)) :: ma(0:lunpkz), mafm(0:lunpck), mb(0:lunpkz), mbfm(0:lunpck), & mc(0:lunpkz), md(0:lunpkz) ! .. ! .. External Functions .. LOGICAL, EXTERNAL :: fmcomp ! .. ! .. External Subroutines .. EXTERNAL errprt, fmi2m, fmipwr, zm2i2m, zmabs, zmadd, zmst2m, zmsub ! .. ! .. Scalars in Common .. REAL (KIND(0.0D0)) :: mbase INTEGER :: jform1, jform2, kdebug, keswch, kflag, krad, kround, kswide, & kw, kwarn, lvltrc, ndig, ntrace ! .. ! .. Common Blocks .. COMMON /fmuser/mbase, ndig, jform1, jform2, krad, kw, ntrace, lvltrc, & kflag, kwarn, kround, kswide, keswch, kdebug ! .. WRITE (kw,90000) ncase = 9 CALL zmst2m('123 + 456 i',ma) CALL zmst2m('789 - 543 i',mb) CALL zmadd(ma,mb,ma) CALL zm2i2m(912,-87,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(0,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMADD ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 10 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) st1 = '.3505154639175257731958762886597938144329896907216495 ' // & '+ .7319587628865979381443298969072164948453608247422680 i' CALL zmst2m(st1,mb) CALL zmadd(ma,mb,ma) st2 = '1.1204269683423045342578231913146610710701578323145698 ' // & '+ 0.2098348690812882036310555606240306541373962229723565 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-49,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMADD ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 11 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) st1 = '.3505154639175257731958762886597938144329896907216495 ' // & '+ .7319587628865979381443298969072164948453608247422680 i' CALL zmst2m(st1,mb) CALL zmsub(ma,mb,ma) st2 = '0.4193960405072529878660706139950734422041784508712709 ' // & '- 1.2540826566919076726576042331904023355533254265121795 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-49,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMSUB ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 12 st1 = '.7699115044247787610619469026548672566371681415929204E3 ' // & '- .5221238938053097345132743362831858407079646017699115E3 i' CALL zmst2m(st1,ma) st1 = '.3505154639175257731958762886597938144329896907216495 ' // & '+ .7319587628865979381443298969072164948453608247422680 i' CALL zmst2m(st1,mb) CALL zmsub(ma,mb,ma) st2 = '769.5609889608612352887510263662074628227351519021987045 ' // & '- 522.8558525681963324514186661800930572028099625946537725 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-47,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMSUB ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF RETURN 90000 FORMAT (/' Testing add and subtract routines.') END SUBROUTINE test2 SUBROUTINE test3(ma,mb,mc,md,mafm,mbfm,st1,st2,ncase,nerror,klog) ! Test multiply, divide and square root. ! .. Parameters .. INTEGER, PARAMETER :: nbits = 64, ndigmx = 256 INTEGER, PARAMETER :: lpack = (ndigmx+1)/2 + 1 INTEGER, PARAMETER :: lpackz = 2*lpack + 1 INTEGER, PARAMETER :: lunpck = (6*ndigmx)/5 + 20 INTEGER, PARAMETER :: lunpkz = 2*lunpck + 1 INTEGER, PARAMETER :: kptimp = lpack + 1 INTEGER, PARAMETER :: kptimu = lunpck + 1 INTEGER, PARAMETER :: ljsums = 8*(lunpck+2) INTEGER, PARAMETER :: lmbuff = ((lunpck+3)*(nbits-1)*301)/2000 + 6 INTEGER, PARAMETER :: lmbufz = 2*lmbuff + 10 INTEGER, PARAMETER :: lmwa = 2*lunpck ! .. ! .. Scalar Arguments .. INTEGER :: klog, ncase, nerror CHARACTER (160) :: st1, st2 ! .. ! .. Array Arguments .. REAL (KIND(0.0D0)) :: ma(0:lunpkz), mafm(0:lunpck), mb(0:lunpkz), mbfm(0:lunpck), & mc(0:lunpkz), md(0:lunpkz) ! .. ! .. External Functions .. LOGICAL, EXTERNAL :: fmcomp ! .. ! .. External Subroutines .. EXTERNAL errprt, fmi2m, fmipwr, zm2i2m, zmabs, zmdiv, zmdivi, zmmpy, & zmmpyi, zmsqr, zmsqrt, zmst2m, zmsub ! .. ! .. Scalars in Common .. REAL (KIND(0.0D0)) :: mbase INTEGER :: jform1, jform2, kdebug, keswch, kflag, krad, kround, kswide, & kw, kwarn, lvltrc, ndig, ntrace ! .. ! .. Common Blocks .. COMMON /fmuser/mbase, ndig, jform1, jform2, krad, kw, ntrace, lvltrc, & kflag, kwarn, kround, kswide, keswch, kdebug ! .. WRITE (kw,90000) ncase = 13 CALL zmst2m('123 + 456 i',ma) CALL zmst2m('789 - 543 i',mb) CALL zmmpy(ma,mb,ma) CALL zm2i2m(344655,292995,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(0,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMMPY ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 14 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) st1 = '.3505154639175257731958762886597938144329896907216495 ' // & '+ .7319587628865979381443298969072164948453608247422680 i' CALL zmst2m(st1,mb) CALL zmmpy(ma,mb,ma) st2 = '0.6520390475321594745005017790347596022260742632971444 ' // & '+ 0.3805309734513274336283185840707964601769911504424779 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMMPY ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 15 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) st1 = '.3505154639175257731958762886597938144329896907216495 ' // & '+ .7319587628865979381443298969072164948453608247422680 i' CALL zmst2m(st1,mb) CALL zmdiv(ma,mb,ma) st2 = '-.1705178497731560089737969128653459210208765017614861 ' // & '- 1.1335073636829696356072949942949842987114804337239972 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-49,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMDIV ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 16 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmmpyi(ma,36,ma) st2 = '27.7168141592920353982300884955752212389380530973451327 ' // & '- 18.7964601769911504424778761061946902654867256637168142 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-48,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMMPYI',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 17 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmdivi(ma,37,ma) st2 = '2.080841903850753408275532169337479071992346328629514E-2 ' // & '- 1.411145658933269552738579287251853623535039464243004E-2 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-52,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMDIVI',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 18 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmsqr(ma,ma) st2 = '0.3201503641632077688150990680554467851828647505677813 ' // & '- 0.8039783851515388832328295089670295246299631921058814 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMSQR ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 19 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmsqrt(ma,ma) st2 = '0.9219999909012323458336720551458583330580388434229845 ' // & '- 0.2831474506279259570386845864488094697732718981999941 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMSQRT',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF RETURN 90000 FORMAT (/' Testing multiply, divide and square root routines.') END SUBROUTINE test3 SUBROUTINE test4(ma,mb,mc,md,mafm,mbfm,st1,st2,ncase,nerror,klog) ! Test exponentials. ! .. Parameters .. INTEGER, PARAMETER :: nbits = 64, ndigmx = 256 INTEGER, PARAMETER :: lpack = (ndigmx+1)/2 + 1 INTEGER, PARAMETER :: lpackz = 2*lpack + 1 INTEGER, PARAMETER :: lunpck = (6*ndigmx)/5 + 20 INTEGER, PARAMETER :: lunpkz = 2*lunpck + 1 INTEGER, PARAMETER :: kptimp = lpack + 1 INTEGER, PARAMETER :: kptimu = lunpck + 1 INTEGER, PARAMETER :: ljsums = 8*(lunpck+2) INTEGER, PARAMETER :: lmbuff = ((lunpck+3)*(nbits-1)*301)/2000 + 6 INTEGER, PARAMETER :: lmbufz = 2*lmbuff + 10 INTEGER, PARAMETER :: lmwa = 2*lunpck ! .. ! .. Scalar Arguments .. INTEGER :: klog, ncase, nerror CHARACTER (160) :: st1, st2 ! .. ! .. Array Arguments .. REAL (KIND(0.0D0)) :: ma(0:lunpkz), mafm(0:lunpck), mb(0:lunpkz), mbfm(0:lunpck), & mc(0:lunpkz), md(0:lunpkz) ! .. ! .. External Functions .. LOGICAL, EXTERNAL :: fmcomp ! .. ! .. External Subroutines .. EXTERNAL errprt, fmi2m, fmipwr, zmabs, zmexp, zmipwr, zmpwr, zmrpwr, & zmst2m, zmsub ! .. ! .. Scalars in Common .. REAL (KIND(0.0D0)) :: mbase INTEGER :: jform1, jform2, kdebug, keswch, kflag, krad, kround, kswide, & kw, kwarn, lvltrc, ndig, ntrace ! .. ! .. Common Blocks .. COMMON /fmuser/mbase, ndig, jform1, jform2, krad, kw, ntrace, lvltrc, & kflag, kwarn, kround, kswide, keswch, kdebug ! .. WRITE (kw,90000) ncase = 20 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmexp(ma,ma) st2 = '1.8718374504057787925867989348073888855260008469310002 ' // & '- 1.0770279996847678711699041910427261417963102075889234 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-49,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMEXP ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 21 st1 = '5.7699115044247787610619469026548672566371681415929204 ' // & '- 4.5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmexp(ma,ma) st2 = '-60.6144766542152809520229386164396710991242264070603612 ' // & '+ 314.7254994809539691403004121118801578835669635535466592 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-47,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMEXP ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 22 st1 = '1.7699115044247787610619469026548672566371681415929204 ' // & '- 1.5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmipwr(ma,45,ma) st2 = '31595668743300099.70429472191424818167262151605608585179 ' // & '- 19209634448276799.67717448173630165852744930837930753788 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-33,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMIPWR',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 23 st1 = '1.7699115044247787610619469026548672566371681415929204 ' // & '- 1.5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmipwr(ma,-122,ma) st2 = '3.1000215641022021714480000129414241564868699479432E-46 ' // & '- 1.1687846789859477815450163510927243367234863123667E-45 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-93,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMIPWR',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 24 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) st1 = '.3505154639175257731958762886597938144329896907216495 ' // & '+ .7319587628865979381443298969072164948453608247422680 i' CALL zmst2m(st1,mb) CALL zmpwr(ma,mb,ma) st2 = '1.4567089343012352449621841355636496276866203747888724 ' // & '- 0.3903177712261966292764255714390622205129978923650749 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-49,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMPWR ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 25 st1 = '.3505154639175257731958762886597938144329896907216495 ' // & '+ .7319587628865979381443298969072164948453608247422680 i' CALL zmst2m(st1,ma) st1 = '2.7699115044247787610619469026548672566371681415929204 ' // & '- 0.5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,mb) CALL zmpwr(ma,mb,ma) st2 = '-1.0053105716678380336247948739245187868180079734997482 ' // & '- 0.0819537653234704467729051473979237153087038930127116 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-49,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMPWR ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 26 st1 = '0.7699115044247787610619469026548672566371681415929204 ' // & '- 0.5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmrpwr(ma,2,7,ma) st2 = '0.9653921326136512316639621651337975772631340364271270 ' // & '- 0.1659768285667051396562270035411852432430188906482848 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMRPWR',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 27 st1 = '0.7699115044247787610619469026548672566371681415929204 ' // & '- 0.5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmrpwr(ma,-19,7,ma) st2 = '-0.0567985880053556315170006800325686036902111276420647 ' // & '+ 1.2154793972711356706410882510363594270389067962568571 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-49,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMRPWR',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF RETURN 90000 FORMAT (/' Testing exponential routines.') END SUBROUTINE test4 SUBROUTINE test5(ma,mb,mc,md,mafm,mbfm,st1,st2,ncase,nerror,klog) ! Test logarithms. ! .. Parameters .. INTEGER, PARAMETER :: nbits = 64, ndigmx = 256 INTEGER, PARAMETER :: lpack = (ndigmx+1)/2 + 1 INTEGER, PARAMETER :: lpackz = 2*lpack + 1 INTEGER, PARAMETER :: lunpck = (6*ndigmx)/5 + 20 INTEGER, PARAMETER :: lunpkz = 2*lunpck + 1 INTEGER, PARAMETER :: kptimp = lpack + 1 INTEGER, PARAMETER :: kptimu = lunpck + 1 INTEGER, PARAMETER :: ljsums = 8*(lunpck+2) INTEGER, PARAMETER :: lmbuff = ((lunpck+3)*(nbits-1)*301)/2000 + 6 INTEGER, PARAMETER :: lmbufz = 2*lmbuff + 10 INTEGER, PARAMETER :: lmwa = 2*lunpck ! .. ! .. Scalar Arguments .. INTEGER :: klog, ncase, nerror CHARACTER (160) :: st1, st2 ! .. ! .. Array Arguments .. REAL (KIND(0.0D0)) :: ma(0:lunpkz), mafm(0:lunpck), mb(0:lunpkz), mbfm(0:lunpck), & mc(0:lunpkz), md(0:lunpkz) ! .. ! .. External Functions .. LOGICAL, EXTERNAL :: fmcomp ! .. ! .. External Subroutines .. EXTERNAL errprt, fmi2m, fmipwr, zmabs, zmlg10, zmln, zmst2m, zmsub ! .. ! .. Scalars in Common .. REAL (KIND(0.0D0)) :: mbase INTEGER :: jform1, jform2, kdebug, keswch, kflag, krad, kround, kswide, & kw, kwarn, lvltrc, ndig, ntrace ! .. ! .. Common Blocks .. COMMON /fmuser/mbase, ndig, jform1, jform2, krad, kw, ntrace, lvltrc, & kflag, kwarn, kround, kswide, keswch, kdebug ! .. WRITE (kw,90000) ncase = 28 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmln(ma,ma) st2 = '-0.0722949652393911311212450699415231782692434885813725 ' // & '- 0.5959180055163009910007765127008371205749515965219804 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMLN ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 29 st1 = '.7699115044247787610619469026548672566371681415929204E28 ' // & '- .5221238938053097345132743362831858407079646017699115E28 i' CALL zmst2m(st1,ma) CALL zmln(ma,ma) st2 = '64.4000876385938880213825156612206746345615981930242708 ' // & '- 0.5959180055163009910007765127008371205749515965219804 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-48,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMLN ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 30 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmlg10(ma,ma) st2 = '-0.0313973044728549715287589498363619677438302809470943 ' // & '- 0.2588039014625211035392823012785304771809982053965284 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMLG10',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 31 st1 = '.7699115044247787610619469026548672566371681415929204E82 ' // & '- .5221238938053097345132743362831858407079646017699115E82 i' CALL zmst2m(st1,ma) CALL zmlg10(ma,ma) st2 = '81.9686026955271450284712410501636380322561697190529057 ' // & '- 0.2588039014625211035392823012785304771809982053965284 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-48,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMLG10',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF RETURN 90000 FORMAT (/' Testing logarithm routines.') END SUBROUTINE test5 SUBROUTINE test6(ma,mb,mc,md,mafm,mbfm,st1,st2,ncase,nerror,klog) ! Test trigonometric functions. ! .. Parameters .. INTEGER, PARAMETER :: nbits = 64, ndigmx = 256 INTEGER, PARAMETER :: lpack = (ndigmx+1)/2 + 1 INTEGER, PARAMETER :: lpackz = 2*lpack + 1 INTEGER, PARAMETER :: lunpck = (6*ndigmx)/5 + 20 INTEGER, PARAMETER :: lunpkz = 2*lunpck + 1 INTEGER, PARAMETER :: kptimp = lpack + 1 INTEGER, PARAMETER :: kptimu = lunpck + 1 INTEGER, PARAMETER :: ljsums = 8*(lunpck+2) INTEGER, PARAMETER :: lmbuff = ((lunpck+3)*(nbits-1)*301)/2000 + 6 INTEGER, PARAMETER :: lmbufz = 2*lmbuff + 10 INTEGER, PARAMETER :: lmwa = 2*lunpck ! .. ! .. Scalar Arguments .. INTEGER :: klog, ncase, nerror CHARACTER (160) :: st1, st2 ! .. ! .. Array Arguments .. REAL (KIND(0.0D0)) :: ma(0:lunpkz), mafm(0:lunpck), mb(0:lunpkz), mbfm(0:lunpck), & mc(0:lunpkz), md(0:lunpkz) ! .. ! .. External Functions .. LOGICAL, EXTERNAL :: fmcomp ! .. ! .. External Subroutines .. EXTERNAL errprt, fmi2m, fmipwr, zmabs, zmcos, zmcssn, zmsin, zmst2m, & zmsub, zmtan ! .. ! .. Scalars in Common .. REAL (KIND(0.0D0)) :: mbase INTEGER :: jform1, jform2, kdebug, keswch, kflag, krad, kround, kswide, & kw, kwarn, lvltrc, ndig, ntrace ! .. ! .. Common Blocks .. COMMON /fmuser/mbase, ndig, jform1, jform2, krad, kw, ntrace, lvltrc, & kflag, kwarn, kround, kswide, keswch, kdebug ! .. WRITE (kw,90000) ncase = 32 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmcos(ma,ma) st2 = '0.8180802525254482451348613286211514555816444253416895 ' // & '+ 0.3801751200076938035500853542125525088505055292851393 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMCOS ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 33 st1 = '34.7699115044247787610619469026548672566371681415929204 ' // & '- 42.5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmcos(ma,ma) st2 = '-1432925478410268113.5816466154230974355002592549420099 ' // & '- 309002816679456015.00151246245263842483282458519462258 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-31,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMCOS ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 34 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmsin(ma,ma) st2 = '0.7931260548991613428648822413402447097755865697557818 ' // & '- 0.3921366045897070762848927655743167937790944353110710 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMSIN ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 35 st1 = '34.7699115044247787610619469026548672566371681415929204 ' // & '- 42.5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmsin(ma,ma) st2 = '-3.090028166794560150015124624526384249047272360765358E17 ' // & '+ 1.432925478410268113581646615423097435166828182950161E18 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-31,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMSIN ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 36 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmtan(ma,ma) st2 = '0.6141156219447569167198437040270236055089243090199979 ' // & '- 0.7647270337230070156308196055474639461102792169274526 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMTAN ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 37 st1 = '35.7699115044247787610619469026548672566371681415929204 ' // & '- 43.5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmtan(ma,ma) st2 = '2.068934241218867332441292427642153175237611151321340E-38 ' // & '- 1.000000000000000000000000000000000000023741659169354 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-49,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMTAN ',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 38 st1 = '0.3505154639175257731958762886597938144329896907216495 ' // & '+ 0.7319587628865979381443298969072164948453608247422680 i' CALL zmst2m(st1,ma) CALL zmcssn(ma,ma,mc) st2 = '1.2022247452809115256533054407001508718694617802593324 ' // & '- 0.2743936538120352873902095801531325075994392065668943 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-49,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMCSSN',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 39 st1 = '0.3505154639175257731958762886597938144329896907216495 ' // & '+ 0.7319587628865979381443298969072164948453608247422680 i' CALL zmst2m(st1,ma) CALL zmcssn(ma,mc,ma) st2 = '0.4395486978082638069281369170831952476351663772871008 ' // & '+ 0.7505035100906417134864779281080728222900154610025883 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMCSSN',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF RETURN 90000 FORMAT (/' Testing trigonometric routines.') END SUBROUTINE test6 SUBROUTINE test7(ma,mb,mc,md,mafm,mbfm,st1,st2,ncase,nerror,klog) ! Test inverse trigonometric functions. ! .. Parameters .. INTEGER, PARAMETER :: nbits = 64, ndigmx = 256 INTEGER, PARAMETER :: lpack = (ndigmx+1)/2 + 1 INTEGER, PARAMETER :: lpackz = 2*lpack + 1 INTEGER, PARAMETER :: lunpck = (6*ndigmx)/5 + 20 INTEGER, PARAMETER :: lunpkz = 2*lunpck + 1 INTEGER, PARAMETER :: kptimp = lpack + 1 INTEGER, PARAMETER :: kptimu = lunpck + 1 INTEGER, PARAMETER :: ljsums = 8*(lunpck+2) INTEGER, PARAMETER :: lmbuff = ((lunpck+3)*(nbits-1)*301)/2000 + 6 INTEGER, PARAMETER :: lmbufz = 2*lmbuff + 10 INTEGER, PARAMETER :: lmwa = 2*lunpck ! .. ! .. Scalar Arguments .. INTEGER :: klog, ncase, nerror CHARACTER (160) :: st1, st2 ! .. ! .. Array Arguments .. REAL (KIND(0.0D0)) :: ma(0:lunpkz), mafm(0:lunpck), mb(0:lunpkz), mbfm(0:lunpck), & mc(0:lunpkz), md(0:lunpkz) ! .. ! .. External Functions .. LOGICAL, EXTERNAL :: fmcomp ! .. ! .. External Subroutines .. EXTERNAL errprt, fmi2m, fmipwr, zmabs, zmacos, zmasin, zmatan, zmst2m, & zmsub ! .. ! .. Scalars in Common .. REAL (KIND(0.0D0)) :: mbase INTEGER :: jform1, jform2, kdebug, keswch, kflag, krad, kround, kswide, & kw, kwarn, lvltrc, ndig, ntrace ! .. ! .. Common Blocks .. COMMON /fmuser/mbase, ndig, jform1, jform2, krad, kw, ntrace, lvltrc, & kflag, kwarn, kround, kswide, keswch, kdebug ! .. WRITE (kw,90000) ncase = 40 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmacos(ma,ma) st2 = '0.8797127900868121872960714368309657795959216549012347 ' // & '+ 0.6342141347945396859119941874681961111936156338608130 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMACOS',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 41 st1 = '.7699115044247787610619469026548672566371681415929204E12 ' // & '- .5221238938053097345132743362831858407079646017699115E12 i' CALL zmst2m(st1,ma) CALL zmacos(ma,ma) st2 = '0.5959180055163009910007767810953294528367807973983794 ' // & '+28.2518733312491023865118844008522768856672089946951468 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-48,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMACOS',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 42 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmasin(ma,ma) st2 = '0.6910835367080844319352502548087856625026630447863182 ' // & '- 0.6342141347945396859119941874681961111936156338608130 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMASIN',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 43 st1 = '.7699115044247787610619469026548672566371681415929204E13 ' // & '- .5221238938053097345132743362831858407079646017699115E13 i' CALL zmst2m(st1,ma) CALL zmasin(ma,ma) st2 = '0.9748783212785956282305451762549693982010148111568094 ' // & '-30.5544584242431480705298759613446206186670533428066404 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-48,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMASIN',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 44 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmatan(ma,ma) st2 = '0.7417952692265900376512911713942700568648670953521258 ' // & '- 0.3162747143126729004878357203292329539837025170484857 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMATAN',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 45 st1 = '.7699115044247787610619469026548672566371681415929204E13 ' // & '- .5221238938053097345132743362831858407079646017699115E13 i' CALL zmst2m(st1,ma) CALL zmatan(ma,ma) st2 = ' 1.570796326794807650905529836436131532596233124329403 ' // & '-6.033484162895927601809954710695221401671437742867605E-14 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-49,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMATAN',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF RETURN 90000 FORMAT (/' Testing inverse trigonometric routines.') END SUBROUTINE test7 SUBROUTINE test8(ma,mb,mc,md,mafm,mbfm,st1,st2,ncase,nerror,klog) ! Test hyperbolic functions. ! .. Parameters .. INTEGER, PARAMETER :: nbits = 64, ndigmx = 256 INTEGER, PARAMETER :: lpack = (ndigmx+1)/2 + 1 INTEGER, PARAMETER :: lpackz = 2*lpack + 1 INTEGER, PARAMETER :: lunpck = (6*ndigmx)/5 + 20 INTEGER, PARAMETER :: lunpkz = 2*lunpck + 1 INTEGER, PARAMETER :: kptimp = lpack + 1 INTEGER, PARAMETER :: kptimu = lunpck + 1 INTEGER, PARAMETER :: ljsums = 8*(lunpck+2) INTEGER, PARAMETER :: lmbuff = ((lunpck+3)*(nbits-1)*301)/2000 + 6 INTEGER, PARAMETER :: lmbufz = 2*lmbuff + 10 INTEGER, PARAMETER :: lmwa = 2*lunpck ! .. ! .. Scalar Arguments .. INTEGER :: klog, ncase, nerror CHARACTER (160) :: st1, st2 ! .. ! .. Array Arguments .. REAL (KIND(0.0D0)) :: ma(0:lunpkz), mafm(0:lunpck), mb(0:lunpkz), mbfm(0:lunpck), & mc(0:lunpkz), md(0:lunpkz) ! .. ! .. External Functions .. LOGICAL, EXTERNAL :: fmcomp ! .. ! .. External Subroutines .. EXTERNAL errprt, fmi2m, fmipwr, zmabs, zmchsh, zmcosh, zmsinh, zmst2m, & zmsub, zmtanh ! .. ! .. Scalars in Common .. REAL (KIND(0.0D0)) :: mbase INTEGER :: jform1, jform2, kdebug, keswch, kflag, krad, kround, kswide, & kw, kwarn, lvltrc, ndig, ntrace ! .. ! .. Common Blocks .. COMMON /fmuser/mbase, ndig, jform1, jform2, krad, kw, ntrace, lvltrc, & kflag, kwarn, kround, kswide, keswch, kdebug ! .. WRITE (kw,90000) ncase = 46 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmcosh(ma,ma) st2 = '1.1365975275870879962259716562608779977957563621412079 ' // & '- 0.4230463404769118342540441830446134405410543954181579 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-49,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMCOSH',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 47 st1 = '34.7699115044247787610619469026548672566371681415929204 ' // & '- 42.5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmcosh(ma,ma) st2 = '69552104658681.7558589320148420094288419217262200765435 ' // & '+ 626163773308016.884007302915197616300902876551542156676 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-35,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMCOSH',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 48 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmsinh(ma,ma) st2 = '0.7352399228186907963608272785465108877302444847897922 ' // & '- 0.6539816592078560369158600079981127012552558121707655 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMSINH',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 49 st1 = '34.7699115044247787610619469026548672566371681415929204 ' // & '- 42.5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmsinh(ma,ma) st2 = '6.955210465868175585893201484192181376093291191637290E 13 ' // & '+ 6.261637733080168840073029151984050820616907795167046E 14 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-35,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMSINH',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 50 st1 = '.7699115044247787610619469026548672566371681415929204 ' // & '- .5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmtanh(ma,ma) st2 = '0.7562684782933185240709480231996041186654551038993505 ' // & '- 0.2938991498221693198532255749292372853685311106820169 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMTANH',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 51 st1 = '35.7699115044247787610619469026548672566371681415929204 ' // & '- 43.5221238938053097345132743362831858407079646017699115 i' CALL zmst2m(st1,ma) CALL zmtanh(ma,ma) st2 = '9.999999999999999999999999999998967653135180689424497E-01 ' // & '+ 1.356718776492102400812550018433337461876455254467192E-31 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMTANH',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 52 st1 = '0.3505154639175257731958762886597938144329896907216495 ' // & '+ 0.7319587628865979381443298969072164948453608247422680 i' CALL zmst2m(st1,ma) CALL zmchsh(ma,ma,mc) st2 = '0.7900326499280864816444807620997665088044412803737969 ' // & '+ 0.2390857359988804105051429301542214823277594407302781 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMCHSH',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF ncase = 53 st1 = '0.3505154639175257731958762886597938144329896907216495 ' // & '+ 0.7319587628865979381443298969072164948453608247422680 i' CALL zmst2m(st1,ma) CALL zmchsh(ma,mc,ma) st2 = '0.2661087555034471983220879532235334422670297141428191 ' // & '+ 0.7098057980612199357870532628105009808447460332437714 i' CALL zmst2m(st2,mc) CALL zmsub(ma,mc,md) CALL zmabs(md,mafm) CALL fmi2m(10,mbfm) CALL fmipwr(mbfm,-50,mbfm) IF ( .NOT. fmcomp(mafm,'LE',mbfm)) THEN CALL errprt('ZMCHSH',ma,'MA',mc,'MC',md,'MD',ncase,nerror,klog) END IF RETURN 90000 FORMAT (/' Testing hyperbolic routines.') END SUBROUTINE test8 SUBROUTINE errprt(nrout,m1,name1,m2,name2,m3,name3,ncase,nerror,klog) ! Print error messages. ! M1 is the value to be tested, as computed by the routine named NROUT. ! M2 is the reference value, usually converted using ZMST2M. ! M3 is ABS(M1-M2), and ERRPRT is called if this is too big. ! NAME1,NAME2,NAME3 are strings identifying which variables in main ! correspond to M1,M2,M3. ! .. Parameters .. INTEGER, PARAMETER :: nbits = 64, ndigmx = 256 INTEGER, PARAMETER :: lpack = (ndigmx+1)/2 + 1 INTEGER, PARAMETER :: lpackz = 2*lpack + 1 INTEGER, PARAMETER :: lunpck = (6*ndigmx)/5 + 20 INTEGER, PARAMETER :: lunpkz = 2*lunpck + 1 INTEGER, PARAMETER :: kptimp = lpack + 1 INTEGER, PARAMETER :: kptimu = lunpck + 1 INTEGER, PARAMETER :: ljsums = 8*(lunpck+2) INTEGER, PARAMETER :: lmbuff = ((lunpck+3)*(nbits-1)*301)/2000 + 6 INTEGER, PARAMETER :: lmbufz = 2*lmbuff + 10 INTEGER, PARAMETER :: lmwa = 2*lunpck ! .. ! .. Scalar Arguments .. INTEGER :: klog, ncase, nerror CHARACTER (2) :: name1, name2, name3 CHARACTER (6) :: nrout ! .. ! .. Array Arguments .. REAL (KIND(0.0D0)) :: m1(0:lunpkz), m2(0:lunpkz), m3(0:lunpkz) ! .. ! .. Local Scalars .. INTEGER :: kwsave ! .. ! .. External Subroutines .. EXTERNAL zmprnt ! .. ! .. Scalars in Common .. REAL (KIND(0.0D0)) :: mbase INTEGER :: jform1, jform2, kdebug, keswch, kflag, krad, kround, kswide, & kw, kwarn, lvltrc, ndig, ntrace ! .. ! .. Common Blocks .. COMMON /fmuser/mbase, ndig, jform1, jform2, krad, kw, ntrace, lvltrc, & kflag, kwarn, kround, kswide, keswch, kdebug ! .. nerror = nerror + 1 WRITE (kw,90000) ncase, nrout WRITE (klog,90000) ncase, nrout ! Temporarily change KW to KLOG so ZMPRNT ! will write to the log file. kwsave = kw kw = klog WRITE (klog,90010) name1 CALL zmprnt(m1) WRITE (klog,90010) name2 CALL zmprnt(m2) WRITE (klog,90010) name3 CALL zmprnt(m3) kw = kwsave RETURN 90000 FORMAT (//' Error in case',I3,'. The routine',' being tested was ',A6) 90010 FORMAT (1X,A2,' =') END SUBROUTINE errprt SHAR_EOF fi # end of overwriting check if test -f 'driver2.f90' then echo shar: will not over-write existing file "'driver2.f90'" else cat << SHAR_EOF > 'driver2.f90' PROGRAM sample ! David M. Smith 9-17-96 ! This is a test program for ZMLIB 1.1, a multiple-precision real ! arithmetic package. A few example ZM calculations are carried ! out using 30 significant digit precision. ! This program uses both ZMLIB.f90 and FMLIB.f90. ! The output is saved in file ZMSAMPLE.LOG. A comparison file, ! ZMSAMPLE.CHK, is provided showing the expected output from 32-bit ! (IEEE arithmetic) machines. When run on other computers, all the ! numerical results should still be the same, but the number of terms ! needed for some of the results might be slightly different. The ! program checks all the results and the last line of the log file ! should be "All results were ok." !----------------------------------------------------------------------- ! These five common blocks contain information that must be saved ! between calls, so they should be declared in the main program. ! The parameter statement defines array sizes and pointers, and ! contains the FMLIB parameters, followed by ZMLIB parameters. !----------------------------------------------------------------------- ! .. Parameters .. INTEGER, PARAMETER :: lhash1 = 0, lhash2 = 256, nbits = 64, ndigmx = 256 INTEGER, PARAMETER :: lpack = (ndigmx+1)/2 + 1 INTEGER, PARAMETER :: lpackz = 2*lpack + 1 INTEGER, PARAMETER :: lunpck = (6*ndigmx)/5 + 20 INTEGER, PARAMETER :: lunpkz = 2*lunpck + 1 INTEGER, PARAMETER :: kptimp = lpack + 1 INTEGER, PARAMETER :: kptimu = lunpck + 1 INTEGER, PARAMETER :: ljsums = 8*(lunpck+2) INTEGER, PARAMETER :: lmbuff = ((lunpck+3)*(nbits-1)*301)/2000 + 6 INTEGER, PARAMETER :: lmbufz = 2*lmbuff + 10 INTEGER, PARAMETER :: lmwa = 2*lunpck ! .. ! .. Local Scalars .. INTEGER :: iter, k, klog, nerror ! Character string used for input and output. CHARACTER (80) :: st1 ! Declare arrays for ZM variables. All are in ! unpacked format. ! .. ! .. Local Arrays .. REAL (KIND(0.0D0)) :: ma(0:lunpkz), mafm(0:lunpck), mb(0:lunpkz), mbfm(0:lunpck), & mc(0:lunpkz), md(0:lunpkz) ! .. ! .. External Functions .. LOGICAL, EXTERNAL :: fmcomp ! .. ! .. External Subroutines .. EXTERNAL fmst2m, zmabs, zmadd, zmaddi, zmdiv, zmdivi, zmeq, zmform, & zmi2m, zmmpy, zmmpyi, zmset, zmst2m, zmsub ! .. ! .. Scalars in Common .. REAL :: alogm2, alogmb, alogmt, alogmx, runkno, spmax REAL (KIND(0.0D0)) :: dlogeb, dlogmb, dlogpi, dlogtn, dlogtp, dlogtw, dpeps, & dpmax, dppi REAL (KIND(0.0D0)) :: maxint, mbase, mblogs, mbse, mbslb, mbsli, & mbspi, mexpab, mexpov, mexpun, munkno, mxbase, mxexp, mxexp2 INTEGER :: intmax, iunkno, jform1, jform2, jformz, jprntz, kaccsw, & kdebug, keswch, kflag, krad, kround, ksub, kswide, kw, kwarn, lvltrc, & ncall, ndg2mx, ndig, ndige, ndiglb, ndigli, ndigpi, ngrd21, ngrd22, & ngrd52, ntrace CHARACTER (1) :: cmchar ! .. ! .. Arrays in Common .. REAL (KIND(0.0D0)) :: mesav(0:lunpck), mlbsav(0:lunpck), mln1(0:lunpck), & mln2(0:lunpck), mln3(0:lunpck), mln4(0:lunpck), mpisav(0:lunpck), & mwa(lmwa) INTEGER :: khasht(lhash1:lhash2), khashv(lhash1:lhash2) CHARACTER (1) :: cmbuff(lmbuff) CHARACTER (6) :: namest(0:50) ! .. ! .. Common Blocks .. COMMON /fm/mwa, ncall, kaccsw, mxexp, mxexp2, mexpun, mexpov, munkno, & iunkno, runkno, mxbase, ndg2mx, spmax, dpmax, maxint, intmax, ksub COMMON /fmbuff/cmbuff, namest, cmchar COMMON /fmsave/ndigpi, ndige, ndiglb, ndigli, mbspi, mbse, mbslb, mbsli, & mpisav, mesav, mlbsav, mln1, mln2, mln3, mln4, mblogs, mexpab, alogmb, & alogm2, alogmx, alogmt, dlogmb, dlogtn, dlogtw, dlogtp, dlogpi, dppi, & dpeps, dlogeb, khasht, khashv, ngrd21, ngrd52, ngrd22 COMMON /fmuser/mbase, ndig, jform1, jform2, krad, kw, ntrace, lvltrc, & kflag, kwarn, kround, kswide, keswch, kdebug COMMON /zmuser/jformz, jprntz ! .. ! Set precision to give at least 30 significant digits ! and initialize both the ZMLIB and FMLIB packages. ! Note that any program using the ZM package MUST call ! ZMSET before using the package. CALL zmset(30) nerror = 0 ! Write output to the standard FM output (unit KW, defined ! in subroutine FMSET), and also to the file ZMSAMPLE.LOG. klog = 18 OPEN (klog,file='ZMSAMPLE.LOG') ! 1. Find a complex root of the equation ! f(x) = x**5 - 3x**4 + x**3 - 4x**2 + x - 6 = 0. ! Newton's method with initial guess x = .56 + 1.06 i. ! This version is not tuned for speed. See the ZMSQRT ! routine for possible ways to increase speed. ! Horner's rule is used to evaluate the function: ! f(x) = ((((x-3)*x+1)*x-4)*x+1)*x-6. ! MA is the previous iterate. ! MB is the current iterate. CALL zmst2m('.56 + 1.06 i',ma) ! Print the first iteration. WRITE (kw,90000) WRITE (klog,90000) CALL zmform('F32.30','F32.30',ma,st1) WRITE (kw,90010) 0, st1(1:69) WRITE (klog,90010) 0, st1(1:69) DO 10 iter = 1, 10 ! MC is f(MA). 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) ! MD is f'(MA). 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) CALL zmdiv(mc,md,mb) CALL zmsub(ma,mb,mb) ! Print each iteration. CALL zmform('F32.30','F32.30',mb,st1) WRITE (kw,90010) iter, st1(1:69) WRITE (klog,90010) iter, st1(1:69) ! Stop iterating if MA and MB agree to over ! 30 places. CALL zmsub(ma,mb,md) CALL zmabs(md,mafm) ! The ABS result is real -- do a real (FM) compare. CALL fmst2m('1.0E-31',mbfm) IF (fmcomp(mafm,'LT',mbfm)) GO TO 20 ! Set MA = MB for the next iteration. CALL zmeq(mb,ma) 10 CONTINUE ! Check the answer. 20 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,90020) WRITE (klog,90020) END IF ! 2. Compute Exp(1.23-2.34i). ! Use the direct Taylor series. See the ZMEXP routine ! for a faster way to get Exp(x). ! MA is x. ! MB is the current term, x**n/n!. ! MC is the current partial sum. CALL zmst2m('1.23-2.34i',ma) CALL zmi2m(1,mb) CALL zmeq(mb,mc) DO 30 k = 1, 100 CALL zmmpy(mb,ma,mb) CALL zmdivi(mb,k,mb) CALL zmadd(mc,mb,mc) ! Test for convergence. KFLAG will be 1 if the result ! of the last add or subtract is the same as one of the ! input arguments. IF (kflag==1) THEN WRITE (kw,90030) k WRITE (klog,90030) k GO TO 40 END IF 30 CONTINUE ! Print the result. 40 CALL zmform('F33.30','F32.30',mc,st1) WRITE (kw,90040) st1(1:70) WRITE (klog,90040) st1(1:70) ! Check the answer. 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,90050) WRITE (klog,90050) END IF IF (nerror==0) THEN WRITE (kw,90060) ' All results were ok.' WRITE (klog,90060) ' All results were ok.' END IF STOP 90000 FORMAT (//' Sample 1. Find a root of f(x) = x**5 - 3x**4 + ', & 'x**3 - 4x**2 + x - 6 = 0.'///' Iteration Newton Approximation') 90010 FORMAT (/I6,4X,A) 90020 FORMAT (/' Error in sample case number 1.'/) 90030 FORMAT (///' Sample 2.',8X,I5,' terms were added to get ', & 'Exp(1.23-2.34i)'/) 90040 FORMAT (' Result= ',A) 90050 FORMAT (/' Error in sample case number 2.'/) 90060 FORMAT (//A/) END PROGRAM sample SHAR_EOF fi # end of overwriting check if test -f 'driver3.f90' then echo shar: will not over-write existing file "'driver3.f90'" else cat << SHAR_EOF > 'driver3.f90' PROGRAM test ! David M. Smith 6-14-96 ! This is a test program for FMLIB 1.1, a multiple-precision real ! arithmetic package. Most of the FM (floating-point) routines ! are tested, and the results are checked to 50 significant digits. ! Most of the IM (integer) routines are tested, with exact results ! required to pass the tests. ! This program uses FMLIB.f90. ! The four common blocks contain information that must be saved ! between calls, so they should be declared in the main program. ! The parameter statement defines various array sizes. ! .. Parameters .. INTEGER, PARAMETER :: lhash1 = 0, lhash2 = 256, nbits = 64, ndigmx = 256 INTEGER, PARAMETER :: lpack = (ndigmx+1)/2 + 1 INTEGER, PARAMETER :: lunpck = (6*ndigmx)/5 + 20 INTEGER, PARAMETER :: ljsums = 8*(lunpck+2) INTEGER, PARAMETER :: lmbuff = ((lunpck+3)*(nbits-1)*301)/2000 + 6 INTEGER, PARAMETER :: lmwa = 2*lunpck ! .. ! .. Local Scalars .. INTEGER :: klog, ncase, nerror ! Character strings used for input and output. CHARACTER (80) :: st1, st2 ! Declare arrays for FM variables. All are in ! unpacked format. ! .. ! .. Local Arrays .. REAL (KIND(0.0D0)) :: ma(0:lunpck), mb(0:lunpck), mc(0:lunpck), md(0:lunpck) ! .. ! .. External Subroutines .. EXTERNAL fmset, test1, test10, test11, test12, test13, test14, test15, & test2, test3, test4, test5, test6, test7, test8, test