C C FM 1.1 David M. Smith 5-19-97 C C C C The FM routines in this package perform floating-point C multiple-precision arithmetic, and the IM routines perform C integer multiple-precision arithmetic. C C C C 1. INITIALIZING THE PACKAGE C C Before calling any routine in the package, several variables in C the common blocks /FMUSER/, /FM/, /FMBUFF/, and /FMSAVE/ must be C initialized. These four common blocks contain information that C is saved between calls, so they should be declared in the main C program. C C Subroutine FMSET initializes these variables to default values and C defines all machine-dependent values in the package. After calling C FMSET once at the start of a program, the user may sometimes want C to reset some of the variables in these common blocks. These C variables are described below. C C C C 2. REPRESENTATION OF FM NUMBERS C C MBASE is the base in which the arithmetic is done. MBASE must be C bigger than one, and less than or equal to the square root of C the largest representable integer. For best efficiency MBASE C should be large, but no more than about 1/4 of the square root C of the largest representable integer. Input and output C conversions are much faster when MBASE is a power of ten. C C NDIG is the number of base MBASE digits that are carried in the C multiple precision numbers. NDIG must be at least two. The C upper limit for NDIG is defined in the PARAMETER statement at C the top of each routine and is restricted only by the amount C of memory available. C C Sometimes it is useful to dynamically vary NDIG during the program. C Use FMEQU to round numbers to lower precision or zero-pad them to C higher precision when changing NDIG. C C It is rare to need to change MBASE during a program. Use FMCONS to C reset some saved constants that depend on MBASE. FMCONS should be C called immediately after changing MBASE. C C There are two representations for a floating multiple precision C number. The unpacked representation used by the routines while C doing the computations is base MBASE and is stored in NDIG+2 words. C A packed representation is available to store the numbers in the C user's program in compressed form. In this format, the NDIG C (base MBASE) digits of the mantissa are packed two per word to C conserve storage. Thus the external, packed form of a number C requires (NDIG+1)/2+2 words. C C This version uses double precision arrays to hold the numbers. C Version 1.0 of FM used integer arrays, which are faster on some C machines. The package can easily be changed to use integer C arrays -- see section 11 on EFFICIENCY below. C C The unpacked format of a floating multiple precision number is as C follows. A number MA is kept in an array with MA(1) containing C the exponent and MA(2) through MA(NDIG+1) containing one digit of C the mantissa, expressed in base MBASE. The array is dimensioned C to start at MA(0), with the approximate number of bits of precision C stored in MA(0). This precision value is intended to be used by FM C functions that need to monitor cancellation error in addition and C subtraction. The cancellation monitor code is usually disabled for C user calls, and FM functions only check for cancellation when they C must. Tracking cancellation causes most routines to run slower, C with addition and subtraction being affected the most. C C The exponent is a power of MBASE and the implied radix point is C immediately before the first digit of the mantissa. Every nonzero C number is normalized so that the second array element (the first C digit of the mantissa) is nonzero. C C In both representations the sign of the number is carried on the C second array element only. Elements 3,4,... are always nonnegative. C The exponent is a signed integer and may be as large in magnitude as C MXEXP (defined in FMSET). C C For MBASE = 10,000 and NDIG = 4, the number -pi would have these C representations: C Word 1 2 3 4 5 C C Unpacked: 1 -3 1415 9265 3590 C Packed: 1 -31415 92653590 C C Word 0 would be 42 in both formats, indicating that the mantissa C has about 42 bits of precision. C C Because of normalization in a large base, the equivalent number C of base 10 significant digits for an FM number may be as small as C LOG10(MBASE)*(NDIG-1) + 1. C C The integer routines use the FMLIB format to represent numbers, C without the number of digits (NDIG) being fixed. Integers in IM C format are essentially variable precision, using the minimum number C of words to represent each value. C C For programs using both FM and IM numbers, FM routines should not C be called with IM numbers, and IM routines should not be called C with FM numbers, since the implied value of NDIG used for an IM C number may not match the explicit NDIG expected by an FM routine. C Use the conversion routines IMFM2I and IMI2FM to change between C the FM and IM formats. C C C C 3. INPUT/OUTPUT ROUTINES C C All versions of the input routines perform free-format conversion C from characters to FM numbers. C C a. Conversion to or from a character array C C FMINP converts from a character*1 array to an FM number. C C FMOUT converts an FM number to base 10 and formats it for output C as an array of type character*1. The output is left C justified in the array, and the format is defined by two C variables in common, so that a separate format definition C does not have to be provided for each output call. C C The user sets JFORM1 and JFORM2 to determine the output format. C C JFORM1 = 0 E format ( .314159M+6 ) C = 1 1PE format ( 3.14159M+5 ) C = 2 F format ( 314159.000 ) C C JFORM2 is the number of significant digits to display (if C JFORM1 = 0 or 1). If JFORM2.EQ.0 then a default number C of digits is chosen. The default is roughly the full C precision of the number. C JFORM2 is the number of digits after the decimal point (if C JFORM1 = 2). See the FMOUT documentation for more details. C C b. Conversion to or from a character string C C FMST2M converts from a character string to an FM number. C C FMFORM converts an FM number to a character string according to C a format provided in each call. The format description C is more like that of a Fortran FORMAT statement, and C integer or fixed-point output is right justified. C C c. Direct read or write C C FMPRNT uses FMOUT to print one FM number. C C FMFPRT uses FMFORM to print one FM number. C C FMWRIT writes FM numbers for later input using FMREAD. C C FMREAD reads FM numbers written by FMWRIT. C C The values given to JFORM1 and JFORM2 can be used to define a C default output format when FMOUT or FMPRNT are called. The C explicit format used in a call to FMFORM or FMFPRT overrides C the settings of JFORM1 and JFORM2. C C KW is the unit number to be used for standard output from C the package, including error and warning messages, and C trace output. C C For multiple precision integers, the corresponding routines C IMINP, IMOUT, IMST2M, IMFORM, IMPRNT, IMFPRT, IMWRIT, and C IMREAD provide similar input and output conversions. For C output of IM numbers, JFORM1 and JFORM2 are ignored and C integer format (JFORM1=2, JFORM2=0) is used. C C For further description of these routines, see sections C 9 and 10 below. C C C C 4. ARITHMETIC TRACING C C NTRACE and LVLTRC control trace printout from the package. C C NTRACE = 0 No printout except warnings and errors. C = 1 The result of each call to one of the routines C is printed in base 10, using FMOUT. C = -1 The result of each call to one of the routines C is printed in internal base MBASE format. C = 2 The input arguments and result of each call to one C of the routines is printed in base 10, using FMOUT. C = -2 The input arguments and result of each call to one C of the routines is printed in base MBASE format. C C LVLTRC defines the call level to which the trace is done. LVLTRC = 1 C means only FM routines called directly by the user are traced, C LVLTRC = 2 also prints traces for FM routines called by other C FM routines called directly by the user, etc. C C In the above description, internal MBASE format means the number is C printed as it appears in the array --- an exponent followed by NDIG C base MBASE digits. C C C C 5. ERROR CONDITIONS C C KFLAG is a condition parameter returned by the package after each C call to one of the routines. Negative values indicate C conditions for which a warning message will be printed C unless KWARN = 0. Positive values indicate conditions C that may be of interest but are not errors. C No warning message is printed if KFLAG is nonnegative. C C KFLAG = 0 Normal operation. C C = 1 One of the operands in FMADD or FMSUB was C insignificant with respect to the other, so C that the result was equal to the argument of C larger magnitude. C = 2 In converting an FM number to a one word integer C in FMM2I, the FM number was not exactly an C integer. The next integer toward zero was C returned. C C = -1 NDIG was less than 2 or more than NDIGMX. C = -2 MBASE was less than 2 or more than MXBASE. C = -3 An exponent was out of range. C = -4 Invalid input argument(s) to an FM routine. C UNKNOWN was returned. C = -5 + or - OVERFLOW was generated as a result from an C FM routine. C = -6 + or - UNDERFLOW was generated as a result from an C FM routine. C = -7 The input string (array) to FMINP was not legal. C = -8 The character array was not large enough in an C input or output routine. C = -9 Precision could not be raised enough to provide all C requested guard digits. Increasing NDIGMX in C all the PARAMETER statements may fix this. C UNKNOWN was returned. C = -10 An FM input argument was too small in magnitude to C convert to the machine's single or double C precision in FMM2SP or FMM2DP. Check that the C definitions of SPMAX and DPMAX in FMSET are C correct for the current machine. C Zero was returned. C C When a negative KFLAG condition is encountered, the value of KWARN C determines the action to be taken. C C KWARN = 0 Execution continues and no message is printed. C = 1 A warning message is printed and execution continues. C = 2 A warning message is printed and execution stops. C C The default setting is KWARN = 1. C C When an overflow or underflow is generated for an operation in which C an input argument was already an overflow or underflow, no additional C message is printed. When an unknown result is generated and an input C argument was already unknown, no additional message is printed. In C these cases the negative KFLAG value is still returned. C C IM routines handle exceptions like OVERFLOW or UNKNOWN in the same C way as FM routines. When using IMMPY, the product of two large C positive integers will return +OVERFLOW. The routine IMMPYM can C be used to obtain a modular result without overflow. The largest C representable IM integer is MBASE**NDIGMX - 1. For example, if C MBASE is 10**7 and NDIGMX is set to 256, integers less than 10**1792 C can be used. C C C C 6. OTHER PARAMETERS C C KRAD = 0 All angles in the trigonometric functions and C inverse functions are measured in degrees. C = 1 All angles are measured in radians. (Default) C C C KROUND = 0 All final results are chopped (rounded toward C zero). Intermediate results are rounded. C = 1 All results are rounded to the nearest FM C number, or to the value with an even last C digit if the result is halfway between two C FM numbers. (Default) C C KSWIDE defines the maximum screen width to be used for C all unit KW output. Default is 80. C C KESWCH controls the action taken in FMINP and other input routines C for strings like 'E7' that have no digits before the exponent C field. Default is for 'E7' to translate like '1.0E+7'. C C CMCHAR defines the exponent letter to be used for FM variable C output. Default is 'M', as in 1.2345M+678. C C KDEBUG = 0 Error checking is not done for valid input arguments C and parameters like NDIG and MBASE upon entry to C each routine. (Default) C = 1 Some error checking is done. (Slower speed) C C See FMSET for additional description of these and other variables C defining various FM conditions. C C C C 7. ARRAY DIMENSIONS C C The dimensions of the arrays in the FM package are defined using C a PARAMETER statement at the top of each routine. The size of C these arrays depends on the values of parameters NDIGMX and NBITS. C NDIGMX is the maximum value the user may set for NDIG. C NBITS is the number of bits used to represent integers for a C given machine. See the EFFICIENCY discussion below. C C The standard version of FMLIB sets NDIGMX = 256, so on a 32-bit C machine using MBASE = 10**7 the maximum precision is about C 7*255+1 = 1786 significant digits. To change dimensions so that C 10,000 significant digit calculation can be done, NDIGMX needs to C be at least 10**4/7 + 5 = 1434. This allows for a few user guard C digits to be defined when the package is initialized using C CALL FMSET(10000). Changing 'NDIGMX=256' to 'NDIGMX=1434' C everywhere in the package and the user's calling program will C define all the new array sizes. C C If NDIG much greater than 256 is to be used and elementary functions C will be needed, they will be faster if array MJSUMS is larger. The C parameter defining the size of MJSUMS is set in the standard version C by LJSUMS = 8*(LUNPCK+2). The 8 means that up to eight concurrent C sums can be used by the elementary functions. The approximate number C needed for best speed is given by the formula C 0.051*Log(MBASE)*NDIG**(1/3) + 1.85 C For example, with MBASE=10**7 and NDIG=1434 this gives 11. Changing C 'LJSUMS = 8*(LUNPCK+2)' to 'LJSUMS =11*(LUNPCK+2)' everywhere in the C package and the user's calling program will give slightly better C speed. C C FM numbers in packed format have dimension 0:LPACK, and those C in unpacked format have dimension 0:LUNPCK. C C C C 8. PORTABILITY C C In FMSET there is some machine-dependent code that attempts to C approximate the largest representable integer value. The current C code works on all machines tested, but if an FM run fails, check C the MAXINT and INTMAX loops in FMSET. Values for SPMAX and DPMAX C are also defined in FMSET that should be set to values near overflow C for single precision and double precision. Setting KDEBUG = 1 may C also identify some errors if a run fails. C C Some compilers object to a function like FMCOMP with side effects C such as changing KFLAG or other common variables. Blocks of code C in FMCOMP and IMCOMP that modify common are identified so they may C be removed or commented out to produce a function without side C effects. This disables trace printing in FMCOMP and IMCOMP, and C error codes are not returned in KFLAG. See FMCOMP and IMCOMP for C further details. C C All variables are explicitly declared in each routine. There is C a commented IMPLICIT NONE statement in each routine that can be C enabled to get more compiler diagnostic information in some testing C or debugging situations. C C C 9. LIST OF ROUTINES C C These are the FM routines that are designed to be called by C the user. All are subroutines except logical function FMCOMP. C MA, MB, MC refer to FM format numbers. C C In each case it is permissible to use the same array more than C once in the calling sequence. The statement MA = MA*MA can C be written CALL FMMPY(MA,MA,MA). C C For each of these routines there is also a version available for C which the argument list is the same but all FM numbers are in packed C format. The routines using packed numbers have the same names except C 'FM' is replaced by 'FP' at the start of each name. C C C FMABS(MA,MB) MB = ABS(MA) C C FMACOS(MA,MB) MB = ACOS(MA) C C FMADD(MA,MB,MC) MC = MA + MB C C FMADDI(MA,IVAL) MA = MA + IVAL Increment an FM number by a one C word integer. Note this call C does not have an "MB" result C like FMDIVI and FMMPYI. C C FMASIN(MA,MB) MB = ASIN(MA) C C FMATAN(MA,MB) MB = ATAN(MA) C C FMATN2(MA,MB,MC) MC = ATAN2(MA,MB) C C FMBIG(MA) MA = Biggest FM number less than overflow. C C FMCHSH(MA,MB,MC) MB = COSH(MA), MC = SINH(MA). Faster than C making two separate calls. C C FMCOMP(MA,LREL,MB) Logical comparison of MA and MB. C LREL is a CHARACTER*2 value identifying C which comparison is made. C Example: IF (FMCOMP(MA,'GE',MB)) ... C C FMCONS Set several saved constants that depend C on MBASE, the base being used. FMCONS C should be called immediately after C changing MBASE. C C FMCOS(MA,MB) MB = COS(MA) C C FMCOSH(MA,MB) MB = COSH(MA) C C FMCSSN(MA,MB,MC) MB = COS(MA), MC = SIN(MA). Faster than C making two separate calls. C C FMDIG(NSTACK,KST) Find a set of precisions to use during C Newton iteration for finding a simple C root starting with about double C precision accuracy. C C FMDIM(MA,MB,MC) MC = DIM(MA,MB) C C FMDIV(MA,MB,MC) MC = MA/MB C C FMDIVI(MA,IVAL,MB) MB = MA/IVAL IVAL is a one word integer. C C FMDP2M(X,MA) MA = X Convert from double precision to FM. C C FMDPM(X,MA) MA = X Convert from double precision to FM. C Much faster than FMDP2M, but MA agrees C with X only to D.P. accuracy. See C the comments in the two routines. C C FMEQ(MA,MB) MB = MA Both have precision NDIG. C This is the version to use for C standard B = A statements. C C FMEQU(MA,MB,NA,NB) MB = MA Version for changing precision. C MA has NA digits (i.e., MA was C computed using NDIG = NA), and MB C will be defined having NB digits. C MB is zero-padded if NB.GT.NA C MB is rounded if NB.LT.NA C C FMEXP(MA,MB) MB = EXP(MA) C C FMFORM(FORM,MA,STRING) MA is converted to a character string C using format FORM and returned in C STRING. FORM can represent I, F, C E, or 1PE formats. Example: C CALL FMFORM('F60.40',MA,STRING) C C FMFPRT(FORM,MA) Print MA on unit KW using FORM format. C C FMI2M(IVAL,MA) MA = IVAL Convert from one word integer C to FM. C C FMINP(LINE,MA,LA,LB) MA = LINE Input conversion. C Convert LINE(LA) through LINE(LB) C from characters to FM. C C FMINT(MA,MB) MB = INT(MA) Integer part of MA. C C FMIPWR(MA,IVAL,MB) MB = MA**IVAL Raise an FM number to a one C word integer power. C C FMLG10(MA,MB) MB = LOG10(MA) C C FMLN(MA,MB) MB = LOG(MA) C C FMLNI(IVAL,MA) MA = LOG(IVAL) Natural log of a one word C integer. C C FMM2DP(MA,X) X = MA Convert from FM to double precision. C C FMM2I(MA,IVAL) IVAL = MA Convert from FM to integer. C C FMM2SP(MA,X) X = MA Convert from FM to single precision. C C FMMAX(MA,MB,MC) MC = MAX(MA,MB) C C FMMIN(MA,MB,MC) MC = MIN(MA,MB) C C FMMOD(MA,MB,MC) MC = MA mod MB C C FMMPY(MA,MB,MC) MC = MA*MB C C FMMPYI(MA,IVAL,MB) MB = MA*IVAL Multiply by a one word integer. C C FMNINT(MA,MB) MB = NINT(MA) Nearest FM integer. C C FMOUT(MA,LINE,LB) LINE = MA Convert from FM to character. C LINE is a character array of C length LB. C C FMPI(MA) MA = pi C C FMPRNT(MA) Print MA on unit KW using current format. C C FMPWR(MA,MB,MC) MC = MA**MB C C FMREAD(KREAD,MA) MA is returned after reading one (possibly C multi-line) FM number on unit KREAD. This C routine reads numbers written by FMWRIT. C C FMRPWR(MA,K,J,MB) MB = MA**(K/J) Rational power. Faster than C FMPWR for functions like the cube root. C C FMSET(NPREC) Set default values and machine-dependent C variables to give at least NPREC base 10 C digits plus three base 10 guard digits. C Must be called to initialize FM package. C C FMSIGN(MA,MB,MC) MC = SIGN(MA,MB) Sign transfer. C C FMSIN(MA,MB) MB = SIN(MA) C C FMSINH(MA,MB) MB = SINH(MA) C C FMSP2M(X,MA) MA = X Convert from single precision to FM. C C FMSQR(MA,MB) MB = MA*MA Faster than FMMPY. C C FMSQRT(MA,MB) MB = SQRT(MA) C C FMST2M(STRING,MA) MA = STRING C Convert from character string to FM. C Often more convenient than FMINP, which C converts an array of CHARACTER*1 values. C Example: CALL FMST2M('123.4',MA). C C FMSUB(MA,MB,MC) MC = MA - MB C C FMTAN(MA,MB) MB = TAN(MA) C C FMTANH(MA,MB) MB = TANH(MA) C C FMULP(MA,MB) MB = One Unit in the Last Place of MA. C C FMWRIT(KWRITE,MA) Write MA on unit KWRITE. C Multi-line numbers will have '&' as the C last nonblank character on all but the last C line. These numbers can then be read C easily using FMREAD. C C C These are the integer routines that are designed to be called by C the user. All are subroutines except logical function IMCOMP. C MA, MB, MC refer to IM format numbers. In each case the version C of the routine to handle packed IM numbers has the same name, C with 'IM' replaced by 'IP'. C C IMABS(MA,MB) MB = ABS(MA) C C IMADD(MA,MB,MC) MC = MA + MB C C IMBIG(MA) MA = Biggest IM number less than overflow. C C IMCOMP(MA,LREL,MB) Logical comparison of MA and MB. C LREL is a CHARACTER*2 value identifying C which comparison is made. C Example: IF (IMCOMP(MA,'GE',MB)) ... C C IMDIM(MA,MB,MC) MC = DIM(MA,MB) C C IMDIV(MA,MB,MC) MC = int(MA/MB) C Use IMDIVR if the remainder is also needed. C C IMDIVI(MA,IVAL,MB) MB = int(MA/IVAL) C IVAL is a one word integer. Use IMDVIR C to get the remainder also. C C IMDIVR(MA,MB,MC,MD) MC = int(MA/MB), MD = MA mod MB C When both the quotient and remainder are C needed, this routine is twice as fast as C calling both IMDIV and IMMOD. C C IMDVIR(MA,IVAL,MB,IREM) MB = int(MA/IVAL), IREM = MA mod IVAL C IVAL and IREM are one word integers. C C IMEQ(MA,MB) MB = MA C C IMFM2I(MAFM,MB) MB = MAFM Convert from real (FM) format C to integer (IM) format. C C IMFORM(FORM,MA,STRING) MA is converted to a character string C using format FORM and returned in C STRING. FORM can represent I, F, C E, or 1PE formats. Example: C CALL IMFORM('I70',MA,STRING) C C IMFPRT(FORM,MA) Print MA on unit KW using FORM format. C C IMGCD(MA,MB,MC) MC = greatest common divisor of MA and MB. C C IMI2FM(MA,MBFM) MBFM = MA Convert from integer (IM) format C to real (FM) format. C C IMI2M(IVAL,MA) MA = IVAL Convert from one word integer C to IM. C C IMINP(LINE,MA,LA,LB) MA = LINE Input conversion. C Convert LINE(LA) through LINE(LB) C from characters to IM. C C IMM2DP(MA,X) X = MA Convert from IM to double precision. C C IMM2I(MA,IVAL) IVAL = MA Convert from IM to one word integer. C C IMMAX(MA,MB,MC) MC = MAX(MA,MB) C C IMMIN(MA,MB,MC) MC = MIN(MA,MB) C C IMMOD(MA,MB,MC) MC = MA mod MB C C IMMPY(MA,MB,MC) MC = MA*MB C C IMMPYI(MA,IVAL,MB) MB = MA*IVAL Multiply by a one word integer. C C IMMPYM(MA,MB,MC,MD) MD = MA*MB mod MC C Slightly faster than calling IMMPY and C IMMOD separately, and it works for cases C where IMMPY would return OVERFLOW. C C IMOUT(MA,LINE,LB) LINE = MA Convert from IM to character. C LINE is a character array of C length LB. C C IMPMOD(MA,MB,MC,MD) MD = MA**MB mod MC C C IMPRNT(MA) Print MA on unit KW. C C IMPWR(MA,MB,MC) MC = MA**MB C C IMREAD(KREAD,MA) MA is returned after reading one (possibly C multi-line) IM number on unit KREAD. This C routine reads numbers written by IMWRIT. C C IMSIGN(MA,MB,MC) MC = SIGN(MA,MB) Sign transfer. C C IMSQR(MA,MB) MB = MA*MA Faster than IMMPY. C C IMST2M(STRING,MA) MA = STRING C Convert from character string to IM. C Often more convenient than IMINP, which C converts an array of CHARACTER*1 values. C Example: CALL IMST2M('12345678901',MA). C C IMSUB(MA,MB,MC) MC = MA - MB C C IMWRIT(KWRITE,MA) Write MA on unit KWRITE. C Multi-line numbers will have '&' as the C last nonblank character on all but the last C line. These numbers can then be read C easily using IMREAD. C C Many of the IM routines call FM routines, but none of the FM C routines call IM routines, so the IM routines can be omitted C if none are called explicitly from a program. C C C C 10. NEW FOR VERSION 1.1 C C Version 1.0 used integer arrays and integer arithmetic internally C to perform the multiple precision operations. Version 1.1 uses C double precision arithmetic and arrays internally. This is usually C faster at higher precisions, and on many machines it is also faster C at lower precisions. Version 1.1 is written so that the arithmetic C used can easily be changed from double precision to integer, or any C other available arithmetic type. This permits the user to make the C best use of a given machine's arithmetic hardware. C See the EFFICIENCY discussion below. C C Several routines have undergone minor modification, but only a few C changes should affect programs that used FM 1.0. Many of the C routines are faster in version 1.1, because code has been added to C take advantage of special cases for individual functions instead of C using general formulas that are more compact. For example, there C are separate routines using series for SINH and COSH instead of C just calling EXP. C C FMEQU was the only routine that required the user to give the value C of the current precision. This was to allow automatic C rounding or zero-padding when changing precision. Since few C user calls change precision, a new routine has been added for C this case. C FMEQ now handles this case and has a simple argument list that C does not include the value of NDIG. C FMEQU is used for changing precision. C See the list of FM routines above for details. C C All variable names beginning with M in the package are now declared C as double precision, so FM common blocks in the user's program need C D.P. declarations, and FM variables (arrays) used in the calling C program need to be D.P. C C /FMUSER/ is a common block holding parameters that define the C arithmetic to be used and other user options. Several C new variables have been added, including screen width to C be used for output. See above for further description. C C /FMSAVE/ is a common block for saving constants to avoid C re-computing them. Several new variables have been added. C C /FMBUFF/ is a common block containing a character array used to C format FM numbers for output. Two new items have been C added. C C New routines: C C All the IM routines are new for version 1.1. C C FMADDI increments an FM number by a small integer. C It runs in O(1) time, on the average. C C FMCHSH returns both SINH(MA) and COSH(MA). C When both are needed, this is almost twice as fast C as making separate calls to FMCOSH and FMSINH. C C FMCSSN returns both SIN(MA) and COS(MA). C When both are needed, this is almost twice as fast C as making separate calls to FMCOS and FMSIN. C C FMFORM uses a format string to convert an FM number to a C character string. C C FMFPRT prints an FM number using a format string. C C FMREAD reads an FM number written using FMWRIT. C C FMRPWR computes an FM number raised to a rational power. For cube C roots and similar rational powers it is usually much faster C than FMPWR. C C FMSQR squares an FM number. It is faster than using FMMPY. C C FMST2M converts character strings to FM format. Since FMINP converts C character arrays, this routine can be more convenient for C easily defining an FM number. C For example, CALL FMST2M('123.4',MA). C C FMWRIT writes an FM number using a format for multi-line numbers C with '&' at the end of all but the last line of a multi-line C number. This allows automatic reading of FM numbers without C needing to know the base, precision or format under which they C were written. C C One extra word has been added to the dimensions of all FM numbers. C Word zero in each array contains a value used to monitor cancellation C error arising from addition or subtraction. This value approximates C the number of bits of precision for an FM value. It allows higher C level FM functions to detect cases where too much cancellation has C occurred. KACCSW is a switch variable in COMMON /FM/ used internally C to enable cancellation error monitoring. C C C 11. EFFICIENCY C C To take advantage of hardware architecture on different machines, the C package has been designed so that the arithmetic used to perform the C multiple precision operations can easily be changed. All variables C that must be changed to get a different arithmetic have names C beginning with 'M' and are declared using DOUBLE PRECISION M.... C C For example, to change the package to use integer arithmetic C internally, make these two changes everywhere in the package: C change 'DOUBLE PRECISION M' to 'INTEGER M', C change 'DINT(' to 'INT('. C On some systems, changing 'DINT(' to '(' may give better speed. C C When changing to a different type of arithmetic, all FM common blocks C and arrays in the user's program must be changed to agree. In a few C places in FM, where a DINT function is not supposed to be changed, it C is spelled 'DINT (' so the global change will not find it. C C This version restricts the base used to be also representable in C integer variables, so using precision above double usually does not C save much time unless integers can also be declared at a higher C precision. Using IEEE Extended would allow a base of around 10**9 C to be chosen, but the delayed digit-normalization method used for C multiplication and division means that a slightly smaller base like C 10**8 would usually run faster. This would usually not be much C faster than using 10**7 with double precision. C C The value of NBITS defined as a parameter in most FM routines C refers to the number of bits used to represent integers in an C M-variable word. Typical values for NBITS are: 24 for IEEE single C precision, 32 for integer, 53 for IEEE double precision. NBITS C controls only array size, so setting it too high is ok, but then C the program will use more memory than necessary. C C For cases where special compiler directives or minor re-writing C of the code may improve speed, several of the most important C loops in FM are identified by comments containing the string C '(Inner Loop)'. C C -------------------------------------------------------------------- C -------------------------------------------------------------------- C C SUBROUTINE FMSET(NPREC) C C Initialize the values in common that must be set before calling C other FM routines. C C Base and precision will be set to give at least NPREC+3 decimal C digits of precision (giving the user three base ten guard digits). C C MBASE is set to a large power of ten. C JFORM1 and JFORM2 are set to 1PE format displaying NPREC C significant digits. C C The trace option is set off. C The mode for angles in trig functions is set to radians. C The rounding mode is set to symmetric rounding. C Warning error message level is set to 1. C Cancellation error monitor is set off. C Screen width for output is set to 80 columns. C The exponent character for FM output is set to 'M'. C Debug error checking is set off. C C KW, the unit number for all FM output, is set to 6. C C The size of all arrays is controlled by defining two parameters: C NDIGMX is the maximum value the user can set NDIG, C NBITS is the number of bits used to represent integers in an C M-variable word. C C IMPLICIT NONE C INTEGER NDIGMX,NBITS,LPACK,LUNPCK,LMWA,LJSUMS,LMBUFF C C Define the array sizes: C PARAMETER ( NDIGMX=256 , NBITS=64 , * LPACK = (NDIGMX+1)/2 + 1 , LUNPCK = (6*NDIGMX)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*(LUNPCK+2), * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C INTEGER NPREC C C Here are all the common blocks used in FM. C C /FMUSER/, /FM/, /FMBUFF/, and /FMSAVE/ should also be declared in the C main program, because some compilers allocate and free space used for C labelled common that is declared only in subprograms. This causes C the saved information to be lost. C C FMUSER contains values that may need to be C changed by the calling program. C INTEGER NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC,KFLAG, * KWARN,KROUND,KSWIDE,KESWCH,KDEBUG DOUBLE PRECISION MBASE C COMMON /FMUSER/ MBASE,NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC, * KFLAG,KWARN,KROUND,KSWIDE,KESWCH,KDEBUG C C FM contains the work array used by the low-level C arithmetic routines, definitions for overflow C and underflow thresholds, and other C machine-dependent values. C DOUBLE PRECISION DPMAX INTEGER NCALL,KACCSW,IUNKNO,NDG2MX,INTMAX,KSUB REAL RUNKNO,SPMAX DOUBLE PRECISION MWA,MXEXP,MXEXP2,MEXPUN,MEXPOV,MUNKNO,MXBASE, * MAXINT C COMMON /FM/ MWA(LMWA),NCALL,KACCSW,MXEXP,MXEXP2,MEXPUN,MEXPOV, * MUNKNO,IUNKNO,RUNKNO,MXBASE,NDG2MX,SPMAX,DPMAX, * MAXINT,INTMAX,KSUB C C FMSAVE contains information about saved constants. C INTEGER LHASH1,LHASH2 PARAMETER (LHASH1=0 , LHASH2=256) DOUBLE PRECISION DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB INTEGER NDIGPI,NDIGE,NDIGLB,NDIGLI,KHASHT,KHASHV, * NGRD21,NGRD52,NGRD22 DOUBLE PRECISION MBSPI,MBSE,MBSLB,MBSLI,MPISAV,MESAV,MLBSAV, * MLN1,MLN2,MLN3,MLN4,MBLOGS,MEXPAB REAL ALOGMB,ALOGMT,ALOGM2,ALOGMX C COMMON /FMSAVE/ NDIGPI,NDIGE,NDIGLB,NDIGLI,MBSPI,MBSE,MBSLB,MBSLI, * MPISAV(0:LUNPCK),MESAV(0:LUNPCK),MLBSAV(0:LUNPCK), * MLN1(0:LUNPCK),MLN2(0:LUNPCK),MLN3(0:LUNPCK), * MLN4(0:LUNPCK),MBLOGS,MEXPAB,ALOGMB,ALOGM2,ALOGMX, * ALOGMT,DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB,KHASHT(LHASH1:LHASH2), * KHASHV(LHASH1:LHASH2),NGRD21,NGRD52,NGRD22 C C MJSUMS is an array that can contain several FM numbers C being used to accumulate concurrent sums in exponential C and trigonometric functions. When NDIGMX = 256, eight is C about the maximum number of sums needed (but this depends C on MBASE). For larger NDIGMX, dimensioning MJSUMS to hold C more than eight FM numbers could increase the speed of the C functions. C DOUBLE PRECISION MJSUMS C COMMON /FMSUMS/ MJSUMS(0:LJSUMS) C C FMWA contains two work arrays similar to MWA. They are C used in routines FMDIVD, FMMPYD, and FMMPYE. C DOUBLE PRECISION MWD,MWE C COMMON /FMWA/ MWD(LMWA),MWE(LMWA) C C CMBUFF is a character array used by FMPRNT for printing C output from FMOUT. This array may also be used C for calls to FMOUT from outside the FM package. C CMCHAR is the letter used before the exponent field C in FMOUT. It is defined in FMSET. C NAMEST is a stack for names of the routines. It is C used for trace printing and error messages. C CHARACTER CMBUFF,CMCHAR CHARACTER *6 NAMEST C COMMON /FMBUFF/ CMBUFF(LMBUFF),NAMEST(0:50),CMCHAR C C FM1 contains scratch arrays for temporary storage of FM C numbers while computing various functions. C DOUBLE PRECISION M01,M02,M03,M04,M05,M06 C COMMON /FM1/ M01(0:LUNPCK),M02(0:LUNPCK),M03(0:LUNPCK), * M04(0:LUNPCK),M05(0:LUNPCK),M06(0:LUNPCK) C C FMPCK contains scratch arrays used to hold input arguments C in unpacked format when the packed versions of functions C are used. C DOUBLE PRECISION MPA,MPB,MPC C COMMON /FMPCK/ MPA(0:LUNPCK),MPB(0:LUNPCK),MPC(0:LUNPCK) C DOUBLE PRECISION ONE,TEMP,TWO,YT C CHARACTER LCHARS(21) DOUBLE PRECISION ML,MLD2,MLM1 INTEGER J,K,KPT,L,LTYPES(21),LVALS(21),NPSAVE C DATA LCHARS/'+','-','0','1','2','3','4','5','6','7','8','9', * '.','E','D','Q','M','e','d','q','m'/ DATA LTYPES/1,1,10*2,3,8*4/ DATA LVALS/1,-1,0,1,2,3,4,5,6,7,8,9,9*0/ C C KW is the unit number for standard output from the C FM package. This includes trace output and error C messages. C KW = 6 C C MAXINT should be set to a very large integer, possibly C the largest representable integer for the current C machine. For most 32-bit machines, MAXINT is set C to 2**53 - 1 = 9.007D+15 when double precision C arithmetic is used for M-variables. Using integer C M-variables usually gives MAXINT = 2**31 - 1 = C 2 147 483 647. C C Setting MAXINT to a smaller number is ok, but this C unnecessarily restricts the permissible range of C MBASE and MXEXP. C C The following code should set MAXINT to the largest C representable number of the form 2**J - 1. C C The FMMSET call keeps some compilers from doing the 110 C loop at the highest precision available and then rounding C to the declared precision. C MAXINT = 3 110 CALL FMMSET(MAXINT,ML,MLD2,MLM1) IF (MLD2.EQ.MAXINT .AND. MLM1.NE.ML) THEN MAXINT = ML GO TO 110 ENDIF C C INTMAX is a large value close to the overflow threshold C for integer variables. It is usually 2**31 - 1 C for machines with 32-bit integer arithmetic. C C WARNING: This loop causes integer overflow to occur, so it C is a likely place for the program to fail when C run on a different machine. The loop below has C been used successfully with Fortran compilers C for many different machines, but even different C versions of the same compiler may give different C results. Check the values of MAXINT and INTMAX C if there are problems installing FM. C INTMAX = 3 120 L = 2*INTMAX + 1 IF (INT(L/2).EQ.INTMAX) THEN INTMAX = L GO TO 120 ENDIF C C DPMAX should be set to a value near the machine's double C precision overflow threshold, so that DPMAX and C 1.0D0/DPMAX are both representable in double C precision. C DPMAX = 1.0D+74 C C SPMAX should be set to a value near the machine's single C precision overflow threshold, so that 1.01*SPMAX C and 1.0/SPMAX are both representable in single C precision. C SPMAX = 1.0E+37 C C NDG2MX is the maximum value for NDIG that can be used C internally. FM routines may raise NDIG above C NDIGMX temporarily, to compute correctly C rounded results. C In the definition of LUNPCK, the '6/5' condition C allows for converting from a large base to the C (smaller) largest power of ten base for output C conversion. C The '+ 20' condition allows for the need to carry C many guard digits when using a small base like 2. C NDG2MX = LUNPCK - 1 C C MXBASE is the maximum value for MBASE. C TEMP = MAXINT MXBASE = INT(MIN(DBLE(INTMAX),SQRT(TEMP))) C C MBASE is the currently used base for arithmetic. C K = INT(LOG10(DBLE(MXBASE)/4)) MBASE = 10**K C C NDIG is the number of digits currently being carried. C NPSAVE = NPREC NDIG = 2 + (NPREC+2)/K IF (NDIG.LT.2 .OR. NDIG.GT.NDIGMX) THEN NDIG = MAX(2,MIN(NDIGMX,NDIG)) WRITE (KW,130) NPREC,NDIG 130 FORMAT(//' Precision out of range when calling FMSET.', * ' NPREC =',I20/' The nearest valid NDIG will be used', * ' instead: NDIG =',I6//) NPSAVE = 0 ENDIF C C KFLAG is the flag for error conditions. C KFLAG = 0 C C NTRACE is the trace switch. Default is no printing. C NTRACE = 0 C C LVLTRC is the trace level. Default is to trace only C routines called directly by the user. C LVLTRC = 1 C C NCALL is the call stack pointer. C NCALL = 0 C C NAMEST is the call stack. C DO 140 J = 0, 50 NAMEST(J) = 'MAIN ' 140 CONTINUE C C Some constants that are often needed are stored with the C maximum precision to which they have been computed in the C currently used base. This speeds up the trig, log, power, C and exponential functions. C C NDIGPI is the number of digits available in the currently C stored value of pi (MPISAV). C NDIGPI = 0 C C MBSPI is the value of MBASE for the currently stored C value of pi. C MBSPI = 0 C C NDIGE is the number of digits available in the currently C stored value of e (MESAV). C NDIGE = 0 C C MBSE is the value of MBASE for the currently stored C value of e. C MBSE = 0 C C NDIGLB is the number of digits available in the currently C stored value of LN(MBASE) (MLBSAV). C NDIGLB = 0 C C MBSLB is the value of MBASE for the currently stored C value of LN(MBASE). C MBSLB = 0 C C NDIGLI is the number of digits available in the currently C stored values of the four logarithms used by FMLNI C MLN1 - MLN4. C NDIGLI = 0 C C MBSLI is the value of MBASE for the currently stored C values of MLN1 - MLN4. C MBSLI = 0 C C MXEXP is the current maximum exponent. C MXEXP2 is the internal maximum exponent. This is used to C define the overflow and underflow thresholds. C C These values are chosen so that FM routines can raise the C overflow/underflow limit temporarily while computing C intermediate results, and so that EXP(INTMAX) is greater C than MXBASE**(MXEXP2+1). C C The overflow threshold is MBASE**(MXEXP+1), and the C underflow threshold is MBASE**(-MXEXP-1). C This means the valid exponents in the first word of an FM C number can range from -MXEXP to MXEXP+1 (inclusive). C MXEXP = INT((DBLE(INTMAX))/(2.0D0*LOG(DBLE(MXBASE))) - 1.0D0) MXEXP2 = INT(2*MXEXP + MXEXP/100) C C KACCSW is a switch used to enable cancellation error C monitoring. Routines where cancellation is C not a problem run faster by skipping the C cancellation monitor calculations. C KACCSW = 0 means no error monitoring, C = 1 means error monitoring is done. C KACCSW = 0 C C MEXPUN is the exponent used as a special symbol for C underflowed results. C MEXPUN = -MXEXP2 - 5*NDIGMX C C MEXPOV is the exponent used as a special symbol for C overflowed results. C MEXPOV = -MEXPUN C C MUNKNO is the exponent used as a special symbol for C unknown FM results (1/0, SQRT(-3.0), ...). C MUNKNO = MEXPOV + 5*NDIGMX C C RUNKNO is returned from FM to real or double conversion C routines when no valid result can be expressed in C real or double precision. On systems that provide C a value for undefined results (e.g., Not A Number) C setting RUNKNO to that value is reasonable. On C other systems set it to a value that is likely to C make any subsequent results obviously wrong that C use it. In either case a KFLAG = -4 condition is C also returned. C RUNKNO = -1.01*SPMAX C C IUNKNO is returned from FM to integer conversion routines C when no valid result can be expressed as a one word C integer. KFLAG = -4 is also set. C IUNKNO = -INT(MXEXP2) C C JFORM1 indicates the format used by FMOUT. C JFORM1 = 1 C C JFORM2 indicates the number of digits used in FMOUT. C JFORM2 = NPSAVE C C KRAD = 1 indicates that trig functions use radians, C = 0 means use degrees. C KRAD = 1 C C KWARN = 0 indicates that no warning message is printed C and execution continues when UNKNOWN or another C exception is produced. C = 1 means print a warning message and continue. C = 2 means print a warning message and stop. C KWARN = 1 C C KROUND = 1 causes all results to be rounded to the C nearest FM number, or to the value with C an even last digit if the result is halfway C between two FM numbers. C = 0 causes all results to be chopped. C KROUND = 1 C C KSWIDE defines the maximum screen width to be used for C all unit KW output. C KSWIDE = 80 C C KESWCH = 1 causes input to FMINP with no digits before C the exponent letter to be treated as if there C were a leading '1'. This is sometimes better C for interactive input: 'E7' converts to C 10.0**7. C = 0 causes a leading zero to be assumed. This C gives compatibility with Fortran: 'E7' C converts to 0.0. C KESWCH = 1 C C CMCHAR defines the exponent letter to be used for C FM variable output from FMOUT, as in 1.2345M+678. C Change it to 'E' for output to be read by a C non-FM program. C CMCHAR = 'M' C C KSUB is an internal flag set during subtraction so that C the addition routine will negate its second argument. C KSUB = 0 C C KDEBUG = 0 Error checking is not done for valid input C arguments and parameters like NDIG and MBASE C upon entry to each routine. C = 1 Error checking is done. C KDEBUG = 0 C C Initialize two hash tables that are used for character C look-up during input conversion. C DO 150 J = LHASH1, LHASH2 KHASHT(J) = 5 KHASHV(J) = 0 150 CONTINUE DO 170 J = 1, 21 KPT = ICHAR(LCHARS(J)) IF (KPT.LT.LHASH1 .OR. KPT.GT.LHASH2) THEN WRITE (KW,160) LCHARS(J),KPT,LHASH1,LHASH2 160 FORMAT(/' Error in input conversion.'/ * ' ICHAR function was out of range for the current', * ' dimensions.'/' ICHAR(''',A,''') gave the value ', * I12,', which is outside the currently'/' dimensioned', * ' bounds of (',I5,':',I5,') for variables KHASHT ', * 'and KHASHV.'/' Re-define the two parameters ', * 'LHASH1 and LHASH2 so the dimensions will'/' contain', * ' all possible output values from ICHAR.'//) ELSE KHASHT(KPT) = LTYPES(J) KHASHV(KPT) = LVALS(J) ENDIF 170 CONTINUE C C DPEPS is the approximate machine precision. C ONE = 1.0D0 TWO = 128.0D0 DPEPS = ONE C 180 DPEPS = DPEPS/TWO CALL FMDBL(ONE,DPEPS,YT) IF (YT.GT.ONE) GO TO 180 DPEPS = DPEPS*TWO TWO = 2.0D0 190 DPEPS = DPEPS/TWO CALL FMDBL(ONE,DPEPS,YT) IF (YT.GT.ONE) GO TO 190 DPEPS = DPEPS*TWO C C FMCONS sets several real and double precision constants. C CALL FMCONS C RETURN END SUBROUTINE FMABS(MA,MB) C C MB = ABS(MA) C C IMPLICIT NONE C INTEGER NDIGMX,NBITS,LPACK,LUNPCK,LMWA,LJSUMS,LMBUFF C PARAMETER ( NDIGMX=256 , NBITS=64 , * LPACK = (NDIGMX+1)/2+1 , LUNPCK = (6*NDIGMX)/5+20, * LMWA = 2*LUNPCK , LJSUMS = 8*(LUNPCK+2), * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000+6 ) C DOUBLE PRECISION MA(0:LUNPCK),MB(0:LUNPCK) C INTEGER NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC,KFLAG, * KWARN,KROUND,KSWIDE,KESWCH,KDEBUG DOUBLE PRECISION MBASE C COMMON /FMUSER/ MBASE,NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC, * KFLAG,KWARN,KROUND,KSWIDE,KESWCH,KDEBUG C DOUBLE PRECISION DPMAX INTEGER NCALL,KACCSW,IUNKNO,NDG2MX,INTMAX,KSUB REAL RUNKNO,SPMAX DOUBLE PRECISION MWA,MXEXP,MXEXP2,MEXPUN,MEXPOV,MUNKNO,MXBASE, * MAXINT C COMMON /FM/ MWA(LMWA),NCALL,KACCSW,MXEXP,MXEXP2,MEXPUN,MEXPOV, * MUNKNO,IUNKNO,RUNKNO,MXBASE,NDG2MX,SPMAX,DPMAX, * MAXINT,INTMAX,KSUB C INTEGER LHASH1,LHASH2 PARAMETER (LHASH1=0 , LHASH2=256) DOUBLE PRECISION DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB INTEGER NDIGPI,NDIGE,NDIGLB,NDIGLI,KHASHT,KHASHV, * NGRD21,NGRD52,NGRD22 DOUBLE PRECISION MBSPI,MBSE,MBSLB,MBSLI,MPISAV,MESAV,MLBSAV, * MLN1,MLN2,MLN3,MLN4,MBLOGS,MEXPAB REAL ALOGMB,ALOGMT,ALOGM2,ALOGMX C COMMON /FMSAVE/ NDIGPI,NDIGE,NDIGLB,NDIGLI,MBSPI,MBSE,MBSLB,MBSLI, * MPISAV(0:LUNPCK),MESAV(0:LUNPCK),MLBSAV(0:LUNPCK), * MLN1(0:LUNPCK),MLN2(0:LUNPCK),MLN3(0:LUNPCK), * MLN4(0:LUNPCK),MBLOGS,MEXPAB,ALOGMB,ALOGM2,ALOGMX, * ALOGMT,DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB,KHASHT(LHASH1:LHASH2), * KHASHV(LHASH1:LHASH2),NGRD21,NGRD52,NGRD22 C CHARACTER CMBUFF,CMCHAR CHARACTER *6 NAMEST C COMMON /FMBUFF/ CMBUFF(LMBUFF),NAMEST(0:50),CMCHAR C DOUBLE PRECISION MD2B INTEGER KWRNSV C NCALL = NCALL + 1 NAMEST(NCALL) = 'FMABS ' IF (NTRACE.NE.0) CALL FMNTR(2,MA,MA,1) C KFLAG = 0 KWRNSV = KWARN KWARN = 0 CALL FMEQ(MA,MB) MB(2) = ABS(MB(2)) KWARN = KWRNSV C IF (KACCSW.EQ.1) THEN MD2B = NINT((NDIG-1)*ALOGM2 + LOG(REAL(ABS(MB(2))+1))/0.69315) MB(0) = MIN(MB(0),MD2B) ENDIF IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1) NCALL = NCALL - 1 RETURN END SUBROUTINE FMACOS(MA,MB) C C MB = ACOS(MA) C C IMPLICIT NONE C INTEGER NDIGMX,NBITS,LPACK,LUNPCK,LMWA,LJSUMS,LMBUFF C PARAMETER ( NDIGMX=256 , NBITS=64 , * LPACK = (NDIGMX+1)/2+1 , LUNPCK = (6*NDIGMX)/5+20, * LMWA = 2*LUNPCK , LJSUMS = 8*(LUNPCK+2), * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000+6 ) C DOUBLE PRECISION MA(0:LUNPCK),MB(0:LUNPCK) C INTEGER NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC,KFLAG, * KWARN,KROUND,KSWIDE,KESWCH,KDEBUG DOUBLE PRECISION MBASE C COMMON /FMUSER/ MBASE,NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC, * KFLAG,KWARN,KROUND,KSWIDE,KESWCH,KDEBUG C DOUBLE PRECISION DPMAX INTEGER NCALL,KACCSW,IUNKNO,NDG2MX,INTMAX,KSUB REAL RUNKNO,SPMAX DOUBLE PRECISION MWA,MXEXP,MXEXP2,MEXPUN,MEXPOV,MUNKNO,MXBASE, * MAXINT C COMMON /FM/ MWA(LMWA),NCALL,KACCSW,MXEXP,MXEXP2,MEXPUN,MEXPOV, * MUNKNO,IUNKNO,RUNKNO,MXBASE,NDG2MX,SPMAX,DPMAX, * MAXINT,INTMAX,KSUB C C Scratch array usage during FMACOS: M01 - M06 C DOUBLE PRECISION M01,M02,M03,M04,M05,M06 C COMMON /FM1/ M01(0:LUNPCK),M02(0:LUNPCK),M03(0:LUNPCK), * M04(0:LUNPCK),M05(0:LUNPCK),M06(0:LUNPCK) C INTEGER LHASH1,LHASH2 PARAMETER (LHASH1=0 , LHASH2=256) DOUBLE PRECISION DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB INTEGER NDIGPI,NDIGE,NDIGLB,NDIGLI,KHASHT,KHASHV, * NGRD21,NGRD52,NGRD22 DOUBLE PRECISION MBSPI,MBSE,MBSLB,MBSLI,MPISAV,MESAV,MLBSAV, * MLN1,MLN2,MLN3,MLN4,MBLOGS,MEXPAB REAL ALOGMB,ALOGMT,ALOGM2,ALOGMX C COMMON /FMSAVE/ NDIGPI,NDIGE,NDIGLB,NDIGLI,MBSPI,MBSE,MBSLB,MBSLI, * MPISAV(0:LUNPCK),MESAV(0:LUNPCK),MLBSAV(0:LUNPCK), * MLN1(0:LUNPCK),MLN2(0:LUNPCK),MLN3(0:LUNPCK), * MLN4(0:LUNPCK),MBLOGS,MEXPAB,ALOGMB,ALOGM2,ALOGMX, * ALOGMT,DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB,KHASHT(LHASH1:LHASH2), * KHASHV(LHASH1:LHASH2),NGRD21,NGRD52,NGRD22 C CHARACTER CMBUFF,CMCHAR CHARACTER *6 NAMEST C COMMON /FMBUFF/ CMBUFF(LMBUFF),NAMEST(0:50),CMCHAR C DOUBLE PRECISION MA2,MACCA,MACMAX,MXSAVE INTEGER K,KASAVE,KOVUN,KRESLT,NDSAVE C IF (MBLOGS.NE.MBASE) CALL FMCONS IF (ABS(MA(1)).GT.MEXPAB .OR. MA(1).GT.0 .OR. MA(2).EQ.0) THEN CALL FMENTR('FMACOS',MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE, * KOVUN) IF (KRESLT.NE.0) RETURN ELSE NCALL = NCALL + 1 NAMEST(NCALL) = 'FMACOS' IF (NTRACE.NE.0) CALL FMNTR(2,MA,MA,1) KOVUN = 0 IF (MA(1).EQ.MEXPOV .OR. MA(1).EQ.MEXPUN) KOVUN = 1 NDSAVE = NDIG IF (NCALL.EQ.1) THEN K = MAX(NGRD52-1,2) NDIG = MAX(NDIG+K,2) IF (NDIG.GT.NDG2MX) THEN KFLAG = -9 CALL FMWARN NDIG = NDSAVE KRESLT = 12 CALL FMRSLT(MA,MA,MB,KRESLT) IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1) NCALL = NCALL - 1 RETURN ENDIF ENDIF KASAVE = KACCSW KACCSW = 0 MXSAVE = MXEXP MXEXP = MXEXP2 ENDIF C MA2 = MA(2) MACCA = MA(0) CALL FMEQ2(MA,MB,NDSAVE,NDIG,0) MB(0) = NINT(NDIG*ALOGM2) C C Use ACOS(X) = ATAN(SQRT(1-X*X)/X) C MB(2) = ABS(MB(2)) CALL FMI2M(1,M05) CALL FMSUB(M05,MB,M03) CALL FMADD(M05,MB,M04) CALL FMMPY(M03,M04,M04) CALL FMSQRT(M04,M04) CALL FMDIV(M04,MB,MB) C CALL FMATAN(MB,MB) C IF (MA2.LT.0) THEN IF (KRAD.EQ.1) THEN CALL FMPI(M05) ELSE CALL FMI2M(180,M05) ENDIF CALL FMSUB(M05,MB,MB) ENDIF C C Round the result and return. C MACMAX = NINT((NDSAVE-1)*ALOGM2 + LOG(REAL(ABS(MB(2))+1))/0.69315) MB(0) = MIN(MB(0),MACCA,MACMAX) CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN) RETURN END SUBROUTINE FMADD(MA,MB,MC) C C MC = MA + MB C C This routine performs the trace printing for addition. C FMADD2 is used to do the arithmetic. C C IMPLICIT NONE C INTEGER NDIGMX,NBITS,LPACK,LUNPCK,LMWA,LJSUMS,LMBUFF C PARAMETER ( NDIGMX=256 , NBITS=64 , * LPACK = (NDIGMX+1)/2+1 , LUNPCK = (6*NDIGMX)/5+20, * LMWA = 2*LUNPCK , LJSUMS = 8*(LUNPCK+2), * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000+6 ) C DOUBLE PRECISION MA(0:LUNPCK),MB(0:LUNPCK),MC(0:LUNPCK) C INTEGER NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC,KFLAG, * KWARN,KROUND,KSWIDE,KESWCH,KDEBUG DOUBLE PRECISION MBASE C COMMON /FMUSER/ MBASE,NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC, * KFLAG,KWARN,KROUND,KSWIDE,KESWCH,KDEBUG C DOUBLE PRECISION DPMAX INTEGER NCALL,KACCSW,IUNKNO,NDG2MX,INTMAX,KSUB REAL RUNKNO,SPMAX DOUBLE PRECISION MWA,MXEXP,MXEXP2,MEXPUN,MEXPOV,MUNKNO,MXBASE, * MAXINT C COMMON /FM/ MWA(LMWA),NCALL,KACCSW,MXEXP,MXEXP2,MEXPUN,MEXPOV, * MUNKNO,IUNKNO,RUNKNO,MXBASE,NDG2MX,SPMAX,DPMAX, * MAXINT,INTMAX,KSUB C CHARACTER CMBUFF,CMCHAR CHARACTER *6 NAMEST C COMMON /FMBUFF/ CMBUFF(LMBUFF),NAMEST(0:50),CMCHAR C NCALL = NCALL + 1 IF (NTRACE.NE.0) THEN NAMEST(NCALL) = 'FMADD ' CALL FMNTR(2,MA,MB,2) C CALL FMADD2(MA,MB,MC) C CALL FMNTR(1,MC,MC,1) ELSE CALL FMADD2(MA,MB,MC) ENDIF NCALL = NCALL - 1 RETURN END SUBROUTINE FMADD2(MA,MB,MC) C C Internal addition routine. MC = MA + MB C C IMPLICIT NONE C INTEGER NDIGMX,NBITS,LPACK,LUNPCK,LMWA,LJSUMS,LMBUFF C PARAMETER ( NDIGMX=256 , NBITS=64 , * LPACK = (NDIGMX+1)/2+1 , LUNPCK = (6*NDIGMX)/5+20, * LMWA = 2*LUNPCK , LJSUMS = 8*(LUNPCK+2), * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000+6 ) C DOUBLE PRECISION MA(0:LUNPCK),MB(0:LUNPCK),MC(0:LUNPCK) C INTEGER NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC,KFLAG, * KWARN,KROUND,KSWIDE,KESWCH,KDEBUG DOUBLE PRECISION MBASE C COMMON /FMUSER/ MBASE,NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC, * KFLAG,KWARN,KROUND,KSWIDE,KESWCH,KDEBUG C DOUBLE PRECISION DPMAX INTEGER NCALL,KACCSW,IUNKNO,NDG2MX,INTMAX,KSUB REAL RUNKNO,SPMAX DOUBLE PRECISION MWA,MXEXP,MXEXP2,MEXPUN,MEXPOV,MUNKNO,MXBASE, * MAXINT C COMMON /FM/ MWA(LMWA),NCALL,KACCSW,MXEXP,MXEXP2,MEXPUN,MEXPOV, * MUNKNO,IUNKNO,RUNKNO,MXBASE,NDG2MX,SPMAX,DPMAX, * MAXINT,INTMAX,KSUB C INTEGER LHASH1,LHASH2 PARAMETER (LHASH1=0 , LHASH2=256) DOUBLE PRECISION DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB INTEGER NDIGPI,NDIGE,NDIGLB,NDIGLI,KHASHT,KHASHV, * NGRD21,NGRD52,NGRD22 DOUBLE PRECISION MBSPI,MBSE,MBSLB,MBSLI,MPISAV,MESAV,MLBSAV, * MLN1,MLN2,MLN3,MLN4,MBLOGS,MEXPAB REAL ALOGMB,ALOGMT,ALOGM2,ALOGMX C COMMON /FMSAVE/ NDIGPI,NDIGE,NDIGLB,NDIGLI,MBSPI,MBSE,MBSLB,MBSLI, * MPISAV(0:LUNPCK),MESAV(0:LUNPCK),MLBSAV(0:LUNPCK), * MLN1(0:LUNPCK),MLN2(0:LUNPCK),MLN3(0:LUNPCK), * MLN4(0:LUNPCK),MBLOGS,MEXPAB,ALOGMB,ALOGM2,ALOGMX, * ALOGMT,DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB,KHASHT(LHASH1:LHASH2), * KHASHV(LHASH1:LHASH2),NGRD21,NGRD52,NGRD22 C CHARACTER CMBUFF,CMCHAR CHARACTER *6 NAMEST C COMMON /FMBUFF/ CMBUFF(LMBUFF),NAMEST(0:50),CMCHAR C DOUBLE PRECISION MA0,MA1,MA2,MB0,MB1,MB2,MB2RD INTEGER J,JCOMP,JSIGN,KRESLT,N1,NGUARD,NMWA REAL B2RDA,B2RDB C IF (MBLOGS.NE.MBASE) CALL FMCONS IF (ABS(MA(1)).GT.MEXPAB .OR. ABS(MB(1)).GT.MEXPAB .OR. * KDEBUG.EQ.1) THEN IF (KSUB.EQ.1) THEN CALL FMARGS('FMSUB ',2,MA,MB,KRESLT) ELSE CALL FMARGS('FMADD ',2,MA,MB,KRESLT) ENDIF IF (KRESLT.NE.0) THEN NCALL = NCALL + 1 IF (KSUB.EQ.1) THEN NAMEST(NCALL) = 'FMSUB ' ELSE NAMEST(NCALL) = 'FMADD ' ENDIF CALL FMRSLT(MA,MB,MC,KRESLT) NCALL = NCALL - 1 RETURN ENDIF ELSE IF (MA(2).EQ.0) THEN MA0 = MIN(MA(0),MB(0)) CALL FMEQ(MB,MC) MC(0) = MA0 KFLAG = 1 IF (KSUB.EQ.1) THEN IF (MC(1).NE.MUNKNO) MC(2) = -MC(2) KFLAG = 0 ENDIF RETURN ENDIF IF (MB(2).EQ.0) THEN MA0 = MIN(MA(0),MB(0)) CALL FMEQ(MA,MC) MC(0) = MA0 KFLAG = 1 RETURN ENDIF ENDIF C MA0 = MA(0) IF (KACCSW.EQ.1) THEN MB0 = MB(0) MA1 = MA(1) MB1 = MB(1) ENDIF KFLAG = 0 N1 = NDIG + 1 C C NGUARD is the number of guard digits used. C IF (NCALL.GT.1) THEN NGUARD = NGRD21 IF (NGUARD.GT.NDIG) NGUARD = NDIG ELSE NGUARD = NGRD52 IF (NGUARD.GT.NDIG) NGUARD = NDIG ENDIF NMWA = N1 + NGUARD C C Save the signs of MA and MB and then work with C positive numbers. C JSIGN is the sign of the result of MA + MB. C JSIGN = 1 MA2 = MA(2) MB2 = MB(2) IF (KSUB.EQ.1) MB2 = -MB2 MA(2) = ABS(MA(2)) MB(2) = ABS(MB(2)) C C See which one is larger in absolute value. C IF (MA(1).GT.MB(1)) THEN JCOMP = 1 GO TO 120 ENDIF IF (MB(1).GT.MA(1)) THEN JCOMP = 3 GO TO 120 ENDIF C DO 110 J = 2, N1 IF (MA(J).GT.MB(J)) THEN JCOMP = 1 GO TO 120 ENDIF IF (MB(J).GT.MA(J)) THEN JCOMP = 3 GO TO 120 ENDIF 110 CONTINUE C JCOMP = 2 C 120 IF (JCOMP.LT.3) THEN IF (MA2.LT.0) JSIGN = -1 IF (MA2*MB2.GT.0) THEN CALL FMADDP(MA,MB,NGUARD,NMWA) ELSE CALL FMADDN(MA,MB,NGUARD,NMWA) ENDIF ELSE IF (MB2.LT.0) JSIGN = -1 IF (MA2*MB2.GT.0) THEN CALL FMADDP(MB,MA,NGUARD,NMWA) ELSE CALL FMADDN(MB,MA,NGUARD,NMWA) ENDIF ENDIF IF (KSUB.EQ.1) MB2 = -MB2 MB(2) = MB2 MA(2) = MA2 C C Transfer to MC and fix the sign of the result. C CALL FMMOVE(MWA,MC) IF (JSIGN.LT.0) MC(2) = -MC(2) C IF (KFLAG.LT.0) THEN IF (KSUB.EQ.1) THEN NAMEST(NCALL) = 'FMSUB ' ELSE NAMEST(NCALL) = 'FMADD ' ENDIF CALL FMWARN ENDIF C IF (KACCSW.EQ.1) THEN B2RDA = LOG(REAL(ABS(MC(2))+1)/REAL(ABS(MA2)+1))/0.69315 + * REAL(MC(1)-MA1)*ALOGM2 + REAL(MA0) B2RDB = LOG(REAL(ABS(MC(2))+1)/REAL(ABS(MB2)+1))/0.69315 + * REAL(MC(1)-MB1)*ALOGM2 + REAL(MB0) MB2RD = NINT(MAX(0.0,MIN(B2RDA,B2RDB, * (NDIG-1)*ALOGM2 + LOG(REAL(ABS(MC(2))+1))/0.69315))) IF (MC(2).EQ.0) THEN MC(0) = 0 ELSE MC(0) = MIN(MAX(MA0,MB0),MB2RD) ENDIF ELSE MC(0) = MA0 ENDIF C RETURN END SUBROUTINE FMADDI(MA,IVAL) C C MA = MA + IVAL C C Increment MA by one word integer IVAL. C C This routine is faster than FMADD when IVAL is small enough so C that it can be added to a single word of MA without often causing C a carry. Otherwise FMI2M and FMADD are used. C C IMPLICIT NONE C INTEGER NDIGMX,NBITS,LPACK,LUNPCK,LMWA,LJSUMS,LMBUFF C PARAMETER ( NDIGMX=256 , NBITS=64 , * LPACK = (NDIGMX+1)/2+1 , LUNPCK = (6*NDIGMX)/5+20, * LMWA = 2*LUNPCK , LJSUMS = 8*(LUNPCK+2), * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000+6 ) C DOUBLE PRECISION MA(0:LUNPCK) INTEGER IVAL C INTEGER NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC,KFLAG, * KWARN,KROUND,KSWIDE,KESWCH,KDEBUG DOUBLE PRECISION MBASE C COMMON /FMUSER/ MBASE,NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC, * KFLAG,KWARN,KROUND,KSWIDE,KESWCH,KDEBUG C DOUBLE PRECISION DPMAX INTEGER NCALL,KACCSW,IUNKNO,NDG2MX,INTMAX,KSUB REAL RUNKNO,SPMAX DOUBLE PRECISION MWA,MXEXP,MXEXP2,MEXPUN,MEXPOV,MUNKNO,MXBASE, * MAXINT C COMMON /FM/ MWA(LMWA),NCALL,KACCSW,MXEXP,MXEXP2,MEXPUN,MEXPOV, * MUNKNO,IUNKNO,RUNKNO,MXBASE,NDG2MX,SPMAX,DPMAX, * MAXINT,INTMAX,KSUB C C Scratch array usage during FMADDI: M01 C DOUBLE PRECISION M01,M02,M03,M04,M05,M06 C COMMON /FM1/ M01(0:LUNPCK),M02(0:LUNPCK),M03(0:LUNPCK), * M04(0:LUNPCK),M05(0:LUNPCK),M06(0:LUNPCK) C INTEGER LHASH1,LHASH2 PARAMETER (LHASH1=0 , LHASH2=256) DOUBLE PRECISION DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB INTEGER NDIGPI,NDIGE,NDIGLB,NDIGLI,KHASHT,KHASHV, * NGRD21,NGRD52,NGRD22 DOUBLE PRECISION MBSPI,MBSE,MBSLB,MBSLI,MPISAV,MESAV,MLBSAV, * MLN1,MLN2,MLN3,MLN4,MBLOGS,MEXPAB REAL ALOGMB,ALOGMT,ALOGM2,ALOGMX C COMMON /FMSAVE/ NDIGPI,NDIGE,NDIGLB,NDIGLI,MBSPI,MBSE,MBSLB,MBSLI, * MPISAV(0:LUNPCK),MESAV(0:LUNPCK),MLBSAV(0:LUNPCK), * MLN1(0:LUNPCK),MLN2(0:LUNPCK),MLN3(0:LUNPCK), * MLN4(0:LUNPCK),MBLOGS,MEXPAB,ALOGMB,ALOGM2,ALOGMX, * ALOGMT,DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB,KHASHT(LHASH1:LHASH2), * KHASHV(LHASH1:LHASH2),NGRD21,NGRD52,NGRD22 C CHARACTER CMBUFF,CMCHAR CHARACTER *6 NAMEST C COMMON /FMBUFF/ CMBUFF(LMBUFF),NAMEST(0:50),CMCHAR C DOUBLE PRECISION MAEXP,MD2B,MKSUM INTEGER KPTMA C NCALL = NCALL + 1 IF (NTRACE.NE.0) THEN NAMEST(NCALL) = 'FMADDI' CALL FMNTR(2,MA,MA,1) CALL FMNTRI(2,IVAL,0) ENDIF KFLAG = 0 C MAEXP = MA(1) IF (MAEXP.LE.0 .OR. MAEXP.GT.NDIG) GO TO 110 KPTMA = INT(MAEXP) + 1 IF (KPTMA.GT.2 .AND. MA(2).LT.0) THEN MKSUM = MA(KPTMA) - IVAL ELSE MKSUM = MA(KPTMA) + IVAL ENDIF C IF (MKSUM.GE.MBASE .OR. MKSUM.LE.(-MBASE)) GO TO 110 IF (MA(2).LT.0) THEN IF (KPTMA.GT.2) THEN IF (MKSUM.GE.0) THEN MA(KPTMA) = MKSUM GO TO 120 ELSE GO TO 110 ENDIF ELSE IF (MKSUM.LT.0) THEN MA(KPTMA) = MKSUM GO TO 120 ELSE GO TO 110 ENDIF ENDIF ELSE IF (KPTMA.GT.2) THEN IF (MKSUM.GE.0) THEN MA(KPTMA) = MKSUM GO TO 120 ELSE GO TO 110 ENDIF ELSE IF (MKSUM.GT.0) THEN MA(KPTMA) = MKSUM GO TO 120 ELSE GO TO 110 ENDIF ENDIF ENDIF C 110 CALL FMI2M(IVAL,M01) CALL FMADD(MA,M01,MA) C 120 IF (KACCSW.EQ.1) THEN MD2B = NINT((NDIG-1)*ALOGM2 + LOG(REAL(ABS(MA(2))+1))/0.69315) MA(0) = MIN(MA(0),MD2B) ENDIF IF (NTRACE.NE.0) THEN CALL FMNTR(1,MA,MA,1) ENDIF NCALL = NCALL - 1 RETURN END SUBROUTINE FMADDN(MA,MB,NGUARD,NMWA) C C Internal addition routine. MWA = MA - MB C The arguments are such that MA.GE.MB.GE.0. C C NGUARD is the number of guard digits being carried. C NMWA is the number of words in MWA that will be used. C C IMPLICIT NONE C INTEGER NDIGMX,NBITS,LPACK,LUNPCK,LMWA,LJSUMS,LMBUFF C PARAMETER ( NDIGMX=256 , NBITS=64 , * LPACK = (NDIGMX+1)/2+1 , LUNPCK = (6*NDIGMX)/5+20, * LMWA = 2*LUNPCK , LJSUMS = 8*(LUNPCK+2), * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000+6 ) C DOUBLE PRECISION MA(0:LUNPCK),MB(0:LUNPCK) INTEGER NGUARD,NMWA C INTEGER NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC,KFLAG, * KWARN,KROUND,KSWIDE,KESWCH,KDEBUG DOUBLE PRECISION MBASE C COMMON /FMUSER/ MBASE,NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC, * KFLAG,KWARN,KROUND,KSWIDE,KESWCH,KDEBUG C DOUBLE PRECISION DPMAX INTEGER NCALL,KACCSW,IUNKNO,NDG2MX,INTMAX,KSUB REAL RUNKNO,SPMAX DOUBLE PRECISION MWA,MXEXP,MXEXP2,MEXPUN,MEXPOV,MUNKNO,MXBASE, * MAXINT C COMMON /FM/ MWA(LMWA),NCALL,KACCSW,MXEXP,MXEXP2,MEXPUN,MEXPOV, * MUNKNO,IUNKNO,RUNKNO,MXBASE,NDG2MX,SPMAX,DPMAX, * MAXINT,INTMAX,KSUB C DOUBLE PRECISION MK,MR INTEGER J,K,KL,KP1,KP2,KPT,KSH,N1,N2,NK,NK1 C N1 = NDIG + 1 C C Check for an insignificant operand. C MK = MA(1) - MB(1) IF (MK.GE.NDIG+2) THEN DO 110 J = 1, N1 MWA(J) = MA(J) 110 CONTINUE MWA(N1+1) = 0 KFLAG = 1 RETURN ENDIF K = INT(MK) IF (NGUARD.LE.1) NMWA = N1 + 2 C C Subtract MB from MA. C KP1 = MIN(N1,K+1) MWA(K+1) = 0 DO 120 J = 1, KP1 MWA(J) = MA(J) 120 CONTINUE KP2 = K + 2 C C (Inner Loop) C DO 130 J = KP2, N1 MWA(J) = MA(J) - MB(J-K) 130 CONTINUE C N2 = NDIG + 2 IF (N2-K.LE.1) N2 = 2 + K NK = MIN(NMWA,N1+K) DO 140 J = N2, NK MWA(J) = -MB(J-K) 140 CONTINUE NK1 = NK + 1 DO 150 J = NK1, NMWA MWA(J) = 0 150 CONTINUE C C Normalize. Fix the sign of any negative digit. C IF (K.GT.0) THEN DO 160 J = NMWA, KP2, -1 IF (MWA(J).LT.0) THEN MWA(J) = MWA(J) + MBASE MWA(J-1) = MWA(J-1) - 1 ENDIF 160 CONTINUE C KPT = KP2 - 1 170 IF (MWA(KPT).LT.0 .AND. KPT.GE.3) THEN MWA(KPT) = MWA(KPT) + MBASE MWA(KPT-1) = MWA(KPT-1) - 1 KPT = KPT - 1 GO TO 170 ENDIF GO TO 190 ENDIF C DO 180 J = N1, 3, -1 IF (MWA(J).LT.0) THEN MWA(J) = MWA(J) + MBASE MWA(J-1) = MWA(J-1) - 1 ENDIF 180 CONTINUE C C Shift left if there are any leading zeros in the mantissa. C 190 DO 200 J = 2, NMWA IF (MWA(J).GT.0) THEN KSH = J - 2 GO TO 210 ENDIF 200 CONTINUE MWA(1) = 0 RETURN C 210 IF (KSH.GT.0) THEN KL = NMWA - KSH DO 220 J = 2, KL MWA(J) = MWA(J+KSH) 220 CONTINUE DO 230 J = KL+1, NMWA MWA(J) = 0 230 CONTINUE MWA(1) = MWA(1) - KSH ENDIF C C Round the result. C MR = 2*MWA(NDIG+2) + 1 IF (MR.GE.MBASE) THEN IF (MR-1.GT.MBASE .AND. MWA(N1).LT.MBASE-1) THEN IF (KROUND.NE.0 .OR. NCALL.GT.1) THEN MWA(N1) = MWA(N1) + 1 MWA(N1+1) = 0 ENDIF ELSE CALL FMRND(MWA,NDIG,NGUARD,0) ENDIF ENDIF C C See if the result is equal to one of the input arguments. C IF (ABS(MA(1)-MB(1)).LT.NDIG) GO TO 250 IF (ABS(MA(1)-MB(1)).GT.NDIG+1) THEN KFLAG = 1 GO TO 250 ENDIF C N2 = NDIG + 4 DO 240 J = 3, N1 IF (MWA(N2-J).NE.MA(N2-J)) GO TO 250 240 CONTINUE IF (MWA(1).NE.MA(1)) GO TO 250 IF (MWA(2).NE.ABS(MA(2))) GO TO 250 KFLAG = 1 C 250 RETURN END SUBROUTINE FMADDP(MA,MB,NGUARD,NMWA) C C Internal addition routine. MWA = MA + MB C The arguments are such that MA.GE.MB.GE.0. C C NMWA is the number of words in MWA that will be used. C C IMPLICIT NONE C INTEGER NDIGMX,NBITS,LPACK,LUNPCK,LMWA,LJSUMS,LMBUFF C PARAMETER ( NDIGMX=256 , NBITS=64 , * LPACK = (NDIGMX+1)/2+1 , LUNPCK = (6*NDIGMX)/5+20, * LMWA = 2*LUNPCK , LJSUMS = 8*(LUNPCK+2), * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000+6 ) C DOUBLE PRECISION MA(0:LUNPCK),MB(0:LUNPCK) INTEGER NGUARD,NMWA C INTEGER NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC,KFLAG, * KWARN,KROUND,KSWIDE,KESWCH,KDEBUG DOUBLE PRECISION MBASE C COMMON /FMUSER/ MBASE,NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC, * KFLAG,KWARN,KROUND,KSWIDE,KESWCH,KDEBUG C DOUBLE PRECISION DPMAX INTEGER NCALL,KACCSW,IUNKNO,NDG2MX,INTMAX,KSUB REAL RUNKNO,SPMAX DOUBLE PRECISION MWA,MXEXP,MXEXP2,MEXPUN,MEXPOV,MUNKNO,MXBASE, * MAXINT C COMMON /FM/ MWA(LMWA),NCALL,KACCSW,MXEXP,MXEXP2,MEXPUN,MEXPOV, * MUNKNO,IUNKNO,RUNKNO,MXBASE,NDG2MX,SPMAX,DPMAX, * MAXINT,INTMAX,KSUB C DOUBLE PRECISION MK,MKT,MR INTEGER J,K,KP,KP2,KPT,KSHIFT,N1,N2,NK C N1 = NDIG + 1 C C Check for an insignificant operand. C MK = MA(1) - MB(1) IF (MK.GE.NDIG+1) THEN MWA(1) = MA(1) + 1 MWA(2) = 0 DO 110 J = 2, N1 MWA(J+1) = MA(J) 110 CONTINUE MWA(N1+2) = 0 KFLAG = 1 RETURN ENDIF K = INT(MK) C C Add MA and MB. C MWA(1) = MA(1) + 1 MWA(2) = 0 DO 120 J = 2, K+1 MWA(J+1) = MA(J) 120 CONTINUE KP2 = K + 2 C C (Inner Loop) C DO 130 J = KP2, N1 MWA(J+1) = MA(J) + MB(J-K) 130 CONTINUE N2 = NDIG + 2 NK = MIN(NMWA,N1+K) DO 140 J = N2, NK MWA(J+1) = MB(J-K) 140 CONTINUE DO 150 J = NK+1, NMWA MWA(J+1) = 0 150 CONTINUE C C Normalize. Fix any digit not less than MBASE. C IF (K.EQ.NDIG) GO TO 220 C IF (K.GT.0) THEN DO 160 J = N1+1, KP2, -1 IF (MWA(J).GE.MBASE) THEN MWA(J) = MWA(J) - MBASE MWA(J-1) = MWA(J-1) + 1 ENDIF 160 CONTINUE C KPT = KP2 - 1 170 IF (MWA(KPT).GE.MBASE .AND. KPT.GE.3) THEN MWA(KPT) = MWA(KPT) - MBASE MWA(KPT-1) = MWA(KPT-1) + 1 KPT = KPT - 1 GO TO 170 ENDIF GO TO 190 ENDIF C DO 180 J = N1+1, 3, -1 IF (MWA(J).GE.MBASE) THEN MWA(J) = MWA(J) - MBASE MWA(J-1) = MWA(J-1) + 1 ENDIF 180 CONTINUE C C Shift right if the leading digit is not less than MBASE. C 190 IF (MWA(2).GE.MBASE) THEN 200 KP = NMWA + 4 DO 210 J = 4, NMWA MWA(KP-J) = MWA(KP-J-1) 210 CONTINUE MKT = DINT(MWA(2)/MBASE) MWA(3) = MWA(2) - MKT*MBASE MWA(2) = MKT MWA(1) = MWA(1) + 1 IF (MWA(2).GE.MBASE) GO TO 200 ENDIF C C Round the result. C 220 KSHIFT = 0 IF (MWA(2).EQ.0) KSHIFT = 1 MR = 2*MWA(NDIG+2+KSHIFT) + 1 IF (MR.GE.MBASE) THEN IF (MR-1.GT.MBASE .AND. MWA(N1+KSHIFT).LT.MBASE-1) THEN IF (KROUND.NE.0 .OR. NCALL.GT.1) THEN MWA(N1+KSHIFT) = MWA(N1+KSHIFT) + 1 MWA(N1+1+KSHIFT) = 0 ENDIF ELSE CALL FMRND(MWA,NDIG,NGUARD,KSHIFT) ENDIF ENDIF C C See if the result is equal to one of the input arguments. C IF (ABS(MA(1)-MB(1)).LT.NDIG) GO TO 240 IF (KSHIFT.EQ.0) GO TO 240 IF (ABS(MA(1)-MB(1)).GT.NDIG+1) THEN KFLAG = 1 GO TO 240 ENDIF C N2 = NDIG + 4 DO 230 J = 3, N1 IF (MWA(N2-J+1).NE.MA(N2-J)) GO TO 240 230 CONTINUE IF (MWA(1).NE.MA(1)+1) GO TO 240 IF (MWA(3).NE.ABS(MA(2))) GO TO 240 KFLAG = 1 C 240 RETURN END SUBROUTINE FMARGS(KROUTN,NARGS,MA,MB,KRESLT) C C Check the input arguments to a routine for special cases. C C KROUTN - Name of the subroutine that was called C NARGS - The number of input arguments (1 or 2) C MA - First input argument C MB - Second input argument (if NARGS is 2) C KRESLT - Result code returned to the calling routine. C C Result codes: C C 0 - Perform the normal operation C 1 - The result is the first input argument C 2 - The result is the second input argument C 3 - The result is -OVERFLOW C 4 - The result is +OVERFLOW C 5 - The result is -UNDERFLOW C 6 - The result is +UNDERFLOW C 7 - The result is -1.0 C 8 - The result is +1.0 C 9 - The result is -pi/2 C 10 - The result is +pi/2 C 11 - The result is 0.0 C 12 - The result is UNKNOWN C 13 - The result is +pi C 14 - The result is -pi/4 C 15 - The result is +pi/4 C C IMPLICIT NONE C INTEGER NDIGMX,NBITS,LPACK,LUNPCK,LMWA,LJSUMS,LMBUFF C PARAMETER ( NDIGMX=256 , NBITS=64 , * LPACK = (NDIGMX+1)/2+1 , LUNPCK = (6*NDIGMX)/5+20, * LMWA = 2*LUNPCK , LJSUMS = 8*(LUNPCK+2), * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000+6 ) C CHARACTER *6 KROUTN DOUBLE PRECISION MA(0:LUNPCK),MB(0:LUNPCK) INTEGER NARGS,KRESLT C INTEGER NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC,KFLAG, * KWARN,KROUND,KSWIDE,KESWCH,KDEBUG DOUBLE PRECISION MBASE C COMMON /FMUSER/ MBASE,NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC, * KFLAG,KWARN,KROUND,KSWIDE,KESWCH,KDEBUG C DOUBLE PRECISION DPMAX INTEGER NCALL,KACCSW,IUNKNO,NDG2MX,INTMAX,KSUB REAL RUNKNO,SPMAX DOUBLE PRECISION MWA,MXEXP,MXEXP2,MEXPUN,MEXPOV,MUNKNO,MXBASE, * MAXINT C COMMON /FM/ MWA(LMWA),NCALL,KACCSW,MXEXP,MXEXP2,MEXPUN,MEXPOV, * MUNKNO,IUNKNO,RUNKNO,MXBASE,NDG2MX,SPMAX,DPMAX, * MAXINT,INTMAX,KSUB C INTEGER LHASH1,LHASH2 PARAMETER (LHASH1=0 , LHASH2=256) DOUBLE PRECISION DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB INTEGER NDIGPI,NDIGE,NDIGLB,NDIGLI,KHASHT,KHASHV, * NGRD21,NGRD52,NGRD22 DOUBLE PRECISION MBSPI,MBSE,MBSLB,MBSLI,MPISAV,MESAV,MLBSAV, * MLN1,MLN2,MLN3,MLN4,MBLOGS,MEXPAB REAL ALOGMB,ALOGMT,ALOGM2,ALOGMX C COMMON /FMSAVE/ NDIGPI,NDIGE,NDIGLB,NDIGLI,MBSPI,MBSE,MBSLB,MBSLI, * MPISAV(0:LUNPCK),MESAV(0:LUNPCK),MLBSAV(0:LUNPCK), * MLN1(0:LUNPCK),MLN2(0:LUNPCK),MLN3(0:LUNPCK), * MLN4(0:LUNPCK),MBLOGS,MEXPAB,ALOGMB,ALOGM2,ALOGMX, * ALOGMT,DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB,KHASHT(LHASH1:LHASH2), * KHASHV(LHASH1:LHASH2),NGRD21,NGRD52,NGRD22 C CHARACTER CMBUFF,CMCHAR CHARACTER *6 NAMEST C COMMON /FMBUFF/ CMBUFF(LMBUFF),NAMEST(0:50),CMCHAR C DOUBLE PRECISION MBS INTEGER J,KWRNSV,NCATMA,NCATMB,NDS INTEGER KADD(15,15),KMPY(15,15),KDIV(15,15),KPWR(15,15), * KSQRT(15),KEXP(15),KLN(15),KSIN(15),KCOS(15),KTAN(15), * KASIN(15),KACOS(15),KATAN(15),KSINH(15),KCOSH(15), * KTANH(15),KLG10(15) C C These tables define the result codes to be returned for C given values of the input argument(s). C C For example, in row 7 column 2 of this DATA statement C KADD(2,7) = 2 means that if the first argument in a call C to FMADD is in category 7 ( -UNDERFLOW ) and the second C argument is in category 2 ( near -OVERFLOW but C representable ) then the result code is 2 ( the value C of the sum is equal to the second input argument). C See routine FMCAT for descriptions of the categories. C DATA KADD/ 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,12,12, 2 3, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0,12, 3 3, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 4, 4 3, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 4, 5 3, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 4, 6 3, 0, 0, 0, 0, 0,12, 1,12, 0, 0, 0, 0, 0, 4, 7 3, 2, 2, 2, 2,12,12, 5,12,12, 2, 2, 2, 2, 4, 8 3, 2, 2, 2, 2, 2, 5, 2, 6, 2, 2, 2, 2, 2, 4, 9 3, 2, 2, 2, 2,12,12, 6,12,12, 2, 2, 2, 2, 4, A 3, 0, 0, 0, 0, 0,12, 1,12, 0, 0, 0, 0, 0, 4, B 3, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 4, C 3, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 4, D 3, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 4, E 12, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 4, F 12,12, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 / C DATA KMPY/ 4, 4, 4, 4,12,12,12,11,12,12,12, 3, 3, 3, 3, 2 4, 0, 0, 0, 0, 0,12,11,12, 0, 0, 1, 0, 0, 3, 3 4, 0, 0, 0, 0, 0,12,11,12, 0, 0, 1, 0, 0, 3, 4 4, 0, 0, 0, 0, 0, 6,11, 5, 0, 0, 1, 0, 0, 3, 5 12, 0, 0, 0, 0, 0, 6,11, 5, 0, 0, 1, 0, 0,12, 6 12, 0, 0, 0, 0, 0, 6,11, 5, 0, 0, 1, 0, 0,12, 7 12,12,12, 6, 6, 6, 6,11, 5, 5, 5, 5,12,12,12, 8 11,11,11,11,11,11,11,11,11,11,11,11,11,11,11, 9 12,12,12, 5, 5, 5, 5,11, 6, 6, 6, 6,12,12,12, A 12, 0, 0, 0, 0, 0, 5,11, 6, 0, 0, 1, 0, 0,12, B 12, 0, 0, 0, 0, 0, 5,11, 6, 0, 0, 1, 0, 0,12, C 3, 2, 2, 2, 2, 2, 5,11, 6, 2, 2, 2, 2, 2, 4, D 3, 0, 0, 0, 0, 0,12,11,12, 0, 0, 1, 0, 0, 4, E 3, 0, 0, 0, 0, 0,12,11,12, 0, 0, 1, 0, 0, 4, F 3, 3, 3, 3,12,12,12,11,12,12,12, 4, 4, 4, 4 / C DATA KDIV/12,12,12, 4, 4, 4, 4,12, 3, 3, 3, 3,12,12,12, 2 12, 0, 0, 0, 0, 0, 4,12, 3, 0, 0, 1, 0, 0,12, 3 12, 0, 0, 0, 0, 0, 4,12, 3, 0, 0, 1, 0, 0,12, 4 6, 0, 0, 0, 0, 0, 4,12, 3, 0, 0, 1, 0, 0, 5, 5 6, 0, 0, 0, 0, 0,12,12,12, 0, 0, 1, 0, 0, 5, 6 6, 0, 0, 0, 0, 0,12,12,12, 0, 0, 1, 0, 0, 5, 7 6, 6, 6, 6,12,12,12,12,12,12,12, 5, 5, 5, 5, 8 11,11,11,11,11,11,11,12,11,11,11,11,11,11,11, 9 5, 5, 5, 5,12,12,12,12,12,12,12, 6, 6, 6, 6, A 5, 0, 0, 0, 0, 0,12,12,12, 0, 0, 1, 0, 0, 6, B 5, 0, 0, 0, 0, 0,12,12,12, 0, 0, 1, 0, 0, 6, C 5, 0, 0, 0, 0, 0, 3,12, 4, 0, 0, 1, 0, 0, 6, D 12, 0, 0, 0, 0, 0, 3,12, 4, 0, 0, 1, 0, 0,12, E 12, 0, 0, 0, 0, 0, 3,12, 4, 0, 0, 1, 0, 0,12, F 12,12,12, 3, 3, 3, 3,12, 4, 4, 4, 4,12,12,12 / C DATA KPWR/12,12, 0, 5,12,12,12, 8,12,12,12, 3, 0,12,12, 2 12,12, 0, 0,12,12,12, 8,12,12,12, 1, 0,12,12, 3 12,12, 0, 0,12,12,12, 8,12,12,12, 1, 0,12,12, 4 12,12, 0, 0,12,12,12, 8,12,12,12, 1, 0,12,12, 5 12,12, 0, 0,12,12,12, 8,12,12,12, 1, 0,12,12, 6 12,12, 0, 0,12,12,12, 8,12,12,12, 1, 0,12,12, 7 12,12, 0, 3,12,12,12, 8,12,12,12, 5, 0,12,12, 8 12,12,12,12,12,12,12,12,11,11,11,11,11,11,11, 9 4, 4, 4, 4,12,12,12, 8,12,12,12, 6, 6, 6, 6, A 4, 4, 0, 0, 0, 8, 8, 8, 8, 0, 0, 1, 0, 6, 6, B 4, 4, 0, 0, 0, 8, 8, 8, 8, 0, 0, 1, 0, 6, 6, C 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, D 6, 6, 0, 0, 0, 8, 8, 8, 8, 8, 0, 1, 0, 4, 4, E 6, 6, 0, 0, 0, 8, 8, 8, 8, 8, 0, 1, 0, 4, 4, F 6, 6, 6, 6,12,12,12, 8,12,12,12, 4, 4, 4, 4 / C DATA KSQRT/12,12,12,12,12,12,12,11,12, 0, 0, 8, 0, 0,12/ DATA KEXP / 6, 6, 0, 0, 0, 8, 8, 8, 8, 8, 0, 0, 0, 4, 4/ DATA KLN /12,12,12,12,12,12,12,12,12, 0, 0,11, 0, 0,12/ DATA KSIN /12,12, 0, 0, 0, 0, 5,11, 6, 0, 0, 0, 0,12,12/ DATA KCOS /12,12, 0, 0, 0, 8, 8, 8, 8, 8, 0, 0, 0,12,12/ DATA KTAN /12,12, 0, 0, 0, 0, 5,11, 6, 0, 0, 0, 0,12,12/ DATA KASIN/12,12,12, 9, 0, 0, 5,11, 6, 0, 0,10,12,12,12/ DATA KACOS/12,12,12,13, 0,10,10,10,10,10, 0,11,12,12,12/ DATA KATAN/ 9, 9, 0,14, 0, 0, 5,11, 6, 0, 0,15, 0,10,10/ DATA KSINH/ 3, 3, 0, 0, 0, 1, 5,11, 6, 1, 0, 0, 0, 4, 4/ DATA KCOSH/ 4, 4, 0, 0, 0, 8, 8, 8, 8, 8, 0, 0, 0, 4, 4/ DATA KTANH/ 7, 7, 0, 0, 0, 1, 5,11, 6, 1, 0, 0, 0, 8, 8/ DATA KLG10/12,12,12,12,12,12,12,12,12, 0, 0,11, 0, 0,12/ C KRESLT = 12 KFLAG = -4 IF (MA(1).EQ.MUNKNO) RETURN IF (NARGS.EQ.2) THEN IF (MB(1).EQ.MUNKNO) RETURN ENDIF IF (MBLOGS.NE.MBASE) CALL FMCONS KFLAG = 0 NAMEST(NCALL) = KROUTN C C Check the validity of parameters if this is a user call. C IF (NCALL.GT.1 .AND. KDEBUG.EQ.0) GO TO 170 C C Check NDIG. C IF (NDIG.LT.2 .OR. NDIG.GT.NDIGMX) THEN KFLAG = -1 CALL FMWARN NDS = NDIG IF (NDIG.LT.2) NDIG = 2 IF (NDIG.GT.NDIGMX) NDIG = NDIGMX WRITE (KW,110) NDS,NDIG 110 FORMAT(' NDIG was',I10,'. It has been changed to',I10,'.') RETURN ENDIF C C Check MBASE. C IF (MBASE.LT.2 .OR. MBASE.GT.MXBASE) THEN KFLAG = -2 CALL FMWARN MBS = MBASE IF (MBASE.LT.2) MBASE = 2 IF (MBASE.GT.MXBASE) MBASE = MXBASE WRITE (KW,120) INT(MBS),INT(MBASE) 120 FORMAT(' MBASE was',I10,'. It has been changed to',I10,'.') CALL FMCONS RETURN ENDIF C C Check exponent range. C IF (MA(1).GT.MXEXP+1 .OR. MA(1).LT.-MXEXP) THEN IF (ABS(MA(1)).NE.MEXPOV .OR. ABS(MA(2)).NE.1) THEN CALL FMIM(0,MA) KFLAG = -3 CALL FMWARN MA(1) = MUNKNO MA(2) = 1 MA(0) = NINT(NDIG*ALOGM2) RETURN ENDIF ENDIF IF (NARGS.EQ.2) THEN IF (MB(1).GT.MXEXP+1 .OR. MB(1).LT.-MXEXP) THEN IF (ABS(MB(1)).NE.MEXPOV .OR. ABS(MB(2)).NE.1) THEN CALL FMIM(0,MB) KFLAG = -3 CALL FMWARN MB(1) = MUNKNO MB(2) = 1 MB(0) = NINT(NDIG*ALOGM2) RETURN ENDIF ENDIF ENDIF C C Check for properly normalized digits in the C input arguments. C IF (ABS(MA(1)-INT(MA(1))).NE.0) KFLAG = 1 IF (MA(2).LE.(-MBASE) .OR. MA(2).GE.MBASE .OR. * ABS(MA(2)-INT(MA(2))).NE.0) KFLAG = 2 IF (KDEBUG.EQ.0) GO TO 140 DO 130 J = 3, NDIG+1 IF (MA(J).LT.0 .OR. MA(J).GE.MBASE .OR. * ABS(MA(J)-INT(MA(J))).NE.0) THEN KFLAG = J GO TO 140 ENDIF 130 CONTINUE 140 IF (KFLAG.NE.0) THEN J = KFLAG MBS = MA(J) CALL FMIM(0,MA) KFLAG = -4 KWRNSV = KWARN IF (KWARN.GE.2) KWARN = 1 CALL FMWARN KWARN = KWRNSV IF (KWARN.GE.1) THEN WRITE (KW,*) ' First invalid array element: MA(', * J,') = ',MBS ENDIF MA(1) = MUNKNO MA(2) = 1 MA(0) = NINT(NDIG*ALOGM2) IF (KWARN.GE.2) THEN STOP ENDIF RETURN ENDIF IF (NARGS.EQ.2) THEN IF (ABS(MB(1)-INT(MB(1))).NE.0) KFLAG = 1 IF (MB(2).LE.(-MBASE) .OR. MB(2).GE.MBASE .OR. * ABS(MB(2)-INT(MB(2))).NE.0) KFLAG = 2 IF (KDEBUG.EQ.0) GO TO 160 DO 150 J = 3, NDIG+1 IF (MB(J).LT.0 .OR. MB(J).GE.MBASE .OR. * ABS(MB(J)-INT(MB(J))).NE.0) THEN KFLAG = J GO TO 160 ENDIF 150 CONTINUE 160 IF (KFLAG.NE.0) THEN J = KFLAG MBS = MB(J) CALL FMIM(0,MB) KFLAG = -4 KWRNSV = KWARN IF (KWARN.GE.2) KWARN = 1 CALL FMWARN KWARN = KWRNSV IF (KWARN.GE.1) THEN WRITE (KW,*) ' First invalid array element: MB(', * J,') = ',MBS ENDIF MB(1) = MUNKNO MB(2) = 1 MB(0) = NINT(NDIG*ALOGM2) IF (KWARN.GE.2) THEN STOP ENDIF RETURN ENDIF ENDIF C C Check for special cases. C 170 CALL FMCAT(MA,NCATMA) NCATMB = 0 IF (NARGS.EQ.2) CALL FMCAT(MB,NCATMB) C IF (KROUTN.EQ.'FMADD ') THEN KRESLT = KADD(NCATMB,NCATMA) GO TO 180 ENDIF C IF (KROUTN.EQ.'FMSUB ') THEN IF (NCATMB.LT.16) NCATMB = 16 - NCATMB KRESLT = KADD(NCATMB,NCATMA) GO TO 180 ENDIF C IF (KROUTN.EQ.'FMMPY ') THEN KRESLT = KMPY(NCATMB,NCATMA) GO TO 180 ENDIF C IF (KROUTN.EQ.'FMDIV ') THEN KRESLT = KDIV(NCATMB,NCATMA) GO TO 180 ENDIF C IF (KROUTN.EQ.'FMPWR ') THEN KRESLT = KPWR(NCATMB,NCATMA) GO TO 180 ENDIF C IF (KROUTN.EQ.'FMSQRT') THEN KRESLT = KSQRT(NCATMA) GO TO 180 ENDIF C IF (KROUTN.EQ.'FMEXP ') THEN KRESLT = KEXP(NCATMA) GO TO 180 ENDIF C IF (KROUTN.EQ.'FMLN ') THEN KRESLT = KLN(NCATMA) GO TO 180 ENDIF C IF (KROUTN.EQ.'FMSIN ') THEN KRESLT = KSIN(NCATMA) GO TO 180 ENDIF C IF (KROUTN.EQ.'FMCOS ') THEN KRESLT = KCOS(NCATMA) GO TO 180 ENDIF C IF (KROUTN.EQ.'FMTAN ') THEN KRESLT = KTAN(NCATMA) GO TO 180 ENDIF C IF (KROUTN.EQ.'FMASIN') THEN KRESLT = KASIN(NCATMA) IF ((NCATMA.EQ.7.OR.NCATMA.EQ.9) .AND. KRAD.EQ.0) KRESLT = 12 GO TO 180 ENDIF C IF (KROUTN.EQ.'FMACOS') THEN KRESLT = KACOS(NCATMA) GO TO 180 ENDIF C IF (KROUTN.EQ.'FMATAN') THEN KRESLT = KATAN(NCATMA) IF ((NCATMA.EQ.7.OR.NCATMA.EQ.9) .AND. KRAD.EQ.0) KRESLT = 12 GO TO 180 ENDIF C IF (KROUTN.EQ.'FMSINH') THEN KRESLT = KSINH(NCATMA) GO TO 180 ENDIF C IF (KROUTN.EQ.'FMCOSH') THEN KRESLT = KCOSH(NCATMA) GO TO 180 ENDIF C IF (KROUTN.EQ.'FMTANH') THEN KRESLT = KTANH(NCATMA) GO TO 180 ENDIF C IF (KROUTN.EQ.'FMLG10') THEN KRESLT = KLG10(NCATMA) GO TO 180 ENDIF C KRESLT = 0 RETURN C 180 IF (KRESLT.EQ.12) THEN KFLAG = -4 CALL FMWARN ENDIF IF (KRESLT.EQ.3 .OR. KRESLT.EQ.4) THEN IF (NCATMA.EQ.1 .OR. NCATMA.EQ.7 .OR. NCATMA.EQ.9 .OR. * NCATMA.EQ.15 .OR. NCATMB.EQ.1 .OR. NCATMB.EQ.7 .OR. * NCATMB.EQ.9 .OR. NCATMB.EQ.15) THEN KFLAG = -5 ELSE KFLAG = -5 CALL FMWARN ENDIF ENDIF IF (KRESLT.EQ.5 .OR. KRESLT.EQ.6) THEN IF (NCATMA.EQ.1 .OR. NCATMA.EQ.7 .OR. NCATMA.EQ.9 .OR. * NCATMA.EQ.15 .OR. NCATMB.EQ.1 .OR. NCATMB.EQ.7 .OR. * NCATMB.EQ.9 .OR. NCATMB.EQ.15) THEN KFLAG = -6 ELSE KFLAG = -6 CALL FMWARN ENDIF ENDIF RETURN END SUBROUTINE FMASIN(MA,MB) C C MB = ARCSIN(MA) C C IMPLICIT NONE C INTEGER NDIGMX,NBITS,LPACK,LUNPCK,LMWA,LJSUMS,LMBUFF C PARAMETER ( NDIGMX=256 , NBITS=64 , * LPACK = (NDIGMX+1)/2+1 , LUNPCK = (6*NDIGMX)/5+20, * LMWA = 2*LUNPCK , LJSUMS = 8*(LUNPCK+2), * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000+6 ) C DOUBLE PRECISION MA(0:LUNPCK),MB(0:LUNPCK) C INTEGER NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC,KFLAG, * KWARN,KROUND,KSWIDE,KESWCH,KDEBUG DOUBLE PRECISION MBASE C COMMON /FMUSER/ MBASE,NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC, * KFLAG,KWARN,KROUND,KSWIDE,KESWCH,KDEBUG C DOUBLE PRECISION DPMAX INTEGER NCALL,KACCSW,IUNKNO,NDG2MX,INTMAX,KSUB REAL RUNKNO,SPMAX DOUBLE PRECISION MWA,MXEXP,MXEXP2,MEXPUN,MEXPOV,MUNKNO,MXBASE, * MAXINT C COMMON /FM/ MWA(LMWA),NCALL,KACCSW,MXEXP,MXEXP2,MEXPUN,MEXPOV, * MUNKNO,IUNKNO,RUNKNO,MXBASE,NDG2MX,SPMAX,DPMAX, * MAXINT,INTMAX,KSUB C C Scratch array usage during FMASIN: M01 - M06 C DOUBLE PRECISION M01,M02,M03,M04,M05,M06 C COMMON /FM1/ M01(0:LUNPCK),M02(0:LUNPCK),M03(0:LUNPCK), * M04(0:LUNPCK),M05(0:LUNPCK),M06(0:LUNPCK) C INTEGER LHASH1,LHASH2 PARAMETER (LHASH1=0 , LHASH2=256) DOUBLE PRECISION DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB INTEGER NDIGPI,NDIGE,NDIGLB,NDIGLI,KHASHT,KHASHV, * NGRD21,NGRD52,NGRD22 DOUBLE PRECISION MBSPI,MBSE,MBSLB,MBSLI,MPISAV,MESAV,MLBSAV, * MLN1,MLN2,MLN3,MLN4,MBLOGS,MEXPAB REAL ALOGMB,ALOGMT,ALOGM2,ALOGMX C COMMON /FMSAVE/ NDIGPI,NDIGE,NDIGLB,NDIGLI,MBSPI,MBSE,MBSLB,MBSLI, * MPISAV(0:LUNPCK),MESAV(0:LUNPCK),MLBSAV(0:LUNPCK), * MLN1(0:LUNPCK),MLN2(0:LUNPCK),MLN3(0:LUNPCK), * MLN4(0:LUNPCK),MBLOGS,MEXPAB,ALOGMB,ALOGM2,ALOGMX, * ALOGMT,DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB,KHASHT(LHASH1:LHASH2), * KHASHV(LHASH1:LHASH2),NGRD21,NGRD52,NGRD22 C CHARACTER CMBUFF,CMCHAR CHARACTER *6 NAMEST C COMMON /FMBUFF/ CMBUFF(LMBUFF),NAMEST(0:50),CMCHAR C DOUBLE PRECISION MACCA,MACMAX,MXSAVE INTEGER K,KASAVE,KOVUN,KRESLT,NDSAVE C IF (MBLOGS.NE.MBASE) CALL FMCONS IF (ABS(MA(1)).GT.MEXPAB .OR. MA(1).GT.0 .OR. MA(2).EQ.0) THEN CALL FMENTR('FMASIN',MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE, * KOVUN) IF (KRESLT.NE.0) RETURN ELSE NCALL = NCALL + 1 NAMEST(NCALL) = 'FMASIN' IF (NTRACE.NE.0) CALL FMNTR(2,MA,MA,1) KOVUN = 0 IF (MA(1).EQ.MEXPOV .OR. MA(1).EQ.MEXPUN) KOVUN = 1 NDSAVE = NDIG IF (NCALL.EQ.1) THEN K = MAX(NGRD52-1,2) NDIG = MAX(NDIG+K,2) IF (NDIG.GT.NDG2MX) THEN KFLAG = -9 CALL FMWARN NDIG = NDSAVE KRESLT = 12 CALL FMRSLT(MA,MA,MB,KRESLT) IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1) NCALL = NCALL - 1 RETURN ENDIF ENDIF KASAVE = KACCSW KACCSW = 0 MXSAVE = MXEXP MXEXP = MXEXP2 ENDIF C MACCA = MA(0) CALL FMEQ2(MA,MB,NDSAVE,NDIG,0) MB(0) = NINT(NDIG*ALOGM2) C C Use ASIN(X) = ATAN(X/SQRT(1-X*X)) C CALL FMI2M(1,M05) CALL FMSUB(M05,MB,M03) CALL FMADD(M05,MB,M04) CALL FMMPY(M03,M04,M04) CALL FMSQRT(M04,M04) CALL FMDIV(MB,M04,MB) C CALL FMATAN(MB,MB) C C Round the result and return. C MACMAX = NINT((NDSAVE-1)*ALOGM2 + LOG(REAL(ABS(MB(2))+1))/0.69315) MB(0) = MIN(MB(0),MACCA,MACMAX) CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN) RETURN END SUBROUTINE FMATAN(MA,MB) C C MB = ARCTAN(MA) C C IMPLICIT NONE C INTEGER NDIGMX,NBITS,LPACK,LUNPCK,LMWA,LJSUMS,LMBUFF C PARAMETER ( NDIGMX=256 , NBITS=64 , * LPACK = (NDIGMX+1)/2+1 , LUNPCK = (6*NDIGMX)/5+20, * LMWA = 2*LUNPCK , LJSUMS = 8*(LUNPCK+2), * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000+6 ) C DOUBLE PRECISION X,XM DOUBLE PRECISION MA(0:LUNPCK),MB(0:LUNPCK) INTEGER NSTACK(19) C INTEGER NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC,KFLAG, * KWARN,KROUND,KSWIDE,KESWCH,KDEBUG DOUBLE PRECISION MBASE C COMMON /FMUSER/ MBASE,NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC, * KFLAG,KWARN,KROUND,KSWIDE,KESWCH,KDEBUG C DOUBLE PRECISION DPMAX INTEGER NCALL,KACCSW,IUNKNO,NDG2MX,INTMAX,KSUB REAL RUNKNO,SPMAX DOUBLE PRECISION MWA,MXEXP,MXEXP2,MEXPUN,MEXPOV,MUNKNO,MXBASE, * MAXINT C COMMON /FM/ MWA(LMWA),NCALL,KACCSW,MXEXP,MXEXP2,MEXPUN,MEXPOV, * MUNKNO,IUNKNO,RUNKNO,MXBASE,NDG2MX,SPMAX,DPMAX, * MAXINT,INTMAX,KSUB C C Scratch array usage during FMATAN: M01 - M06 C DOUBLE PRECISION M01,M02,M03,M04,M05,M06 C COMMON /FM1/ M01(0:LUNPCK),M02(0:LUNPCK),M03(0:LUNPCK), * M04(0:LUNPCK),M05(0:LUNPCK),M06(0:LUNPCK) C INTEGER LHASH1,LHASH2 PARAMETER (LHASH1=0 , LHASH2=256) DOUBLE PRECISION DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB INTEGER NDIGPI,NDIGE,NDIGLB,NDIGLI,KHASHT,KHASHV, * NGRD21,NGRD52,NGRD22 DOUBLE PRECISION MBSPI,MBSE,MBSLB,MBSLI,MPISAV,MESAV,MLBSAV, * MLN1,MLN2,MLN3,MLN4,MBLOGS,MEXPAB REAL ALOGMB,ALOGMT,ALOGM2,ALOGMX C COMMON /FMSAVE/ NDIGPI,NDIGE,NDIGLB,NDIGLI,MBSPI,MBSE,MBSLB,MBSLI, * MPISAV(0:LUNPCK),MESAV(0:LUNPCK),MLBSAV(0:LUNPCK), * MLN1(0:LUNPCK),MLN2(0:LUNPCK),MLN3(0:LUNPCK), * MLN4(0:LUNPCK),MBLOGS,MEXPAB,ALOGMB,ALOGM2,ALOGMX, * ALOGMT,DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB,KHASHT(LHASH1:LHASH2), * KHASHV(LHASH1:LHASH2),NGRD21,NGRD52,NGRD22 C CHARACTER CMBUFF,CMCHAR CHARACTER *6 NAMEST C COMMON /FMBUFF/ CMBUFF(LMBUFF),NAMEST(0:50),CMCHAR C DOUBLE PRECISION MA1,MA2,MACCA,MACMAX,MXSAVE INTEGER J,K,KASAVE,KOVUN,KRESLT,KRSAVE,KST,KWRNSV,NDSAV1,NDSAVE, * NDSV C IF (MBLOGS.NE.MBASE) CALL FMCONS IF (ABS(MA(1)).GT.MEXPAB .OR. MA(2).EQ.0) THEN CALL FMENTR('FMATAN',MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE, * KOVUN) IF (KRESLT.NE.0) RETURN ELSE NCALL = NCALL + 1 NAMEST(NCALL) = 'FMATAN' IF (NTRACE.NE.0) CALL FMNTR(2,MA,MA,1) KOVUN = 0 IF (MA(1).EQ.MEXPOV .OR. MA(1).EQ.MEXPUN) KOVUN = 1 NDSAVE = NDIG IF (NCALL.EQ.1) THEN K = MAX(NGRD52-1,2) NDIG = MAX(NDIG+K,2) IF (NDIG.GT.NDG2MX) THEN KFLAG = -9 CALL FMWARN NDIG = NDSAVE KRESLT = 12 CALL FMRSLT(MA,MA,MB,KRESLT) IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1) NCALL = NCALL - 1 RETURN ENDIF ENDIF KASAVE = KACCSW KACCSW = 0 MXSAVE = MXEXP MXEXP = MXEXP2 ENDIF C MACCA = MA(0) CALL FMEQ2(MA,M05,NDSAVE,NDIG,0) M05(0) = NINT(NDIG*ALOGM2) C C If MA.GE.1 work with 1/MA. C MA1 = MA(1) MA2 = MA(2) M05(2) = ABS(M05(2)) IF (MA1.GE.1) THEN CALL FMI2M(1,MB) CALL FMDIV(MB,M05,M05) ENDIF C KRSAVE = KRAD KRAD = 1 KWRNSV = KWARN C X = M05(1) XM = MXBASE C C In case pi has not been computed at the current precision C and will be needed here, get it to full precision first C to avoid repeated calls at increasing precision during C Newton iteration. C IF (MA1.GE.1 .OR. KRSAVE.EQ.0) THEN IF (MBSPI.NE.MBASE .OR. NDIGPI.LT.NDIG) THEN NDSV = NDIG NDIG = MIN(NDIG+2,NDG2MX) NCALL = NCALL + 1 NAMEST(NCALL) = 'NOEQ ' CALL FMPI(MPISAV) NCALL = NCALL - 1 NDIG = NDSV ENDIF ENDIF C C If the argument is small, use the Taylor series, C otherwise use Newton iteration. C IF (X*DLOGMB.LT.-5.0D0*LOG(XM)) THEN KWARN = 0 CALL FMEQ(M05,MB) IF (MB(1).LE.-NDIG) GO TO 130 CALL FMSQR(M05,M06) J = 3 NDSAV1 = NDIG C 110 CALL FMMPY(M05,M06,M05) IF (M05(1).NE.MUNKNO) M05(2) = -M05(2) CALL FMDIVI(M05,J,M03) NDIG = NDSAV1 CALL FMADD(MB,M03,MB) IF (KFLAG.NE.0) THEN KFLAG = 0 GO TO 130 ENDIF NDIG = NDSAV1 - INT((MB(1)-M03(1))) IF (NDIG.LT.2) NDIG = 2 J = J + 2 GO TO 110 ELSE C CALL FMM2DP(M05,X) X = ATAN(X) CALL FMDPM(X,MB) CALL FMDIG(NSTACK,KST) C C Newton iteration. C DO 120 J = 1, KST NDIG = NSTACK(J) CALL FMSIN(MB,M06) CALL FMSQR(M06,M03) CALL FMI2M(1,M04) CALL FMSUB(M04,M03,M03) CALL FMSQRT(M03,M04) CALL FMDIV(M06,M04,M04) CALL FMSUB(M04,M05,M04) CALL FMMPY(M03,M04,M04) CALL FMSUB(MB,M04,MB) 120 CONTINUE MB(0) = NINT(NDIG*ALOGM2) ENDIF C C If MA.GE.1 use pi/2 - ATAN(1/MA) C 130 IF (MA1.GE.1) THEN CALL FMDIVI(MPISAV,2,M06) CALL FMSUB(M06,MB,MB) ENDIF C C Convert to degrees if necessary, round and return. C KRAD = KRSAVE IF (KRAD.EQ.0) THEN CALL FMMPYI(MB,180,MB) CALL FMDIV(MB,MPISAV,MB) ENDIF IF (MB(1).NE.MUNKNO .AND. MA2.LT.0) MB(2) = -MB(2) C IF (KFLAG.EQ.1) KFLAG = 0 KWARN = KWRNSV MACMAX = NINT((NDSAVE-1)*ALOGM2 + LOG(REAL(ABS(MB(2))+1))/0.69315) MB(0) = MIN(MB(0),MACCA,MACMAX) CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN) RETURN END SUBROUTINE FMATN2(MA,MB,MC) C C MC = ATAN2(MA,MB) C C MC is returned as the angle between -pi and pi (or -180 and 180 if C degree mode is selected) for which TAN(MC) = MA/MB. MC is an angle C for the point (MB,MA) in polar coordinates. C C IMPLICIT NONE C INTEGER NDIGMX,NBITS,LPACK,LUNPCK,LMWA,LJSUMS,LMBUFF C PARAMETER ( NDIGMX=256 , NBITS=64 , * LPACK = (NDIGMX+1)/2+1 , LUNPCK = (6*NDIGMX)/5+20, * LMWA = 2*LUNPCK , LJSUMS = 8*(LUNPCK+2), * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000+6 ) C DOUBLE PRECISION MA(0:LUNPCK),MB(0:LUNPCK),MC(0:LUNPCK) C INTEGER NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC,KFLAG, * KWARN,KROUND,KSWIDE,KESWCH,KDEBUG DOUBLE PRECISION MBASE C COMMON /FMUSER/ MBASE,NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC, * KFLAG,KWARN,KROUND,KSWIDE,KESWCH,KDEBUG C DOUBLE PRECISION DPMAX INTEGER NCALL,KACCSW,IUNKNO,NDG2MX,INTMAX,KSUB REAL RUNKNO,SPMAX DOUBLE PRECISION MWA,MXEXP,MXEXP2,MEXPUN,MEXPOV,MUNKNO,MXBASE, * MAXINT C COMMON /FM/ MWA(LMWA),NCALL,KACCSW,MXEXP,MXEXP2,MEXPUN,MEXPOV, * MUNKNO,IUNKNO,RUNKNO,MXBASE,NDG2MX,SPMAX,DPMAX, * MAXINT,INTMAX,KSUB C C Scratch array usage during FMATN2: M01 - M06 C DOUBLE PRECISION M01,M02,M03,M04,M05,M06 C COMMON /FM1/ M01(0:LUNPCK),M02(0:LUNPCK),M03(0:LUNPCK), * M04(0:LUNPCK),M05(0:LUNPCK),M06(0:LUNPCK) C INTEGER LHASH1,LHASH2 PARAMETER (LHASH1=0 , LHASH2=256) DOUBLE PRECISION DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB INTEGER NDIGPI,NDIGE,NDIGLB,NDIGLI,KHASHT,KHASHV, * NGRD21,NGRD52,NGRD22 DOUBLE PRECISION MBSPI,MBSE,MBSLB,MBSLI,MPISAV,MESAV,MLBSAV, * MLN1,MLN2,MLN3,MLN4,MBLOGS,MEXPAB REAL ALOGMB,ALOGMT,ALOGM2,ALOGMX C COMMON /FMSAVE/ NDIGPI,NDIGE,NDIGLB,NDIGLI,MBSPI,MBSE,MBSLB,MBSLI, * MPISAV(0:LUNPCK),MESAV(0:LUNPCK),MLBSAV(0:LUNPCK), * MLN1(0:LUNPCK),MLN2(0:LUNPCK),MLN3(0:LUNPCK), * MLN4(0:LUNPCK),MBLOGS,MEXPAB,ALOGMB,ALOGM2,ALOGMX, * ALOGMT,DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB,KHASHT(LHASH1:LHASH2), * KHASHV(LHASH1:LHASH2),NGRD21,NGRD52,NGRD22 C CHARACTER CMBUFF,CMCHAR CHARACTER *6 NAMEST C COMMON /FMBUFF/ CMBUFF(LMBUFF),NAMEST(0:50),CMCHAR C DOUBLE PRECISION MACCA,MACCB,MACMAX,MXEXP1,MXSAVE INTEGER JQUAD,K,KASAVE,KOVUN,KRESLT,KWRNSV,NDSAVE C IF (MBLOGS.NE.MBASE) CALL FMCONS IF (ABS(MA(1)).GT.MEXPAB .OR. ABS(MB(1)).GT.MEXPAB) THEN CALL FMENTR('FMATN2',MA,MB,2,MC,KRESLT,NDSAVE,MXSAVE,KASAVE, * KOVUN) IF (KRESLT.NE.0) RETURN ELSE NCALL = NCALL + 1 NAMEST(NCALL) = 'FMATN2' IF (NTRACE.NE.0) CALL FMNTR(2,MA,MB,2) KOVUN = 0 IF (MA(1).EQ.MEXPOV .OR. MA(1).EQ.MEXPUN) KOVUN = 1 NDSAVE = NDIG IF (NCALL.EQ.1) THEN K = MAX(NGRD52-1,2) NDIG = MAX(NDIG+K,2) IF (NDIG.GT.NDG2MX) THEN KFLAG = -9 CALL FMWARN NDIG = NDSAVE KRESLT = 12 CALL FMRSLT(MA,MB,MC,KRESLT) IF (NTRACE.NE.0) CALL FMNTR(1,MC,MC,1) NCALL = NCALL - 1 RETURN ENDIF ENDIF KASAVE = KACCSW KACCSW = 0 MXSAVE = MXEXP MXEXP = MXEXP2 ENDIF C KWRNSV = KWARN KWARN = 0 C MACCA = MA(0) MACCB = MB(0) CALL FMEQ2(MA,M01,NDSAVE,NDIG,0) M01(0) = NINT(NDIG*ALOGM2) CALL FMEQ2(MB,M02,NDSAVE,NDIG,0) M02(0) = NINT(NDIG*ALOGM2) C C Check for special cases. C IF (MA(1).EQ.MUNKNO .OR. MB(1).EQ.MUNKNO .OR. * (MA(2).EQ.0 .AND. MB(2).EQ.0)) THEN CALL FMIM(0,MC) MC(1) = MUNKNO MC(2) = 1 MC(0) = NINT(NDIG*ALOGM2) KFLAG = -4 GO TO 110 ENDIF C IF (MB(2).EQ.0 .AND. MA(2).GT.0) THEN IF (KRAD.EQ.0) THEN CALL FMI2M(90,MC) ELSE CALL FMPI(MC) CALL FMDIVI(MC,2,MC) ENDIF GO TO 110 ENDIF C IF (MB(2).EQ.0 .AND. MA(2).LT.0) THEN IF (KRAD.EQ.0) THEN CALL FMI2M(-90,MC) ELSE CALL FMPI(MC) CALL FMDIVI(MC,-2,MC) ENDIF GO TO 110 ENDIF C MXEXP1 = INT(MXEXP2/2.01D0) IF (MA(1).EQ.MEXPOV .AND. MB(1).LT.MXEXP1-NDIG-2) THEN IF (KRAD.EQ.0) THEN CALL FMI2M(90,MC) ELSE CALL FMPI(MC) CALL FMDIVI(MC,2,MC) ENDIF IF (M01(2).LT.0) MC(2) = -MC(2) GO TO 110 ENDIF C IF (MA(1).EQ.MEXPUN .AND. (-MB(1)).LT.MXEXP1-NDIG-2 .AND. * MB(2).LT.0) THEN IF (KRAD.EQ.0) THEN CALL FMI2M(180,MC) ELSE CALL FMPI(MC) ENDIF IF (M01(2).LT.0) MC(2) = -MC(2) GO TO 110 ENDIF C IF (MB(1).EQ.MEXPOV .AND. MA(1).LT.MXEXP1-NDIG-2 .AND. * MB(2).LT.0) THEN IF (KRAD.EQ.0) THEN CALL FMI2M(180,MC) ELSE CALL FMPI(MC) ENDIF IF (M01(2).LT.0) MC(2) = -MC(2) GO TO 110 ENDIF C IF (MB(1).EQ.MEXPUN .AND. MA(2).EQ.0) THEN IF (MB(2).LT.0) THEN IF (KRAD.EQ.0) THEN CALL FMI2M(180,MC) ELSE CALL FMPI(MC) ENDIF ELSE CALL FMI2M(0,MC) ENDIF GO TO 110 ENDIF C IF (MB(1).EQ.MEXPUN .AND. (-MA(1)).LT.MXEXP1-NDIG-2) THEN IF (KRAD.EQ.0) THEN CALL FMI2M(90,MC) ELSE CALL FMPI(MC) CALL FMDIVI(MC,2,MC) ENDIF IF (M01(2).LT.0) MC(2) = -MC(2) GO TO 110 ENDIF C C Determine the quadrant for the result, then use FMATAN. C IF (MA(2).GE.0 .AND. MB(2).GT.0) JQUAD = 1 IF (MA(2).GE.0 .AND. MB(2).LT.0) JQUAD = 2 IF (MA(2).LT.0 .AND. MB(2).LT.0) JQUAD = 3 IF (MA(2).LT.0 .AND. MB(2).GT.0) JQUAD = 4 C CALL FMDIV(M01,M02,MC) MC(2) = ABS(MC(2)) CALL FMATAN(MC,MC) C IF (JQUAD.EQ.2 .OR. JQUAD.EQ.3) THEN IF (KRAD.EQ.0) THEN CALL FMI2M(180,M05) CALL FMSUB(M05,MC,MC) ELSE CALL FMPI(M05) CALL FMSUB(M05,MC,MC) ENDIF ENDIF C IF ((JQUAD.EQ.3 .OR. JQUAD.EQ.4) .AND. MC(1).NE.MUNKNO) * MC(2) = -MC(2) C C Round the result and return. C 110 IF (KFLAG.EQ.1) KFLAG = 0 KWARN = KWRNSV MACMAX = NINT((NDSAVE-1)*ALOGM2 + LOG(REAL(ABS(MC(2))+1))/0.69315) MC(0) = MIN(MC(0),MACCA,MACCB,MACMAX) CALL FMEXIT(MC,MC,NDSAVE,MXSAVE,KASAVE,KOVUN) RETURN END SUBROUTINE FMBIG(MA) C C MA = The biggest representable FM number using the current base C and precision. C The smallest positive number is then 1.0/MA. C Because of rounding, 1.0/(1.0/MA) will then overflow. C C IMPLICIT NONE C INTEGER NDIGMX,NBITS,LPACK,LUNPCK,LMWA,LJSUMS,LMBUFF C PARAMETER ( NDIGMX=256 , NBITS=64 , * LPACK = (NDIGMX+1)/2+1 , LUNPCK = (6*NDIGMX)/5+20, * LMWA = 2*LUNPCK , LJSUMS = 8*(LUNPCK+2), * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000+6 ) C DOUBLE PRECISION MA(0:LUNPCK) C INTEGER NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC,KFLAG, * KWARN,KROUND,KSWIDE,KESWCH,KDEBUG DOUBLE PRECISION MBASE C COMMON /FMUSER/ MBASE,NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC, * KFLAG,KWARN,KROUND,KSWIDE,KESWCH,KDEBUG C DOUBLE PRECISION DPMAX INTEGER NCALL,KACCSW,IUNKNO,NDG2MX,INTMAX,KSUB REAL RUNKNO,SPMAX DOUBLE PRECISION MWA,MXEXP,MXEXP2,MEXPUN,MEXPOV,MUNKNO,MXBASE, * MAXINT C COMMON /FM/ MWA(LMWA),NCALL,KACCSW,MXEXP,MXEXP2,MEXPUN,MEXPOV, * MUNKNO,IUNKNO,RUNKNO,MXBASE,NDG2MX,SPMAX,DPMAX, * MAXINT,INTMAX,KSUB C INTEGER LHASH1,LHASH2 PARAMETER (LHASH1=0 , LHASH2=256) DOUBLE PRECISION DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB INTEGER NDIGPI,NDIGE,NDIGLB,NDIGLI,KHASHT,KHASHV, * NGRD21,NGRD52,NGRD22 DOUBLE PRECISION MBSPI,MBSE,MBSLB,MBSLI,MPISAV,MESAV,MLBSAV, * MLN1,MLN2,MLN3,MLN4,MBLOGS,MEXPAB REAL ALOGMB,ALOGMT,ALOGM2,ALOGMX C COMMON /FMSAVE/ NDIGPI,NDIGE,NDIGLB,NDIGLI,MBSPI,MBSE,MBSLB,MBSLI, * MPISAV(0:LUNPCK),MESAV(0:LUNPCK),MLBSAV(0:LUNPCK), * MLN1(0:LUNPCK),MLN2(0:LUNPCK),MLN3(0:LUNPCK), * MLN4(0:LUNPCK),MBLOGS,MEXPAB,ALOGMB,ALOGM2,ALOGMX, * ALOGMT,DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB,KHASHT(LHASH1:LHASH2), * KHASHV(LHASH1:LHASH2),NGRD21,NGRD52,NGRD22 C CHARACTER CMBUFF,CMCHAR CHARACTER *6 NAMEST C COMMON /FMBUFF/ CMBUFF(LMBUFF),NAMEST(0:50),CMCHAR C INTEGER J,N1 C NCALL = NCALL + 1 NAMEST(NCALL) = 'FMBIG ' C IF (MBLOGS.NE.MBASE) CALL FMCONS KFLAG = 0 N1 = NDIG + 1 DO 110 J = 2, N1 MA(J) = MBASE - 1 110 CONTINUE MA(1) = MXEXP + 1 MA(0) = NINT(NDIG*ALOGM2) C IF (NTRACE.NE.0) CALL FMNTR(1,MA,MA,1) NCALL = NCALL - 1 RETURN END SUBROUTINE FMCAT(MA,NCAT) C C NCAT is returned as the category of MA. This is used by the various C arithmetic routines to handle special cases such as: C 'number greater than 1' + 'underflowed result' is the first argument, C 'overflowed result' / 'overflowed result' is 'unknown'. C C NCAT range C C 1. -OV OV stands for overflowed results. C 2. (-OV , -OVTH) ( MA(1) .GE. MAXEXP+2 ) C 3. (-OVTH , -1) C 4. -1 OVTH stands for a representable C 5. (-1 , -UNTH) number near the overflow C 6. (-UNTH , -UN) threshold. C 7. -UN ( MA(1) .GE. MAXEXP-NDIG+1 ) C 8. 0 C 9. +UN UN stands for underflowed results. C 10. (+UN , +UNTH) ( MA(1) .LE. -MAXEXP-1 ) C 11. (+UNTH , +1) C 12. +1 UNTH stands for a representable C 13. (+1 , +OVTH) number near the underflow C 14. (+OVTH , +OV) threshold. C 15. +OV ( MA(1) .LE. -MAXEXP+NDIG-1 ) C 16. UNKNOWN C C IMPLICIT NONE C INTEGER NDIGMX,NBITS,LPACK,LUNPCK,LMWA,LJSUMS,LMBUFF C PARAMETER ( NDIGMX=256 , NBITS=64 , * LPACK = (NDIGMX+1)/2+1 , LUNPCK = (6*NDIGMX)/5+20, * LMWA = 2*LUNPCK , LJSUMS = 8*(LUNPCK+2), * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000+6 ) C DOUBLE PRECISION MA(0:LUNPCK) INTEGER NCAT C INTEGER NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC,KFLAG, * KWARN,KROUND,KSWIDE,KESWCH,KDEBUG DOUBLE PRECISION MBASE C COMMON /FMUSER/ MBASE,NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC, * KFLAG,KWARN,KROUND,KSWIDE,KESWCH,KDEBUG C DOUBLE PRECISION DPMAX INTEGER NCALL,KACCSW,IUNKNO,NDG2MX,INTMAX,KSUB REAL RUNKNO,SPMAX DOUBLE PRECISION MWA,MXEXP,MXEXP2,MEXPUN,MEXPOV,MUNKNO,MXBASE, * MAXINT C COMMON /FM/ MWA(LMWA),NCALL,KACCSW,MXEXP,MXEXP2,MEXPUN,MEXPOV, * MUNKNO,IUNKNO,RUNKNO,MXBASE,NDG2MX,SPMAX,DPMAX, * MAXINT,INTMAX,KSUB C DOUBLE PRECISION MA2,MXEXP1 INTEGER J,NLAST C C Check for special symbols. C NCAT = 16 IF (MA(1).EQ.MUNKNO) RETURN C IF (MA(1).EQ.MEXPOV) THEN NCAT = 15 IF (MA(2).LT.0) NCAT = 1 RETURN ENDIF C IF (MA(1).EQ.MEXPUN) THEN NCAT = 9 IF (MA(2).LT.0) NCAT = 7 RETURN ENDIF C IF (MA(2).EQ.0) THEN NCAT = 8 RETURN ENDIF C C Check for +1 or -1. C MA2 = ABS(MA(2)) IF (MA(1).EQ.1 .AND. MA2.EQ.1) THEN NLAST = NDIG + 1 IF (NLAST.GE.3) THEN DO 110 J = 3, NLAST IF (MA(J).NE.0) GO TO 120 110 CONTINUE ENDIF NCAT = 12 IF (MA(2).LT.0) NCAT = 4 RETURN ENDIF C 120 MXEXP1 = INT(MXEXP2/2.01D0) IF (MA(1).GE.MXEXP1-NDIG+1) THEN NCAT = 14 IF (MA(2).LT.0) NCAT = 2 RETURN ENDIF C IF (MA(1).GE.1) THEN NCAT = 13 IF (MA(2).LT.0) NCAT = 3 RETURN ENDIF C IF (MA(1).GE.-MXEXP1+NDIG) THEN NCAT = 11 IF (MA(2).LT.0) NCAT = 5 RETURN ENDIF C IF (MA(1).GE.-MXEXP2) THEN NCAT = 10 IF (MA(2).LT.0) NCAT = 6 RETURN ENDIF C RETURN END SUBROUTINE FMCHSH(MA,MB,MC) C C MB = COSH(MA), MC = SINH(MA) C C If both the hyperbolic sine and cosine are needed, this routine C is faster than calling both FMCOSH and FMSINH. C C MB and MC must be distinct arrays. C C IMPLICIT NONE C INTEGER NDIGMX,NBITS,LPACK,LUNPCK,LMWA,LJSUMS,LMBUFF C PARAMETER ( NDIGMX=256 , NBITS=64 , * LPACK = (NDIGMX+1)/2+1 , LUNPCK = (6*NDIGMX)/5+20, * LMWA = 2*LUNPCK , LJSUMS = 8*(LUNPCK+2), * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000+6 ) C DOUBLE PRECISION MA(0:LUNPCK),MB(0:LUNPCK),MC(0:LUNPCK) C INTEGER NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC,KFLAG, * KWARN,KROUND,KSWIDE,KESWCH,KDEBUG DOUBLE PRECISION MBASE C COMMON /FMUSER/ MBASE,NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC, * KFLAG,KWARN,KROUND,KSWIDE,KESWCH,KDEBUG C DOUBLE PRECISION DPMAX INTEGER NCALL,KACCSW,IUNKNO,NDG2MX,INTMAX,KSUB REAL RUNKNO,SPMAX DOUBLE PRECISION MWA,MXEXP,MXEXP2,MEXPUN,MEXPOV,MUNKNO,MXBASE, * MAXINT C COMMON /FM/ MWA(LMWA),NCALL,KACCSW,MXEXP,MXEXP2,MEXPUN,MEXPOV, * MUNKNO,IUNKNO,RUNKNO,MXBASE,NDG2MX,SPMAX,DPMAX, * MAXINT,INTMAX,KSUB C INTEGER LHASH1,LHASH2 PARAMETER (LHASH1=0 , LHASH2=256) DOUBLE PRECISION DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB INTEGER NDIGPI,NDIGE,NDIGLB,NDIGLI,KHASHT,KHASHV, * NGRD21,NGRD52,NGRD22 DOUBLE PRECISION MBSPI,MBSE,MBSLB,MBSLI,MPISAV,MESAV,MLBSAV, * MLN1,MLN2,MLN3,MLN4,MBLOGS,MEXPAB REAL ALOGMB,ALOGMT,ALOGM2,ALOGMX C COMMON /FMSAVE/ NDIGPI,NDIGE,NDIGLB,NDIGLI,MBSPI,MBSE,MBSLB,MBSLI, * MPISAV(0:LUNPCK),MESAV(0:LUNPCK),MLBSAV(0:LUNPCK), * MLN1(0:LUNPCK),MLN2(0:LUNPCK),MLN3(0:LUNPCK), * MLN4(0:LUNPCK),MBLOGS,MEXPAB,ALOGMB,ALOGM2,ALOGMX, * ALOGMT,DLOGMB,DLOGTN,DLOGTW,DLOGTP,DLOGPI,DPPI, * DPEPS,DLOGEB,KHASHT(LHASH1:LHASH2), * KHASHV(LHASH1:LHASH2),NGRD21,NGRD52,NGRD22 C C Scratch array usage during FMCHSH: M01 - M04 C DOUBLE PRECISION M01,M02,M03,M04,M05,M06 C COMMON /FM1/ M01(0:LUNPCK),M02(0:LUNPCK),M03(0:LUNPCK), * M04(0:LUNPCK),M05(0:LUNPCK),M06(0:LUNPCK) C CHARACTER CMBUFF,CMCHAR CHARACTER *6 NAMEST C COMMON /FMBUFF/ CMBUFF(LMBUFF),NAMEST(0:50),CMCHAR C DOUBLE PRECISION MA2,MACCA,MACMAX,MXSAVE INTEGER K,KASAVE,KOVUN,KRESLT,KWRNSV,NCSAVE,NDSAVE C IF (MBLOGS.NE.MBASE) CALL FMCONS MACCA = MA(0) MA2 = MA(2) IF (ABS(MA(1)).GT.MEXPAB) THEN NCSAVE = NCALL CALL FMENTR('FMCHSH',MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE, * KOVUN) IF (MA(1).EQ.MUNKNO) KOVUN = 2 NCALL = NCSAVE + 1 CALL FMEQ(MA,M04) M04(0) = NINT(NDIG*ALOGM2) M04(2) = ABS(M04(2)) CALL FMCOSH(M04,MB) CALL FMSINH(M04,MC) GO TO 110 ELSE NCALL = NCALL + 1 NAMEST(NCALL) = 'FMCHSH' IF (NTRACE.NE.0) CALL FMNTR(2,MA,MA,1) KOVUN = 0 IF (MA(1).EQ.MEXPOV .OR. MA(1).EQ.MEXPUN) KOVUN = 1 NDSAVE = NDIG IF (NCALL.EQ.1) THEN K = MAX(NGRD52,2) NDIG = MAX(NDIG+K,2) IF (NDIG.GT.NDG2MX) THEN NCALL = NCALL - 1 NDIG = NDSAVE CALL FMEQ(MA,M04) CALL FMCOSH(M04,MB) CALL FMSINH(M04,MC) KFLAG = -9 RETURN ENDIF ENDIF KASAVE = KACCSW KACCSW = 0 MXSAVE = MXEXP MXEXP = MXEXP2 ENDIF C CALL FMEQ2(MA,M04,NDSAVE,NDIG,0) M04(0) = NINT(NDIG*ALOGM2) M04(2) = ABS(M04(2)) C K = 1 IF (M04(1).EQ.0 .AND. M04(2).NE.0) THEN IF (MBASE/M04(2).GE.100) K = 2 ENDIF IF (M04(1).GE.0 .AND. M04(2).NE.0 .AND. K.EQ.1) THEN CALL FMCOSH(M04,MB) IF (MB(1).GT.NDIG) THEN CALL FMEQ(MB,MC) GO TO 110 ENDIF CALL FMSQR(MB,M03) CALL FMI2M(-1,M02) CALL FMADD(M03,M02,M03) CALL FMSQRT(M03,MC) ELSE CALL FMSINH(M04,MC) CALL FMSQR(MC,M03) CALL FMI2M(1,M02) CALL FMADD(M03,M02,M03) CALL FMSQRT(M03,MB) ENDIF C C Round and return. C 110 MACMAX = NINT((NDSAVE-1)*ALOGM2 + LOG(REAL(ABS(MC(2))+1))/0.69315) MC(0) = MIN(MC(0),MACCA,MACMAX) IF (MA2.LT.0 .AND. MC(1).NE.MUNKNO) MC(2) = -MC(2) CALL FMEQ2(MC,MC,NDIG,NDSAVE,1) MACMAX = NINT((NDSAVE-1)*ALOGM2 + LOG(REAL(ABS(MB(2))+1))/0.69315) MB(0) = MIN(MB(0),MACCA,MACMAX) IF (KOVUN.EQ.2) THEN KWRNSV = KWARN KWARN = 0 ENDIF CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN) IF (KOVUN.EQ.2) THEN KWARN = KWRNSV ENDIF IF (NTRACE.NE.0) THEN IF (ABS(NTRACE).GE.1 .AND. NCALL+1.LE.LVLTRC) THEN IF (NTRACE.LT.0) THEN CALL FMNTRJ(MC,NDIG) ELSE CALL FMPRNT(MC) ENDIF ENDIF ENDIF RETURN END FUNCTION FMCOMP(MA,LREL,MB) C C Logical comparison of FM numbers MA and MB. C C LREL is a CHARACTER *2 description of the comparison to be done: C LREL = 'EQ' returns FMCOMP = .TRUE. if MA.EQ.MB C = 'NE', 'GE', 'GT', 'LE', 'LT' also work like a logical IF. C C For comparisons involving 'UNKNOWN' or two identical special symbols C such as +OVERFLOW,'EQ',+OVERFLOW, FMCOMP is returned FALSE and a C KFLAG = -4 error condition is returned. C C Some compilers object to functions with side effects such as C changing KFLAG or other common variables. Blocks of code that C modify common are identified by: C C DELETE START C ... C C DELETE STOP C These may be removed or commented out to produce a function without C side effects. This disables trace printing in FMCOMP, and error C codes are not returned in KFLAG. C C IMPLICIT NONE C INTEGER NDIGMX,NBITS,LPACK,LUNPCK,LMWA,LJSUMS,LMBUFF C PARAMETER ( NDIGMX=256 , NBITS=64 , * LPACK = (NDIGMX+1)/2+1 , LUNPCK = (6*NDIGMX)/5+20, * LMWA = 2*LUNPCK , LJSUMS = 8*(LUNPCK+2), * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000+6 ) C LOGICAL FMCOMP CHARACTER *2 JREL,LREL DOUBLE PRECISION MA(0:LUNPCK),MB(0:LUNPCK) C INTEGER NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC,KFLAG, * KWARN,KROUND,KSWIDE,KESWCH,KDEBUG DOUBLE PRECISION MBASE C COMMON /FMUSER/ MBASE,NDIG,JFORM1,JFORM2,KRAD,KW,NTRACE,LVLTRC, * KFLAG,KWARN,KROUND,KSWIDE,KESWCH,KDEBUG C DOUBLE PRECISION DPMAX INTEGER NCALL,KACCSW,IUNKNO,NDG2MX,INTMAX,KSUB REAL RUNKNO,SPMAX DOUBLE PRECISION MWA,MXEXP,MXEXP2,MEXPUN,MEXPOV,MUNKNO,MXBASE, *