Diffusion COSOlver (improved_sugama) Browse

FM_User_Manual.txt
cd1286e2d410improved_sugama

FM_User_Manual.txt
View Options



			=========================================== Contents ===========================================


			Installing the FM package

			Troubleshooting

			User's guide for the FM package

			Functions available

			Converting a Fortran program to use the FM package

			FMZM Interface Notes

			Initialization and internal names

			Conversion functions

			Inquiry functions

			Multiple precision versions of Fortran operations and functions

			Formatting functions

			Integer functions

			Special functions

			Operations defined using subroutines

			Array operations

			Array functions

			FM.f95 Notes -- low-level routines




			====================================================================================================
			================================== Installing the FM package ===================================




			The files for version 1.3:

			1. FMSAVE.f95 Module for FM internal global variables

			2. FM.f95 Subroutine library for multiple-precision operations

			3. FMZM90.f95 Modules for interfaces and definitions of derived-types

			4. TestFM.f95 Test program for the FM routines

			5. SampleFM.f95 Small sample program using FM


			The first three files form the FM library routines that are used by an application program.
			They need to be compiled once, then any application program using FM can be compiled and
			linked to these library object files.

			TestFM.f95 is a large program that does thousands of tests and calls all the library routines.
			It checks each operation and should be run once after compiling the FM library routines, to
			make sure FM is properly installed.

			SampleFM.f95 is a small program with several examples showing how a typical user's program
			would use FM for some multiple precision real, integer, and complex calculations. It should
			be helpful as a model for getting started with FM.




			Here are two example sets of compiler/linker commands for building the programs using the
			gfortran compiler (free -- click on the download link "Binaries for Windows, Linux, MacOS
			and much more" from this page: gcc.gnu.org/wiki/GFortran).

			The first three files are compiled as object code libraries FMSAVE.o, FM.o, FMZM90.o, and then
			each program that uses FM is compiled and linked to those three libraries. Some compilers name
			these object files FMSAVE.obj, etc.

			Most compilers also produce files FMVALS.mod, FMZM.mod, etc., containing module information from
			the first three files.


			1. For Windows, after installing the compiler, run in a PowerShell command prompt window:


			gfortran fmsave.f95 -c -O3

			gfortran FM.f95 -c -O3

			gfortran FMZM90.f95 -c -O3


			gfortran TestFM.f95 -c -O3

			gfortran fmsave.o FM.o FMZM90.o TestFM.o -o TestFM.exe

			./TestFM


			gfortran SampleFM.f95 -c -O3

			gfortran fmsave.o FM.o FMZM90.o SampleFM.o -o SampleFM.exe

			./SampleFM


			2. For a Mac, after installing the compiler, run in a Terminal window:


			gfortran fmsave.f95 -c -O3

			gfortran FM.f95 -c -O3

			gfortran FMZM90.f95 -c -O3


			gfortran TestFM.f95 -c -O3

			gfortran fmsave.o FM.o FMZM90.o TestFM.o -o TestFM

			./TestFM


			gfortran SampleFM.f95 -c -O3

			gfortran fmsave.o FM.o FMZM90.o SampleFM.o -o SampleFM

			./SampleFM



			The compiler options used in the examples were:
			-c compile to object code -- don't make executable
			-O3 optimization level 3
			-o output file name for the executable program


			Most other compilers use options very similar to these. The g95 and NAG compilers use the same
			commands ("g95 fmsave.f95 -c -O3" and "nagfor fmsave.f95 -c -O3" respectively), and Lahey's
			compiler uses
			lf95 fmsave.f95 -ap -c -o1
			...
			lf95 fmsave.obj FM.obj FMZM90.obj TestFM.obj -out TestFM.exe
			Older versions of Lahey's compiler needed the -ap option (see below) for all files -- I have
			not tested the latest version.

			Some compilers have options that might improve the speed of the program beyond the basic -O
			optimization. For example, "gfortran fm.f95 -c -O3 -funroll-loops". This runs slightly
			faster when precision is above 100 digits.

			Some programs need 64-bit integers, and the easiest way to get them is often by using
			command-line compiler directives to change all integer constants and variables in a program.

			Prior to June, 2015, the FM package was not compatible with 64-bit integers, since a few routines
			assumed that any integer value could be represented exactly when converted to double precision.
			But 64-bit integers can have more than 16 decimal digits, causing errors when converted to double
			precision. With the current version, FM can be used with either 32-bit or 64-bit integers.

			With the gfortran compiler the command "gfortran prog.f95 -fdefault-integer-8 -o prog"
			will make 64-bit integers the default. When using 64-bit integers, all FM files and all
			files from the user's program should be compiled with the -fdefault-integer-8 option.

			From the user's point of view, the only difference in results that come from the previous version
			and this one is that using 64-bit integers allows the range between FM's underflow and overflow
			thresholds to be greater.

			In previous versions of FM, exp(1.0e+8) = 1.5500e+43429448 did not overflow, while exp(1.0e+9)
			overflowed. Using 64-bit integers with this version, exp(1.0e+15) = 6.7244e+434294481903251
			does not overflow, while exp(1.0e+16) overflows.

			FM has been run using many different compilers, both free and commercial. Most are used in a
			similar way to gfortran, although many can also be used with development environments that can
			make the process of compiling, linking, and executing the files even easier.

			Two other files define optional multiprecision operations for exact rational arithmetic
			(fm_rational.f95) and interval arithmetic (fm_interval.f95). To use these, compile and
			link as with the examples above. See the FM web page for testing and sample programs
			for each of these, analogous to TestFM.f95 and SampleFM.f95.





			----------------------------------------------------------------------------------------------------
			--------------------------------------- Troubleshooting ----------------------------------------



			After compiling and running the programs TestFM and SampleFM, each should say "no errors were
			found" at the end. If there were problems compiling the programs or some errors were found
			when they ran, read the rest of this section for possible fixes.


			1. After downloading the files, if the compiler gives many error messages or it appears to see
			no code in the file at all, check that the lines in the file have the proper end-of-line
			characters for your system.

			For a PC, this means each line should end with both a carriage return <cr> character (ascii 13)
			and a line feed <lf> character (ascii 10). If the file appears to be one huge line when viewed
			in an editor, one of these two characters is probably missing and should be added to each line.

			To use FM on a Unix system, lines end with <lf>, and for a Mac system they might end with <cr>,
			but should end in <lf> because the Terminal window is really Unix. On these systems, failing
			to fix the end-of-line characters might mean the file seems to have twice the expected number
			of lines, with a blank line between each line of code when viewed in an editor, or the compiler
			might give an error message like "Unexpected end of file".

			For a Mac with <cr> lines, a file like fmsave.f95 can be fixed by changing its name to
			fmsave0.f95 and then giving this command in Terminal:

			tr \\r \\n < fmsave0.f95 > fmsave.f95

			Then the new file fmsave.f95 should have the proper format.

			Some good text editors will recognize a foreign end-of-line format and automatically fix each
			file the first time it is opened. Many recent Fortran compilers also automatically convert
			input source files to the proper type.

			2. The compiler gives an "out of memory" error message or crashes during compile of one or
			more of the files.

			It might be necessary to break the file into smaller pieces or split it into separate files
			for each routine or module. This could be caused by lack of system memory, lack of virtual
			memory, or a bug (memory leak) in the compiler.

			Some compilers have an option (e.g., -split) to do this automatically.

			This problem is less likely on more recent compilers and computers. With several recent
			compilers, each file can compile (barely) on machines with 1/2 Gb of memory, and should
			compile easily with 1 Gb or more.

			3. Most of the routines compile, but a few fail with error messages like
			"symbol 120 is not the label of a branch target statement".
			However, looking at the code shows there is a label 120 in that routine.

			This might happen in the larger routines. Some older compilers may require additional options
			be enabled (e.g., to force 32-bit branches or addresses to be used). Check in the compiler
			manual and try turning on any options that mention "long branches", "32-bit addresses", etc.

			4. All files compile, but the TestFM program reports a few errors when it runs. There are other
			possibilities, but one thing to check is whether the compiler has any options controlling
			arithmetic precision of intermediate results.

			Because the FM numbers are stored as integer values in double precision arrays, any sloppy
			rounding can cause problems. In one case, a compiler optimized an expression by leaving the
			result of a division in an 80-bit register and then used that result later in the calculation.
			Rounding the division back to double precision would have fixed the error, but using the
			inaccurate extended precision value caused the final result to be off by one when it was
			returned to an integer value.

			This compiler had an option (-ap) to force intermediate results to not be left in registers,
			and that fixed the problem.

			Another way to check to see if this is the problem is to create a version of FM that uses
			integer arrays instead of double precision. See the section titled "EFFICIENCY" below to see
			how to make this change. On most current machines, there is not a large speed penalty for
			using integer arrays as long as the precision is under 100 significant digits (i.e., NDIG < 15
			or so with MBASE = 10**7).

			5. TestFM or SampleFM gives an error message beginning something like this:

			Element ( 1 ) of a multiple precision one-dimensional array is undefined in an expression.

			This could happen with older compilers that support Fortran-90 but not Fortran-95. Version 1.3
			of FM needs Fortran-95 or later. Try upgrading to the latest version of your compiler, or
			try the latest version of the free gfortran compiler to check to see if this is causing that
			error message.





			====================================================================================================
			=============================== User's guide for the FM package ================================




			The various lists of available multiple precision operations and routines have been collected here,
			along with some general advice on using the package.

			See the program SampleFM.f95 for some examples of initializing and using the package.





			----------------------------------------------------------------------------------------------------
			------------------------------------ Functions available ---------------------------------------



			The FM package provides 5 types of multiple precision operations:

			1. type(FM) - floating-point

			2. type(IM) - integer

			3. type(ZM) - complex

			4. type(FM_RATIONAL) - rational

			5. type(FM_INTERVAL) - interval

			The type definitions and interfaces for the first three are in file fmzm90.f95, and the other
			two are in files fm_rational.f95 and fm_interval.f95.

			Some multiple precision functions take input arguments that are intrinsic types. In the table
			below, the types of arguments allowed for each function are abbreviated as:
			FM, IM, ZM, RAT, IVL, for the 5 types above, and int, sp, dp, spz, dpz, str for the intrinsic
			Fortran types integer, single precision real, double precision real, single precision complex,
			double precision complex, character string.

			Further description of many of these functions can be found later in this file.


			ABS FM IM ZM RAT IVL
			ACOS FM ZM IVL
			ACOSH FM ZM IVL
			AIMAG ZM
			AINT FM ZM IVL
			ANINT FM ZM IVL
			ASIN FM ZM IVL
			ASINH FM ZM IVL
			ATAN FM ZM IVL
			ATANH FM ZM IVL
			ATAN2 FM IVL
			BERNOULLI int
			BESSEL_J0 FM IVL
			BESSEL_J1 FM IVL
			BESSEL_JN FM IVL
			BESSEL_Y0 FM IVL
			BESSEL_Y1 FM IVL
			BESSEL_YN FM IVL
			BETA FM IVL
			BINOMIAL FM IM IVL int
			BTEST IM
			CEILING FM IM ZM RAT IVL
			CMPLX FM IM
			CONJG ZM
			COS FM ZM IVL
			COSH FM ZM IVL
			COS_INTEGRAL FM IVL
			COSH_INTEGRAL FM IVL
			DBLE FM IM ZM IVL
			DIGITS FM IM ZM IVL
			DIM FM IM RAT IVL
			DINT FM ZM IVL
			DOT_PRODUCT FM IM ZM RAT IVL
			EPSILON FM IVL
			ERF FM IVL
			ERFC FM IVL
			ERFC_SCALED FM
			EXP FM ZM IVL
			EXPONENT FM IVL
			EXP_INTEGRAL_EI FM IVL
			EXP_INTEGRAL_EN FM IVL
			FACTORIAL FM IM IVL int
			FLOOR FM IM ZM RAT IVL
			FRACTION FM ZM IVL
			FRESNEL_C FM IVL
			FRESNEL_S FM IVL
			GAMMA FM IVL
			GCD IM
			HUGE FM IM ZM IVL
			HYPOT FM
			INCOMPLETE_BETA FM IVL
			INCOMPLETE_GAMMA1 FM IVL
			INCOMPLETE_GAMMA2 FM IVL
			INT FM IM ZM RAT IVL
			IS_OVERFLOW FM IM ZM IVL
			IS_UNDERFLOW FM IM ZM IVL
			IS_UNKNOWN FM IM ZM RAT IVL
			LEFT_ENDPOINT IVL
			LOG FM ZM IVL
			LOG10 FM ZM IVL
			LOG_ERFC FM IVL
			LOG_GAMMA FM IVL
			LOG_INTEGRAL FM IVL
			MATMUL FM IM ZM RAT IVL
			MAX FM IM RAT IVL
			MAXEXPONENT FM IVL
			MAXLOC FM IM ZM RAT
			MAXVAL FM IM ZM RAT IVL
			MIN FM IM RAT IVL
			MINEXPONENT FM IVL
			MINLOC FM IM ZM RAT
			MINVAL FM IM ZM RAT IVL
			MOD FM IM RAT IVL
			MODULO FM IM RAT IVL
			MULTIPLY_MOD IM
			NEAREST FM IVL
			NINT FM IM ZM RAT IVL
			NORM2 FM
			POCHHAMMER FM IVL
			POLYGAMMA FM IVL
			POWER_MOD IM
			PRECISION FM ZM IVL
			PRODUCT FM IM ZM RAT IVL
			PSI FM IVL
			RADIX FM IM ZM IVL
			RANGE FM IM ZM IVL
			RATIONAL_DENOMINATOR RAT
			RATIONAL_NUMERATOR RAT
			REAL FM IM ZM IVL
			RIGHT_ENDPOINT IVL
			RRSPACING FM IVL
			SCALE FM ZM IVL
			SETEXPONENT FM IVL
			SIGN FM IM IVL
			SIN FM ZM IVL
			SINH FM ZM IVL
			SIN_INTEGRAL FM IVL
			SINH_INTEGRAL FM IVL
			SPACING FM IVL
			SQRT FM ZM IVL
			SUM FM IM ZM RAT IVL
			TAN FM ZM IVL
			TANH FM ZM IVL
			TO_DP FM IM ZM
			TO_DPZ FM IM ZM
			TO_FM FM IM ZM RAT IVL int sp dp spz dpz str
			TO_FM_INTERVAL FM IM int sp dp str
			TO_FM_RATIONAL IM int str
			TO_IM FM IM ZM RAT IVL int sp dp spz dpz str
			TO_INT FM IM ZM
			TO_SP FM IM ZM
			TO_SPZ FM IM ZM
			TO_ZM FM IM ZM RAT IVL int sp dp spz dpz str
			TINY FM IM ZM IVL
			TRANSPOSE FM IM ZM RAT





			----------------------------------------------------------------------------------------------------
			------------------------- Converting a program to use the FM package ---------------------------



			0. If you want to write a program from scratch that uses FM, instead of converting an existing
			double precision version, consider writing a d.p. version first anyway. It is very useful to
			have a working d.p. version to compare the FM results and quickly locate large errors that
			might be caused by mistakes in the conversion.


			1. Before any variable declarations in each PROGRAM, SUBROUTINE, MODULE, or FUNCTION that will
			use multiple precision variables, insert

			USE FMZM

			This module contains all the rules needed by the compiler for doing the multiple precision
			operations.


			2. In all routines using multiple precision variables,
			change real or double precision declarations to TYPE (FM)
			change complex or complex d.p. declarations to TYPE (ZM)
			if any integers need to be multiple precision, declare as TYPE (IM)

			For example, if the original main program had these declarations,

			REAL (KIND(1.0D0)) :: X, Y, A(50)
			COMPLEX (KIND(1.0D0)) :: C, Z(20)

			change them to this for the FM version.

			TYPE (FM) :: X, Y, A(50)
			TYPE (ZM) :: C, Z(20)

			Local multiple precision variables in a subprogram are those that are not input arguments,
			module variables, or function names. These local variables should be saved, because a compiler
			might create a new instance of the local variable each time the routine is called. That would
			not cause an error, but it would leak memory and possibly cause the program to fail later for
			lack of memory. So if the original code was part of a subroutine,

			SUBROUTINE SOLVE(A,Z,C)
			REAL (KIND(1.0D0)) :: X, Y, A(50)
			COMPLEX (KIND(1.0D0)) :: C, Z(20)

			then the only local variables were X and Y, so the new declaration becomes

			TYPE (FM) :: X, Y, A(50)
			TYPE (ZM) :: C, Z(20)
			SAVE :: X, Y

			or equivalently,

			TYPE (FM), SAVE :: X, Y
			TYPE (FM) :: A(50)
			TYPE (ZM) :: C, Z(20)

			Using SAVE for variables in the main program is ok, but not required.


			3. Variables that were initialized in the declarations of the original program must
			be initialized separately as FM variables.

			DOUBLE PRECISION :: X = 1.2D0

			becomes

			TYPE (FM), SAVE :: X
			...
			X = TO_FM('1.2')

			If this is in a subroutine and the value might change during one call and need to be
			remembered in a subsequent call, something like this can be done:

			TYPE (FM), SAVE :: X
			LOGICAL, SAVE :: FIRST_CALL = .TRUE.
			...
			IF (FIRST_CALL) THEN
			X = TO_FM('1.2')
			FIRST_CALL = .FALSE.
			ENDIF


			4. At the beginning of the main program, call FM_SET to set the FM precision.
			For example, to get 50 significant digits,

			CALL FM_SET(50)

			Since increasing FM's internal precision level by one gives several extra base 10 significant
			digits, this call will actually set the user's precision to slightly more than 50 digits.


			5. At the beginning of every function subprogram that returns a multiple precision value
			insert

			CALL FM_ENTER_USER_FUNCTION(your_function_name)

			where "your_function_name" is replaced by the actual name of the function.

			Before each RETURN statement and before the END FUNCTION statement if it is not
			immediately preceded by a RETURN, insert

			CALL FM_EXIT_USER_FUNCTION(your_function_name)

			These are needed so that FM does not discard temporary variables too soon when they
			are created by the FMZM interface routines. The function name itself is used as a
			variable name to return the function value, so it is one of these temporary variables.

			If a function subprogram uses multiple precision variables but returns a value that is
			not multiple precision, then the function name is not a temporary multiple precision
			object, but FM still needs to be notified not to discard temporaries within the function
			until that function is done. In this case, use

			CALL FM_ENTER_USER_ROUTINE

			and

			CALL FM_EXIT_USER_ROUTINE


			6. Check constants that are now part of multiple precision expressions and convert them.

			X = Y/3 need not be converted (since integers are exact in binary), but

			X = Y/3.7 should become X = Y/TO_FM('3.7')
			Since 3.7 is not represented exactly in the machine's single or double
			precision, leaving the statement as X = Y/3.7 would give X accurate only
			the machine's single precision, even though X and Y are multiple precision.

			Also, constants in routine argument lists that now refer to multiple precision variables
			must be converted.

			CALL SUB(A,B,2.6D0,X)

			becomes

			C = TO_FM('2.6')
			CALL SUB(A,B,C,X)

			Arithmetic expressions in subroutine calls should also be removed if they will use
			multiple precision in the new version, since they also create temporary FM variables.

			CALL SUB(A,B,C+3*D,X)

			becomes

			T = C+3*D
			CALL SUB(A,B,T,X)

			If many such argument expressions occur in a program, an alternate way to convert them is
			by leaving the temporary object in the argument list, like

			CALL SUB(A,B,TO_FM('2.6'),X)
			CALL SUB(A,B,C+3*D,X)

			and then inserting calls to FM_ENTER_USER_ROUTINE and FM_EXIT_USER_ROUTINE into subroutine SUB
			as described in section 5 above.


			7. If a routine uses allocatable type FM, ZM, or IM arrays and allocates and deallocates with
			each call, then a call to FM_DEALLOCATE should be made before actually deallocating each array,
			to avoid leaking memory in FM's internal variable database. For example:

			DEALLOCATE(T)

			becomes:

			CALL FM_DEALLOCATE(T)
			DEALLOCATE(T)


			8. Multiple precision variables in WRITE statements can be handled in several ways:

			(a) If we use higher precision arithmetic for the calculations, but we only need to see the
			final output at double precision, the simplest option is to convert the multiple precision
			variables back to double for printing. Then no changes to formats are needed.

			WRITE (*,"(' Step size = ',F15.6,' tolerance = ',E15.7)"),H,T

			becomes

			WRITE (*,"(' Step size = ',F15.6,' tolerance = ',E15.7)"),TO_DP(H),TO_DP(T)

			now that H and T are TYPE (FM) variables.

			If the FM automatic tracing option is on (see NTRACE below), some Fortran-95 compilers
			might generate a "recursive write" error message here, since another write statement would
			be executed during the TO_DP call. A fix is to turn the tracing off before this write
			statement. This can also happen if an FM error message is written during the TO_DP call.

			(b) Format the writes for multiple precision.

			WRITE (*,"(' Step size = ',F15.6,' tolerance = ',E15.7)"),H,T

			becomes

			WRITE (*,"(' Step size = ',A,' tolerance = ',A)"), &
			TRIM(FM_FORMAT('F15.6',H)),TRIM(FM_FORMAT('E15.7',T))

			FM_FORMAT is a formatting function used when the number of digits being shown
			is small enough to fit on one line.

			Often after converting to multiple precision, we want to see more digits, so here
			F15.6 and E15.7 might become F35.20 and E35.25 in the FM version.

			(c) Subroutine FM_FORM does similar formatting, but we supply a character string for
			the formatted result. After declaring the strings at the top of the routine, as with

			CHARACTER(80) :: ST1,ST2

			the WRITE above could become

			CALL FM_FORM('F15.6',H,ST1)
			CALL FM_FORM('E15.7',T,ST2)
			WRITE (*,"(' Step size = ',A,' tolerance = ',A)") TRIM(ST1),TRIM(ST2)

			FM_FORM must be used instead of FM_FORMAT when there are more than 200 characters
			in the formatted string. These longer numbers usually need to be broken into several
			lines.

			FM_FORM should also be used when the FM trace option is on, since some compilers may
			generate an error message about a "recursive I/O reference" if a trace write executes
			from within another write statement via FM_FORMAT.

			(d) To use the current FM default format and handle any line breaks automatically, subroutine
			FM_PRINT can be used. Calling FM_SET to set precision at the beginning of the program also
			initializes this format. For example, FM_PRINT displays 50 significant digits after
			CALL FM_SET(50). See the discussion of FM's settings for JFORM1, JFORM2, and KSWIDE for
			changing the default format. The FM numbers will print on separate lines.

			WRITE (,) ' Step size = '
			CALL FM_PRINT(H)
			WRITE (,) ' Tolerance = '
			CALL FM_PRINT(T)


			9. Multiple precision variables in READ statements can be done with FM's free-format input:

			(a) Read the line as a character string then convert using TO_FM.

			READ (,) A,B,C

			becomes this (with ST1 declared at the top of the routine as CHARACTER with
			length large enough to hold each input data line).

			READ (*,"(A)") ST1
			CALL FMSCAN(ST1,1,JA,JB)
			A = TO_FM(ST1(JA:JB))
			J1 = JB
			CALL FMSCAN(ST1,J1,JA,JB)
			B = TO_FM(ST1(JA:JB))
			J1 = JB
			CALL FMSCAN(ST1,J1,JA,JB)
			C = TO_FM(ST1(JA:JB))

			Where FMSCAN is defined by:

			SUBROUTINE FMSCAN(STRING,JSTART,JA,JB)
			! Scan STRING from position JSTART and return JA as the next non-blank and
			! JB as the next blank after JA.
			CHARACTER (*) :: STRING
			JA = 0
			JB = 0
			DO J = JSTART, LEN(STRING)
			IF (JA == 0) THEN
			IF (STRING(J:J) /= ' ') JA = J
			ELSE
			IF (STRING(J:J) == ' ') THEN
			JB = J
			RETURN
			ENDIF
			ENDIF
			ENDDO
			END SUBROUTINE FMSCAN

			This assumes all three numbers are on one line. If they could appear on two or three
			lines, more code would be needed to check for that.

			It also assumes that blanks separate the numbers. If input records use commas to separate
			numbers, repeat counts on input items, or slashes, then either the code above can be made
			more elaborate to handle those cases, or the data file can be edited so the simpler code
			works.

			(b) Formatted reads can be converted directly to calls to TO_FM without scanning to find
			where each number appears on the line.

			READ (*,"(F20.15,E26.16,E20.10)") A,B,C

			becomes this

			READ (*,"(A)") ST1
			A = TO_FM(ST1( 1:20))
			B = TO_FM(ST1(21:46))
			C = TO_FM(ST1(47:66))

			(c) Declare double precision variables so the original read statement still works, then
			convert to multiple precision.

			READ (,) A,B,C

			becomes this

			READ (,) A_DP,B_DP,C_DP
			A = A_DP
			B = B_DP
			C = C_DP

			where A_DP,B_DP,C_DP are double precision and A,B,C are multiple precision.

			A possible drawback to this method is that the values are read as double precision, so
			after conversion to FM they are usually accurate only to double precision. As with the
			TO_FM example in section 6 above, a number such as 3.7 is not exactly representable in
			binary, so if that is the value being read for A in this example, reading it as a string
			in (b) and then converting the string gives full multiple precision accuracy for 3.7,
			but reading it as in (c) gives double precision accuracy.



			See the page from the link "Automatic conversion of Fortran programs to use FM" on the main FM
			web page for a program that tries to automate most of these changes needed to convert a normal
			Fortran program to a multiple-precision version using FM.



			----------------------------------------------------------------------------------------------------
			------------------------------------ FMZM Interface Notes --------------------------------------



			The FMZM module extends the definition of the basic Fortran arithmetic and function operations
			so they also apply to multiple precision numbers.

			There are three multiple precision data types:
			FM (multiple precision real)
			IM (multiple precision integer)
			ZM (multiple precision complex)

			A routine using any of these types needs this statement at the top:
			USE FMZM

			For some examples and general advice about using these multiple-precision data types, see the
			program SampleFM.f95.

			Most of the functions defined in the FMZM module are multiple precision versions of standard
			Fortran functions. In addition, there are functions for direct conversion, formatting, and some
			mathematical special functions.

			An attempt to use a multiple precision variable that has not been defined will be detected by
			the routines in FMZM and an error message printed.



			----------------------------- Initialization and internal names --------------------------------



			Initialization: The default precision for the multiple-precision numbers is about 50
			significant digits.

			To set precision to a different value, put this

			CALL FM_SET(N)

			in the main program before any multiple precision operations are done, with N
			replaced by the number of decimal digits of accuracy to be used.

			Routine names: For each multiple precision operation there are several routines with related
			names that perform variations of that operation. For example, the addition
			operation has these forms:

			Using the FMZM interface module to perform real (floating-point) multiple
			precision addition, declare the variables as FM derived types with

			USE FMZM
			TYPE ( FM ) A,B,C

			and then after values are assigned to A and B, doing a multiple precision
			addition looks the same as if the variables were real or double.

			C = A + B

			Normally, using the interface module avoids the need to know the name of the
			FM routine being called. For some operations, usually those that are not
			numerical Fortran functions (such as formatting a number), a direct call may be
			needed. For the addition above there is no reason to write it as a call in the
			user's program, but it could be done as

			CALL FM_ADD(A,B,C)

			Routines with names starting with FM_ in the FMZM module (file FMZM90.f95) then
			call the low-level arithmetic routines in file FM.f95. The low-level routines
			are not usually called directly by the user's program, since they do not operate
			on the derived type variables that the user sees, but on the internal components
			of the types.

			The low-level routines in FM.f95 usually have similar names to those in
			FMZM90.f95, but with no underscore after the first two letters. In this case
			the low-level routine is named FMADD.

			For a few routines that don't have multiple precision arguments, like FM_SET and
			FM_SETVAR, the corresponding low-level names FMSET AND FMSETVAR are also
			available, and either form can be used.




			------------------------------------ Conversion functions --------------------------------------


			TO_FM is a function for converting other types of numbers to type FM. Note that TO_FM(3.12)
			converts the REAL constant to FM, but it is accurate only to single precision, since the number
			3.12 cannot be represented exactly in binary and has already been rounded to single precision.
			Similarly, TO_FM(3.12D0) agrees with 3.12 to double precision accuracy, and TO_FM('3.12') or
			TO_FM(312)/TO_FM(100) agrees to full FM accuracy.

			TO_IM converts to type IM, and TO_ZM converts to type ZM.

			Functions are also supplied for converting the three multiple precision types to the other
			numeric data types:
			TO_INT converts to machine precision integer
			TO_SP converts to single precision
			TO_DP converts to double precision
			TO_SPZ converts to single precision complex
			TO_DPZ converts to double precision complex

			WARNING: When multiple precision type declarations are inserted in an existing program, take
			care in converting functions like DBLE(X), where X has been declared as a multiple
			precision type. If X was single precision in the original program, then replacing
			the DBLE(X) by TO_DP(X) in the new version could lose accuracy. For this reason, the
			Fortran type-conversion functions defined in the module assume that results should
			be multiple precision whenever inputs are. Examples:
			DBLE(TO_FM('1.23E+123456')) is type FM
			REAL(TO_FM('1.23E+123456')) is type FM
			REAL(TO_ZM('3.12+4.56i')) is type FM = TO_FM('3.12')
			INT(TO_FM('1.23')) is type IM = TO_IM(1)
			INT(TO_IM('1E+23')) is type IM
			CMPLX(TO_FM('1.23'),TO_FM('4.56')) is type ZM




			------------------------------------- Inquiry functions ----------------------------------------


			IS_OVERFLOW, IS_UNDERFLOW, and IS_UNKNOWN are logical functions for checking whether a multiple
			precision number is in one of the exception categories. Testing to see if a type FM number is in
			the +overflow category by directly using an IF can be tricky. When X is +overflow, the statement

			IF (X == TO_FM(' +OVERFLOW ')) THEN

			will return false, since the comparison routine cannot be sure that two different overflowed
			results would have been equal if the overflow threshold had been higher. Instead, use

			IF (IS_OVERFLOW(X)) THEN

			which will be true if X is + or - overflow.




			------------------------ Multiple precision operations and functions ---------------------------


			For each of the operations =, == , /= , < , <= , > , >= , +, -, , /, and *, the FMZM
			interface module defines all mixed mode variations involving one of the three multiple precision
			derived types and another argument having one of these types:

			{ integer, real, double, complex, complex double, FM, IM, ZM }

			So mixed mode expressions such as

			X = 12
			X = X + 1
			IF (ABS(X) > 1.0D-23) THEN

			are handled correctly.

			Not all the named functions are defined for all three multiple precision derived types, so the
			list below shows which can be used. The labels "real", "integer", and "complex" refer to types
			FM, IM, and ZM respectively, "string" means the function accepts character strings (e.g.,
			TO_FM('3.45')), and "other" means the function can accept any of the machine precision data
			types integer, real, double, complex, or complex double. For functions that accept two or more
			arguments, like ATAN2 or MAX, all the arguments must be of the same type.

			TO_ZM also has a 2-argument form: TO_ZM(2,3) for getting 2 + 3*i.
			CMPLX can be used for that, as in CMPLX( TO_FM(2) , TO_FM(3) ), and so can the string form,
			TO_ZM(' 2 + 3 i '), but the 2-argument form is sometimes more convenient. The 2-argument form
			is available for machine precision integer, single and double precision real pairs. For others,
			such as X and Y being type(fm), just use CMPLX(X,Y).




			-------------- Multiple precision versions of Fortran operations and functions -----------------


			=
			+
			-
			*
			/
			**
			==
			/=
			<
			<=
			>
			>=
			ABS real integer complex
			ACOS real complex
			ACOSH real complex
			AIMAG complex
			AINT real complex
			ANINT real complex
			ASIN real complex
			ASINH real complex
			ATAN real complex
			ATANH real complex
			ATAN2 real
			BTEST integer
			CEILING real integer complex
			CMPLX real integer
			CONJG complex
			COS real complex
			COSH real complex
			DBLE real integer complex
			DIGITS real integer complex
			DIM real integer
			DINT real complex
			EPSILON real
			EXP real complex
			EXPONENT real
			FLOOR real integer complex
			FRACTION real complex
			HUGE real integer complex
			HYPOT real
			INT real integer complex
			LOG real complex
			LOG10 real complex
			MAX real integer
			MAXEXPONENT real
			MIN real integer
			MINEXPONENT real
			MOD real integer
			MODULO real integer
			NEAREST real
			NINT real integer complex
			NORM2 real
			PRECISION real complex
			RADIX real integer complex
			RANGE real integer complex
			REAL real integer complex
			RRSPACING real
			SCALE real complex
			SETEXPONENT real
			SIGN real integer
			SIN real complex
			SINH real complex
			SPACING real
			SQRT real complex
			TAN real complex
			TANH real complex
			TINY real integer complex


			----------------------------- Conversion and inquiry functions ---------------------------------


			TO_FM real integer complex string other
			TO_IM real integer complex string other
			TO_ZM real integer complex string other
			TO_INT real integer complex
			TO_SP real integer complex
			TO_DP real integer complex
			TO_SPZ real integer complex
			TO_DPZ real integer complex
			IS_OVERFLOW real integer complex
			IS_UNDERFLOW real integer complex
			IS_UNKNOWN real integer complex


			----------------------------------- Formatting functions ---------------------------------------


			FM_FORMAT real
			IM_FORMAT integer
			ZM_FORMAT complex


			------------------------------------ Integer functions -----------------------------------------


			BINOMIAL integer
			FACTORIAL integer
			GCD integer
			MULTIPLY_MOD integer
			POWER_MOD integer


			------------------------------------ Special functions -----------------------------------------


			BERNOULLI(N) real
			BESSEL_J0(X) real
			BESSEL_J1(X) real
			BESSEL_JN(N,X) real
			BESSEL_JN(N1,N2,X) real
			BESSEL_Y0(X) real
			BESSEL_Y1(X) real
			BESSEL_YN(N,X) real
			BESSEL_YN(N1,N2,X) real
			BETA(A,B) real
			BINOMIAL(A,B) real
			COS_INTEGRAL(X) real
			COSH_INTEGRAL(X) real
			ERF(X) real
			ERFC(X) real
			ERFC_SCALED(X) real
			EXP_INTEGRAL_EI(X) real
			EXP_INTEGRAL_EN(N,X) real
			FACTORIAL(X) real
			FRESNEL_C(X) real
			FRESNEL_S(X) real
			GAMMA(X) real
			INCOMPLETE_BETA(X,A,B) real
			INCOMPLETE_GAMMA1(A,X) real
			INCOMPLETE_GAMMA2(A,X) real
			LOG_ERFC(X) real
			LOG_GAMMA(X) real
			LOG_INTEGRAL(X) real
			POCHHAMMER(X,N) real
			POLYGAMMA(N,X) real
			PSI(X) real
			SIN_INTEGRAL(X) real
			SINH_INTEGRAL(X) real


			Several of these functions are described in more detail below.


			----------------- Subroutines that do not correspond to any function above ---------------------


			1. Type (FM). MA, MB, MC refer to type (FM) numbers.

			These are subroutines instead of functions, so they are invoked as with
			CALL FM_COS_SIN(MA,MB,MC)

			FM_COS_SIN(MA,MB,MC) MB = COS(MA), MC = SIN(MA)
			Faster than making two separate calls.

			FM_COSH_SINH(MA,MB,MC) MB = COSH(MA), MC = SINH(MA)
			Faster than making two separate calls.

			FM_EULER(MA) MA = Euler's constant ( 0.5772156649... )

			FM_EQU(MA,MB,NA,NB) MB = MA where precision is being changed.
			MA is defined with NDIG = NA digits and
			MB will be defined having NB digits.
			MB is rounded if NB < NA
			MB is zero-padded if NB > NA

			FM_FLAG(K) K = KFLAG get the value of the FM condition flag -- stored in
			the internal FM variable KFLAG in module FMVALS.

			FM_FORM(FORM,MA,STRING) MA is converted to a character string using format FORM and
			returned in STRING. FORM can represent I, F, E, or ES formats.
			Example:
			CALL FMFORM('F60.40',MA,STRING)

			FM_FPRINT(FORM,MA) Print MA on unit KW using FORM format.

			FM_PI(MA) MA = pi

			FM_PRINT(MA) Print MA on unit KW using the current default format.

			FM_RANDOM_NUMBER(X) X is returned as a double precision random number, uniformly
			distributed on the open interval (0,1). It is a high-quality,
			long-period generator based on 49-digit prime numbers.
			Note that X is double precision, unlike the similar Fortran
			intrinsic random number routine, which returns a single-precision
			result. A default initial seed is used if FM_RANDOM_NUMBER is
			called without calling FM_RANDOM_SEED_PUT first.

			FM_RANDOM_SEED_GET(SEED) returns the seven integers SEED(1) through SEED(7) as the current
			seed for the FM_RANDOM_NUMBER generator.

			FM_RANDOM_SEED_PUT(SEED) initializes the FM_RANDOM_NUMBER generator using the seven integers
			SEED(1) through SEED(7). These get and put functions are slower
			than FM_RANDOM_NUMBER, so FM_RANDOM_NUMBER should be called many
			times between FM_RANDOM_SEED_PUT calls. Also, some generators that
			used a 9-digit modulus have failed randomness tests when used with
			only a few numbers being generated between calls to re-start with
			a new seed.

			FM_RANDOM_SEED_SIZE(SIZE) returns integer SIZE as the size of the SEED array used by the
			FM_RANDOM_NUMBER generator. Currently, SIZE = 7.

			FM_RATIONAL_POWER(MA,K,J,MB)
			MB = MA**(K/J) Rational power.
			Faster than MB = MA**(TO_FM(K)/J) for functions like the cube root.

			FM_READ(KREAD,MA) MA is returned after reading one (possibly multi-line) FM number
			on unit KREAD. This routine reads numbers written by FM_WRITE.

			FM_SET(NPREC) Set the internal FM variables so that the precision is at least
			NPREC base 10 digits plus three base 10 guard digits.

			FM_SETVAR(STRING) Define a new value for one of the internal FM variables in module
			FMVALS that controls one of the FM options. STRING has the form
			variable = value.
			Example: To change the screen width for FM output:
			CALL FM_SETVAR(' KSWIDE = 120 ')
			The variables that can be changed and the options they control are
			listed in sections 2 through 6 of the "FM.f95 Notes" section below.
			Only one variable can be set per call. The variable name in STRING
			must have no embedded blanks. The value part of STRING can be in
			any numerical format, except in the case of variable CMCHAR, which
			is character type. To set CMCHAR to 'E', don't use quotes in STRING:
			CALL FM_SETVAR(' CMCHAR = E ')

			FM_ULP(MA,MB) MB = One Unit in the Last Place of MA. For positive MA this is the
			same as the Fortran function SPACING, but MB < 0 if MA < 0.
			Examples: If MBASE = 10 and NDIG = 30, then ulp(1.0) =
			1.0E-29, ulp(-4.5E+67) = -1.0E+38.

			FM_VARS Write the current values of the internal FM variables on unit KW.

			FM_WRITE(KWRITE,MA) Write MA on unit KWRITE.
			Multi-line numbers will have '&' as the last nonblank character on
			all but the last line. These numbers can then be read easily using
			FM_READ.


			2. Type (IM). MA, MB, MC refer to type (IM) numbers.

			IM_DIVR(MA,MB,MC,MD) MC = int(MA/MB), MD = MA mod MB
			When both the quotient and remainder are needed, this routine
			is twice as fast as doing MC = MA/MB and MD = MOD(MA,MB)
			separately.

			IM_DVIR(MA,IVAL,MB,IREM) MB = int(MA/IVAL), IREM = MA mod IVAL
			IVAL and IREM are one word integers. Faster than doing separately.

			IM_FORM(FORM,MA,STRING) MA is converted to a character string using format FORM and
			returned in STRING. FORM can represent I, F, E, or ES formats.
			Example: CALL IMFORM('I70',MA,STRING)

			IM_FPRINT(FORM,MA) Print MA on unit KW using FORM format.

			IM_PRINT(MA) Print MA on unit KW using the current default format.

			IM_READ(KREAD,MA) MA is returned after reading one (possibly multi-line) IM number
			on unit KREAD. This routine reads numbers written by IM_WRITE.

			IM_WRITE(KWRITE,MA) Write MA on unit KWRITE. Multi-line numbers will have '&' as the
			last nonblank character on all but the last line.
			These numbers can then be read easily using IM_READ.


			3. Type (ZM). MA, MB, MC refer to type (ZM) numbers. MBFM is type (FM).

			ZM_ARG(MA,MBFM) MBFM = complex argument of MA. MBFM is the (real) angle in the
			interval ( -pi , pi ] from the positive real axis to the
			point (x,y) when MA = x + y*i.

			ZM_COS_SIN(MA,MB,MC) MB = COS(MA), MC = SIN(MA).
			Faster than 2 calls.

			ZM_COSH_SINH(MA,MB,MC) MB = COSH(MA), MC = SINH(MA).
			Faster than 2 calls.

			ZM_FORM(FORM1,FORM2,MA,STRING)
			STRING = MA
			MA is converted to a character string using format FORM1 for the
			real part and FORM2 for the imaginary part. The result is returned
			in STRING. FORM1 and FORM2 can represent I, F, E, or ES formats.
			Example:
			CALL ZMFORM('F20.10','F15.10',MA,STRING)

			ZM_FPRINT(FORM1,FORM2,MA) Print MA on unit KW using formats FORM1 and FORM2.

			ZM_PRINT(MA) Print MA on unit KW using the current default format.

			ZM_READ(KREAD,MA) MA is returned after reading one (possibly multi-line) ZM number
			on unit KREAD. This routine reads numbers written by ZMWRITE.

			ZM_RATIONAL_POWER(MA,IVAL,JVAL,MB)
			MB = MA ** (IVAL/JVAL)
			Faster than MB = MA**(TO_FM(K)/J) for functions like the cube root.

			ZM_WRITE(KWRITE,MA) Write MA on unit KWRITE. Multi-line numbers are formatted for
			automatic reading with ZMREAD.


			-------------------------------- More detail about functions -----------------------------------


			Some other functions are defined that do not correspond to machine precision intrinsic
			functions. These include formatting functions, integer modular functions and GCD, and some
			mathematical special functions.
			N, K below are machine precision integers, J1, J2, J3 are TYPE (IM), FMT, FMTR, FMTI are
			character strings, A,B,X are TYPE (FM), and Z is TYPE (ZM).
			The three formatting functions return a character string containing the formatted number, the
			three TYPE (IM) functions return a TYPE (IM) result, and the 12 special functions return
			TYPE (FM) results.


			Formatting functions:

			FM_FORMAT(FMT,A) Put A into FMT (real) format
			IM_FORMAT(FMT,J1) Put J1 into FMT (integer) format
			ZM_FORMAT(FMTR,FMTI,Z) Put Z into (complex) format, FMTR for the real
			part and FMTI for the imaginary part

			Examples:
			ST = FM_FORMAT('F65.60',A)
			WRITE (,) ' A = ',TRIM(ST)
			ST = FM_FORMAT('E75.60',B)
			WRITE (,) ' B = ',ST(1:75)
			ST = IM_FORMAT('I50',J1)
			WRITE (,) ' J1 = ',ST(1:50)
			ST = ZM_FORMAT('F35.30','F30.25',Z)
			WRITE (,) ' Z = ',ST(1:70)

			These functions are used for one-line output. The returned character strings are of
			length 200.

			For higher precision numbers, the output can be broken onto multiple lines automatically by
			calling subroutines FM_PRINT, IM_PRINT, ZM_PRINT, or the line breaks can be done by hand after
			calling one of the subroutines FM_FORM, IM_FORM, ZM_FORM.

			For ZM_FORMAT the length of the output is 5 more than the sum of the two field widths.

			For more examples using the various formatting subroutines and functions, see the program
			SampleFM.f95.


			Integer functions:

			BINOMIAL(N,K) Binomial coefficient N choose K. Returns the exact result as a
			type IM value.
			BINOMIAL(J1,J2) Binomial coefficient J1 choose J2. Like factorial below, the result
			might be too large unless min(J2,J1-J2) is fairly small,
			FACTORIAL(N) N! Returns the exact result as a type IM value.
			FACTORIAL(J1) J1! Note that the factorial function grows so rapidly that if type IM
			variable J1 is larger than the largest machine precision integer,
			then J1! has over 10 billion digits and the calculation would
			likely fail due to memory or time constraints. This version is
			provided for convenience, and will return UNKNOWN if J1 cannot
			be represented as a machine precision integer.
			GCD(J1,J2) Greatest Common Divisor of J1 and J2.
			MULTIPLY_MOD(J1,J2,J3) J1 * J2 mod J3
			POWER_MOD(J1,J2,J3) J1 ** J2 mod J3


			Special functions:

			BERNOULLI(N) Nth Bernoulli number
			BESSEL_JN(N,X) Bessel function of the first kind J_n(x)
			BESSEL_J0(X) J_0(x)
			BESSEL_J1(X) J_1(x)
			BESSEL_JN(N1,N2,X) Returns the array (/ J_n1(x) , ..., J_n2(x) /)
			BESSEL_YN(N,X) Bessel function of the second kind Y_n(x)
			BESSEL_Y0(X) Y_0(x)
			BESSEL_Y1(X) Y_1(x)
			BESSEL_YN(N1,N2,X) Returns the array (/ Y_n1(x) , ..., Y_n2(x) /)
			BETA(A,B) Integral (0 to 1) t*(a-1) (1-t)**(b-1) dt
			BINOMIAL(A,B) Binomial Coefficient a! / ( b! (a-b)! )
			COS_INTEGRAL(X) Cosine Integral Ci(x)
			COSH_INTEGRAL(X) Hyperbolic Cosine Integral Chi(x)
			ERF(X) Error function Erf(x)
			ERFC(X) Complimentary error function Erfc(x)
			ERFC_SCALED(X) Exp(x*2) Erfc(x)
			EXP_INTEGRAL_EI(X) Exponential Integral Ei(x)
			EXP_INTEGRAL_EN(N,X) Exponential Integral E_n(x)
			FACTORIAL(X) x! = Gamma(x+1)
			FRESNEL_C(X) Fresnel Cosine Integral C(x)
			FRESNEL_S(X) Fresnel Sine Integral S(x)
			GAMMA(X) Integral (0 to infinity) t*(x-1) exp(-t) dt
			INCOMPLETE_BETA(X,A,B) Integral (0 to x) t*(a-1) (1-t)**(b-1) dt
			INCOMPLETE_GAMMA1(A,X) Integral (0 to x) t*(a-1) exp(-t) dt
			INCOMPLETE_GAMMA2(A,X) Integral (x to infinity) t*(a-1) exp(-t) dt
			LOG_ERFC(X) Ln( Erfc(x) )
			LOG_GAMMA(X) Ln( Gamma(x) )
			LOG_INTEGRAL(X) Logarithmic Integral Li(x)
			POLYGAMMA(N,X) Nth derivative of Psi(x)
			POCHHAMMER(X,N) x(x+1)(x+2)...(x+n-1)
			PSI(X) Derivative of Ln(Gamma(x))
			SIN_INTEGRAL(X) Sine Integral Si(x)
			SINH_INTEGRAL(X) Hyperbolic Sine Integral Shi(x)


			Array operations:

			Arithmetic operations and functions on arrays of dimension (rank) one or two are supported for
			each of the three multiple-precision types. Binary operations (+-*/) require both arguments to
			have the same rank and shape.

			Examples:
			TYPE (FM), SAVE, DIMENSION(10) :: A, B
			TYPE (FM), SAVE, DIMENSION(3,3) :: C
			TYPE (IM), SAVE, DIMENSION(10) :: J, K
			TYPE (IM), SAVE, DIMENSION(3,3) :: L
			...
			A = 0 ! Set the whole array to zero
			J = J * K ! Set J(i) = J(i) * K(i) for i = 1, ..., 10
			B = A - K ! Mixed-mode operations are ok
			C = 7.3D0 * C - ( C + 2*L )/3

			Array functions:

			DOT_PRODUCT(X,Y) Dot product of rank 1 vectors of the same type.
			Note that when X and Y are complex, the result is not just
			the sum of the products of the corresponding array elements,
			as it is for types FM and IM. For ZM the formula is the sum of
			conjg(X(j)) * Y(j).
			IS_OVERFLOW(X) Returns true if any element is + or - overflow.
			IS_UNDERFLOW(X) Returns true if any element is + or - underflow.
			IS_UNKNOWN(X) Returns true if any element is unknown.
			MATMUL(X,Y) Matrix multiplication of arrays of the same type
			Cases for valid argument shapes:
			(1) (n,m) * (m,k) --> (n,k)
			(2) (m) * (m,k) --> (k)
			(3) (n,m) * (m) --> (n)
			MAXLOC(X) Location of the maximum value in the array
			MAXVAL(X) Maximum value in the array
			MINLOC(X) Location of the minimum value in the array
			MINVAL(X) Minimum value in the array
			PRODUCT(X) Product of all values in the array
			SUM(X) Sum of all values in the array
			TRANSPOSE(X) Matrix transposition. If X is a rank 2 array with shape (n,m), then
			Y = TRANSPOSE(X) has shape (m,n) with Y(i,j) = X(j,i).
			TO_FM(X) Rank 1 or 2 arrays are converted to similar type (fm) arrays.
			TO_IM(X) Rank 1 or 2 arrays are converted to similar type (im) arrays.
			TO_ZM(X) Rank 1 or 2 arrays are converted to similar type (zm) arrays.
			TO_INT(X) Rank 1 or 2 arrays are converted to similar integer arrays.
			TO_SP(X) Rank 1 or 2 arrays are converted to similar single precision arrays.
			TO_DP(X) Rank 1 or 2 arrays are converted to similar double precision arrays.
			TO_SPZ(X) Rank 1 or 2 arrays are converted to similar single complex arrays.
			TO_DPZ(X) Rank 1 or 2 arrays are converted to similar double complex arrays.

			The arithmetic array functions DOT_PRODUCT, MATMUL, PRODUCT, and SUM work like the other
			functions in the FM package in that they raise precision and compute the sums and/or products
			at the higher precision, then round the final result back to the user's precision to provide
			a more accurate result.

			Fortran's optional [,mask] argument for these functions is not provided.

			Many of the 1-argument functions can be used with array arguments, with the result being an
			array of the same size and shape where the function has been applied to each element.

			Examples:
			TYPE (FM), SAVE, DIMENSION(10) :: A, B, C
			...
			A = ABS(B) ! Set A(i) = ABS(B(i)) for i = 1, ..., 10
			C = SQRT(A+4+BB) ! Set C(i) = SQRT(A(i)+4+B(i)B(i)) for i = 1, ..., 10

			Functions that can have array arguments. As above, "real", "integer", and "complex" refer
			to types FM, IM, and ZM respectively.

			ABS real integer complex
			ACOS real complex
			ACOSH real complex
			AIMAG complex
			AINT real complex
			ANINT real complex
			ASIN real complex
			ASINH real complex
			ATAN real complex
			ATANH real complex
			CEILING real integer complex
			CONJG complex
			COS real complex
			COSH real complex
			EXP real complex
			FLOOR real integer complex
			FRACTION real complex
			INT real integer complex
			LOG real complex
			LOG10 real complex
			NINT real integer complex
			SIN real complex
			SINH real complex
			SQRT real complex
			TAN real complex
			TANH real complex
			COS_INTEGRAL real
			COSH_INTEGRAL real
			ERF real
			ERFC real
			ERFC_SCALED real
			EXP_INTEGRAL_EI real
			FACTORIAL real integer machine-precision integer
			FRESNEL_C real
			FRESNEL_S real
			GAMMA real
			LOG_ERFC real
			LOG_GAMMA real
			LOG_INTEGRAL real
			PSI real
			SIN_INTEGRAL real
			SINH_INTEGRAL real





			----------------------------------------------------------------------------------------------------
			---------------------------------------- FM.f95 Notes ------------------------------------------



			The routines in the file FM.f95 perform the low-level multiple precision arithmetic and functions
			on three kinds of numbers.

			FM routines handle floating-point real multiple precision numbers,
			IM routines handle integer multiple precision numbers, and
			ZM routines handle floating-point complex multiple precision numbers.

			References to FM numbers below mean the low-level array form of the number used by the routines
			in FM.f95, and not the derived type (fm) numbers handled by the FMZM module. Logically, both
			may refer to the same multiple precision number, but the syntax for dealing with the two types
			of objects is different. The same is true of references to IM numbers and ZM numbers below.

			These are the basic routines for the FM package, and the expectation is that the user will not
			call these routines directly. The typical usage is for a program to declare multiple precision
			variables with the three derived types defined in module FMZM in file FMZM90.f95. Then that
			module provides the interface between the user's program and the routines in this file. See the
			documentation in the sections above for advice on using the FMZM module. The information below
			is intended as a technical reference on the inner workings of FM, and most FM users should not
			need to study it.


			1. INITIALIZING THE PACKAGE

			The variables that contain values to be shared by the different routines are located in module
			FMVALS in file FMSAVE.f95. Variables that are described below for controlling various features
			of the FM package are found in this module. They are initialized to default values assuming
			32-bit integers and 64-bit double precision representation of the arrays holding multiple
			precision numbers. The base and number of digits to be used are initialized to give slightly
			more than 50 decimal digits. Subroutine FMVARS can be used to get a list of these variables
			and their values.

			The intent of module FMVALS is to hide the FM internal variables from the user's program, so that
			no name conflicts can occur. Subroutine FM_SETVAR can be used to change the variables listed
			below to new values. It is not always safe to try to change these variables directly by putting
			USE FMVALS into the calling program and then changing them by hand. Some of the saved constants
			depend upon others, so that changing one variable may cause errors if others depending on that
			one are not also changed. FM_SETVAR automatically updates any others that depend upon the one
			being changed.

			Subroutine FMSET also initializes these variables. It tries to compute the best value for each,
			and it checks several of the default values set in FMVALS to see that they are reasonable for a
			given machine. FMSET can also be called to set or change the current precision level for the
			multiple precision numbers.

			Calling FMSET is optional starting in version 1.2 of the FM package. In previous versions one
			call was required before any other routine in the package could be used.

			The routine ZMSET from version 1.1 is no longer needed, and the complex operations are
			automatically initialized in FMVALS. It has been left in the package for compatibility with
			version 1.1.


			2. REPRESENTATION OF FM NUMBERS

			MBASE is the base in which the arithmetic is done. MBASE must be bigger than one, and less than
			or equal to the square root of the largest representable integer. For best efficiency
			MBASE should be large, but no more than about 1/4 of the square root of the largest
			representable integer. Input and output conversions are much faster when MBASE is a
			power of ten.

			NDIG is the number of base MBASE digits that are carried in the multiple precision numbers.
			NDIG must be at least two. The upper limit for NDIG is restricted only by the amount of
			memory available.

			Sometimes it is useful to dynamically vary NDIG during the program. Routine FMEQU should be used
			to round numbers to lower precision or zero-pad them to higher precision when changing NDIG.

			The default value of MBASE is a large power of ten. FMSET also sets MBASE to a large power of
			ten. For an application where another base is used, such as simulating a given machine's base
			two arithmetic, use subroutine FM_SETVAR (equivalent to routine FMSETVAR in FM.f95) to change
			MBASE, so that the other internal values depending on MBASE will be changed accordingly.

			There are two representations for a floating point multiple precision number. The unpacked
			representation used by the routines while doing the computations is base MBASE and is stored
			in NDIG+3 words. A packed representation is available to store the numbers in compressed form.
			In this format, the NDIG (base MBASE) digits of the mantissa are packed two per word to conserve
			storage. Thus the external, packed form of a number requires (NDIG+1)/2+3 words.

			This version uses double precision arrays to hold the numbers. Version 1.0 of FM used integer
			arrays, which are faster on some machines. The package can be changed to use integer arrays ---
			see section 10 on EFFICIENCY below.

			The unpacked format of a floating multiple precision number is as follows. A number MA refers
			to elements of an array with the multiple precision number stored as follows:
			1 Sign of the number
			2 Accuracy
			3 Exponent of the number
			4 First digit of the number
			...
			NDIG+3 Last digit of the number.

			The accuracy is the approximate number of bits of precision of the number. This precision value
			is intended to be used by FM functions that need to monitor cancellation error in addition and
			subtraction. The cancellation monitor code is usually disabled for user calls, and FM functions
			only check for cancellation when they must. Tracking cancellation causes most routines to run
			slower, with addition and subtraction being affected the most.

			The exponent is a power of MBASE and the implied radix point is immediately before the first
			digit of the mantissa. The exponent is a signed integer. The overflow threshold is
			MBASE(MXEXP+1), and the underflow threshold is MBASE(-MXEXP-1). This means the valid
			exponents for an FM number can range from -MXEXP to MXEXP+1 (inclusive).
			Every nonzero number is normalized so that the first digit of the mantissa is nonzero.

			For MBASE = 10,000 and NDIG = 4, if MA is the number -pi, it would have these representations:

			Word 1 2 3 4 5 6 7

			Unpacked: -1 42 1 3 1415 9265 3590
			Packed: -1 42 1 31415 92653590

			The mantissa has about 42 bits of precision, and the number represented is
			(-1)(100001)(.0003141592653590).

			Because of the normalization of the digits with a large base, the equivalent number of base 10
			significant digits for an FM number may be as small as LOG10(MBASE)*(NDIG-1) + 1. In the -pi
			example above, this is 4*3 + 1 = 13.

			In version 1.2 and before, each variable like MA was a fixed-size array. Now the array values
			for all multiple precision numbers are stored together in one dynamic array, MWK. A variable
			name like MA is just an integer index giving the location within MWK where the digits are put.
			This allows the memory management of the package to be more flexible than before, to handle
			large arrays of multiple precision numbers at low precision (30 to 50 s.d.) as well as some
			at high precision (millions of s.d.).

			The integer routines use the FM format to represent numbers, without the number of digits (NDIG)
			being fixed. Integers in IM format are essentially variable precision, using the minimum number
			of words to represent each value.

			The unpacked format is the default. As machines' memories have gotten bigger, few applications
			need the packed format. A program that uses packed format numbers should not use the FMZM module
			or the multiple precision derived types defined in FMZM. Packed numbers are treated as temporary
			values by the routines in FMZM, so mixing packed numbers with derived type operations from FMZM
			means the packed numbers could be deleted before the main program is finished with them.

			For programs using both FM and IM numbers, FM routines should not be called with IM numbers, and
			IM routines should not be called with FM numbers, since the implied value of NDIG used for an IM
			number may not match the explicit NDIG expected by an FM routine. Use the conversion routines
			IMFM2I and IMI2FM to change between the FM and IM formats.

			The format for complex FM numbers (called ZM numbers below) is very similar to that for real FM
			numbers. Each ZM number consists of two FM numbers representing the real and imaginary parts of
			a complex number. If MA is a ZM number, then the real part is MA(1) and the imaginary part is
			MA(2). As with FM, there are packed and unpacked formats for the numbers.


			3. INPUT/OUTPUT ROUTINES

			All versions of the input routines perform free-format conversion from characters to FM numbers.

			a. Conversion to or from a character array

			FMINP converts from a character(1) array to an FM number.

			FMOUT converts an FM number to base 10 and formats it for output as an array of type
			character(1). The output is left justified in the array, and the format is defined
			by two variables in module FMVALS, so that a separate format definition does not have
			to be provided for each output call.

			Variables JFORM1 and JFORM2 define a default output format.

			JFORM1 = 0 E format ( .314159M+6 )
			= 1 ES format ( 3.14159M+5 )
			= 2 F format ( 314159.000 )

			JFORM2 is the number of significant digits to display (if JFORM1 = 0 or 1).
			If JFORM2 = 0 then a default number of digits is chosen. The default is roughly
			the full precision of the number.
			JFORM2 is the number of digits after the decimal point (if JFORM1 = 2).
			See the FMOUT documentation for more details.

			b. Conversion to or from a character string

			FMST2M converts from a character string to an FM number.

			FMFORM converts an FM number to a character string according to a format provided in each
			call. The format description is more like that of a Fortran FORMAT statement, and
			integer or fixed-point output is right justified.

			c. Direct read or write

			FMPRINT uses FMOUT to print one FM number.

			FMFPRINT uses FMFORM to print one FM number.

			FMWRITE writes FM numbers for later input using FMREAD.

			FMREAD reads FM numbers written by FMWRITE.

			The values given to JFORM1 and JFORM2 can be used to define a default output format when FMOUT
			or FMPRINT are called. The explicit format used in a call to FMFORM or FMFPRINT overrides the
			settings of JFORM1 and JFORM2.

			KW is the unit number to be used for standard output from the package, including error and
			warning messages, and trace output.

			For multiple precision integers, the corresponding routines IMINP, IMOUT, IMST2M, IMFORM,
			IMPRINT, IMFPRINT, IMWRITE, and IMREAD provide similar input and output conversions. For output
			of IM numbers, JFORM1 and JFORM2 are ignored and integer format (JFORM1=2, JFORM2=0) is used.

			For ZM numbers, the corresponding routines ZMINP, ZMOUT, ZMST2M, ZMFORM, ZMPRINT, ZMFPRINT,
			ZMWRITE, and ZMREAD provide similar input and output conversions.

			For the output format of ZM numbers, JFORM1 and JFORM2 determine the default format for the
			individual parts of a complex number as with FM numbers.

			JFORMZ determines the combined output format of the real and imaginary parts.

			JFORMZ = 1 normal setting : 1.23 - 4.56 i
			= 2 use capital I : 1.23 - 4.56 I
			= 3 parenthesis format: ( 1.23 , -4.56 )

			JPRNTZ controls whether to print real and imaginary parts on one line whenever possible.

			JPRNTZ = 1 print both parts as a single string :
			1.23456789M+321 - 9.87654321M-123 i
			= 2 print on separate lines without the 'i' :
			1.23456789M+321
			-9.87654321M-123

			For further description of these routines, see section 9 below.


			4. ARITHMETIC TRACING

			Variables NTRACE and LVLTRC control trace printout from the package.

			NTRACE = 0 No output except warnings and errors. (Default)
			= 1 The result of each call to one of the routines is printed in base 10, using FMOUT.
			= -1 The result of each call to one of the routines is printed in internal base MBASE
			format.
			= 2 The input arguments and result of each call to one of the routines is printed in
			base 10, using FMOUT.
			= -2 The input arguments and result of each call to one of the routines is printed in
			base MBASE format.

			LVLTRC defines the call level to which the trace is done. LVLTRC = 1 means only FM routines
			called directly by the user are traced, LVLTRC = 2 also prints traces for FM routines
			called by other FM routines called directly by the user, etc. Default is 1.

			In the above description, internal MBASE format means the number is printed as it appears in the
			array --- the sign, accuracy, exponent, then the NDIG base MBASE digits.


			5. ERROR CONDITIONS

			KFLAG is a condition value returned by the package after each call to one of the routines.
			Negative values indicate conditions for which a warning message will be printed unless
			KWARN = 0.
			Positive values indicate conditions that may be of interest but are not errors. No warning
			message is printed if KFLAG is nonnegative.

			Subroutine FMFLAG is provided to give the user access to the current condition code. For
			example, to set the user's local variable LFLAG to FM's internal KFLAG value:
			CALL FMFLAG(LFLAG)

			KFLAG = 0 Normal operation.

			= 1 One of the operands in FMADD or FMSUB was insignificant with respect to the
			other. This means that in the default (symmetric) rounding mode the result
			is equal to the argument of larger magnitude. KFLAG = 1 is still returned
			with the other three rounding modes (see KROUND below), but the result may
			not be equal to either input argument.
			= 2 In converting an FM number to a one word integer in FMM2I, the FM number was
			not exactly an integer. The next integer toward zero was returned.

			= -1 NDIG was less than 2.
			= -2 MBASE was less than 2 or more than MXBASE.
			= -3 An exponent was out of range.
			= -4 Invalid input argument(s) to an FM routine. UNKNOWN was returned.
			= -5 + or - OVERFLOW was generated as a result from an FM routine.
			= -6 + or - UNDERFLOW was generated as a result from an FM routine.
			= -7 The input string (array) to FMINP was not legal.
			= -8 The character array was not large enough in an input or output routine.
			= -9 Precision could not be raised enough to provide all requested guard digits.
			This means the program has run out of memory.
			UNKNOWN was returned.
			= -10 An FM input argument was too small in magnitude to convert to the machine's
			single or double precision in FMM2SP or FMM2DP. Check that the definitions
			of SPMAX and DPMAX in file FMSAVE.f95 are correct for the current machine.
			Zero was returned.
			= -11 Array MBERN is not dimensioned large enough for the requested number of
			Bernoulli numbers.
			= -12 Array MJSUMS is not dimensioned large enough for the number of coefficients
			needed in the reflection formula in FMPGAM.

			When a negative KFLAG condition is encountered, the value of KWARN determines the action to
			be taken.

			KWARN = 0 Execution continues and no message is printed.
			= 1 A warning message is printed and execution continues.
			= 2 A warning message is printed and execution stops.

			The default setting is KWARN = 1.

			When an overflow or underflow is generated for an operation in which an input argument was
			already an overflow or underflow, no additional message is printed. When an unknown result
			is generated and an input argument was already unknown, no additional message is printed.
			In these cases the negative KFLAG value is still returned.

			IM routines handle exceptions like OVERFLOW or UNKNOWN in the same way as FM routines, but there
			are some differences because the number of digits carried for IM numbers is not fixed. For
			example, in computing the product of two large integers FM will try to allocate more space rather
			than returning +OVERFLOW. If this allocation fails, FM will write an error message indicating it
			could not get more memory, and the program will stop. The routines IMMPY_MOD and IMPOWER_MOD can
			be used to obtain modular products and powers without as much chance of running out of memory.


			6. OTHER OPTIONS

			KRAD = 0 All angles in the real trigonometric functions and inverse functions are measured
			in degrees.
			= 1 All angles are measured in radians. (Default)

			KROUND = -1 All results are rounded toward minus infinity.
			= 0 All results are rounded toward zero (chopped).
			= 1 All results are rounded to the nearest FM number, or to the value with an even last
			digit if the result is exactly halfway between two FM numbers. (Default)
			= 2 All results are rounded toward plus infinity.

			KSWIDE defines the maximum screen width to be used for all unit KW output. Default is 80.

			KESWCH controls the action taken in FMINP and other input routines for strings like 'E7' that
			have no digits before the exponent field. This is sometimes a convenient abbreviation
			when doing interactive keyboard input.
			KESWCH = 1 causes 'E7' to translate like '1.0E+7'. (Default)
			KESWCH = 0 causes 'E7' to translate like '0.0E+7' and give 0.

			CMCHAR defines the exponent letter to be used for FM variable output.
			Default is 'M', as in 1.2345M+678.
			Change it to 'E' for output to be read by a non-FM program.

			KDEBUG = 0 No error checking is done to see if input arguments are valid and parameters like
			NDIG and MBASE are correct upon entry to each routine. (Default)
			= 1 Some error checking is done. (Slower speed)

			See module FMVALS in file FMSAVE.f95 for additional description of these and other variables
			defining various FM conditions.


			7. ARRAY DIMENSIONS

			Before version 1.3 the multiple-precision numbers were each stored in a separate array. Now all
			these arrays have been combined into one, MWK. Each multiple-precision number is represented by
			a single integer that is used to access the number in the MWK array. This makes the memory usage
			within the package much more flexible.


			8. PORTABILITY

			In FMSET several variables are set to machine-dependent values, and many of the variables
			initialized in module FMVALS in file FMSAVE.f95 are checked to see that they have reasonable
			values. FMSET will print warning messages on unit KW for any of the FMVALS variables that
			seem to be poorly initialized.

			If an FM run fails, call FMVARS to get a list of all the FMVALS variables printed on unit KW.
			Setting KDEBUG = 1 at the start may also identify some errors.

			In the routines for special functions, several constants are used that require the machine's
			integer word size to be at least 32 bits.


			9. LIST OF ROUTINES

			First are the routines that deal with multiple precision real numbers. All of these are
			subroutines except logical function FMCOMPARE.

			MA, MB, MC refer to FM format numbers (i.e., integers as opposed to the derived types with
			integer components that are defined in file FMZM90.f95)

			In Fortran-90 and later versions of the Fortran standard, it is potentially unsafe to use the
			same variable more than once in the calling sequence. The operation MA = MA + MB should not be
			written as
			CALL FMADD(MA,MB,MA)
			since the code for the subroutine will not know that the first and third arguments are the same,
			and some code optimizations under the assumption that all three arguments are different could
			cause errors.

			One solution is to use a third array and then put the result back in MA:
			CALL FMADD(MA,MB,MC)
			CALL FMEQ(MC,MA)

			When the first call is doing one of the "fast" operations like addition, the extra call to move
			the result back to MA can cause a noticeable loss in efficiency. To avoid this, separate
			routines are provided for the basic arithmetic operations when the result is to be returned in
			the same array as one of the inputs.

			A routine name with a suffix of "_R1" returns the result in the first input array, and a suffix
			of "_R2" returns the result in the second input array. The example above would then be:
			CALL FMADD_R1(MA,MB)

			These routines each have one less argument than the original version, since the output is
			re-directed to one of the inputs. The result array should not be the same as any input array
			when the original version of the routine is used.

			The routines that can be used this way are listed below. For others, like
			CALL FMEXP(MA,MA)
			the relative cost of doing an extra copy is small. This one should become
			CALL FMEXP(MA,MB)
			CALL FMEQ(MB,MA)

			When the derived-type interface is used, as in
			TYPE (FM), SAVE :: A, B
			...
			A = A + B
			there is no problem putting the result back into A, since the interface routine creates a
			temporary scratch array for the result of A + B.

			For each of these routines there is also a version available for which the argument list is
			the same but all FM numbers are in packed format. The routines using packed numbers have the
			same names except 'FM' is replaced by 'FP' at the start of each name.

			Some of the routine names were restricted to 6 characters in earlier versions of FM. The old
			names have been retained for compatibility, but new names that are longer and more readable
			have been added. For example, the old routine FMCSSN can now also be called as FMCOS_SIN.
			Both old and new names are listed below.


			FMABS(MA,MB) MB = ABS(MA)

			FMACOS(MA,MB) MB = ACOS(MA)

			FMACOSH(MA,MB) MB = ACOSH(MA)

			FMADD(MA,MB,MC) MC = MA + MB

			FMADD_R1(MA,MB) MA = MA + MB

			FMADD_R2(MA,MB) MB = MA + MB

			FMADDI(MA,IVAL) MA = MA + IVAL Increment an FM number by a one word integer.
			Note this call does not have an "MB" result
			like FMDIVI and FMMPYI.

			FMASIN(MA,MB) MB = ASIN(MA)

			FMASINH(MA,MB) MB = ASINH(MA)

			FMATAN(MA,MB) MB = ATAN(MA)

			FMATANH(MA,MB) MB = ATANH(MA)

			FMATAN2(MA,MB,MC) MC = ATAN2(MA,MB) < old name: FMATN2 >

			FMBIG(MA) MA = Biggest FM number less than overflow.

			FMCHANGEBASE(MA,MB,NEW_MBASE,NEW_NDIG)
			MB is returned with the base NEW_MBASE and precision NEW_NDIG
			representation of MA, where MA is given in the current base (MBASE)
			and precision (NDIG). This routine is primarily meant to be used
			for input and output conversion when a base is being used that is
			not a power of ten.

			FMCOMPARE(MA,LREL,MB) Logical comparison of MA and MB. < old name: FMCOMP >
			LREL is a character(2) value identifying which of the six comparisons
			is to be made.
			Example: IF (FMCOMPARE(MA,'>=',MB)) ...
			Also can be: IF (FMCOMPARE(MA,'GE',MB)) ...
			character(1) is ok: IF (FMCOMPARE(MA,'>',MB)) ...

			FMCONS Set several saved constants that depend on MBASE, the base being used.
			FMCONS should be called immediately after changing MBASE.

			FMCOS(MA,MB) MB = COS(MA)

			FMCOS_SIN(MA,MB,MC) MB = COS(MA), MC = SIN(MA). < old name: FMCSSN >
			Faster than making two separate calls.

			FMCOSH(MA,MB) MB = COSH(MA)

			FMCOSH_SINH(MA,MB,MC) MB = COSH(MA), MC = SINH(MA). < old name: FMCHSH >
			Faster than making two separate calls.

			FMDIG(NSTACK,KST) Find a set of precisions to use during Newton iteration for finding a
			simple root starting with about double precision accuracy.

			FMDIM(MA,MB,MC) MC = DIM(MA,MB)

			FMDIV(MA,MB,MC) MC = MA / MB

			FMDIV_R1(MA,MB) MA = MA / MB

			FMDIV_R2(MA,MB) MB = MA / MB

			FMDIVI(MA,IVAL,MB) MB = MA/IVAL IVAL is a one word integer.

			FMDIVI_R1(MA,IVAL) MA = MA/IVAL

			FMDP2M(X,MA) MA = X Convert from double precision to FM.

			FMDPM(X,MA) MA = X Convert from double precision to FM.
			Faster than FMDP2M, but MA agrees with X only to D.P.
			accuracy. See the comments in the two routines.

			FMEQ(MA,MB) MB = MA Both have precision NDIG.
			This is the version to use for standard B = A statements.

			FMEQU(MA,MB,NA,NB) MB = MA Version for changing precision.
			MA has NA digits (i.e., MA was computed using NDIG = NA), and
			MB will be defined having NB digits.
			MB is rounded if NB < NA
			MB is zero-padded if NB > NA

			FMEXP(MA,MB) MB = EXP(MA)

			FMFLAG(K) K = KFLAG get the value of the FM condition flag -- stored in the
			internal FM variable KFLAG in module FMVALS.

			FMFORM(FORM,MA,STRING) MA is converted to a character string using format FORM and returned in
			STRING. FORM can represent I, F, E, or ES formats. Example:
			CALL FMFORM('F60.40',MA,STRING)

			FMFPRINT(FORM,MA) Print MA on unit KW using FORM format. < old name: FMFPRT >

			FMHYPOT(MA,MB,MC) MA = SQRT(MA2 + MB2)

			FMI2M(IVAL,MA) MA = IVAL Convert from one word integer to FM.

			FMINP(LINE,MA,LA,LB) MA = LINE Input conversion.
			Convert LINE(LA) through LINE(LB) from characters to FM.

			FMINT(MA,MB) MB = INT(MA) Integer part of MA.

			FMIPOWER(MA,IVAL,MB) MB = MA**IVAL Raise an FM number to a one word integer power.
			< old name: FMIPWR >

			FMLOG10(MA,MB) MB = LOG10(MA) < old name: FMLG10 >

			FMLN(MA,MB) MB = LOG(MA)

			FMLNI(IVAL,MA) MA = LOG(IVAL) Natural log of a one word integer.

			FMM2DP(MA,X) X = MA Convert from FM to double precision.

			FMM2I(MA,IVAL) IVAL = MA Convert from FM to integer.

			FMM2SP(MA,X) X = MA Convert from FM to single precision.

			FMMAX(MA,MB,MC) MC = MAX(MA,MB)

			FMMIN(MA,MB,MC) MC = MIN(MA,MB)

			FMMOD(MA,MB,MC) MC = MA mod MB

			FMMPY(MA,MB,MC) MC = MA * MB

			FMMPY_R1(MA,MB) MA = MA * MB

			FMMPY_R2(MA,MB) MB = MA * MB

			FMMPYI(MA,IVAL,MB) MB = MA*IVAL Multiply by a one word integer.

			FMMPYI_R1(MA,IVAL) MA = MA*IVAL

			FMNINT(MA,MB) MB = NINT(MA) Nearest FM integer.

			FMNORM2(MA,N,MB) MB = SQRT( MA(1)2 + MA(2)2 + ... + MA(N)**2 )

			FMOUT(MA,LINE,LB) LINE = MA Convert from FM to character.
			LINE is a character array of length LB.

			FMPI(MA) MA = pi

			FMPRINT(MA) Print MA on unit KW using current format. < old name: FMPRNT >

			FMPOWER(MA,MB,MC) MC = MA**MB < old name: FMPWR >

			FM_RANDOM_NUMBER(X) X is returned as a double precision random number, uniformly
			distributed on the open interval (0,1). It is a high-quality,
			long-period generator based on 49-digit prime numbers.
			Note that X is double precision, unlike the similar Fortran intrinsic
			random number routine, which returns a single-precision result.
			A default initial seed is used if FM_RANDOM_NUMBER is called without
			calling FM_RANDOM_SEED_PUT first. See the comments in section 11 below
			and also those in the routine for more details.

			FM_RANDOM_SEED_GET(SEED) returns the seven integers SEED(1) through SEED(7) as the current seed
			for the FM_RANDOM_NUMBER generator.

			FM_RANDOM_SEED_PUT(SEED) initializes the FM_RANDOM_NUMBER generator using the seven integers
			SEED(1) through SEED(7). These get and put functions are slower than
			FM_RANDOM_NUMBER, so FM_RANDOM_NUMBER should be called many times
			between FM_RANDOM_SEED_PUT calls. Also, some generators that used a
			9-digit modulus have failed randomness tests when used with only a few
			numbers being generated between calls to re-start with a new seed.

			FM_RANDOM_SEED_SIZE(SIZE) returns integer SIZE as the size of the SEED array used by the
			FM_RANDOM_NUMBER generator. Currently, SIZE = 7.

			FMRATIONAL_POWER(MA,K,J,MB)
			MB = MA**(K/J) Rational power. < old name: FMRPWR >
			Faster than FMPOWER for functions like the cube root.

			FMREAD(KREAD,MA) MA is returned after reading one (possibly multi-line) FM number
			on unit KREAD. This routine reads numbers written by FMWRITE.

			FMSET(NPREC) Set the internal FM variables so that the precision is at least NPREC
			base 10 digits plus three base 10 guard digits.

			FMSETVAR(STRING) Define a new value for one of the internal FM variables in module
			FMVALS that controls one of the FM options. STRING has the form
			variable = value.
			Example: To change the screen width for FM output:
			CALL FMSETVAR(' KSWIDE = 120 ')
			The variables that can be changed and the options they control are
			listed in sections 2 through 6 above. Only one variable can be set
			per call. The variable name in STRING must have no embedded blanks.
			The value part of STRING can be in any numerical format, except in
			the case of variable CMCHAR, which is character type. To set CMCHAR
			to 'E', don't use any quotes in STRING:
			CALL FMSETVAR(' CMCHAR = E ')

			FMSIGN(MA,MB,MC) MC = SIGN(MA,MB) Returns the absolute value of MA times the sign
			of MB.

			FMSIN(MA,MB) MB = SIN(MA)

			FMSINH(MA,MB) MB = SINH(MA)

			FMSP2M(X,MA) MA = X Convert from single precision to FM.

			FMSQR(MA,MB) MB = MA * MA Faster than FMMPY.

			FMSQR_R1(MA) MA = MA * MA

			FMSQRT(MA,MB) MB = SQRT(MA)

			FMSQRT_R1(MA) MA = SQRT(MA)

			FMST2M(STRING,MA) MA = STRING
			Convert from character string to FM. STRING may be in any
			numerical format. FMST2M is often more convenient than FMINP,
			which converts an array of character(1) values. Example:
			CALL FMST2M('123.4',MA)

			FMSUB(MA,MB,MC) MC = MA - MB

			FMSUB_R1(MA,MB) MA = MA - MB

			FMSUB_R2(MA,MB) MB = MA - MB

			FMTAN(MA,MB) MB = TAN(MA)

			FMTANH(MA,MB) MB = TANH(MA)

			FMTINY(MA) MA = Smallest positive FM number greater than underflow.

			FMULP(MA,MB) MB = One Unit in the Last Place of MA. For positive MA this is the
			same as the Fortran function SPACING, but MB < 0 if MA < 0.
			Examples: If MBASE = 10 and NDIG = 30, then ulp(1.0) = 1.0E-29,
			ulp(-4.5E+67) = -1.0E+38.

			FMVARS Write the current values of the internal FM variables on unit KW.

			FMWRITE(KWRITE,MA) Write MA on unit KWRITE. < old name: FMWRIT >
			Multi-line numbers will have '&' as the last nonblank character on all
			but the last line. These numbers can then be read easily using FMREAD.



			These are the available mathematical special functions.

			FMBERNOULLI(N,MA) MA = B(N) Nth Bernoulli number

			FMBESJ(N,MA,MB) MB = J(N,MA) Bessel function of the first kind

			FMBESJ2(N1,N2,MA,MB) MB = (/ J(N1,MA) , ..., J(N2,MA) /) returns an array

			FMBESY(N,MA,MB) MB = Y(N,MA) Bessel function of the second kind

			FMBESY2(N1,N2,MA,MB) MB = (/ Y(N1,MA) , ..., Y(N2,MA) /) returns an array

			FMBETA(MA,MB,MC) MC = Beta(MA,MB)

			FMC(MA,MB) MB = C(MA) Fresnel Cosine Integral

			FMCHI(MA,MB) MB = Chi(MA) Hyperbolic Cosine Integral

			FMCI(MA,MB) MB = Ci(MA) Cosine Integral

			FMCOMB(MA,MB,MC) MC = Combination MA choose MB (Binomial coefficient)

			FMEI(MA,MB) MB = Ei(MA) Exponential Integral

			FMEN(N,MA,MB) MB = E(N,MA) Exponential Integral E_n

			FMERF(MA,MB) MB = Erf(MA) Error function

			FMERFC(MA,MB) MB = Erfc(MA) Complimentary Error function

			FMERFC_SCALED(MA,MB) MB = Exp(x*2) Erfc(MA)

			FMEULER(MA) MA = Euler's constant ( 0.5772156649... ) < old name: FMEULR >

			FMFACT(MA,MB) MB = MA Factorial (Gamma(MA+1))

			FMGAM(MA,MB) MB = Gamma(MA)

			FMIBTA(MX,MA,MB,MC) MC = Incomplete Beta(MX,MA,MB)

			FMIGM1(MA,MB,MC) MC = Incomplete Gamma(MA,MB). Lower case Gamma(a,x)

			FMIGM2(MA,MB,MC) MC = Incomplete Gamma(MA,MB). Upper case Gamma(a,x)

			FMLERC(MA,MB) MB = Ln(Erfc(MA)) Log Erfc

			FMLI(MA,MB) MB = Li(MA) Logarithmic Integral

			FMLNGM(MA,MB) MB = Ln(Gamma(MA))

			FMPGAM(N,MA,MB) MB = Polygamma(N,MA) (Nth derivative of Psi)

			FMPOCH(MA,N,MB) MB = MA(MA+1)(MA+2)...(MA+N-1) (Pochhammer)

			FMPSI(MA,MB) MB = Psi(MA) (Derivative of Ln(Gamma(MA))

			FMS(MA,MB) MB = S(MA) Fresnel Sine Integral

			FMSHI(MA,MB) MB = Shi(MA) Hyperbolic Sine Integral

			FMSI(MA,MB) MB = Si(MA) Sine Integral


			These are the routines that deal with multiple precision integer numbers.
			All are subroutines except logical function IMCOMPARE. MA, MB, MC refer to IM format numbers.
			In each case the version of the routine to handle packed IM numbers has the same name, with
			'IM' replaced by 'IP'.

			IMABS(MA,MB) MB = ABS(MA)

			IMADD(MA,MB,MC) MC = MA + MB

			IMBIG(MA) MA = 10(106)
			Larger IM numbers can be obtained, but setting MA to the largest
			possible value would leave no room for any other numbers.

			IMCOMPARE(MA,LREL,MB) Logical comparison of MA and MB. < old name: IMCOMP >
			LREL is a character(2) value identifying which of the six comparisons
			is to be made.
			Example: IF (IMCOMPARE(MA,'GE',MB)) ...
			Also can be: IF (IMCOMPARE(MA,'>=',MB))
			character(1) is ok: IF (IMCOMPARE(MA,'>',MB)) ...

			IMDIM(MA,MB,MC) MC = DIM(MA,MB)

			IMDIV(MA,MB,MC) MC = int(MA/MB)
			Use IMDIVR if the remainder is also needed.

			IMDIVI(MA,IVAL,MB) MB = int(MA/IVAL)
			IVAL is a one word integer. Use IMDVIR to get the remainder also.

			IMDIVR(MA,MB,MC,MD) MC = int(MA/MB), MD = MA mod MB
			When both the quotient and remainder are needed, this routine is
			twice as fast as calling both IMDIV and IMMOD.

			IMDVIR(MA,IVAL,MB,IREM) MB = int(MA/IVAL), IREM = MA mod IVAL
			IVAL and IREM are one word integers.

			IMEQ(MA,MB) MB = MA

			IMFM2I(MAFM,MB) MB = MAFM Convert from real (FM) format to integer (IM) format.

			IMFORM(FORM,MA,STRING) MA is converted to a character string using format FORM and
			returned in STRING. FORM can represent I, F, E, or ES formats.
			Example: CALL IMFORM('I70',MA,STRING)

			IMFPRINT(FORM,MA) Print MA on unit KW using FORM format. < old name: IMFPRT >

			IMGCD(MA,MB,MC) MC = greatest common divisor of MA and MB.

			IMI2FM(MA,MBFM) MBFM = MA Convert from integer (IM) format to real (FM) format.

			IMI2M(IVAL,MA) MA = IVAL Convert from one word integer to IM.

			IMINP(LINE,MA,LA,LB) MA = LINE Input conversion.
			Convert LINE(LA) through LINE(LB) from characters to IM.

			IMM2DP(MA,X) X = MA Convert from IM to double precision.

			IMM2I(MA,IVAL) IVAL = MA Convert from IM to one word integer.

			IMM2SP(MA,X) X = MA Convert from IM to single precision.

			IMMAX(MA,MB,MC) MC = MAX(MA,MB)

			IMMIN(MA,MB,MC) MC = MIN(MA,MB)

			IMMOD(MA,MB,MC) MC = MA mod MB

			IMMPY(MA,MB,MC) MC = MA*MB

			IMMPYI(MA,IVAL,MB) MB = MA*IVAL Multiply by a one word integer.

			IMMPY_MOD(MA,MB,MC,MD) MD = MA*MB mod MC < old name: IMMPYM >
			Slightly faster than calling IMMPY and IMMOD separately.

			IMOUT(MA,LINE,LB) LINE = MA Convert from IM to character.
			LINE is a character array of length LB.

			IMPOWER(MA,MB,MC) MC = MA**MB < old name: IMPWR >

			IMPOWER_MOD(MA,MB,MC,MD) MD = MA**MB mod MC < old name: IMPMOD >

			IMPRINT(MA) Print MA on unit KW. < old name: IMPRNT >

			IMREAD(KREAD,MA) MA is returned after reading one (possibly multi-line)
			IM number on unit KREAD.
			This routine reads numbers written by IMWRITE.

			IMSIGN(MA,MB,MC) MC = SIGN(MA,MB) Returns the absolute value of MA times the
			sign of MB.

			IMSQR(MA,MB) MB = MA*MA Faster than IMMPY.

			IMST2M(STRING,MA) MA = STRING
			Convert from character string to IM.
			IMST2M is often more convenient than IMINP, which converts
			an array of character(1) values. Example:
			CALL IMST2M('12345678901',MA)

			IMSUB(MA,MB,MC) MC = MA - MB

			IMWRITE(KWRITE,MA) Write MA on unit KWRITE.
			Multi-line numbers will have '&' as the last nonblank character on all
			but the last line. These numbers can then be read easily using IMREAD.


			These are the routines that deal with multiple precision complex numbers.
			All are subroutines, and in each case the version of the routine to handle packed ZM numbers has
			the same name, with 'ZM' replaced by 'ZP'.

			MA, MB, MC refer to ZM format complex numbers.
			MAFM, MBFM, MCFM refer to FM format real numbers.
			INTEG is a Fortran INTEGER variable.
			ZVAL is a Fortran COMPLEX variable.

			ZMABS(MA,MBFM) MBFM = ABS(MA) Result is real.

			ZMACOS(MA,MB) MB = ACOS(MA)

			ZMACOSH(MA,MB) MB = ACOSH(MA)

			ZMADD(MA,MB,MC) MC = MA + MB

			ZMADDI(MA,INTEG) MA = MA + INTEG Increment an ZM number by a one word integer.
			Note this call does not have an "MB" result
			like ZMDIVI and ZMMPYI.

			ZMARG(MA,MBFM) MBFM = Argument(MA) Result is real.

			ZMASIN(MA,MB) MB = ASIN(MA)

			ZMASINH(MA,MB) MB = ASINH(MA)

			ZMATAN(MA,MB) MB = ATAN(MA)

			ZMATANH(MA,MB) MB = ATANH(MA)

			ZMCOMPLEX(MAFM,MBFM,MC) MC = CMPLX(MAFM,MBFM) < old name: ZMCMPX >

			ZMCONJUGATE(MA,MB) MB = CONJG(MA) < old name: ZMCONJ >

			ZMCOS(MA,MB) MB = COS(MA)

			ZMCOS_SIN(MA,MB,MC) MB = COS(MA), MC = SIN(MA). < old name: ZMCSSN >
			Faster than 2 calls.

			ZMCOSH(MA,MB) MB = COSH(MA)

			ZMCOSH_SINH(MA,MB,MC) MB = COSH(MA), MC = SINH(MA). < old name: ZMCHSH >
			Faster than 2 calls.

			ZMDIV(MA,MB,MC) MC = MA / MB

			ZMDIVI(MA,INTEG,MB) MB = MA / INTEG

			ZMEQ(MA,MB) MB = MA

			ZMEQU(MA,MB,NDA,NDB) MB = MA Version for changing precision.
			(NDA and NDB are as in FMEQU)

			ZMEXP(MA,MB) MB = EXP(MA)

			ZMFORM(FORM1,FORM2,MA,STRING)
			STRING = MA
			MA is converted to a character string using format FORM1 for the real
			part and FORM2 for the imaginary part. The result is returned in
			STRING. FORM1 and FORM2 can represent I, F, E, or ES formats.
			Example:
			CALL ZMFORM('F20.10','F15.10',MA,STRING)

			ZMFPRINT(FORM1,FORM2,MA) Print MA on unit KW using formats FORM1 and FORM2.
			< old name: ZMFPRT >

			ZMI2M(INTEG,MA) MA = CMPLX(INTEG,0)

			ZM2I2M(INTEG1,INTEG2,MA) MA = CMPLX(INTEG1,INTEG2)

			ZMIMAG(MA,MBFM) MBFM = IMAG(MA) Imaginary part.

			ZMINP(LINE,MA,LA,LB) MA = LINE Input conversion.
			Convert LINE(LA) through LINE(LB) from characters to ZM.
			LINE is a character array of length at least LB.

			ZMINT(MA,MB) MB = INT(MA) Integer part of both Real and Imaginary parts of MA.

			ZMIPOWER(MA,INTEG,MB) MB = MA ** INTEG Integer power function. < old name: ZMIPWR >

			ZMLOG10(MA,MB) MB = LOG10(MA) < old name: ZMLG10 >

			ZMLN(MA,MB) MB = LOG(MA)

			ZMM2I(MA,INTEG) INTEG = INT(REAL(MA))

			ZMM2Z(MA,ZVAL) ZVAL = MA

			ZMMPY(MA,MB,MC) MC = MA * MB

			ZMMPYI(MA,INTEG,MB) MB = MA * INTEG

			ZMNINT(MA,MB) MB = NINT(MA) Nearest integer of both Real and Imaginary.

			ZMOUT(MA,LINE,LB,LAST1,LAST2)
			LINE = MA
			Convert from FM to character.
			LINE is the returned character(1) array.
			LB is the dimensioned size of LINE.
			LAST1 is returned as the position in LINE of the last character
			of REAL(MA)
			LAST2 is returned as the position in LINE of the last character
			of AIMAG(MA)

			ZMPOWER(MA,MB,MC) MC = MA ** MB < old name: ZMPWR >

			ZMPRINT(MA) Print MA on unit KW using current format. < old name: ZMPRNT >

			ZMRATIONAL_POWER(MA,IVAL,JVAL,MB)
			MB = MA ** (IVAL/JVAL) < old name: ZMRPWR >

			ZMREAD(KREAD,MA) MA is returned after reading one (possibly multi-line) ZM number on
			unit KREAD. This routine reads numbers written by ZMWRITE.

			ZMREAL(MA,MBFM) MBFM = REAL(MA) Real part.

			ZMSET(NPREC) Set precision to the equivalent of a few more than NPREC base 10
			digits. This is now the same as FMSET, but is retained for
			compatibility with earlier versions of the package.

			ZMSIN(MA,MB) MB = SIN(MA)

			ZMSINH(MA,MB) MB = SINH(MA)

			ZMSQR(MA,MB) MB = MA*MA Faster than ZMMPY.

			ZMSQRT(MA,MB) MB = SQRT(MA)

			ZMST2M(STRING,MA) MA = STRING
			Convert from character string to ZM. ZMST2M is often more
			convenient than ZMINP, which converts an array of character(1)
			values. Example:
			CALL ZMST2M('123.4+5.67i',MA).

			ZMSUB(MA,MB,MC) MC = MA - MB

			ZMTAN(MA,MB) MB = TAN(MA)

			ZMTANH(MA,MB) MB = TANH(MA)

			ZMWRITE(KWRITE,MA) Write MA on unit KWRITE. Multi-line numbers are formatted for
			automatic reading with ZMREAD. < old name: ZMWRIT >

			ZMZ2M(ZVAL,MA) MA = ZVAL


			10. EFFICIENCY

			To take advantage of hardware architecture on different machines, the package has been designed
			so that the arithmetic used to perform the multiple precision operations can easily be changed.
			All variables that must be changed to get a different arithmetic have names beginning with 'M'
			and are declared using REAL (KIND(1.0D0)) ...

			For example, to change the package to use integer arithmetic internally, make these two changes
			everywhere in the FM.f95 file.
			Change 'REAL (KIND(1.0D0))' to 'INTEGER'.
			Change 'AINT (' to '('. Note the blank between AINT and (.

			In many places in FM, an AINT function is not supposed to be changed. These are written 'AINT(',
			with no embedded blank, so they will not be changed by the global change above.

			The first of these changes must also be made throughout the file FMSAVE.f95.
			Change 'REAL (KIND(1.0D0))' to 'INTEGER'.

			This version of FM restricts the base used to be also representable in integer variables, so
			using precision above double usually does not save much time unless integers can also be declared
			at a higher precision. Using IEEE Extended would allow a base of around 10**9 to be chosen, but
			the delayed digit-normalization method used for multiplication and division means that a slightly
			smaller base like 10**8 would probably run faster. This would usually not be much faster than
			using the usual base 10**7 with double precision.

			For cases where special compiler directives or minor re-writing of the code may improve speed,
			several of the most important loops in FM are identified by comments containing the string
			'(Inner Loop)'.


			11. NEW FOR VERSION 1.3

			This is the second edition of version 1.3. The first edition appeared in ACM Transactions on
			Mathematical Software (2-2011). Since then several additions have been made.
			(a) New Fortran-08 functions are available in FMZM:
			ACOSH(X), ASINH(X), ATANH(X) for real and complex X
			ATAN(X,Y) can be used in place of ATAN2(X,Y)
			BESSEL_J0(X), BESSEL_J1(X), BESSEL_JN(N,X), BESSEL_JN(N1,N2,X)
			BESSEL_Y0(X), BESSEL_Y1(X), BESSEL_YN(N,X), BESSEL_YN(N1,N2,X)
			The older FM names, BESSEL_J(N,X) and BESSEL_Y(N,X) are still available.
			ERFC_SCALED(X) for exp(x*2) erfc(x)
			The older FM function LOG_ERFC(X) is also still available for avoiding underflow in erfc.
			HYPOT(X,Y) for sqrt(x2 + y2)
			NORM2(A) for sqrt( a(1)2 + a(2)2 + ... + a(n)**2 )
			This could previously have been done with array operations as SQRT(DOT_PRODUCT(A,A)).
			(b) Many of the elementary and special functions are now faster, after some code-tuning was
			done and a few new methods were added.

			The routines for the exponential integral function and related mathematical special functions
			are new in version 1.3. These routines are:
			FMBESJ, FMBESY, FMC, FMCHI, FMCI, FMEI, FMEN, FMERF, FMERFC, FMLERC, FMLI, FMS, FMSHI, FMSI.

			Some of the routines were moved between files FM.f95 and FMZM90.f95 so that now all routines
			using the module FMZM (in file FMZM90.f95) for multiple precision derived types and operator
			overloading are located in FMZM90.f95. This means that programs not using derived types can
			skip compiling and/or linking FMZM90.f95.

			The array function DOTPRODUCT in FMZM has been re-named DOT_PRODUCT to agree with the Fortran
			standard. For type ZM complex arguments its definition has been changed to agree with the
			Fortran intrinsic function. When X and Y are complex, DOT_PRODUCT(X,Y) is not just the sum of
			the products of the corresponding array elements, as it is for types FM and IM. For type ZM,
			the formula is the sum of conjg(X(j)) * Y(j). This definition is used so that the complex dot
			product will be an inner product in the mathematical sense.

			New routines have been added to module FMZM to provide array syntax for the three multiple
			precision derived types. This means statements like V = 1 and A = B + C now work when these
			variables are vectors or matrices of multiple precision numbers.

			One routine from FM 1.2 has been split into three routines in version 1.3. The routine
			FM_RANDOM_SEED from FM 1.2 has become three subroutines, so that the optional arguments and
			the need for an explicit interface can be avoided. See the three routines starting with
			FM_RANDOM_SEED in the list above. The same multiplicative congruential generator as before
			is used, but the shuffling of those values has been removed, so that saving seeds and
			re-starting the generator now works more like the standard Fortran random function.

			Multiple precision variables were separate fixed-size arrays in previous versions. Now they are
			single integers that serve as index values to a single large array (MWK, defined in file
			FMSAVE.f95) where the actual values are stored. This often improves both efficiency and memory
			utilization, since many compilers implemented the derived type operations using copy in and copy
			out of the arguments for a given operation. Copying entire arrays was slower, and there were
			often memory leaks when the compiler automatically created temporary derived type objects while
			evaluating derived type expressions. The static arrays in previous versions also meant that
			memory was wasted when only a few kinds of operations were used at high precision. Now the
			space needed by any unused operations never gets allocated.

			Some new error checking is now done for the derived type multiple precision variables. Attempting
			to use an undefined variable will cause an error message to be printed.

			Much higher precision can be attained in version 1.3, since machines are faster and have more
			memory. To support higher precision, a routine for FFT-based multiplication has been included,
			and when precision gets high enough, the algorithms for multiplication, division, squares, square
			roots, etc., will switch to the FFT routine.

			Binary splitting algorithms are used for the mathematical constants at high precision. At the
			time version 1.3 was released, computing a million digits of e, pi, or the logarithm of a small
			integer took a few seconds, while a million digits of Euler's constant took a few minutes.

			Perfect rounding is now done all the time. In version 1.2 perfect rounding was an option, but
			the default rounding could round the wrong direction once every few million operations, when the
			exact result was very close to halfway between two adjacent representable numbers.


			================================================================================
			================================================================================

cosolver/basis_transformation/FMlib/FM_User_Manual.txtcd1286e2d410improved_sugama

FM_User_Manual.txtView Options

cosolver/basis_transformation/FMlib/FM_User_Manual.txt
cd1286e2d410improved_sugama

FM_User_Manual.txt
View Options