**************************************************************************MAPMVersion 4.9.5December 10, 2007Michael C. Ringringx004@tc.umn.eduLatest release will be available athttp://tc.umn.edu/~ringx004**************************************************************************** ** Copyright (C) 1999 - 2007 Michael C. Ring ** ** This software is Freeware. ** ** Permission to use, copy, and distribute this software and its ** documentation for any purpose with or without fee is hereby granted, ** provided that the above copyright notice appear in all copies and ** that both that copyright notice and this permission notice appear ** in supporting documentation. ** ** Permission to modify the software is granted. Permission to distribute ** the modified code is granted. Modifications are to be distributed by ** using the file 'license.txt' as a template to modify the file header. ** 'license.txt' is available in the official MAPM distribution. ** ** To distribute modified source code, insert the file 'license.txt' ** at the top of all modified source code files and edit accordingly. ** ** This software is provided "as is" without express or implied warranty. ** ****************************************************************************---------------------------------------Mike's Arbitrary Precision Math Library---------------------------------------Mike's Arbitrary Precision Math Library is a set of functions thatallow the user to perform math to any level of accuracy that isdesired. The inspiration for this library was Lloyd Zusman's similarAPM package that was released in ~1988. I borrowed some of his ideasin my implementation, creating a new data type (MAPM) and how the datatype was used by the user. However, there were a few things I wantedmy library to do that the original library did not :1) Round a value to any desired precision. This is very handy whenmultiplying for many iterations. Since multiplication guarantees anexact result, the number of digits will grow without bound. I wanteda way to trim the number of significant digits that were retained.2) A natural support for floating point values. From most of the otherlibraries I looked at, they seem to have a preference for integeronly type math manipulations. (However, this library will also dointeger only math if you desire).And if a library can only do integers, it can't do ...3) Trig functions and other common C math library functions. This librarywill perform the following functions to any desired precision level :SQRT, CBRT, SIN, COS, TAN, ARC-SIN, ARC-COS, ARC-TAN, ARC-TAN2, LOG,LOG10, EXP, POW, SINH, COSH, TANH, ARC-SINH, ARC-COSH, ARC-TANH, andalso FACTORIAL. The full 'math.h' is not duplicated, though I thinkthese are most of the important ones. My definition of what's importantis what I've actually used in a real application.**************************************************************************NOTE:There is a COMPILER BUG in Microsoft's Visual C++ 7.x (VS.NET 2003) whichprevents a C++ MAPM application from compiling.This only affects C++ applications. C applications are OK.The compiler bug creates an error C2676 similar to this:<path>...\mapm-49\M_APM.H(###) : error C2676: binary operator '-': 'MAPM'doesn't define this operator or a conversion to a suitable type for thepredefined operator.To work around this bug, go to http://www.microsoft.comIn the upper right corner of web page, search for "814455".The results of the search will point you to an article on how to workaround the problem.**************************************************************************NOTE: MAPM Library History can now be found in 'history.txt'**************************************************************************ANOTHER NOTE: For the Windows/DOS distribution, the filename conventionwill always be in 8.3 format. This is because there are some users whostill use 16 bit DOS ....(I really wasn't sure what to call this library. 'Arbitrary Precision Math'is a defacto standard for what this does, but that name was already taken,so I just put an 'M' in front of it ...)**************************************************************************MAPM LIBRARY NUMERICAL LIMITATIONS:A general floating point number is of the form:Sn.nnnnnnnnE+yyyy ex: -4.318384739357974916E+6215Sn.nnnnnnnnE-yyyy ex: +8.208237913789131096523645193E-12873'S' is the sign, + or -.In MAPM, a number (n.nnnn...) may contain up to ( INT_MAX - 1 ) digits.For example, an MAPM number with a 16 bit integer may contain 2^15 or 32,767digits. An MAPM number with a 32 bit integer may contain 2^31 or 2,147,483,647digits. All MAPM numbers are stored in RAM, there is no "data on disk" option.So to represent the maximum number of digits on a 32 bit CPU will requiregreater than 2 Gig of RAM.If you have a CPU with 64 bit ints, then the limitation is 2^63 or9,223,372,036,854,775,807. It should be a very long time before computershave this much RAM in them.For the exponent (yyyy), the limitations are also INT_MAX and INT_MIN.For a 16 bit CPU, the largest number you can represent is approx0.9999999....E+32767. (H)For a 16 bit CPU, the smallest number you can represent (other than 0)is approx 0.1E-32767. (L)For a 32 bit CPU, the largest number you can represent is approx0.9999999....E+2147483647. (H)For a 32 bit CPU, the smallest number you can represent (other than 0)is approx 0.1E-2147483647. (L)The limitations for negative numbers are the same as positive numbers.Real Number Axis+------------------------+ --- +--------------------------+| | | | |-H -L 0.0 +L +HMAPM can represent real numbers from -H to -L, 0.0, +L to +H.The number of significant digits is typically only limited by available RAM.In MAPM, numerical overflows and underflows are not handled very well(actually not at all). There really isn't a clean and portable way todetect integer overflow and underflow. Per K&R C, the result of integeroverflow/underflow is implementation dependent. In assembly language, whenyou add two numbers, you have access to a "carry flag" to see if an overflowoccurred. C has no analogous operator to a carry flag.It is up to the user to decide if the results are valid for a given operation.In a 32 bit environment, the limit is so large that this is likely not anissue for most typical applications. However, it doesn't take much to overflowa 16 bit int so care should taken in a 16 bit environment.The reaction to an integer overflow/underflow is unknown at run-time:o) adding 2 large positive numbers may silently roll over to anegative number.o) in some embedded applications an integer overflow/underflow triggersa hardware exception.Since I don't have control over where this library could possibly run,I chose to ignore the problem for now. If anyone has some suggestions(that's portable), please let me know.**************************************************************************KNOWN BUGS : (Other than integer overflow discussed above....) None**************************************************************************IF YOU ARE IN A HURRY ...UNIX: (assumes gcc compiler)run make (build library + 4 executables)run make -f makefile.osx (for MAC OSX)--- OR ---run: mklib (this will create the library, lib_mapm.a)run: mkdemo (this will create 4 executables, 'calc', 'validate','primenum', and 'cpp_demo')DOS / Win NT/9x (in a DOS window for NT/9x):see the file 'README.DOS' for instructions.**************************************************************************calc: This is a command line version of an RPN calculator. If you arefamiliar with RPN calculators, the use of this program will bequite obvious. The default is 30 decimal places but this can bechanged with the '-d' command line option. This is not aninteractive program, it just computes an expression from the commandline. Run 'calc' with no arguments to get a list of the operators.compute : (345.2 * 87.33) - (11.88 / 3.21E-2)calc 345.2 87.33 x 11.88 3.21E-2 / -result: 29776.22254205607476635514018692compute PI to 70 decimal places : (6 * arcsin(0.5))calc -d70 6 .5 as xresult :3.1415926535897932384626433832795028841971693993751058209749445923078164calc -d70 -1 ac (arccos(-1) for fastest possible way)validate : This program will compare the MAPM math functions to the Cstandard library functions, like sqrt, sin, exp, etc. Thisshould give the user a good idea that the library is operatingcorrectly. The program will also perform high precision mathusing known quantities like PI, log(2), etc.'validate' also attempts to obtain as much code coverage of thelibrary as is practical. I used 'gcov' (available with the gccdistribution) to test the code coverage. 100% coverage is notobtained, compromises must be made in order to have the programrun in a reasonable amount of time.primenum: This program will generate the first 10 prime numbers startingwith the number entered as an argument.Example: primenum 1234567890 will output (actually 1 per line):this took ~3 sec on my Linux PC, a 350MHz PII.1234567891, 1234567907, 1234567913, 1234567927, 1234567949,1234567967, 1234567981, 1234568021, 1234568029, 1234568047**************************************************************************To use the library, simply include 'm_apm.h' and link your programwith the library, libmapm.a (unix) or mapm.lib (dos).Note: If your system creates libraries with a '.a' extension, then thelibrary will be named libmapm.a. (This conforms to typical unix namingconventions).Note: If your system creates libraries with a '.lib' extension, then thelibrary will be named mapm.lib.For unix builds, you also may need to specify the math library (-lm) whenlinking. The reason is some of the MAPM functions use an iterative algorithm.When you use an iterative solution, you have to supply an initial guess. Iuse the standard math library to generate this initial guess. I debatedwhether this library should be stand-alone, i.e. generate it's own initialguesses with some algorithm. In the end, I decided using the standard mathlibrary would not be a big inconvienence and also it was just too temptinghaving an immediate 15 digits of precision. When you prime the iterativeroutines with 15 accurate digits, the MAPM functions converge faster.See the file 'algorithms.used' to see a description of the algorithms Iused in the library. Some I derived on my own, others I borrowed from peoplesmarter than me. Version 2 of the library supports a 'fast' multiplicationalgorithm. The algorithm used is described in the algorithms.used file. Aconsiderable amount of time went into finding fast algorithms for thelibrary. However, some could possibly be even better. If anyone has amore efficient algorithm for any of these functions, I would like to herefrom you.See the file 'function.ref' to see a description of all the functions inthe library and the calling conventions, prototypes, etc.See the file 'struct.ref' which documents how I store the numbers internallyin the MAPM data structure. This will not be needed for normal use, but itwill be very useful if you need to change/add to the existing library.**************************************************************************USING MAPM IN A MULTI-THREADED APPLICATION :Note that the default MAPM library is NOT thread safe. MAPM internal datastructures could get corrupted if multiple MAPM functions are active at thesame time. The user should guarantee that only one thread is performingMAPM functions. This can usually be achieved by a call to the operatingsystem to obtain a 'semaphore', 'mutex', or 'critical code section' so theoperating system will guarantee that only one MAPM thread will be activeat a time.The necessary function wrappers for thread safe operation can be found inthe sub-directory 'multi_thread' (unix) or 'multithd' (Win/Dos). For now,only Microsoft Visual C++ 6.0 is supported.**************************************************************************QUICK TEMPLATE FOR NORMAL USE :The MAPM math is done on a new data type called "M_APM". This isactually a pointer to a structure, but the contents of the structureshould never be manipulated: all operations on MAPM entities are donethrough a functional interface.The MAPM routines will automatically allocate enough space in theirresults to hold the proper number of digits.The caller must initialize all MAPM values before the routines canoperate on them (including the values intended to contain results ofcalculations). Once this initialization is done, the user never needsto worry about resizing the MAPM values, as this is handled inside theMAPM routines and is totally invisible to the user.The result of a MAPM operation cannot be one of the other MAPM operands.If you want this to be the case, you must put the result into atemporary variable and then assign (m_apm_copy) it to the appropriate operand.All of the MAPM math functions begin with m_apm_*.There are some predefined constants for your use. In case it's not obvious,these should never appear as a 'result' parameter in a function call.The following constants are available : (declared in m_apm.h)MM_Zero MM_One MM_Two MM_ThreeMM_Four MM_Five MM_TenMM_PI MM_HALF_PI MM_2_PI MM_EMM_LOG_E_BASE_10 MM_LOG_10_BASE_E MM_LOG_2_BASE_E MM_LOG_3_BASE_EThe non-integer constants above (PI, log(2), etc) are accurate to 128decimal places. The file mapmcnst.c contains these constants. I'veincluded 512 digit constants in this file also that are commented out.If you need more than 512 digits, you can simply use the 'calc' programto generate more precise constants (or create more precise constants atrun-time in your app). The number of significant digits in the constantsshould be 6-8 more than the value specified in the #define.Basic plan of attack:(1) get your 'numbers' into M_APM format.(2) do your high precision math.(3) get the M_APM numbers back into a format you can use.-------- (1) --------#include "m_apm.h" /* must be included before any M_APM 'declares' */M_APM area_mapm; /* declare your variables */M_APM tmp_mapm;M_APM array_mapm[10]; /* can use normal array notation */M_APM array2d_mapm[10][10];area_mapm = m_apm_init() /* init your variables */tmp_mapm = m_apm_init();for (i=0; i < M; i++) /* must init every element of the array */array_mapm[i] = m_apm_init();for (i=0; i < M; i++)for (j=0; j < N; j++)array2d_mapm[i][j] = m_apm_init();/** there are 3 ways to convert your number into an M_APM number* (see the file function.ref)** a) literal string (exponential notation OK)* b) long variable* c) double variable*/m_apm_set_string(tmp_mapm, "5.3286E-7");m_apm_set_long(array_mapm[6], -872253L);m_apm_set_double(array2d_mapm[3][7], -529.4486711);-------- (2) --------do your math ...m_apm_add(cc_mapm, aa_mapm, MM_PI);m_apm_divide(bb_mapm, DECIMAL_PLACES, aa_mapm, MM_LOG_2_BASE_E);m_apm_sin(bb_mapm, DECIMAL_PLACES, aa_mapm);whatever ...-------- (3) --------There are 5 total functions for converting an M_APM number into somethinguseful. (See the file 'function.ref' for full function descriptions)For these 5 functions, M_APM number -> string is the conversion======================= METHOD 1 ==== : floating point values (m_apm_to_string)===================the format will be in scientific (exponential) notationoutput string = [-]n.nnnnnE+x or ...E-xwhere 'n' are digits and the exponent will be always be present,including E+0it's easy to convert this to a double:double dtmp = atof(out_buffer);======================= METHOD 2 ==== : floating point values (m_apm_to_fixpt_string)=================== (m_apm_to_fixpt_stringex)(m_apm_to_fixpt_stringexp)the format will be in fixed point notationoutput string = [-]mmm.nnnnnnwhere 'm' & 'n' are digits.======================= METHOD 3 ==== : integer values (m_apm_to_integer_string)===================the format will simply be digits with a possible leading '-' sign.output string = [-]nnnnnnwhere 'n' are digits.it's easy to convert this to a long :long mtmp = atol(out_buffer);... or an int :int itmp = atoi(out_buffer);Note that if the M_APM number has a fractional portion, the fractionwill be truncated and only the integer portion will be output.char out_buffer[1024];m_apm_to_string(out_buffer, DECIMAL_PLACES, mapm_number);m_apm_to_fixpt_string(out_buffer, DECIMAL_PLACES, mapm_number);m_apm_to_integer_string(out_buffer, mapm_number);*************************************************************************************************************************************** NOTES on the fixed point formatting functions *************************************************************Assume you have the following code:--> m_apm_set_string(aa_mapm, "2.0E18");--> m_apm_sqrt(bb_mapm, 40, aa_mapm);--> m_apm_to_string(buffer, 40, bb_mapm);--> fprintf(stdout,"[%s]\n",buffer);--> m_apm_to_fixpt_string(buffer, 40, bb_mapm);--> fprintf(stdout,"[%s]\n",buffer);It is desired to compute the sqrt(2.0E+18) to40 significant digits. You then want the resultoutput with 40 decimal places. But the outputfrom above is :[1.4142135623730950488016887242096980785697E+9][1414213562.3730950488016887242096980785697000000000]Why are there 9 '0' in the fixed point formatted string??The sqrt calculation computed 40 significant digits relativeto the number in EXPONENTIAL format. When the number isoutput in exponential format, the 40 digits are as expectedwith an exponent of 'E+9'.The same number formatted as fixed point appears to be anerror. Remember, we computed 40 significant digits. However,the result has an exponent of '+9'. So, 9 of the digits areneeded *before* the decimal point. In our calculation, only31 digits of precision remain from our original 40. We thenasked the fixed point formatting function to format 40 digits.Only 31 are left so 9 zeros are used as pad at the end tofulfill the 40 places asked for.Keep this in mind if you truly desire more accurate resultsin fixed point formatting and your result contains a largepositive exponent.**************************************************************************MAPM C++ WRAPPER CLASS:Orion Sky Lawlor (olawlor@acm.org) has added a very nice C++ wrapperclass to m_apm.h. This C++ class will have no effect if you just usea normal C compiler. The library will operate as before with no userimpacts.For now, I recommend compiling the library as 'C'. The library willcompile as a C++ library, but then it is likely that you would notbe able to use the library in a straight C application. Since the C++wrapper class works very nicely as is, there is no pressing need to compilethe library as C++.See the file 'cpp_function.ref' to see a description of how to usethe MAPM class.To compile and link the C++ demo program: (assuming the library isalready built)UNIX:g++ cpp_demo.cpp lib_mapm.a -s -o cpp_demo -lmGCC for DOS: (gxx (or g++) is the C++ compiler)gxx cpp_demo.cpp lib_mapm.a -s -o cpp_demo.exe -lmUsing the C++ wrapper allows you to do things like:// Compute the factorial of the integer nMAPM factorial(MAPM n){MAPM i;MAPM product = 1;for (i=2; i <= n; i++)product *= i;return product;}The syntax is the same as if you were just writing normal code, but allthe computations will be performed with the high precision math library,using the new 'datatype' MAPM.The default precision of the computations is as follows:Addition, subtraction, and multiplication will maintain ALL significantdigits.All other operations (divide, sin, etc) will use the following rules:1) if the operation uses only one input value [y = sin(x)], the result 'y'will be the same precision as 'x', with a minimum of 30 digits if 'x' isless than 30 digits.2) if the operation uses two input values [z = atan2(y,x)], the result 'z'will be the max digits of 'x' or 'y' with a minimum of 30.The default precision is 30 digits. You can change the precision atany time with the function 'm_apm_cpp_precision'. (See function.ref)----> m_apm_cpp_precision(80);will result in all operations being accurate to a minimum of 80 significantdigits. If any operand contains more than the minimum number of digits, thenthe result will contain the max number of digits of the operands.NOTE!: Some real life use with the C++ wrapper has revealed a certaintendency for a program to become quite slow after many iterations(like in a for/while loop). After a little debug, the reasonbecame clear. Remember that multiplication will maintain ALLsignificant digits :20 digit number x 20 digit number = 40 digits40 digit number x 40 digit number = 80 digits80 digit number x 80 digit number = 160 digitsetc.So after numerous iterations, the number of significant digitswas growing without bound. The easy way to fix the problem isto simply *round* your result after a multiply or some othercomplex operation. For example:#define MAX_DIGITS 256p1 = (p0 * sin(b1) * exp(1 + u1)) / sqrt(1 + b1);p1 = p1.round(MAX_DIGITS);If you 'round' as shown here, your program will likely benearly as fast as a straight 'C' version.NOTE #2!Reference the following code snippet:...MAPM pi1, pi2;char obuf[256];m_apm_cpp_precision(62);pi1 = 2 * asin("1");pi2 = 2 * asin(1.0);pi1.toString(obuf, 60); printf("PI1 = [%s] \n",obuf);pi2.toString(obuf, 60); printf("PI2 = [%s] \n",obuf);...On my system, the output is :PI1 = [3.141592653589793238462643383279502884197169399375105820974945E+0]PI2 = [3.141592653589790000000000000000000000000000000000000000000000E+0]PI2 only has 15 significant digits! This is due to how the secondasin is called. It is called with a 'double' as the argument, hencethe compiler will use the default double asin function from thestandard library. This is likely not the intent but this would beeasy to miss if this was a complex calculation and we didn't knowthe 'right' answer.In order to force the use of the overloaded MAPM functions, call theMAPM functions with a quoted string as the argument (if the argumentis a constant and not a variable).This would also work (though it seems less elegant ...) :MAPM t = 1;pi2 = 2 * asin(t);-----------If you have any questions or problems with the C++ wrapper, please letme know. I am not very C++ proficient, but I'd still like to know about anyproblems. Orion Sky Lawlor (olawlor@acm.org) is the one who implementedthe MAPM class, so he'll have to resolve any real hardcore problems, if youhave any.**************************************************************************TESTING :Testing the library was interesting. How do I know it's right? Since Itest the library against the standard C runtime math library (see below)I have high confidence that the MAPM library is giving correct results.The logic here is the basic algorithms are independent of the number ofdigits calculated, more digits just takes longer.The MAPM library has been tested in the following environments :Linux i486 / gcc 2.7.2.3, gcc 2.95.2Linux i686 / gcc 2.91.66, gcc 2.95.2, gcc 2.95.3, gcc 3.0.4, gcc 3.2.3Linux i686 / Intel Linux C/C++ Compiler Verison 7.0 / 8.1FreeBSD 4.1 / gcc 2.95.1FreeBSD 4.8 / gcc 2.95.4FreeBSD 5.x / gcc 2.95.2, gcc 2.95.3Redhat Linux 8.2 / gcc 3.2HP-UX 9.x /gcc 2.5.8HP-UX 10.x / gcc 2.95.2Sun 5.5.1 (output from uname), gcc 3.1.1Sun Solaris 2.6 / gcc 2.95.1, gcc 2.95.3, gcc 3.2.3MAC OSX / gcc ?DOS 5.0 using GCC 2.8.1 for DOSDOS 5.0 using GCC 2.95.2 for DOSDOS ??? using Borland Turbo C++ 3.0WIN NT+SP5 using Borland C++ 5.02 IDE, 5.2 & 5.5 command line.WIN 2000 using National Instruments LabWindows CVI 6.0WIN98 using Borland C++ 5.5 command line.WIN98 & NT & 2000 & XP using LCC-WIN32 Ver 3.2, 3.3WIN98 & NT using Watcom C 11.xWIN95 & NT using Open Watcom 1.0WIN95 & WIN98 using MINGW-32 version mingw-1.0.1-20010726WIN95 & WIN2000 using DEV-CPP 5.0 Beta 8, 4.9.8.0MINGW-32 with gcc 3.2 (mingw special 20020817-1)MINGW-32 with gcc 3.2.3 (mingw special 20030504-1)WINXP & MINGW-32 with gcc 3.4.5WINXP & Digital Mars Compiler 8.49WIN?? using Metrowerks CodeWarrior Pro 7.0 for WindowsDOS 5.0 using Microsoft C 5.1 and 8.00c (16 bit EXEs)WIN98 & NT using Microsoft Visual C++ 6.0WIN98 & NT using Microsoft Visual C++ 7.x (VS.NET 2003 except forknown compiler bug C2676)As a general rule, when calculating a quantity to a given number of decimalplaces, I calculated 4-6 extra digits and then rounded the result to what wasasked for. I decided to be conservative and give a correct answer rather thanto be faster and possibly have the last 2-3 digits in error. Also, some ofthe functions call other functions (calculating arc-cos will call cos, logwill call exp, etc.) so I had to calculate a few extra digits in each iterationto guarantee the loops terminated correctly.1) I debugged the 4 basic math operations. I threw numerous test cases ateach of the operations until I knew they were correct.Also note that the math.h type functions all call the 4 basic operationsnumerous times. So if all the math.h functions work, it is highlyprobable the 4 basic math operations work also.2) 'math.h' type functions.SQRT: Not real hard to check. Single stepping through the iterativeloop showed it was always converging to the sqrt.CBRT: Similar to sqrt, single stepping through the iterative loopshowed it was always converging to the cube root.EXP: I wrote a separate algorithm which expanded the Taylor seriesmanually and compared the results against the library.POW: Straightforward since this just calls 'EXP'.LOG: I wrote a separate algorithm which expanded the Taylor seriesmanually and compared the results against the library. Thistook a long time to execute since the normal series convergesVERY slowly for the log function. This is why the LOG functionuses an iterative algorithm.LOG10: Straightforward since this just calls 'LOG'.SIN/COS: I wrote a separate algorithm which expanded the Taylor seriesmanually and compared the results against the library.TAN: Straightforward since this just calls 'SIN' and 'COS'.ARC-x: Single stepping through the iterative loop showed the arcfamily of functions were always converging. Also used theseto compute PI. The value of PI is now known to many, manydigits. I computed PI to 1000+ digits by computing:6 * arcsin(0.5) and 4 * arctan(1) and 3 * arccos(0.5)and compared the output to the published 'real' values of PI.The arc family of functions exercise considerable portionsof the library.HYPERBOLIC: The hyperbolic functions just use exp, log, and the 4 basicmath operations. All of these functions simply use otherexisting functions in the library.FINALLY: Run the program 'validate'. This program compares the first13-14 significant digits of the standard C library againstthe MAPM library. If this program passes, you can feelconfident that the MAPM library is giving correct results.**************************************************************************