C ALGORITHM 693, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 17, NO. 2, PP. 273-283. JUNE, 1991. C C FM 1.0 David M Smith 1-06-90 C C C The FM package performs floating point multiple precision arithmetic. C C Before calling any routine in the package, several variables in the C common blocks /FMUSER/, /FM/, and /FMSAVE/ must be initialized. C These three common blocks contain information which must be saved C between calls, so they should be declared in the main 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 common block /FMUSER/. C C C JBASE is the base in which the arithmetic is done. C JBASE must be bigger than one, and less than or equal C to the square root of the largest representable integer. C For best efficiency JBASE should be about 1/4 to 1/2 as big as C the square root of the largest representable integer. C Input and output conversion are much faster when JBASE is a C power of ten. C C NDIG is the number of base JBASE digits which are carried in the C multiple precision numbers. NDIG must be at least two. C The upper limit for NDIG is defined in the PARAMETER statement C at the top of each routine and is restricted only by the amount C of memory available. C 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 JBASE and is stored in NDIG+1 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 JBASE) 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+1 words. C C The unpacked format of a floating multiple precision number is as C follows. The number is kept in an integer array with the first C element containing the exponent and each of the succeeding NDIG C locations containing one digit of the mantissa, expressed in base C JBASE. The exponent is a power of JBASE and the implied radix point C is 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 JBASE = 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 Because of normalization, the equivalent number of base 10 C significant digits for an FM number may be as small as C LOG10(JBASE)*(NDIG-1) + 1. C C C Subroutine FMOUT performs conversion of an FM number to base 10 and C formats it for output as a character array. 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 = Number of significant digits to display (if JFORM1 = 0, 1). C If JFORM2.EQ.0 then a default number of digits is chosen. C The default is roughly the full precision of the number. C JFORM2 = Number of digits after the decimal point (if JFORM1 = 2). C See the FMOUT documentation for more details. C C C KW is the unit number to be used for all output from the package, C including error and warning messages, and trace output. C 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 JBASE 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 JBASE format. C 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 JBASE format means the number is C printed as it appears in the array as an exponent followed by NDIG C base JBASE digits. C C C KFLAG is a condition parameter returned by the package. Negative C values indicate conditions for which a warning message will C be printed unless KWARN = 0. Positive values indicate C conditions which 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 MXNDIG. C = -2 JBASE 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 to FMINP was not legal. C = -8 The output array for FMOUT was not large enough. C = -9 Precision could not be raised enough to provide all C requested guard digits. UNKNOWN was returned. C = -10 An FM input argument was too small in magnitude to C convert in FMM2SP or FMM2DP. Zero was C returned. C C When a negative KFLAG condition is encountered the routine calls C FMWARN, which uses the value of KWARN as follows. 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 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 C KRAD = 0 Causes all angles in the trigonometric functions and C inverse functions to be measured in degrees. C = 1 Causes all angles to be measured in radians. C C C KROUND = 0 Causes all final results to be chopped (rounded toward C zero). Intermediate results are rounded. C = 1 Causes all results to be rounded to the nearest FM C number, or to the value with an even last digit if C the result is halfway between two FM numbers. C C C Here is a list of the routines in FM which are designed to be called C by the user. All are subroutines except logical function FMCOMP. C MA, MB, MC refer to FM format numbers. C C FMABS(MA,MB) MB = ABS(MA) C FMACOS(MA,MB) MB = ACOS(MA) C FMADD(MA,MB,MC) MC = MA + MB C FMASIN(MA,MB) MB = ASIN(MA) C FMATAN(MA,MB) MB = ATAN(MA) C FMATN2(MA,MB,MC) MC = ATAN2(MA,MB) C FMBIG(MA) MA = Biggest FM number less than overflow. C FMCOMP(MA,LREL,MB) Logical comparison of MA and MB. LREL is a C CHARACTER*2 value identifying which C comparison is made. C Example: IF (FMCOMP(MA,'GE',MB)) ... C FMCOS(MA,MB) MB = COS(MA) C FMCOSH(MA,MB) MB = COSH(MA) C FMDIG(NSTACK,KST) Find a set of precisions to use during Newton C iteration for finding a simple root C starting with about double precision C accuracy. C FMDIM(MA,MB,MC) MC = DIM(MA,MB) C FMDIV(MA,MB,MC) MC = MA/MB C FMDIVI(MA,INT,MB) MB = MA/INT for one word integer INT. C FMDP2M(X,MA) MA = X conversion from double precision to FM. C FMEQU(MA,MB,NDA,NDB) MB = MA where MA has NDA digits and MB has C NDB digits. C FMEXP(MA,MB) MB = EXP(MA) C FMI2M(INT,MA) MA = INT conversion from one word integer to FM. C FMINP(LINE,MA,LA,LB) MA = LINE input conversion of LINE(LA) through C LINE(LB) from characters to FM. C FMINT(MA,MB) MB = INT(MA) integer part of MA. C FMIPWR(MA,INT,MB) MB = MA**INT raise FM number to a one word C integer power. C FMLG10(MA,MB) MB = LOG10(MA) C FMLN(MA,MB) MB = LOG(MA) C FMLNI(INT,MA) MA = LOG(INT) natural log of one word integer. C FMM2DP(MA,X) X = MA conversion from FM to double precision. C FMM2I(MA,INT) INT = MA conversion from FM to one word integer. C FMM2SP(MA,X) X = MA conversion from FM to single precision. C FMMAX(MA,MB,MC) MC = MAX(MA,MB) C FMMIN(MA,MB,MC) MC = MIN(MA,MB) C FMMOD(MA,MB,MC) MC = MA mod MB C FMMPY(MA,MB,MC) MC = MA*MB C FMMPYI(MA,INT,MB) MB = MA*INT multiplication by one word integer. C FMNINT(MA,MB) MB = NINT(MA) nearest integer. C FMOUT(MA,LINE,LB) LINE = MA conversion from FM to character. LB C is the size of array LINE. C FMPI(MA) MA = pi C FMPRNT(MA) Print MA using current format. C FMPWR(MA,MB,MC) MC = MA**MB 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 FMSIGN(MA,MB,MC) MC = SIGN(MA,MB) sign transfer. C FMSIN(MA,MB) MB = SIN(MA) C FMSINH(MA,MB) MB = SINH(MA) C FMSP2M(X,MA) MA = X conversion from single precision to FM. C FMSQRT(MA,MB) MB = SQRT(MA) C FMSUB(MA,MB,MC) MC = MA - MB C FMTAN(MA,MB) MB = TAN(MA) C FMTANH(MA,MB) MB = TANH(MA) C FMULP(MA,MB) MB = 1 Unit in the Last Place of 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 packed versions have the same names except 'FM' is C replaced by 'FP' at the start of each name. C C C NOTES ON 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 MXNDIG and NBITS. C MXNDIG is the maximum value the user may set for NDIG. C NBITS is the number of bits used to represent integers. C C FM numbers in packed format have size LPACK, and those in unpacked C format have size LUNPCK. C C C PORTABILITY NOTES. C C In routines FMEQU and FMSUB there is code which attempts to C determine if two input arrays refer to identical memory locations. C Some optimizing compilers assume the arrays must be distinct and C may remove code which would then be redundant. This code removal C could cause errors, so the tests are done in a way which should C keep the compiler from removing code. The current version works C correctly on all compilers tested. Computing SIN(1.0) in radian C mode should reveal whether other compilers handle it correctly. C If there is a problem, SIN(1) gives 0.999... instead of 0.841.... C To fix such a problem, MB can be copied to a local temporary array C and then negated in FMSUB before calling FMADD2. For FMEQU C simply set KSAME = 0 after the section which tries to determine if C MA and MB are the same array. This makes both routines run slower. C A simpler fix which often works is to re-compile at a lower C optimization (e.g., OPT=0). C C In FMSET there is machine-dependent code which attempts to C approximate the largest one word integer value. The current code C works on all machines tested, but if an FM run fails, check the C MAXINT loop in FMSET in addition to the three routines mentioned C above. Values for SPMAX and DPMAX are also defined which should C be set to values near overflow for single precision and double C precision. C C Using the CFT77 compiler on a Cray X-MP computer there are C problems using a large value for JBASE when the default 46-bit C integer arithmetic mode is used. In particular, FMSET chooses C too large a JBASE value since some of the arithmetic in the C MAXINT loop is done with 64-bit precision. Setting JBASE = 10**6 C or less may be ok, but the preferred solution is to select the C 64-bit integer compiler option. Then JBASE = 10**9 can be used. C C -------------------------------------------------------------------- C C SUBROUTINE FMSET(NPREC) C C Initialize the values in common which 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 JBASE is set to the largest permissible power of ten. C JFORM1 and JFORM2 are set to 1PE format displaying NPREC significant C digits. C C The trace option is set off, and the mode for angles in trig C functions is set to radians. The rounding mode is set to symmetric C rounding. 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 MXNDIG is the maximum value the user can set NDIG, C NBITS is the number of bits per integer word. C PARAMETER ( MXNDIG=256 , NBITS=32 , C C Define the array sizes: C * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C C Here are all the common blocks used in FM. C C /FMUSER/, /FM/, and /FMSAVE/ should also be declared in the main C program, because some compilers allocate and free space used for C labelled common which is declared only in subprograms. This causes C the saved information to be lost. C C FMUSER contains values which may need to be changed by C the calling program. C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C C FM contains the work array used by the low-level C arithmetic routines, definitions for overflow and C underflow thresholds, etc. C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C C FMSAVE contains information about saved constants. C COMMON /FMSAVE/ NDIGPI,NJBPI,NDIGE,NJBE,NDIGLB,NJBLB,NDIGLI,NJBLI, * MPISAV(LUNPCK),MESAV(LUNPCK),MLBSAV(LUNPCK), * MLN1(LUNPCK),MLN2(LUNPCK),MLN3(LUNPCK), * MLN4(LUNPCK) C C MJSUMS is an array which can contain several FM numbers C being used to accumulate the concurrent sums in FMEXP2 C and FMSIN2. When MXNDIG is 256 eight is about the maximum C number of sums needed (but this depends on JBASE). For C larger MXNDIG dimensioning MJSUMS to hold more than eight C FM numbers could speed the functions up. C COMMON /FMSUMS/ MJSUMS(LJSUMS) C C MBUFF is a character array used by FMPRNT for printing C output from FMOUT. This array may also be used for calls C to FMOUT from outside the FM package. C CHARACTER MBUFF C COMMON /FMBUFF/ MBUFF(LMBUFF) C C FM1 contains scratch arrays for temporary storage of FM C numbers while computing various functions. C COMMON /FM1/ M01(LUNPCK),M02(LUNPCK),M03(LUNPCK),M04(LUNPCK), * M05(LUNPCK),M06(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 COMMON /FMPCK/ MPA(LUNPCK),MPB(LUNPCK) C DOUBLE PRECISION TEMP C C KW is the unit number for all output from the FM package. C This includes trace output and error 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 to C 2**31 - 1 = 2 147 483 647. C C Setting MAXINT to a smaller number is ok, but this C unnecessarily restricts the permissible range of C JBASE and MXEXP. Too small a value of MAXINT will C also slow the elementary functions like SIN, EXP, C etc., since MXBASE = SQRT(MAXINT) is used to C determine how many terms can be combined when C summing series. C C The following code should set MAXINT to the largest C representable number of the form 2**J - 1. 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 on a Fortran-77 compiler C on each of these machines: C IBM 3090, CDC 176, CRAY XMP, MAGNUSON, MIPS 800, C SUN 4/260, IBM PC, MACINTOSH, MACINTOSH II, C COMPAQ 386/20, TRS-80/16. C However, even different versions of the same C compiler may react differently, so check the value C of MAXINT if there are problems installing FM on C a new machine. C MAXINT = 3 110 L = 2*MAXINT + 1 IF (INT(L/2).EQ.MAXINT) THEN MAXINT = L GO TO 110 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 MXNDG2 is the maximum value for NDIG which can be used C internally. Several of the FM routines may raise C NDIG above MXNDIG temporarily, in order to C compute correctly 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 MXNDG2 = LUNPCK - 1 C C MXBASE is the maximum value for JBASE. C TEMP = MAXINT MXBASE = SQRT(TEMP) C C JBASE is the currently used base for arithmetic. C K = LOG10(REAL(MXBASE)) JBASE = 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.MXNDIG) THEN NDIG = MAX(2,MIN(MXNDIG,NDIG)) WRITE (KW,120) NPREC,NDIG 120 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 Some constants which are often needed are stored in the C maximum precision which they have been computed with 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 NJBPI is the value of JBASE for the currently stored C value of pi. C NJBPI = 0 C C NDIGE is the number of digits available in the currently C stored value of e (MESAV). C NDIGE = 0 C C NJBE is the value of JBASE for the currently stored C value of e. C NJBE = 0 C C NDIGLB is the number of digits available in the currently C stored value of LN(JBASE) (MLBSAV). C NDIGLB = 0 C C NJBLB is the value of JBASE for the currently stored C value of LN(JBASE). C NJBLB = 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 NJBLI is the value of JBASE for the currently stored C values of MLN1 - MLN4. C NJBLI = 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(MAXINT) is greater C than MXBASE**(MXEXP2+1). C C The overflow threshold is JBASE**(MXEXP+1), and the C underflow threshold is JBASE**(-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 TEMP = MXBASE MXEXP = (TEMP*TEMP)/(2.0*LOG(TEMP)) - 1.0 MXEXP2 = 2*MXEXP + MXEXP/100 C C KARGSW is a switch used by some of the elementary function C routines to disable argument checking while doing C calculations where no exceptions can occur. C See FMARGS for a description of argument checking. C KARGSW = 0 is the normal setting, C = 1 means argument checking is disabled. C KARGSW = 0 C C KEXPUN is the exponent used as a special symbol for C underflowed results. C KEXPUN = -MXEXP2 - 5*MXNDIG C C KEXPOV is the exponent used as a special symbol for C overflowed results. C KEXPOV = -KEXPUN C C KUNKNO is the exponent used as a special symbol for C unknown FM results (1/0, SQRT(-3.0), etc). C KUNKNO = KEXPOV + 5*MXNDIG 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 which 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 which is likely to C make any subsequent results obviously wrong which 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 = -MXBASE*MXBASE 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 RETURN END SUBROUTINE FMABS(MA,MB) C C MB = ABS(MA) C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C NCALL = NCALL + 1 KSTACK(NCALL) = 1 IF (NTRACE.NE.0) CALL FMNTR(2,MA,MA,1) C KFLAG = 0 KWSV = KWARN KWARN = 0 CALL FMEQU(MA,MB,NDIG,NDIG) MB(2) = ABS(MB(2)) KWARN = KWSV C 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 PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C C Scratch array usage during FMACOS: M01 - M06 C COMMON /FM1/ M01(LUNPCK),M02(LUNPCK),M03(LUNPCK),M04(LUNPCK), * M05(LUNPCK),M06(LUNPCK) C CALL FMENTR(2,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN) IF (KRESLT.NE.0) RETURN C MA2 = MA(2) CALL FMEQU(MA,MB,NDSAVE,NDIG) 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 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 PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C NCALL = NCALL + 1 KSTACK(NCALL) = 3 IF (NTRACE.NE.0) CALL FMNTR(2,MA,MB,2) C CALL FMADD2(MA,MB,MC) C IF (NTRACE.NE.0) CALL FMNTR(1,MC,MC,1) NCALL = NCALL - 1 RETURN END SUBROUTINE FMADD2(MA,MB,MC) C C Internal addition routine. MC = MA + MB C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C IF (KARGSW.NE.1) THEN CALL FMARGS(3,2,MA,MB,KRESLT) IF (KRESLT.NE.0) THEN CALL FMRSLT(MA,MB,MC,KRESLT) RETURN ENDIF ELSE IF (MA(2).EQ.0) THEN CALL FMEQU(MB,MC,NDIG,NDIG) KFLAG = 1 RETURN ENDIF IF (MB(2).EQ.0) THEN CALL FMEQU(MA,MC,NDIG,NDIG) KFLAG = 1 RETURN ENDIF ENDIF C KFLAG = 0 N1 = NDIG + 1 C C NGUARD is the number of guard digits used. C IF (NCALL.GT.1) THEN NGUARD = 1 ELSE B = JBASE NGUARD = 5.0*LOG(10.0)/LOG(B) + 2.0 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) 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 NLAST = NDIG + 1 C DO 110 J = 2, NLAST IF (ABS(MA(J)).GT.ABS(MB(J))) THEN JCOMP = 1 GO TO 120 ENDIF IF (ABS(MB(J)).GT.ABS(MA(J))) THEN JCOMP = 3 GO TO 120 ENDIF 110 CONTINUE C JCOMP = 2 C 120 IF (JCOMP.LT.3) THEN KBIGMA = 1 IF (MA2.LT.0) JSIGN = -1 IF (MA2*MB2.GT.0) THEN CALL FMADDP(MA,MB,NMWA) ELSE CALL FMADDN(MA,MB,NGUARD,NMWA) ENDIF ELSE KBIGMA = 0 IF (MB2.LT.0) JSIGN = -1 IF (MA2*MB2.GT.0) THEN CALL FMADDP(MB,MA,NMWA) ELSE CALL FMADDN(MB,MA,NGUARD,NMWA) ENDIF ENDIF MA(2) = MA2 MB(2) = MB2 C C C Round the result. C CALL FMRND(NDIG,NGUARD,0) C C See if the result is equal to one of the input arguments. C K = ABS(MA(1)-MB(1)) IF (K.LT.NDIG) GO TO 150 IF (K.GT.NDIG+1) THEN KFLAG = 1 GO TO 150 ENDIF C N2 = NDIG + 4 IF (KBIGMA.EQ.1) THEN DO 130 J = 3, N1 IF (MWA(N2-J).NE.MA(N2-J)) GO TO 150 130 CONTINUE IF (MWA(1).NE.MA(1)) GO TO 150 IF (MWA(2).NE.ABS(MA(2))) GO TO 150 ELSE DO 140 J = 3, N1 IF (MWA(N2-J).NE.MB(N2-J)) GO TO 150 140 CONTINUE IF (MWA(1).NE.MB(1)) GO TO 150 IF (MWA(2).NE.ABS(MB(2))) GO TO 150 ENDIF KFLAG = 1 C C Transfer to MC and fix the sign of the result. C 150 CALL FMMOVE(MC) IF (JSIGN.LT.0) MC(2) = -MC(2) C 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 which will be used. C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C N1 = NDIG + 1 C C Check for an insignificant operand. C K = MA(1) - MB(1) IF (K.GE.NDIG+2) THEN DO 110 J = 1, N1 MWA(J) = MA(J) 110 CONTINUE MWA(N1+1) = 0 RETURN ENDIF 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 IF (KP2.LE.N1) THEN DO 130 J = KP2, N1 MWA(J) = MA(J) - MB(J-K) 130 CONTINUE ENDIF N2 = NDIG + 2 IF (N2-K.LE.1) N2 = 2 + K NK = MIN(NMWA,N1+K) IF (N2.LE.NK) THEN DO 140 J = N2, NK MWA(J) = -MB(J-K) 140 CONTINUE ENDIF NK1 = NK + 1 IF (NK1.LE.NMWA) THEN DO 150 J = NK1, NMWA MWA(J) = 0 150 CONTINUE ENDIF C C Normalize. Fix the sign of any negative digit. C IF (K.GT.0) THEN KB = NMWA - KP2 + 1 K2 = NMWA + 1 DO 160 J = 1, KB IF (MWA(K2-J).LT.0) THEN MWA(K2-J) = MWA(K2-J) + JBASE MWA(NMWA-J) = MWA(NMWA-J) - 1 ENDIF 160 CONTINUE KPT = KP2 170 KPT = KPT - 1 IF (MWA(KPT).LT.0 .AND. KPT.GE.3) THEN MWA(KPT) = MWA(KPT) + JBASE MWA(KPT-1) = MWA(KPT-1) - 1 GO TO 170 ENDIF GO TO 190 ENDIF C IF (K.EQ.0) THEN KP = N1 + 3 KPM1 = KP - 1 DO 180 J = 3, N1 IF (MWA(KP-J).LT.0) THEN MWA(KP-J) = MWA(KP-J) + JBASE MWA(KPM1-J) = MWA(KPM1-J) - 1 ENDIF 180 CONTINUE ENDIF C C Shift left if there are any leading zeros in the mantissa. C 190 DO 200 J = 2, NMWA KSH = J - 2 IF (MWA(J).GT.0) GO TO 210 200 CONTINUE MWA(1) = 0 RETURN C 210 KL = NMWA - KSH IF (KSH.GT.0) THEN DO 220 J = 2, KL L = J + KSH MWA(J) = MWA(L) 220 CONTINUE KL = KL + 1 DO 230 J = KL, NMWA MWA(J) = 0 230 CONTINUE MWA(1) = MWA(1) - KSH ENDIF C RETURN END SUBROUTINE FMADDP(MA,MB,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 which will be used. C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C N1 = NDIG + 1 C C Check for an insignificant operand. C K = MA(1) - MB(1) IF (K.GE.NDIG+1) THEN DO 110 J = 1, N1 MWA(J) = MA(J) 110 CONTINUE MWA(N1+1) = 0 RETURN ENDIF C C Add MA and MB. C KP1 = K + 1 DO 120 J = 1, KP1 MWA(J) = MA(J) 120 CONTINUE KP2 = K + 2 IF (KP2.LE.N1) THEN DO 130 J = KP2, N1 MWA(J) = MA(J) + MB(J-K) 130 CONTINUE ENDIF N2 = NDIG + 2 NK = MIN(NMWA,N1+K) IF (N2.LE.NK) THEN DO 140 J = N2, NK MWA(J) = MB(J-K) 140 CONTINUE ENDIF NK1 = NK + 1 IF (NK1.LE.NMWA) THEN DO 150 J = NK1, NMWA MWA(J) = 0 150 CONTINUE ENDIF C C Normalize. Fix any digit not less than JBASE. C IF (K.EQ.NDIG) RETURN C IF (K.GT.0) THEN KB = N1 - KP2 + 1 N2 = N1 + 1 DO 160 J = 1, KB IF (MWA(N2-J).GE.JBASE) THEN MWA(N2-J) = MWA(N2-J) - JBASE MWA(N1-J) = MWA(N1-J) + 1 ENDIF 160 CONTINUE KPT = KP2 170 KPT = KPT - 1 IF (MWA(KPT).GE.JBASE .AND. KPT.GE.3) THEN MWA(KPT) = MWA(KPT) - JBASE MWA(KPT-1) = MWA(KPT-1) + 1 GO TO 170 ENDIF GO TO 190 ENDIF C IF (K.EQ.0) THEN KP = N1 + 3 KPM1 = KP - 1 DO 180 J = 3, N1 IF (MWA(KP-J).GE.JBASE) THEN MWA(KP-J) = MWA(KP-J) - JBASE MWA(KPM1-J) = MWA(KPM1-J) + 1 ENDIF 180 CONTINUE ENDIF C C Shift right if the leading digit is not less than JBASE. C 190 IF (MWA(2).GE.JBASE) THEN 200 KP = NMWA + 4 DO 210 J = 4, NMWA MWA(KP-J) = MWA(KP-J-1) 210 CONTINUE KT = MWA(2)/JBASE MWA(3) = MWA(2) - KT*JBASE MWA(2) = KT MWA(1) = MWA(1) + 1 IF (MWA(2).GE.JBASE) GO TO 200 ENDIF C RETURN END SUBROUTINE FMARGS(KROUTN,NARGS,MA,MB,KRESLT) C C Check the input arguments to a routine for special cases. C C KROUTN - ID number of the subroutine which 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 PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK),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 COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX 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.KUNKNO) RETURN IF (NARGS.EQ.2) THEN IF (MB(1).EQ.KUNKNO) RETURN ENDIF KFLAG = 0 C C Check the validity of parameters if this is a user call. C IF (NCALL.GT.1) GO TO 130 C C Check NDIG. C IF (NDIG.LT.2 .OR. NDIG.GT.MXNDIG) THEN KFLAG = -1 CALL FMWARN NDS = NDIG IF (NDIG.LT.2) NDIG = 2 IF (NDIG.GT.MXNDIG) NDIG = MXNDIG WRITE (KW,110) NDS,NDIG 110 FORMAT(' NDIG WAS',I10,'. IT HAS BEEN CHANGED TO',I10,'.') RETURN ENDIF C C Check JBASE. C IF (JBASE.LT.2 .OR. JBASE.GT.MXBASE) THEN KFLAG = -2 CALL FMWARN JBS = JBASE IF (JBASE.LT.2) JBASE = 2 IF (JBASE.GT.MXBASE) JBASE = MXBASE WRITE (KW,120) JBS,JBASE 120 FORMAT(' JBASE WAS',I10,'. IT HAS BEEN CHANGED TO',I10,'.') 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.KEXPOV .OR. ABS(MA(2)).NE.1) THEN CALL FMIM(0,MA) KFLAG = -3 CALL FMWARN MA(1) = KUNKNO MA(2) = 1 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.KEXPOV .OR. ABS(MB(2)).NE.1) THEN CALL FMIM(0,MB) KFLAG = -3 CALL FMWARN MB(1) = KUNKNO MB(2) = 1 RETURN ENDIF ENDIF ENDIF C C Check for special cases. C 130 CALL FMCAT(MA,NCATMA) NCATMB = 0 IF (NARGS.EQ.2) CALL FMCAT(MB,NCATMB) C IF (KROUTN.EQ.3) THEN KRESLT = KADD(NCATMB,NCATMA) GO TO 140 ENDIF C IF (KROUTN.EQ.29) THEN KRESLT = KMPY(NCATMB,NCATMA) GO TO 140 ENDIF C IF (KROUTN.EQ.12) THEN KRESLT = KDIV(NCATMB,NCATMA) GO TO 140 ENDIF C IF (KROUTN.EQ.34) THEN KRESLT = KPWR(NCATMB,NCATMA) GO TO 140 ENDIF C IF (KROUTN.EQ.39) THEN KRESLT = KSQRT(NCATMA) GO TO 140 ENDIF C IF (KROUTN.EQ.15) THEN KRESLT = KEXP(NCATMA) GO TO 140 ENDIF C IF (KROUTN.EQ.21) THEN KRESLT = KLN(NCATMA) GO TO 140 ENDIF C IF (KROUTN.EQ.36) THEN KRESLT = KSIN(NCATMA) GO TO 140 ENDIF C IF (KROUTN.EQ.9) THEN KRESLT = KCOS(NCATMA) GO TO 140 ENDIF C IF (KROUTN.EQ.41) THEN KRESLT = KTAN(NCATMA) GO TO 140 ENDIF C IF (KROUTN.EQ.4) THEN KRESLT = KASIN(NCATMA) IF ((NCATMA.EQ.7.OR.NCATMA.EQ.9) .AND. KRAD.EQ.0) KRESLT = 12 GO TO 140 ENDIF C IF (KROUTN.EQ.2) THEN KRESLT = KACOS(NCATMA) GO TO 140 ENDIF C IF (KROUTN.EQ.5) THEN KRESLT = KATAN(NCATMA) IF ((NCATMA.EQ.7.OR.NCATMA.EQ.9) .AND. KRAD.EQ.0) KRESLT = 12 GO TO 140 ENDIF C IF (KROUTN.EQ.37) THEN KRESLT = KSINH(NCATMA) GO TO 140 ENDIF C IF (KROUTN.EQ.10) THEN KRESLT = KCOSH(NCATMA) GO TO 140 ENDIF C IF (KROUTN.EQ.42) THEN KRESLT = KTANH(NCATMA) GO TO 140 ENDIF C IF (KROUTN.EQ.20) THEN KRESLT = KLG10(NCATMA) GO TO 140 ENDIF C KRESLT = 0 RETURN C 140 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 PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C C Scratch array usage during FMASIN: M01 - M06 C COMMON /FM1/ M01(LUNPCK),M02(LUNPCK),M03(LUNPCK),M04(LUNPCK), * M05(LUNPCK),M06(LUNPCK) C CALL FMENTR(4,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN) IF (KRESLT.NE.0) RETURN C CALL FMEQU(MA,MB,NDSAVE,NDIG) 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 CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN) RETURN END SUBROUTINE FMATAN(MA,MB) C C MB = ARCTAN(MA) C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DOUBLE PRECISION X,XB,XM DIMENSION MA(LUNPCK),MB(LUNPCK),NSTACK(19) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C C Scratch array usage during FMATAN: M01 - M06 C COMMON /FM1/ M01(LUNPCK),M02(LUNPCK),M03(LUNPCK),M04(LUNPCK), * M05(LUNPCK),M06(LUNPCK) C COMMON /FMSAVE/ NDIGPI,NJBPI,NDIGE,NJBE,NDIGLB,NJBLB,NDIGLI,NJBLI, * MPISAV(LUNPCK),MESAV(LUNPCK),MLBSAV(LUNPCK), * MLN1(LUNPCK),MLN2(LUNPCK),MLN3(LUNPCK), * MLN4(LUNPCK) C CALL FMENTR(5,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN) IF (KRESLT.NE.0) RETURN C CALL FMEQU(MA,M05,NDSAVE,NDIG) 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 KWSV = KWARN C X = M05(1) XB = JBASE 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 (NJBPI.NE.JBASE .OR. NDIGPI.LT.NDIG) THEN NDSV = NDIG NDIG = MIN(NDIG+2,MXNDG2) CALL FMPI2(MPISAV) NJBPI = JBASE NDIGPI = NDIG NDIG = NDSV ENDIF ENDIF C C If the argument is small, use the Taylor series, C otherwise use Newton iteration. C IF (X*LOG(XB).LT.-5.0*LOG(XM)) THEN KWARN = 0 CALL FMEQU(M05,MB,NDIG,NDIG) CALL FMMPY(M05,M05,M06) J = 3 NDSAV1 = NDIG C 110 CALL FMMPY(M05,M06,M05) M05(2) = -M05(2) CALL FMDIVI(M05,J,M03) NDIG = NDSAV1 CALL FMADD(MB,M03,MB) IF (KFLAG.EQ.1) THEN KFLAG = 0 GO TO 130 ENDIF NDIG = NDSAV1 - (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 FMDP2M(X,MB) CALL FMDIG(NSTACK,KST) C C Newton iteration. C DO 120 J = 1, KST NDIG = NSTACK(J) CALL FMSIN(MB,M06) CALL FMMPY(M06,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 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 (MA2.LT.0) MB(2) = -MB(2) C IF (KFLAG.EQ.1) KFLAG = 0 KWARN = KWSV 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 PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C C Scratch array usage during FMATN2: M01 - M06 C COMMON /FM1/ M01(LUNPCK),M02(LUNPCK),M03(LUNPCK),M04(LUNPCK), * M05(LUNPCK),M06(LUNPCK) C CALL FMENTR(6,MA,MB,2,MC,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN) IF (KRESLT.NE.0) RETURN C KARGSW = 0 KWSV = KWARN KWARN = 0 CALL FMEQU(MA,M01,NDSAVE,NDIG) CALL FMEQU(MB,M02,NDSAVE,NDIG) C C Check for special cases. C IF (MA(1).EQ.KUNKNO .OR. MB(1).EQ.KUNKNO .OR. * (MA(2).EQ.0 .AND. MB(2).EQ.0)) THEN CALL FMIM(0,MC) MC(1) = KUNKNO MC(2) = 1 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 IF (MA(1).EQ.KEXPOV .AND. MB(1).LT.MXSAVE-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.KEXPUN .AND. (-MB(1)).LT.MXSAVE-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.KEXPOV .AND. MA(1).LT.MXSAVE-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.KEXPUN .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.KEXPUN .AND. (-MA(1)).LT.MXSAVE-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.KUNKNO) * MC(2) = -MC(2) C C Round the result and return. C 110 IF (KFLAG.EQ.1) KFLAG = 0 KWARN = KWSV 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 PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C NCALL = NCALL + 1 KSTACK(NCALL) = 7 C KFLAG = 0 N1 = NDIG + 1 DO 110 J = 2, N1 MA(J) = JBASE - 1 110 CONTINUE MA(1) = MXEXP + 1 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 PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C C Check for special symbols. C NCAT = 16 IF (MA(1).EQ.KUNKNO) RETURN C IF (MA(1).EQ.KEXPOV) THEN NCAT = 15 IF (MA(2).LT.0) NCAT = 1 RETURN ENDIF C IF (MA(1).EQ.KEXPUN) 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 IF (MA(1).GE.MXEXP-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.-MXEXP+NDIG) THEN NCAT = 11 IF (MA(2).LT.0) NCAT = 5 RETURN ENDIF C IF (MA(1).GE.-MXEXP) THEN NCAT = 10 IF (MA(2).LT.0) NCAT = 6 RETURN ENDIF C RETURN END FUNCTION FMCOMP(MA,LREL,MB) LOGICAL FMCOMP CHARACTER *2 LREL 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 PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C NCALL = NCALL + 1 KSTACK(NCALL) = 8 C IF (NCALL.LE.LVLTRC .AND. ABS(NTRACE).GE.2) THEN WRITE (KW,110) 110 FORMAT(' INPUT TO FMCOMP') C IF (NTRACE.GT.0) THEN CALL FMPRNT(MA) WRITE (KW,120) LREL 120 FORMAT(7X,'.',A2,'.') CALL FMPRNT(MB) ELSE CALL FMNTRJ(MA,NDIG) WRITE (KW,120) LREL CALL FMNTRJ(MB,NDIG) ENDIF ENDIF C C JCOMP will be 1 if MA.GT.MB C 2 if MA.EQ.MB C 3 if MA.LT.MB C C Check for special cases. C IF (LREL.NE.'EQ' .AND. LREL.NE.'NE' .AND. LREL.NE.'LT' .AND. * LREL.NE.'GT' .AND. LREL.NE.'LE' .AND. LREL.NE.'GE') THEN KFLAG = -4 FMCOMP = .FALSE. IF (KWARN.LE.0) GO TO 170 WRITE (KW,130) LREL 130 FORMAT(/' ERROR OF TYPE KFLAG = -4 IN FM PACKAGE IN ROUTINE', * ' FMCOMP'//1X,A,' IS NOT ONE OF THE SIX RECOGNIZED', * ' COMPARISONS.'//' .FALSE. HAS BEEN RETURNED.'/) GO TO 170 ENDIF C IF (MA(1).EQ.KUNKNO .OR. MB(1).EQ.KUNKNO) THEN KFLAG = -4 FMCOMP = .FALSE. GO TO 170 ENDIF C IF (ABS(MA(1)).EQ.KEXPOV .AND. MA(1).EQ.MB(1) .AND. * MA(2).EQ.MB(2)) THEN KFLAG = -4 FMCOMP = .FALSE. IF (KWARN.LE.0) GO TO 170 WRITE (KW,140) 140 FORMAT(/' ERROR OF TYPE KFLAG = -4 IN FM PACKAGE IN ROUTINE', * ' FMCOMP'//' TWO NUMBERS IN THE SAME OVERFLOW OR', * ' UNDERLOW CATEGORY CANNOT BE COMPARED.'// * ' .FALSE. HAS BEEN RETURNED.'/) GO TO 170 ENDIF C C Check for zero. C KFLAG = 0 IF (MA(2).EQ.0) THEN JCOMP = 2 IF (MB(2).LT.0) JCOMP = 1 IF (MB(2).GT.0) JCOMP = 3 GO TO 160 ENDIF IF (MB(2).EQ.0) THEN JCOMP = 1 IF (MA(2).LT.0) JCOMP = 3 GO TO 160 ENDIF C Check for opposite signs. C IF (MA(2).GT.0 .AND. MB(2).LT.0) THEN JCOMP = 1 GO TO 160 ENDIF IF (MB(2).GT.0 .AND. MA(2).LT.0) THEN JCOMP = 3 GO TO 160 ENDIF C C See which one is larger in absolute value. C IF (MA(1).GT.MB(1)) THEN JCOMP = 1 GO TO 160 ENDIF IF (MB(1).GT.MA(1)) THEN JCOMP = 3 GO TO 160 ENDIF NLAST = NDIG + 1 C DO 150 J = 2, NLAST IF (ABS(MA(J)).GT.ABS(MB(J))) THEN JCOMP = 1 GO TO 160 ENDIF IF (ABS(MB(J)).GT.ABS(MA(J))) THEN JCOMP = 3 GO TO 160 ENDIF 150 CONTINUE C JCOMP = 2 C C Now match the JCOMP value to the requested comparison. C 160 IF (JCOMP.EQ.1 .AND. MA(2).LT.0) THEN JCOMP = 3 ELSE IF (JCOMP.EQ.3 .AND. MB(2).LT.0) THEN JCOMP = 1 ENDIF C FMCOMP = .FALSE. IF (JCOMP.EQ.1 .AND. (LREL.EQ.'GT' .OR. LREL.EQ.'GE' .OR. * LREL.EQ.'NE')) FMCOMP = .TRUE. C IF (JCOMP.EQ.2 .AND. (LREL.EQ.'EQ' .OR. LREL.EQ.'GE' .OR. * LREL.EQ.'LE')) FMCOMP = .TRUE. C IF (JCOMP.EQ.3 .AND. (LREL.EQ.'NE' .OR. LREL.EQ.'LT' .OR. * LREL.EQ.'LE')) FMCOMP = .TRUE. C 170 IF (NTRACE.NE.0) THEN IF (NCALL.LE.LVLTRC .AND. ABS(NTRACE).GE.1) THEN IF (KFLAG.EQ.0) THEN WRITE (KW,180) NCALL,JBASE,NDIG 180 FORMAT(' FMCOMP',15X,'CALL LEVEL =',I2,5X,'JBASE =', * I10,5X,'NDIG =',I6) ELSE WRITE (KW,190) NCALL,JBASE,NDIG,KFLAG 190 FORMAT(' FMCOMP',6X,'CALL LEVEL =',I2,4X,'JBASE =', * I10,4X,'NDIG =',I6,4X,'KFLAG =',I3) ENDIF IF (FMCOMP) THEN WRITE (KW,200) 200 FORMAT(7X,'.TRUE.') ELSE WRITE (KW,210) 210 FORMAT(7X,'.FALSE.') ENDIF ENDIF ENDIF NCALL = NCALL - 1 RETURN END SUBROUTINE FMCOS(MA,MB) C C MB = COS(MA) C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C COMMON /FMSAVE/ NDIGPI,NJBPI,NDIGE,NJBE,NDIGLB,NJBLB,NDIGLI,NJBLI, * MPISAV(LUNPCK),MESAV(LUNPCK),MLBSAV(LUNPCK), * MLN1(LUNPCK),MLN2(LUNPCK),MLN3(LUNPCK), * MLN4(LUNPCK) C C Scratch array usage during FMCOS: M01 - M04 C COMMON /FM1/ M01(LUNPCK),M02(LUNPCK),M03(LUNPCK),M04(LUNPCK), * M05(LUNPCK),M06(LUNPCK) C CALL FMENTR(9,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN) IF (KRESLT.NE.0) RETURN C CALL FMEQU(MA,MB,NDSAVE,NDIG) C C Reduce the argument, convert to radians if the input is C in degrees, and evaluate the function. C CALL FMRDC(MB,MB,JSIN,JCOS,JSWAP) IF (MB(1).EQ.KUNKNO) GO TO 110 IF (KRAD.EQ.0) THEN IF (NJBPI.NE.JBASE .OR. NDIGPI.LT.NDIG) THEN NDSV = NDIG NDIG = MIN(NDIG+2,MXNDG2) CALL FMPI2(MPISAV) NJBPI = JBASE NDIGPI = NDIG NDIG = NDSV ENDIF CALL FMMPY(MB,MPISAV,MB) CALL FMDIVI(MB,180,MB) ENDIF IF (MB(1).NE.KUNKNO) THEN IF (JSWAP.EQ.0) THEN CALL FMCOS2(MB,MB) ELSE CALL FMSIN2(MB,MB) ENDIF ENDIF C C Append the sign, round, and return. C IF (JCOS.EQ.-1) MB(2) = -MB(2) 110 CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN) RETURN END SUBROUTINE FMCOS2(MA,MB) C C Internal routine for computing MB = COS(MA) where 0.LE.MA.LE.1. C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C C Scratch array usage during FMCOS2: M01 - M03 C COMMON /FM1/ M01(LUNPCK),M02(LUNPCK),M03(LUNPCK),M04(LUNPCK), * M05(LUNPCK),M06(LUNPCK) C C Use COS(X) = SQRT(1-SIN(X)**2). C CALL FMSIN2(MA,MB) CALL FMMPY(MB,MB,MB) CALL FMI2M(1,M02) CALL FMSUB(M02,MB,M02) CALL FMSQRT(M02,MB) RETURN END SUBROUTINE FMCOSH(MA,MB) C C MB = COSH(MA) C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C C Scratch array usage during FMCOSH: M01 - M03 C COMMON /FM1/ M01(LUNPCK),M02(LUNPCK),M03(LUNPCK),M04(LUNPCK), * M05(LUNPCK),M06(LUNPCK) C CALL FMENTR(10,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN) IF (KRESLT.NE.0) RETURN C KWRNSV = KWARN KWARN = 0 KARGSW = 0 CALL FMEQU(MA,MB,NDSAVE,NDIG) CALL FMEXP(MB,MB) IF (MB(1).EQ.KEXPOV) THEN GO TO 110 ELSE IF (MB(1).EQ.KEXPUN) THEN MB(1) = KEXPOV GO TO 110 ENDIF CALL FMI2M(1,M01) CALL FMDIV(M01,MB,M01) CALL FMADD(MB,M01,MB) CALL FMDIVI(MB,2,MB) C C Round and return. C 110 KWARN = KWRNSV CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN) RETURN END SUBROUTINE FMDIG(NSTACK,KST) C C Compute the number of intermediate digits to be used in Newton C iteration. This assumes that a starting approximation which is C accurate to double precision is used, and the root is simple. C C KST is the number of iterations needed for final accuracy NDIG. C NSTACK(J) holds the value of NDIG to be used for the C Jth iteration. C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DOUBLE PRECISION XBASE,ONE,Y,YT DIMENSION NSTACK(19) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C C Compute the approximate machine precision. This will be C the initial accuracy obtained by using double precision C to compute the first approximation before starting Newton C iteration. C XBASE = JBASE ONE = 1.0 Y = ONE NE = 1 C 110 Y = Y/XBASE YT = ONE + Y IF (YT.GT.ONE) THEN NE = NE + 1 GO TO 110 ENDIF C C Fill the intermediate digit stack (backwards). C KST = 1 ND = NDIG NSTACK(1) = ND IF (ND.LT.2*NE .OR. ND.LE.2) RETURN C 120 Y = ND C C The 1.9 accounts for the fact that the number of correct C digits approximately doubles at each iteration. C NDT = Y/1.9 IF (2*NDT.LE.ND) NDT = NDT + 1 ND = NDT KST = KST + 1 NSTACK(KST) = ND IF (ND.GE.2*NE .AND. ND.GT.2) GO TO 120 C C Reverse the stack. C L = KST/2 DO 130 J = 1, L JT = NSTACK(J) NSTACK(J) = NSTACK(KST+1-J) NSTACK(KST+1-J) = JT 130 CONTINUE C RETURN END SUBROUTINE FMDIM(MA,MB,MC) C C MC = DIM(MA,MB) C C Positive difference. MC = MA - MB if MA.GE.MB, C = 0 otherwise. C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C LOGICAL FMCOMP DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C C Scratch array usage during FMDIM: M01 - M02 C COMMON /FM1/ M01(LUNPCK),M02(LUNPCK),M03(LUNPCK),M04(LUNPCK), * M05(LUNPCK),M06(LUNPCK) C CALL FMENTR(11,MA,MB,2,MC,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN) IF (KRESLT.NE.0) RETURN KARGSW = 0 KWSV = KWARN KWARN = 0 MXEXP = MXSAVE C CALL FMEQU(MA,M01,NDSAVE,NDIG) CALL FMEQU(MB,M02,NDSAVE,NDIG) C IF (FMCOMP(M01,'LT',M02)) THEN CALL FMI2M(0,MC) ELSE CALL FMSUB(M01,M02,MC) ENDIF C IF (KFLAG.EQ.1) KFLAG = 0 KWARN = KWSV CALL FMEXIT(MC,MC,NDSAVE,MXSAVE,KASAVE,KOVUN) RETURN END SUBROUTINE FMDIV(MA,MB,MC) C C MC = MA / MB C C This routine performs the trace printing for division. C FMDIV2 is used to do the arithmetic. C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C NCALL = NCALL + 1 KSTACK(NCALL) = 12 IF (NTRACE.NE.0) CALL FMNTR(2,MA,MB,2) C CALL FMDIV2(MA,MB,MC) C IF (NTRACE.NE.0) CALL FMNTR(1,MC,MC,1) NCALL = NCALL - 1 RETURN END SUBROUTINE FMDIV2(MA,MB,MC) C C Internal division routine. MC = MA / MB C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C IF (KARGSW.NE.1) THEN CALL FMARGS(12,2,MA,MB,KRESLT) IF (KRESLT.NE.0) THEN CALL FMRSLT(MA,MB,MC,KRESLT) RETURN ENDIF ELSE IF (MB(2).EQ.0) THEN CALL FMIM(0,MC) MC(1) = KUNKNO MC(2) = 1 KFLAG = -4 CALL FMWARN RETURN ENDIF IF (MA(2).EQ.0) THEN CALL FMIM(0,MC) RETURN ENDIF ENDIF KFLAG = 0 C C NGUARD is the number of guard digits used. C IF (NCALL.GT.1) THEN NGUARD = 1 ELSE B = JBASE NGUARD = 5.0*LOG(10.0)/LOG(B) + 2.0 ENDIF N1 = NDIG + 1 NG = N1 + NGUARD C C Copy MA into the working array. C DO 110 J = 3, N1 MWA(J+1) = MA(J) 110 CONTINUE MWA(1) = MA(1) - MB(1) + 1 MWA(2) = 0 N3 = NDIG + 3 NL = N1 + NGUARD + 2 DO 120 J = N3, NL MWA(J) = 0 120 CONTINUE C C Save the sign of MA and MB and then work only with C positive numbers. C MA2 = MA(2) MB1 = MB(1) MB2 = MB(2) MA(2) = ABS(MA(2)) MWA(3) = MA(2) MB(1) = 0 MB(2) = ABS(MB(2)) C C NMBWDS is the number of words of MB used to compute C the trial quotient digit LT. C NMBWDS = 2 IF (JBASE.LT.100) NMBWDS = 4 C C XB is an approximation of MB used in selecting C trial quotients. C XBASE = JBASE XB = 0 JL = NMBWDS + 1 DO 130 J = 2, JL IF (J.LE.N1) THEN XB = XB*XBASE + MB(J) ELSE XB = XB*XBASE ENDIF 130 CONTINUE IF (JL+1.LE.N1) XB = XB + MB(JL+1)/XBASE XBLOG = LOG(XBASE) C C LMAX determines when to normalize all of MWA. C NGAP determines how many digits must be normalized for C each trial quotient. C IF (MXBASE**2/(JBASE-1)**2.GE.NL) THEN LMAX = NL*(JBASE-1) ELSE LMAX = MXBASE**2/(JBASE-1) ENDIF X = LMAX NGAP = LOG(X)/XBLOG + 2.0 C C Count the trailing zero digits of MB. C DO 140 J = 1, NDIG IF (MB(NDIG+2-J).NE.0) THEN NZDMB = J - 1 GO TO 150 ENDIF 140 CONTINUE C C MAXMWA is an upper bound on the size of values in MWA C divided by JBASE-1. It is used to determine whether C normalization can be postponed. C 150 MAXMWA = 0 C C KPTMWA points to the next digit in the quotient. C KPTMWA = 2 C C This is the start of the division loop. C C LT is the trial quotient digit. C Divide the first (NMBWDS+1) digits of MWA by MB so that C LT has a high probability of being the correct value. C 160 KL = KPTMWA + NMBWDS XMWA = 0 DO 170 J = KPTMWA, KL IF (J.LE.NL) THEN XMWA = XMWA*XBASE + MWA(J) ELSE XMWA = XMWA*XBASE ENDIF 170 CONTINUE C LT = XMWA/XB IF (LT.GE.JBASE) LT = JBASE - 1 C C Subtract LT*MB from MWA. C KA = KPTMWA + 1 KB = MIN(KA+NDIG-1-NZDMB,NL) JB = KA - 2 IF (LT.GT.0) THEN DO 180 J = KA, KB MWA(J) = MWA(J) - LT*MB(J-JB) 180 CONTINUE ENDIF C C Normalize enough digits in MWA so that LQ, the correct C quotient digit, can be determined. C LQ = LT MAXMWA = MAXMWA + LT C C NORMPT points to the last digit currently normalized. C NORMPT = MIN(KA+NGAP+1,KB) C C If the correct LQ has been found, then the current section C of MWA will be nonnegative and less than MB. C If MWA.LT.0 then too large a trial quotient was used and C a '+ MB' correction step is needed. C If MWA.GE.MB then too small a trial quotient was used and C a '- MB' correction step is needed. C 190 J = NORMPT 200 IF (MWA(J).LT.0) THEN KT = (-MWA(J)-1)/JBASE + 1 MWA(J-1) = MWA(J-1) - KT MWA(J) = MWA(J) + KT*JBASE ENDIF J = J - 1 IF (J.GE.KA) GO TO 200 IF (NORMPT.EQ.KB) MAXMWA = 0 C C NONZPT points to the first nonzero digit in MWA after the C subtraction and normalization have been done. C IF (MWA(KA-1).NE.0) THEN NONZPT = KA - 1 ELSE IF (MWA(KA).NE.0) THEN NONZPT = KA ELSE NONZPT = 0 K1 = KA + 1 IF (K1.LE.NORMPT) THEN DO 210 J = K1, NORMPT IF (MWA(J).NE.0) THEN NONZPT = J GO TO 220 ENDIF 210 CONTINUE ENDIF ENDIF C C KDIFF points to the first digit in MWA which differs C from the corresponding digit in MB. C 220 IF (MWA(KA-1).NE.0) THEN KDIFF = KA - 1 ELSE IF (MWA(KA).NE.MB(2)) THEN KDIFF = KA ELSE KDIFF = 0 K1 = KA + 1 IF (K1.LE.NORMPT) THEN DO 230 J = K1, NORMPT IF (MWA(J).NE.MB(J-JB)) THEN KDIFF = J GO TO 240 ENDIF 230 CONTINUE ENDIF ENDIF C C Protect against cases like MWA = 2 0 0 0 ... C and MB = 1 9 9 9 .... Since MWA may not be completely C normalized, MWA(KDIFF) might still be equal to C MB(KDIFF-JB). C 240 IF (KDIFF.NE.0 .AND. KDIFF.LT.NL .AND. KDIFF+1-JB.LE.N1) THEN IF (MWA(KDIFF+1).EQ.0 .AND. NORMPT.LT.KB .AND. * MB(KDIFF+1-JB).EQ.JBASE-1) KDIFF = 0 ENDIF C C If NONZPT is zero then not enough digits were normalized C to tell whether MWA is negative. C IF (NONZPT.EQ.0 .OR. NORMPT-NONZPT.LT.NGAP) THEN IF (NORMPT.EQ.KB) GO TO 250 NORMPT = KB GO TO 190 ENDIF C C See if the current section of MWA is negative. C 250 IF (NONZPT.LE.0) GO TO 270 IF (MWA(NONZPT).LT.0) THEN C C At the end of the division LT could be off by more than C one because part of MB has shifted off the end of the C active part of MWA. Skip the correction step for the C last guard digit. C IF (KA.EQ.KB .AND. KB.EQ.NL) GO TO 310 C C Correct LQ by adding MB. C LQ = LQ - 1 MAXMWA = MAXMWA - 1 C C This is the only place where values greater than JBASE C can be introduced into MWA. Normalize them here. C K = KB KC = 0 260 MWA(K) = MWA(K) + MB(K-JB) + KC IF (MWA(K).GE.JBASE) THEN KC = 1 MWA(K) = MWA(K) - JBASE ELSE KC = 0 ENDIF K = K - 1 IF (K.GE.KA) GO TO 260 IF (KC.EQ.1) MWA(KA-1) = MWA(KA-1) + 1 GO TO 190 ENDIF C C If KDIFF is zero then not enough digits were normalized C to tell if MWA.GE.MB. C 270 IF (KDIFF.EQ.0) THEN IF (NORMPT.EQ.KB) GO TO 280 NORMPT = KB GO TO 190 ENDIF C C See if the current section of MWA is greater than MB. C IF (MWA(KDIFF).LT.MB(KDIFF-JB)) GO TO 310 C C If NORMPT-KDIFF.GE.NGAP then not enough digits were C normalized to tell if MWA.GE.MB. C IF (NORMPT-KDIFF.LT.NGAP .AND. NORMPT.LT.KB) THEN NORMPT = KB GO TO 190 ENDIF C C Correct LQ by subtracting MB. C 280 IF (KA.EQ.KB .AND. KB.EQ.NL) GO TO 310 LQ = LQ + 1 C C Due to part of MB shifting off the end of the active C part of MWA, LQ could be equal to JBASE here. This C means the correct value of LQ is JBASE-1 for all C subsequent digits. C IF (LQ.GE.JBASE) THEN NG1 = NG + 1 DO 290 K = KPTMWA, NG1 MWA(K) = JBASE - 1 290 CONTINUE GO TO 330 ENDIF C MAXMWA = MAXMWA + 1 DO 300 K = KA, KB MWA(K) = MWA(K) - MB(K-JB) 300 CONTINUE GO TO 190 C C Here the correct value for LQ has been found. C 310 MWA(KPTMWA) = LQ C C See if MWA must be normalized. C IF (MAXMWA+JBASE-1.GE.LMAX) THEN J = KB 320 IF (MWA(J).LT.0) THEN KT = (-MWA(J)-1)/JBASE + 1 MWA(J-1) = MWA(J-1) - KT MWA(J) = MWA(J) + KT*JBASE ENDIF J = J - 1 IF (J.GE.KA) GO TO 320 MAXMWA = 0 ENDIF C KPTMWA = KPTMWA + 1 IF (KPTMWA.LE.NG) GO TO 160 IF (MWA(2).EQ.0 .AND. KPTMWA.LE.NG+1) GO TO 160 C C Round, affix the sign, and return. C 330 MA(2) = MA2 MB(1) = MB1 MB(2) = MB2 IF (MWA(2).EQ.0) THEN CALL FMRND(NDIG,NGUARD,1) ELSE CALL FMRND(NDIG,NGUARD,0) ENDIF CALL FMMOVE(MC) IF (MA2*MB2.LT.0) MC(2) = -MC(2) RETURN END SUBROUTINE FMDIVI(MA,INT,MB) C C MB = MA / INT C C Divide FM number MA by one word integer INT. C C This routine is faster than FMDIV when the divisor is less than C MXBASE. When INT is not less than MXBASE, FMDIV2 is used. In C this case if INT is known to be a product of two integers less than C MXBASE it is usually faster to make two calls to FMDIVI with half- C word factors than one call with their product. C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C C Scratch array usage during FMDIVI: M01 C COMMON /FM1/ M01(LUNPCK),M02(LUNPCK),M03(LUNPCK),M04(LUNPCK), * M05(LUNPCK),M06(LUNPCK) C NCALL = NCALL + 1 KSTACK(NCALL) = 13 IF (NTRACE.NE.0) THEN CALL FMNTR(2,MA,MA,1) CALL FMNTRI(2,INT,0) ENDIF KFLAG = 0 N1 = NDIG + 1 C C Check for special cases. C IF (MA(1).EQ.KUNKNO .OR. INT.EQ.0) THEN MA1 = MA(1) CALL FMIM(0,MB) MB(1) = KUNKNO MB(2) = 1 KFLAG = -4 IF (MA1.NE.KUNKNO) CALL FMWARN IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1) NCALL = NCALL - 1 RETURN ENDIF C IF (ABS(INT).EQ.1) THEN DO 110 J = 1, N1 MB(J) = MA(J) 110 CONTINUE MB(2) = MA(2)*INT IF (MA(1).EQ.KEXPOV) KFLAG = -5 IF (MA(1).EQ.KEXPUN) KFLAG = -6 IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1) NCALL = NCALL - 1 RETURN ENDIF C IF (MA(2).EQ.0) THEN DO 120 J = 1, N1 MB(J) = 0 120 CONTINUE IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1) NCALL = NCALL - 1 RETURN ENDIF C IF (MA(1).EQ.KEXPUN) THEN MA2 = MA(2) CALL FMIM(0,MB) MB(1) = KEXPUN MB(2) = 1 IF ((MA2.LT.0 .AND. INT.GT.0) .OR. * (MA2.GT.0 .AND. INT.LT.0)) MB(2) = -1 KFLAG = -6 IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1) NCALL = NCALL - 1 RETURN ENDIF C IF (MA(1).EQ.KEXPOV) THEN CALL FMIM(0,MB) MB(1) = KUNKNO MB(2) = 1 KFLAG = -4 CALL FMWARN IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1) NCALL = NCALL - 1 RETURN ENDIF C C NGUARD is the number of guard digits used. C IF (NCALL.GT.1) THEN NGUARD = 1 ELSE B = JBASE NGUARD = 5.0*LOG(10.0)/LOG(B) + 2.0 ENDIF N1 = NDIG + 1 C C If ABS(INT).GE.MXBASE use FMDIV. C IF (ABS(INT).GT.MXBASE) THEN CALL FMIM(INT,M01) CALL FMDIV2(MA,M01,MB) IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1) NCALL = NCALL - 1 RETURN ENDIF C C Work with positive numbers. C MA2 = MA(2) MA(2) = ABS(MA(2)) INTP = ABS(INT) C C Find the first significant digit of the quotient. C KT = 0 DO 130 J = 2, N1 KPT = J KT = KT*JBASE + MA(J) IF (KT.GE.INTP) GO TO 150 130 CONTINUE C 140 KPT = KPT + 1 KT = KT*JBASE IF (KT.LT.INTP) GO TO 140 C C Do the rest of the division. C 150 KA = KPT + 1 MWA(1) = MA(1) + 2 - KPT MWA(2) = KT/INTP MODINT = MOD(KT,INTP) KPTWA = 2 IF (KA.LE.N1) THEN KL = 3 - KA DO 160 J = KA, N1 KT = MODINT*JBASE + MA(J) MWA(KL+J) = KT/INTP MODINT = MOD(KT,INTP) 160 CONTINUE KPTWA = KL + N1 ENDIF C KA = KPTWA + 1 KB = N1 + NGUARD IF (KA.LE.KB) THEN DO 170 J = KA, KB KT = MODINT*JBASE MWA(J) = KT/INTP MODINT = MOD(KT,INTP) 170 CONTINUE ENDIF C C Round the result, put the sign on MB and return. C MA(2) = MA2 CALL FMRND(NDIG,NGUARD,0) CALL FMMOVE(MB) IF ((MA2.LT.0 .AND. INT.GT.0) .OR. (MA2.GT.0 .AND. INT.LT.0)) * MB(2) = -MB(2) IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1) NCALL = NCALL - 1 RETURN END SUBROUTINE FMDM(X,MA) C C Internal routine for converting double precision to multiple C precision. C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DOUBLE PRECISION ONE,X,XBASE,Y,YT DIMENSION MA(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C KFLAG = 0 N1 = NDIG + 1 ONE = 1.0 XBASE = JBASE Y = ONE K = 0 NE = 2 C C Find the approximate machine precision. C NE-1 is the number of words at the current precision and C base roughly equal to machine precision. C DO 110 J = 2, NDIG Y = Y/XBASE NE = NE + 1 YT = ONE + Y IF (YT.LE.ONE) GO TO 120 110 CONTINUE 120 Y = X IF (X.LT.0.0) Y = -X C IF (X.EQ.0.0) THEN DO 130 J = 1, N1 MA(J) = 0 130 CONTINUE GO TO 240 ENDIF C C Get the exponent. C IF (Y.GT.ONE) THEN IF (Y/XBASE.LT.Y) THEN 140 K = K + 1 Y = Y/XBASE IF (Y.GT.ONE) GO TO 140 IF (Y.LT.ONE) THEN MA(1) = K GO TO 180 ENDIF GO TO 160 ELSE CALL FMIM(0,MA) MA(1) = KUNKNO MA(2) = 1 KFLAG = -4 CALL FMWARN RETURN ENDIF ENDIF C IF (Y.LT.ONE) THEN IF (Y*XBASE.GT.Y) THEN 150 K = K - 1 Y = Y*XBASE IF (Y.LT.ONE) GO TO 150 IF (Y.GT.ONE) THEN K = K + 1 Y = Y/XBASE MA(1) = K GO TO 180 ENDIF ELSE CALL FMIM(0,MA) MA(1) = KUNKNO MA(2) = 1 KFLAG = -4 CALL FMWARN RETURN ENDIF ENDIF C 160 MA(1) = K + 1 MA(2) = 1 DO 170 J = 3, N1 MA(J) = 0 170 CONTINUE GO TO 240 C C Build the rest of the number. C 180 DO 190 J = 2, NE Y = Y*XBASE K = Y YT = K Y = Y - YT MA(J) = K IF (J.GE.N1) GO TO 210 190 CONTINUE K = NE + 1 DO 200 J = K, N1 MA(J) = 0 200 CONTINUE C C Normalize. C 210 IF (ABS(MA(2)).GE.JBASE) THEN K = N1 + 1 DO 220 J = 3, N1 K = K - 1 MA(K) = MA(K-1) 220 CONTINUE MN = MA(2)/JBASE MA(3) = MA(2) - MN*JBASE MA(2) = MN MA(1) = MA(1) + 1 GO TO 240 ENDIF C IF (MA(2).EQ.0) THEN DO 230 J = 2, NDIG MA(J) = MA(J+1) 230 CONTINUE MA(1) = MA(1) - 1 MA(N1) = 0 ENDIF C 240 IF (X.LT.0.0) MA(2) = -MA(2) RETURN END SUBROUTINE FMDP2M(X,MA) C C MA = X C C Convert a double precision floating point number to FM format. C C In general the relative accuracy of the number returned is only C the relative accuracy of a machine precision number. This may be C true even if X can be represented exactly in the machine floating C point number system. C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DOUBLE PRECISION X DIMENSION MA(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C NCALL = NCALL + 1 KSTACK(NCALL) = 14 IF (NTRACE.NE.0) CALL FMNTRR(2,X,1) C CALL FMDM(X,MA) C IF (NTRACE.NE.0) CALL FMNTR(1,MA,MA,1) NCALL = NCALL - 1 RETURN END SUBROUTINE FMENTR(NROUTN,MA,MB,NARGS,MC,KRESLT,NDSAVE,MXSAVE, * KASAVE,KOVUN) C C Do the argument checking and increasing of precision, overflow C threshold, etc., upon entry to an FM routine. C C NROUTN - routine number of calling routine C MA - first input argument C MB - second input argument (optional) C NARGS - number of input arguments C MC - result argument C KRESLT - returned nonzero if the input arguments give the result C immediately (e.g., MA*0 or OVERFLOW*MB) C NDSAVE - saves the value of NDIG after NDIG is increased C MXSAVE - saves the value of MXEXP C KASAVE - saves the value of KARGSW C KOVUN - returned nonzero if an input argument is (+ or -) overflow C or underflow. C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C NCALL = NCALL + 1 KSTACK(NCALL) = NROUTN IF (NTRACE.NE.0) CALL FMNTR(2,MA,MB,NARGS) CALL FMARGS(NROUTN,NARGS,MA,MB,KRESLT) C KOVUN = 0 IF (MA(1).EQ.KEXPOV .OR. MA(1).EQ.KEXPUN) KOVUN = 1 IF (NARGS.EQ.2) THEN IF (MB(1).EQ.KEXPOV .OR. MB(1).EQ.KEXPUN) KOVUN = 1 ENDIF C C Increase the working precision. C NDSAVE = NDIG IF (NCALL.EQ.1) THEN B = JBASE K = 5.0*LOG(10.0)/LOG(B) + 2.0 NDIG = MAX(NDIG+K,2) IF (NDIG.GT.MXNDG2) THEN KFLAG = -9 CALL FMWARN KRESLT = 12 ENDIF ENDIF C IF (KRESLT.NE.0) THEN IF (KRESLT.EQ.9 .OR. KRESLT.EQ.10 .OR. KRESLT.GE.13) THEN IF (KRAD.EQ.1) THEN CALL FMPI(MC) ELSE CALL FMI2M(180,MC) ENDIF IF (KRESLT.LE.10) CALL FMDIVI(MC,2,MC) IF (KRESLT.GE.14) CALL FMDIVI(MC,4,MC) CALL FMEQU(MC,MC,NDIG,NDSAVE) NDIG = NDSAVE IF (KRESLT.EQ.9 .OR. KRESLT.EQ.14) MC(2) = -MC(2) IF (NTRACE.NE.0) CALL FMNTR(1,MC,MC,1) NCALL = NCALL - 1 RETURN ENDIF C NDIG = NDSAVE CALL FMRSLT(MA,MB,MC,KRESLT) IF (NTRACE.NE.0) CALL FMNTR(1,MC,MC,1) NCALL = NCALL - 1 RETURN ENDIF C C Turn off argument checking for low-level routines while C a high-level call is in progress. C KASAVE = KARGSW KARGSW = 1 C C Extend the overflow/underflow threshold. C MXSAVE = MXEXP MXEXP = MXEXP2 RETURN END SUBROUTINE FMEQU(MA,MB,NDA,NDB) C C Set MB (having NDB digits) equal to MA (having NDA digits). C C If MB has less precision than MA the result is rounded to NDB digits. C C If MB has more precision the result has zero digits padded on the C right. C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C C See if MA and MB are the same array. If so, the loop C which copies MA to MB can be skipped. C C The MIN calls try to obscure the test so that an C optimizing compiler will not remove the K = ... and C MA(3) = ... lines and sometimes change an input argument. C KSAME = 0 K = MIN(MA(3),MXBASE) MB(3) = K + 1 IF (MA(3).EQ.MB(3)) KSAME = 1 MA(3) = MIN(K,MXBASE+MXEXP) C C Check for special symbols. C KFLAG = 0 IF (ABS(MA(1)).GE.KEXPOV) THEN DO 110 J = 2, NDB MB(J+1) = 0 110 CONTINUE MB(1) = MA(1) MB(2) = MA(2) GO TO 240 ENDIF C IF (NDB.EQ.NDA) GO TO 190 C IF (NDB.GT.NDA) GO TO 210 C C Round to NDB digits. C NDG = NDB N1 = NDB + 1 IF (KSAME.EQ.0) THEN DO 120 J = 1, N1 MB(J) = MA(J) 120 CONTINUE ENDIF IF (NDG.LT.1 .OR. (KROUND.EQ.0 .AND. NCALL.LE.1)) GO TO 240 C L = NDB + 2 IF (2*(MA(L)+1).LT.JBASE) GO TO 240 IF (MOD(JBASE,2).EQ.0) THEN IF (2*MA(L).LT.JBASE) GO TO 240 IF (2*MA(L).EQ.JBASE) THEN IF (L.LE.NDA) THEN DO 130 J = L, NDA IF (MA(J+1).GT.0) GO TO 150 130 CONTINUE ENDIF C C Round to even. C IF (MOD(MB(N1),2).EQ.0) GO TO 240 ENDIF ELSE IF (2*MA(L)+1.EQ.JBASE) THEN IF (L.LE.NDA) THEN DO 140 J = L, NDA IF (2*(MA(J+1)+1).LT.JBASE) GO TO 240 IF (2*MA(J+1).GT.JBASE) GO TO 150 140 CONTINUE GO TO 240 ENDIF ENDIF ENDIF C 150 MB(NDG+1) = MB(NDG+1) + 1 MB(NDG+2) = 0 C C Check whether there was a carry in the rounded digit. C MB2 = MB(2) MB(2) = ABS(MB(2)) KB = NDG + 1 IF (KB.GE.3) THEN K = KB + 1 DO 160 J = 3, KB K = K - 1 IF (MB(K).LT.JBASE) GO TO 180 MB(K-1) = MB(K-1) + MB(K)/JBASE MB(K) = MOD(MB(K),JBASE) 160 CONTINUE ENDIF C C If there is a carry in the first digit then the exponent C must be adjusted and the number shifted right. C IF (MB(2).LT.JBASE) GO TO 180 IF (KB.GE.4) THEN K = KB + 1 DO 170 J = 4, KB K = K - 1 MB(K) = MB(K-1) 170 CONTINUE ENDIF C IF (KB.GE.3) MB(3) = MOD(MB(2),JBASE) MB(2) = MB(2)/JBASE MB(1) = MB(1) + 1 C 180 IF (MB2.LT.0) MB(2) = -MB(2) GO TO 240 C C MA and MB have the same precision. C 190 IF (KSAME.EQ.0) THEN NDA1 = NDA + 1 DO 200 J = 1, NDA1 MB(J) = MA(J) 200 CONTINUE ENDIF GO TO 240 C C Extend to NDB digits by padding with zeros. C 210 IF (KSAME.EQ.0) THEN NDA1 = NDA + 1 DO 220 J = 1, NDA1 MB(J) = MA(J) 220 CONTINUE ENDIF NDA2 = NDA + 2 NDB2 = NDB + 1 DO 230 J = NDA2, NDB2 MB(J) = 0 230 CONTINUE C C Check for overflow or underflow. C 240 IF (ABS(MB(1)).GT.MXEXP) THEN IF (MB(1).NE.KUNKNO .OR. MB(2).NE.1) THEN KWRNSV = KWARN KWARN = 0 NCALL = NCALL + 1 CALL FMTRAP(MB) NCALL = NCALL - 1 KWARN = KWRNSV ENDIF IF (MB(1).EQ.KUNKNO) KFLAG = -4 ENDIF C RETURN END SUBROUTINE FMEXIT(MT,MC,NDSAVE,MXSAVE,KASAVE,KOVUN) C C Upon exit from an FM routine the result MT (having precision NDIG) C is rounded and returned in MC (having precision NDSAVE). C The values of NDIG, MXEXP, and KARGSW are restored. C KOVUN is nonzero if one of the routine's input arguments was overflow C or underflow. C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MT(LUNPCK),MC(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C KWSV = KWARN KWARN = 0 MXEXP = MXSAVE KFSAVE = KFLAG CALL FMEQU(MT,MC,NDIG,NDSAVE) IF (KFLAG.NE.-5 .AND. KFLAG.NE.-6) KFLAG = KFSAVE NDIG = NDSAVE KWARN = KWSV IF (KFLAG.EQ.1) KFLAG = 0 IF ((MC(1).EQ.KUNKNO .AND. KFLAG.NE.-9) * .OR. (MC(1).EQ.KEXPUN .AND. KOVUN.EQ.0) * .OR. (MC(1).EQ.KEXPOV .AND. KOVUN.EQ.0)) CALL FMWARN IF (NTRACE.NE.0) CALL FMNTR(1,MC,MC,1) NCALL = NCALL - 1 KARGSW = KASAVE RETURN END SUBROUTINE FMEXP(MA,MB) C C MB = EXP(MA) C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C COMMON /FMSAVE/ NDIGPI,NJBPI,NDIGE,NJBE,NDIGLB,NJBLB,NDIGLI,NJBLI, * MPISAV(LUNPCK),MESAV(LUNPCK),MLBSAV(LUNPCK), * MLN1(LUNPCK),MLN2(LUNPCK),MLN3(LUNPCK), * MLN4(LUNPCK) C C Scratch array usage during FMEXP: M01 - M03 C COMMON /FM1/ M01(LUNPCK),M02(LUNPCK),M03(LUNPCK),M04(LUNPCK), * M05(LUNPCK),M06(LUNPCK) C CALL FMENTR(15,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN) IF (KRESLT.NE.0) RETURN C MA1 = MA(1) MA2 = MA(2) CALL FMEQU(MA,MB,NDSAVE,NDIG) C C Check for obvious underflow or overflow. C XOV is LN(LN(slightly above overflow)) C XMA is LN(LN(EXP(MA))) approximately. C XOV = LOG(1.01*REAL(MXEXP)) + LOG(LOG(REAL(JBASE))) XMA = LOG(REAL(MAX(ABS(MA2),1))) - LOG(REAL(JBASE)) + * MA1*LOG(REAL(JBASE)) C 110 IF (XMA.GE.XOV) THEN CALL FMIM(0,MB) IF (MA2.GT.0) THEN KFLAG = -5 MB(1) = KEXPOV MB(2) = 1 ELSE KFLAG = -6 MB(1) = KEXPUN MB(2) = 1 ENDIF NDIG = NDSAVE MXEXP = MXSAVE KARGSW = KASAVE CALL FMWARN IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1) NCALL = NCALL - 1 RETURN ENDIF C C Split MA into integer and fraction parts. C Work with a positive argument. C M02 = integer part of ABS(MA) C MB = fraction part of ABS(MA) C MB(2) = ABS(MB(2)) CALL FMINT(MB,M02) CALL FMSUB(MB,M02,MB) C C If the integer part is not zero, use FMIPWR to compute C E**(M02). If M02 is too large to represent as a one word C integer, the definition of MXEXP insures that E**(M02) C overflows or underflows. C KWSV = KWARN KWARN = 0 CALL FMM2I(M02,KT) KWARN = KWSV IF (KFLAG.NE.0) THEN XMA = XOV GO TO 110 ENDIF IF (KT.GT.0) THEN C C Compute IEXTRA, the number of extra digits required in C order to get EXP(KT) correct to the current precision. C IEXTRA = (LN(KT) - 1)/LN(JBASE) C XT = KT XBASE = JBASE IEXTRA = LOG(XT)/LOG(XBASE) + 0.5 IF (IEXTRA.GT.0 .AND. NDIG+IEXTRA.LE.MXNDG2) THEN DO 120 J = 1, IEXTRA MB(NDIG+1+J) = 0 120 CONTINUE ENDIF NDIG = NDIG + IEXTRA IF (NDIG.GT.MXNDG2) THEN KFLAG = -9 CALL FMWARN MB(1) = KUNKNO MB(2) = 1 DO 130 J = 2, NDSAVE MB(J+1) = 0 130 CONTINUE CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN) RETURN ENDIF C C Check whether the current precision of e is large C enough. C IF (NJBE.NE.JBASE .OR. NDIG.GT.NDIGE) THEN NDSV = NDIG NDIG = MIN(NDIG+2,MXNDG2) CALL FMI2M(1,MESAV) CALL FMEXP2(MESAV,MESAV) NJBE = JBASE NDIGE = NDIG NDIG = NDSV ENDIF C ENDIF C C Now do the fraction part of MA and combine the results. C KWSV = KWARN KWARN = 0 IF (MB(2).NE.0 .AND. KT.GT.0) THEN CALL FMEXP2(MB,MB) CALL FMIPWR(MESAV,KT,M03) CALL FMMPY(MB,M03,MB) ELSE IF (MB(2).NE.0 .AND. KT.EQ.0) THEN CALL FMEXP2(MB,MB) ELSE IF (MB(2).EQ.0 .AND. KT.GT.0) THEN CALL FMIPWR(MESAV,KT,MB) ELSE CALL FMI2M(1,MB) ENDIF C C Invert if MA was negative. C IF (MA2.LT.0) THEN CALL FMI2M(1,M02) CALL FMDIV(M02,MB,MB) ENDIF KWARN = KWSV C C Round the result and return. C CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN) RETURN END SUBROUTINE FMEXP2(MA,MB) C C MB = EXP(MA) C C Internal exponential routine (called with 0.LT.MA.LE.1). C C Upon return M03 (in COMMON /FM1/) contains an accurate value for C EXP(MA)-1 which can be used to avoid loss of significance in other C routines if MA is close to zero. C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK),MB(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C C Scratch array usage during FMEXP2: M01 - M03 C COMMON /FM1/ M01(LUNPCK),M02(LUNPCK),M03(LUNPCK),M04(LUNPCK), * M05(LUNPCK),M06(LUNPCK) C COMMON /FMSUMS/ MJSUMS(LJSUMS) C C LJSUMS = 8*LUNPCK allows for up to eight concurrent sums. C Increasing this value will begin to improve the speed of C EXP when the base is large and precision exceeds about C 1,500 decimal digits. C NDSAVE = NDIG IF (MA(1).EQ.1) THEN C C Here the special case EXP(1.0) is computed. C Use the direct series e = 1/0! + 1/1! + 1/2! + ... C Do as much of the work as possible using small integers C to minimize the number of FM calls. C Reduce NDIG while computing each term in the C sum as the terms get smaller. C B = JBASE T = NDIG XN = T*LOG(B)/LOG(T) K = LOG(XN)/LOG(B) NDIG = MAX(NDIG+K,2) IF (NDIG.GT.MXNDG2) THEN KFLAG = -9 CALL FMWARN MB(1) = KUNKNO MB(2) = 1 DO 110 J = 2, NDSAVE MB(J+1) = 0 110 CONTINUE NDIG = NDSAVE RETURN ENDIF NDSAV1 = NDIG C CALL FMI2M(2,MB) CALL FMI2M(1,M02) J = 2 NBIG = MXBASE C 120 NTOP = 1 NBOT = J 130 IF (NBOT.GT.NBIG/(J+1)) GO TO 140 J = J + 1 NTOP = J*NTOP + 1 NBOT = J*NBOT GO TO 130 C 140 CALL FMDIVI(M02,NBOT,M02) IF (NTOP.GT.1) THEN CALL FMMPYI(M02,NTOP,M03) NDIG = NDSAV1 CALL FMADD(MB,M03,MB) NDIG = NDSAV1 - (MB(1)-M03(1)) ELSE NDIG = NDSAV1 CALL FMADD(MB,M02,MB) NDIG = NDSAV1 - (MB(1)-M02(1)) ENDIF IF (NDIG.LT.2) NDIG = 2 IF (KFLAG.NE.1) THEN J = J + 1 GO TO 120 ENDIF NDIG = NDSAVE CALL FMI2M(-1,M02) CALL FMADD(MB,M02,M03) KFLAG = 0 RETURN ENDIF C C Here is the general case. Compute EXP(MA) where C 0 .LT. MA .LT. 1. C C Use the direct series C EXP(X) = 1 + X + X**2/2! + X**3/3! + ... C C The argument will be halved K2 times before the series C is summed. The series will be added as J2 concurrent C series. The approximately optimal values of K2 and J2 C are now computed to try to minimize the time required. C N2 is the approximate number of terms of the series which C will be needed, and L2 guard digits will be carried. C B = JBASE K = 5.0*LOG(10.0)/LOG(B) + 2.0 T = MAX(NDIG-K,2) ALOGB = LOG(B) ALOG2 = LOG(2.0) ALOGT = LOG(T) J2 = 0.081*ALOGB*T**0.3333 + 1.85 J2 = MAX(1,MIN(J2,LJSUMS/MXNDG2)) K2 = 0.7*SQRT(T*ALOGB/REAL(J2)) - 0.57*ALOGT + 0.5 C L = -(REAL(MA(1))*ALOGB+LOG(REAL(MA(2))/B + * REAL(MA(3))/(B*B)))/ALOG2 - 0.3 K2 = K2 - L IF (K2.LT.0) THEN K2 = 0 J2 = .43*SQRT(T*ALOGB/(ALOGT+REAL(L)*ALOG2)) + .33 ENDIF IF (J2.LE.1) J2 = 1 C N2 = T*ALOGB/(ALOGT+REAL(L)*ALOG2) L2 = LOG(REAL(N2)+2.0**K2)/ALOGB NDIG = NDIG + L2 IF (NDIG.GT.MXNDG2) THEN KFLAG = -9 CALL FMWARN MB(1) = KUNKNO MB(2) = 1 DO 150 J = 2, NDSAVE MB(J+1) = 0 150 CONTINUE NDIG = NDSAVE RETURN ENDIF NDSAV1 = NDIG C C Halve the argument K2 times. C CALL FMEQU(MA,M02,NDSAVE,NDIG) KTWO = 1 MAXVAL = MXBASE/2 IF (K2.GT.0) THEN DO 160 J = 1, K2 KTWO = 2*KTWO IF (KTWO.GT.MAXVAL) THEN CALL FMDIVI(M02,KTWO,M02) KTWO = 1 ENDIF 160 CONTINUE IF (KTWO.GT.1) CALL FMDIVI(M02,KTWO,M02) ENDIF C C Sum the series X + X**2/2! + X**3/3! + .... C Split into J2 concurrent sums and reduce NDIG while C computing each term in the sum as the terms get smaller. C CALL FMEQU(M02,MB,NDIG,NDIG) NTERM = 1 DO 170 J = 1, J2 CALL FMDIVI(MB,NTERM,MB) NTERM = NTERM + 1 KPT = (J-1)*(NDIG+1) + 1 CALL FMEQU(MB,MJSUMS(KPT),NDIG,NDIG) 170 CONTINUE CALL FMIPWR(M02,J2,M03) C 180 CALL FMMPY(MB,M03,MB) DO 190 J = 1, J2 CALL FMDIVI(MB,NTERM,MB) KPT = (J-1)*(NDSAV1+1) + 1 NDIG = NDSAV1 CALL FMADD(MJSUMS(KPT),MB,MJSUMS(KPT)) IF (KFLAG.EQ.1) GO TO 200 NDIG = NDSAV1 - (MJSUMS(KPT)-MB(1)) IF (NDIG.LT.2) NDIG = 2 NTERM = NTERM + 1 190 CONTINUE GO TO 180 C C Next put the J2 separate sums back together. C 200 KFLAG = 0 KPT = (J2-1)*(NDIG+1) + 1 CALL FMEQU(MJSUMS(KPT),M03,NDIG,NDIG) IF (J2.GE.2) THEN DO 210 J = 2, J2 CALL FMMPY(M02,M03,M03) KPT = (J2-J)*(NDIG+1) + 1 CALL FMADD(M03,MJSUMS(KPT),M03) 210 CONTINUE ENDIF C C Reverse the effect of halving the argument to C compute EXP(MA). C NDIG = NDSAV1 IF (K2.GT.0) THEN CALL FMI2M(2,MB) DO 220 J = 1, K2 CALL FMADD(M03,MB,M02) CALL FMMPY(M03,M02,M03) 220 CONTINUE ENDIF C CALL FMI2M(1,MB) CALL FMADD(M03,MB,MB) CALL FMEQU(MB,MB,NDSAV1,NDSAVE) NDIG = NDSAVE C RETURN END SUBROUTINE FMI2M(INTEG,MA) C C MA = INTEG C C Convert an integer to FM format. C C The conversion is exact if INTEG is less than JBASE**NDIG, C otherwise the result is an approximation. C C This routine performs the trace printing for the conversion. C FMIM is used to do the arithmetic. C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C NCALL = NCALL + 1 KSTACK(NCALL) = 16 IF (NTRACE.NE.0) CALL FMNTRI(2,INTEG,1) C CALL FMIM(INTEG,MA) C IF (NTRACE.NE.0) CALL FMNTR(1,MA,MA,1) NCALL = NCALL - 1 RETURN END SUBROUTINE FMIM(INTEG,MA) C C MA = INTEG. Internal integer conversion routine. C C The conversion is exact if INTEG is less than JBASE**NDIG. C Otherwise FMDM is used to get an approximation. C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C DIMENSION MA(LUNPCK) DOUBLE PRECISION X C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C KFLAG = 0 N1 = NDIG + 1 C C Check for INTEG equal to zero. C INT = ABS(INTEG) IF (INT.EQ.0) THEN DO 110 J = 1, N1 MA(J) = 0 110 CONTINUE GO TO 150 ENDIF C C Compute and store the digits, right to left. C MA(1) = 0 J = NDIG + 1 C 120 K = INT/JBASE L = INT - K*JBASE MA(1) = MA(1) + 1 MA(J) = L IF (K.GT.0) THEN INT = K J = J - 1 IF (J.GE.2) GO TO 120 C C Here INTEG cannot be expressed exactly. C X = INTEG CALL FMDM(X,MA) GO TO 150 ENDIF C C Normalize MA. C KB = N1 - J + 2 JM2 = J - 2 DO 130 J = 2, KB MA(J) = MA(J+JM2) 130 CONTINUE KB1 = KB + 1 IF (KB1.LE.N1) THEN DO 140 J = KB1, N1 MA(J) = 0 140 CONTINUE ENDIF C IF (INTEG.LT.0) MA(2) = -MA(2) C 150 RETURN END SUBROUTINE FMINP(LINE,MA,NPT,LB) C C Convert an A1 character string to floating point multiple precision C format. C C LINE is an A1 character array of length LB to be converted C to FM format and returned in MA. C NPT is a pointer telling the routine where in the array to begin C the conversion. This allows more than one number to be stored C in an array and converted in place. C LB is a pointer to the last character of the field for that number. C C The input number may be in integer or any real format. C The convention is made that if no digits appear before 'E' then 1.0 C is assumed. For example 'E6' is converted as '1.0E+6'. C In exponential format the 'E' may also be 'D', 'Q', or 'M'. C C So that FMINP will convert any output from FMOUT, LINE is tested C to see if the input is one of the special symbols +OVERFLOW, C -OVERFLOW, +UNDERFLOW, -UNDERFLOW, or UNKNOWN. C For user input the abbreviations OVFL, UNFL, UNKN may be used. C PARAMETER ( MXNDIG=256 , NBITS=32 , * LPACK = (MXNDIG+1)/2 + 1 , LUNPCK = (6*MXNDIG)/5 + 20, * LMWA = 2*LUNPCK , LJSUMS = 8*LUNPCK, * LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6 ) C CHARACTER LINE(LB),KBLANK,KMINUS,KOVFL(4),KUNFL(4),KUNKN(4) DIMENSION JTRANS(8,4) DIMENSION MA(LUNPCK) C COMMON /FMUSER/ NDIG,JBASE,JFORM1,JFORM2,KRAD, * KW,NTRACE,LVLTRC,KFLAG,KWARN,KROUND C DOUBLE PRECISION DPMAX C COMMON /FM/ MWA(LMWA),NCALL,MXEXP,MXEXP2,KARGSW,KEXPUN,KEXPOV, * KUNKNO,IUNKNO,RUNKNO,MXBASE,MXNDG2,KSTACK(19), * MAXINT,SPMAX,DPMAX C C Scratch array usage during FMINP: M01 - M05 C COMMON /FM1/ M01(LUNPCK),M02(LUNPCK),M03(LUNPCK),M04(LUNPCK), * M05(LUNPCK),M06(LUNPCK) C C Simulate a finite-state automaton to scan the input line C and build the number. States of the machine: C C 1. Initial entry to the subroutine C 2. Sign of the number C 3. Scanning digits before a decimal point C 4. Decimal point C 5. Scanning digits after a decimal point C 6. E, D, Q, or M - precision indicator before the exponent C 7. Sign of the exponent C 8. Scanning exponent C 9. Syntax error C C Character types recognized by the machine: C C 1. Sign (+,-) C 2. Numeral (0,1,...,9) C 3. Decimal point (.) C 4. Precision indicator (E,D,Q,M) C 5. Illegal character for number C C All blanks are ignored. The analysis of the number proceeds as C follows: If the simulated machine is in state JSTATE and a character C of type JTYPE is encountered the new state of the machine is given by C JTRANS(JSTATE,JTYPE). C C In this DATA statement note the array is loaded by columns. C C State 1 2 3 4 5 6 7 8 C Type DATA JTRANS/ 2, 9, 9, 9, 9, 7, 9, 9, N 3, 3, 3, 5, 5, 8, 8, 8, D 4, 4, 4, 9, 9, 9, 9, 9, P 6, 6, 6, 6, 6, 9, 9, 9 / C DATA KBLANK/' '/,KMINUS/'-'/,KOVFL/'O','V','F','L'/, * KUNFL/'U','N','F','L'/, KUNKN/'U','N','K','N'/ C NCALL = NCALL + 1 KSTACK(NCALL) = 17 IF (ABS(NTRACE).GE.2 .AND. NCALL.LE.LVLTRC) THEN WRITE (KW,110) (LINE(L),L=NPT,LB) 110 FORMAT(' INPUT TO FMINP'/(14X,65A1)) ENDIF NDSAVE = NDIG KASAVE = KARGSW KARGSW = 1 KRSAVE = KROUND KROUND = 1 KFLAG = 0 C C Since arithmetic tracing is not usually desired during C I/O conversion, disable tracing during this routine. C NTRSAV = NTRACE NTRACE = 0 C C Check for special symbols. C KMN = 1 KOF = 1 KUF = 1 KUK = 1 DO 120 J = NPT, LB IF (LINE(J).EQ.KMINUS) KMN = -1 IF (LINE(J).EQ.KOVFL(KOF)) THEN KOF = KOF + 1 IF (KOF.EQ.5) THEN CALL FMIM(0,MA) MA(1) = KEXPOV MA(2) = KMN GO TO 250 ENDIF ENDIF IF (LINE(J).EQ.KUNFL(KUF)) THEN KUF = KUF + 1 IF (KUF.EQ.5) THEN CALL FMIM(0,MA) MA(1) = KEXPUN MA(2) = KMN GO TO 250 ENDIF ENDIF IF (LINE(J).EQ.KUNKN(KUK)) THEN KUK = KUK + 1 IF (KUK.EQ.5) THEN CALL FMIM(0,MA) MA(1) = KUNKNO MA(2) = 1 GO TO 250 ENDIF ENDIF 120 CONTINUE C C Increase the working precision. C IF (NCALL.EQ.1) THEN B = JBASE K = 5.0*LOG(10.0)/LOG(B) + 2.0 NDIG = MAX(NDIG+K,2) IF (NDIG.GT.MXNDG2) THEN KFLAG = -9 CALL FMWARN MA(1) = KUNKNO MA(2) = 1 DO 130 J = 2, NDSAVE MA(J+1) = 0 130 CONTINUE GO TO 250 ENDIF ENDIF NDSAV1 = NDIG KSTART = NPT KSTOP = LB JSTATE = 1 KSIGN = 1 N1 = NDIG + 1 DO 140 J = 1, N1 M02(J) = 0 M03(J) = 0 M04(J) = 0 M05(J) = 0 140 CONTINUE C C If JBASE is a power of ten then call FMINP2 for a much C faster input conversion. C KPOWER = LOG10(REAL(JBASE)) + 0.5 IF (JBASE.EQ.10**KPOWER) THEN CALL FMINP2(MA,LINE,KSTART,KSTOP,JTRANS,KBLANK,KPOWER) GO TO 240 ENDIF C N2 = 0 KSIGNX = 1 KF1 = 0 KF2 = 0 KEXP = 0 KTENF1 = 1 KTENF2 = 1 KTENEX = 1 K10PWR = 0 C C Large is a threshold used in order to do as much of the C conversion as possible in one-word integer arithmetic. C LARGE = (MXBASE*MXBASE - 10)/10 C C KDFLAG will be 1 if any digits are found before 'E'. C KDFLAG = 0 C C Scan the number. C DO 210 J = KSTART, KSTOP IF (LINE(J).EQ.KBLANK) GO TO 210 CALL FMINPT(LINE(J),KTYPE,KVAL) IF (KTYPE.GE.5) GO TO 260 C JSTATE = JTRANS(JSTATE,KTYPE) C GO TO (260,150,160,210,170,180,190,200,260),JSTATE C C State 2. Sign of the number. C 150 KSIGN = KVAL GO TO 210 C C State 3. Digits before a decimal point. C 160 KDFLAG = 1 KF1 = 10*KF1 + KVAL KTENF1 = 10*KTENF1 IF (KTENF1.GT.LARGE) THEN IF (KTENF1.NE.K10PWR .AND. M03(2).NE.0) THEN CALL FMI2M(KTENF1,MA) K10PWR = KTENF1 ENDIF IF (M03(2).EQ.0) THEN CALL FMI2M(KF1,M03) ELSE NDIG = MAX(2,MIN(M03(1)+MA(1),NDSAV1)) CALL FMMPY(M03,MA,M03) NDIG = NDSAV1 CALL FMI2M(KF1,M02) NDIG = MAX(2,MIN(MAX(M03(1),M02(1))+1,NDSAV1)) IF (KF1.NE.0) CALL FMADD(M03,M02,M03) NDIG = NDSAV1 ENDIF KF1 = 0 KTENF1 = 1 ENDIF GO TO 210 C C State 5. Digits after a decimal point. C 170 KDFLAG = 1 N2 = N2 + 1 KF2 = 10*KF2 + KVAL KTENF2 = 10*KTENF2 IF (KTENF2.GT.LARGE) THEN IF (KTENF2.NE.K10PWR .AND. M04(2).NE.0) THEN CALL FMI2M(KTENF2,MA) K10PWR = KTENF2 ENDIF IF (M04(2).EQ.0) THEN CALL FMI2M(KF2,M04) ELSE NDIG = MAX(2,MIN(M04(1)+MA(1),NDSAV1)) CALL FMMPY(M04,MA,M04) NDIG = NDSAV1 CALL FMI2M(KF2,M02) NDIG = MAX(2,MIN(MAX(M04(1),M02(1))+1,NDSAV1)) IF (KF2.NE.0) CALL FMADD(M04,M02,M04) NDIG = NDSAV1 ENDIF KF2 = 0 KTENF2 = 1 ENDIF GO TO 210 C C State 6. Precision indicator. C 180 IF (KDFLAG.EQ.0) CALL FMI2M(1,M03) GO TO 210 C C State 7. Sign of the exponent. C 190 KSIGNX = KVAL GO TO 210 C C State 8. Digits of the exponent. C 200 KEXP = 10*KEXP + KVAL KTENEX = 10*KTENEX IF (KTENEX.GT.LARGE) THEN IF (KTENEX.NE.K10PWR .AND. M05(2).NE.0) THEN CALL FMI2M(KTENEX,MA) K10PWR = KTENEX ENDIF IF (M05(2).EQ.0) THEN CALL FMI2M(KEXP,M05) ELSE NDIG = MAX(2,MIN(M05(1)+MA(1),NDSAV1)) CALL FMMPY(M05,MA,M05) NDIG = NDSAV1 CALL FMI2M(KEXP,M02) NDIG = MAX(2,MIN(MAX(M05(1),M02(1))+1,NDSAV1)) IF (KEXP.NE.0) CALL FMADD(M05,M02,M05) NDIG = NDSAV1 ENDIF KEXP = 0 KTENEX = 1 ENDIF C 210 CONTINUE C C Form the number and return. C MA = KSIGN*(M03 + M04/10.0**N2)*10.0**M05 C IF (KTENF1.GT.1) THEN IF (KTENF1.NE.K10PWR