1**************************************************************************
2  
3				   MAPM
4
5			       Version 4.9.5
6
7			     December 10, 2007
8
9			      Michael C. Ring
10
11			    ringx004@tc.umn.edu
12
13		    Latest release will be available at
14		        http://tc.umn.edu/~ringx004
15
16***************************************************************************
17*									  *
18*  Copyright (C) 1999 - 2007   Michael C. Ring                            *
19*									  *
20*  This software is Freeware.						  *
21*									  *
22*  Permission to use, copy, and distribute this software and its          *
23*  documentation for any purpose with or without fee is hereby granted,   *
24*  provided that the above copyright notice appear in all copies and      *
25*  that both that copyright notice and this permission notice appear      *
26*  in supporting documentation.                                           *
27*									  *
28*  Permission to modify the software is granted. Permission to distribute *
29*  the modified code is granted. Modifications are to be distributed by   *
30*  using the file 'license.txt' as a template to modify the file header.  *
31*  'license.txt' is available in the official MAPM distribution.          *
32*									  *
33*  To distribute modified source code, insert the file 'license.txt'      *
34*  at the top of all modified source code files and edit accordingly.     *
35*									  *
36*  This software is provided "as is" without express or implied warranty. *
37*									  *
38***************************************************************************
39
40		---------------------------------------
41		Mike's Arbitrary Precision Math Library
42		---------------------------------------
43
44Mike's Arbitrary Precision Math Library is a set of functions that
45allow the user to perform math to any level of accuracy that is
46desired. The inspiration for this library was Lloyd Zusman's similar
47APM package that was released in ~1988. I borrowed some of his ideas
48in my implementation, creating a new data type (MAPM) and how the data
49type was used by the user. However, there were a few things I wanted
50my library to do that the original library did not :
51
521) Round a value to any desired precision. This is very handy when
53   multiplying for many iterations. Since multiplication guarantees an
54   exact result, the number of digits will grow without bound. I wanted
55   a way to trim the number of significant digits that were retained.
56
572) A natural support for floating point values. From most of the other
58   libraries I looked at, they seem to have a preference for integer
59   only type math manipulations. (However, this library will also do
60   integer only math if you desire).
61
62   And if a library can only do integers, it can't do ...
63
643) Trig functions and other common C math library functions. This library
65   will perform the following functions to any desired precision level :
66   SQRT, CBRT, SIN, COS, TAN, ARC-SIN, ARC-COS, ARC-TAN, ARC-TAN2, LOG,
67   LOG10, EXP, POW, SINH, COSH, TANH, ARC-SINH, ARC-COSH, ARC-TANH, and
68   also FACTORIAL.  The full 'math.h' is not duplicated, though I think
69   these are most of the important ones. My definition of what's important
70   is what I've actually used in a real application.
71
72**************************************************************************
73
74NOTE:
75
76There is a COMPILER BUG in Microsoft's Visual C++ 7.x (VS.NET 2003) which
77prevents a C++ MAPM application from compiling.
78
79This only affects C++ applications. C applications are OK.
80
81The compiler bug creates an error C2676 similar to this:
82
83<path>...\mapm-49\M_APM.H(###) : error C2676: binary operator '-': 'MAPM'
84doesn't define this operator or a conversion to a suitable type for the
85predefined operator.
86
87To work around this bug, go to http://www.microsoft.com
88
89In the upper right corner of web page, search for "814455".
90
91The results of the search will point you to an article on how to work
92around the problem.
93
94**************************************************************************
95
96NOTE: MAPM Library History can now be found in 'history.txt'
97
98**************************************************************************
99
100ANOTHER NOTE: For the Windows/DOS distribution, the filename convention
101will always be in 8.3 format. This is because there are some users who
102still use 16 bit DOS ....
103
104(I really wasn't sure what to call this library. 'Arbitrary Precision Math'
105is a defacto standard for what this does, but that name was already taken,
106so I just put an 'M' in front of it ...)
107
108**************************************************************************
109
110MAPM LIBRARY NUMERICAL LIMITATIONS:
111
112A general floating point number is of the form:
113
114Sn.nnnnnnnnE+yyyy        ex: -4.318384739357974916E+6215
115Sn.nnnnnnnnE-yyyy        ex: +8.208237913789131096523645193E-12873
116
117'S' is the sign, + or -.
118
119In MAPM, a number (n.nnnn...) may contain up to ( INT_MAX - 1 ) digits.
120
121For example, an MAPM number with a 16 bit integer may contain 2^15 or 32,767
122digits. An MAPM number with a 32 bit integer may contain 2^31 or 2,147,483,647
123digits. All MAPM numbers are stored in RAM, there is no "data on disk" option.
124So to represent the maximum number of digits on a 32 bit CPU will require
125greater than 2 Gig of RAM.
126
127If you have a CPU with 64 bit ints, then the limitation is 2^63 or
1289,223,372,036,854,775,807. It should be a very long time before computers
129have this much RAM in them.
130
131For the exponent (yyyy), the limitations are also INT_MAX and INT_MIN.
132
133For a 16 bit CPU, the largest number you can represent is approx
1340.9999999....E+32767.    (H)
135
136For a 16 bit CPU, the smallest number you can represent (other than 0)
137is approx 0.1E-32767.   (L)
138
139For a 32 bit CPU, the largest number you can represent is approx
1400.9999999....E+2147483647.   (H)
141
142For a 32 bit CPU, the smallest number you can represent (other than 0)
143is approx 0.1E-2147483647.  (L)
144
145The limitations for negative numbers are the same as positive numbers.
146
147
148                            Real Number Axis
149
150     +------------------------+    ---    +--------------------------+
151     |                        |     |     |                          |
152    -H                       -L    0.0   +L                         +H
153
154
155
156MAPM can represent real numbers from -H to -L, 0.0, +L to +H.
157
158The number of significant digits is typically only limited by available RAM.
159
160In MAPM, numerical overflows and underflows are not handled very well
161(actually not at all). There really isn't a clean and portable way to
162detect integer overflow and underflow. Per K&R C, the result of integer
163overflow/underflow is implementation dependent. In assembly language, when
164you add two numbers, you have access to a "carry flag" to see if an overflow
165occurred. C has no analogous operator to a carry flag.
166
167It is up to the user to decide if the results are valid for a given operation.
168In a 32 bit environment, the limit is so large that this is likely not an
169issue for most typical applications. However, it doesn't take much to overflow
170a 16 bit int so care should taken in a 16 bit environment.
171
172The reaction to an integer overflow/underflow is unknown at run-time:
173
174    o) adding 2 large positive numbers may silently roll over to a
175       negative number.
176    o) in some embedded applications an integer overflow/underflow triggers
177       a hardware exception.
178
179Since I don't have control over where this library could possibly run,
180I chose to ignore the problem for now. If anyone has some suggestions
181(that's portable), please let me know.
182
183**************************************************************************
184
185KNOWN BUGS : (Other than integer overflow discussed above....) None
186
187**************************************************************************
188
189IF YOU ARE IN A HURRY ...
190
191UNIX:  (assumes gcc compiler)
192
193run    make          (build library + 4 executables)
194
195run    make -f makefile.osx      (for MAC OSX)
196
197--- OR ---
198
199run:   mklib         (this will create the library, lib_mapm.a)
200
201run:   mkdemo        (this will create 4 executables, 'calc', 'validate',
202                      'primenum', and 'cpp_demo')
203
204
205DOS / Win NT/9x  (in a DOS window for NT/9x):
206
207see the file 'README.DOS' for instructions.
208
209
210**************************************************************************
211
212calc:  This is a command line version of an RPN calculator. If you are
213       familiar with RPN calculators, the use of this program will be
214       quite obvious. The default is 30 decimal places but this can be
215       changed with the '-d' command line option. This is not an
216       interactive program, it just computes an expression from the command
217       line. Run 'calc' with no arguments to get a list of the operators.
218
219       compute : (345.2 * 87.33) - (11.88 / 3.21E-2)
220
221       calc 345.2 87.33 x 11.88 3.21E-2 / -
222       result: 29776.22254205607476635514018692
223
224
225       compute PI to 70 decimal places :  (6 * arcsin(0.5))
226
227       calc -d70 6 .5 as x
228       result :
2293.1415926535897932384626433832795028841971693993751058209749445923078164
230
231       calc -d70 -1 ac             (arccos(-1) for fastest possible way)
232
233validate : This program will compare the MAPM math functions to the C
234	   standard library functions, like sqrt, sin, exp, etc.  This
235	   should give the user a good idea that the library is operating
236	   correctly. The program will also perform high precision math
237	   using known quantities like PI, log(2), etc.
238
239	   'validate' also attempts to obtain as much code coverage of the
240	   library as is practical. I used 'gcov' (available with the gcc
241	   distribution) to test the code coverage. 100% coverage is not
242	   obtained, compromises must be made in order to have the program
243	   run in a reasonable amount of time.
244
245primenum:  This program will generate the first 10 prime numbers starting
246	   with the number entered as an argument.
247
248	   Example:  primenum 1234567890 will output (actually 1 per line):
249		     this took ~3 sec on my Linux PC, a 350MHz PII.
250
251	   1234567891, 1234567907, 1234567913, 1234567927, 1234567949,
252	   1234567967, 1234567981, 1234568021, 1234568029, 1234568047
253
254
255**************************************************************************
256
257To use the library, simply include 'm_apm.h' and link your program
258with the library, libmapm.a (unix) or mapm.lib (dos).
259
260Note: If your system creates libraries with a '.a' extension, then the
261library will be named libmapm.a. (This conforms to typical unix naming
262conventions).
263
264Note: If your system creates libraries with a '.lib' extension, then the
265library will be named mapm.lib.
266
267For unix builds, you also may need to specify the math library (-lm) when
268linking. The reason is some of the MAPM functions use an iterative algorithm.
269When you use an iterative solution, you have to supply an initial guess. I
270use the standard math library to generate this initial guess. I debated
271whether this library should be stand-alone, i.e. generate it's own initial
272guesses with some algorithm. In the end, I decided using the standard math
273library would not be a big inconvienence and also it was just too tempting
274having an immediate 15 digits of precision. When you prime the iterative
275routines with 15 accurate digits, the MAPM functions converge faster.
276
277See the file 'algorithms.used' to see a description of the algorithms I
278used in the library. Some I derived on my own, others I borrowed from people
279smarter than me. Version 2 of the library supports a 'fast' multiplication
280algorithm. The algorithm used is described in the algorithms.used file. A
281considerable amount of time went into finding fast algorithms for the
282library. However, some could possibly be even better. If anyone has a
283more efficient algorithm for any of these functions, I would like to here
284from you.
285
286See the file 'function.ref' to see a description of all the functions in
287the library and the calling conventions, prototypes, etc.
288
289See the file 'struct.ref' which documents how I store the numbers internally
290in the MAPM data structure. This will not be needed for normal use, but it
291will be very useful if you need to change/add to the existing library.
292
293**************************************************************************
294
295USING MAPM IN A MULTI-THREADED APPLICATION :
296
297Note that the default MAPM library is NOT thread safe. MAPM internal data
298structures could get corrupted if multiple MAPM functions are active at the
299same time. The user should guarantee that only one thread is performing
300MAPM functions. This can usually be achieved by a call to the operating
301system to obtain a 'semaphore', 'mutex', or 'critical code section' so the
302operating system will guarantee that only one MAPM thread will be active
303at a time.
304
305The necessary function wrappers for thread safe operation can be found in
306the sub-directory 'multi_thread' (unix) or 'multithd' (Win/Dos). For now,
307only Microsoft Visual C++ 6.0 is supported.
308
309**************************************************************************
310
311QUICK TEMPLATE FOR NORMAL USE :
312
313The MAPM math is done on a new data type called "M_APM".  This is
314actually a pointer to a structure, but the contents of the structure
315should never be manipulated: all operations on MAPM entities are done
316through a functional interface.
317
318The MAPM routines will automatically allocate enough space in their
319results to hold the proper number of digits.
320
321The caller must initialize all MAPM values before the routines can
322operate on them (including the values intended to contain results of
323calculations).  Once this initialization is done, the user never needs
324to worry about resizing the MAPM values, as this is handled inside the
325MAPM routines and is totally invisible to the user.
326
327The result of a MAPM operation cannot be one of the other MAPM operands.
328If you want this to be the case, you must put the result into a
329temporary variable and then assign (m_apm_copy) it to the appropriate operand.
330
331All of the MAPM math functions begin with m_apm_*.
332
333There are some predefined constants for your use. In case it's not obvious,
334these should never appear as a 'result' parameter in a function call.
335
336The following constants are available : (declared in m_apm.h)
337
338MM_Zero             MM_One              MM_Two             MM_Three
339MM_Four             MM_Five             MM_Ten
340MM_PI               MM_HALF_PI          MM_2_PI            MM_E
341MM_LOG_E_BASE_10    MM_LOG_10_BASE_E    MM_LOG_2_BASE_E    MM_LOG_3_BASE_E
342
343
344The non-integer constants above (PI, log(2), etc) are accurate to 128
345decimal places. The file mapmcnst.c contains these constants. I've
346included 512 digit constants in this file also that are commented out.
347If you need more than 512 digits, you can simply use the 'calc' program
348to generate more precise constants (or create more precise constants at
349run-time in your app).  The number of significant digits in the constants
350should be 6-8 more than the value specified in the #define.
351
352
353Basic plan of attack:
354
355(1)  get your 'numbers' into M_APM format.
356(2)  do your high precision math.
357(3)  get the M_APM numbers back into a format you can use.
358
359
360--------  (1)  --------
361
362#include "m_apm.h"        /* must be included before any M_APM 'declares' */
363
364M_APM    area_mapm;                /* declare your variables */
365M_APM    tmp_mapm;
366M_APM    array_mapm[10];           /* can use normal array notation */
367M_APM    array2d_mapm[10][10];
368
369area_mapm = m_apm_init()           /* init your variables */
370tmp_mapm  = m_apm_init();
371
372for (i=0; i < M; i++)              /* must init every element of the array */
373  array_mapm[i] = m_apm_init();
374
375for (i=0; i < M; i++)
376   for (j=0; j < N; j++)
377     array2d_mapm[i][j] = m_apm_init();
378
379/*
380 *   there are 3 ways to convert your number into an M_APM number
381 *   (see the file function.ref)
382 *
383 *   a) literal string   (exponential notation OK)
384 *   b) long variable
385 *   c) double variable
386 */
387
388m_apm_set_string(tmp_mapm, "5.3286E-7");
389m_apm_set_long(array_mapm[6], -872253L);
390m_apm_set_double(array2d_mapm[3][7], -529.4486711);
391
392
393--------  (2)  --------
394
395do your math ...
396
397m_apm_add(cc_mapm, aa_mapm, MM_PI);
398m_apm_divide(bb_mapm, DECIMAL_PLACES, aa_mapm, MM_LOG_2_BASE_E);
399m_apm_sin(bb_mapm, DECIMAL_PLACES, aa_mapm);
400whatever ...
401
402
403--------  (3)  --------
404
405There are 5 total functions for converting an M_APM number into something
406useful.  (See the file 'function.ref' for full function descriptions)
407
408For these 5 functions, M_APM number -> string is the conversion
409
410===================
411====  METHOD 1 ==== : floating point values	(m_apm_to_string)
412===================
413
414the format will be in scientific (exponential) notation
415
416output string = [-]n.nnnnnE+x   or ...E-x
417
418where 'n' are digits and the exponent will be always be present,
419including E+0
420
421it's easy to convert this to a double:
422
423double dtmp = atof(out_buffer);
424
425
426===================
427====  METHOD 2 ==== : floating point values	(m_apm_to_fixpt_string)
428===================				(m_apm_to_fixpt_stringex)
429						(m_apm_to_fixpt_stringexp)
430the format will be in fixed point notation
431
432output string = [-]mmm.nnnnnn
433
434where 'm' & 'n' are digits.
435
436
437===================
438====  METHOD 3 ==== : integer values		(m_apm_to_integer_string)
439===================
440
441the format will simply be digits with a possible leading '-' sign.
442
443output string = [-]nnnnnn
444
445where 'n' are digits.
446
447it's easy to convert this to a long :
448long mtmp = atol(out_buffer);
449
450... or an int :
451int itmp = atoi(out_buffer);
452
453Note that if the M_APM number has a fractional portion, the fraction
454will be truncated and only the integer portion will be output.
455
456
457char  out_buffer[1024];
458
459m_apm_to_string(out_buffer, DECIMAL_PLACES, mapm_number);
460
461m_apm_to_fixpt_string(out_buffer, DECIMAL_PLACES, mapm_number);
462
463m_apm_to_integer_string(out_buffer, mapm_number);
464
465
466**************************************************************************
467
468*********************************************************
469****  NOTES on the fixed point formatting functions  ****
470*********************************************************
471
472Assume you have the following code:
473
474
475-->   m_apm_set_string(aa_mapm, "2.0E18");
476-->   m_apm_sqrt(bb_mapm, 40, aa_mapm);
477
478-->   m_apm_to_string(buffer, 40, bb_mapm);
479-->   fprintf(stdout,"[%s]\n",buffer);
480
481-->   m_apm_to_fixpt_string(buffer, 40, bb_mapm);
482-->   fprintf(stdout,"[%s]\n",buffer);
483
484
485It is desired to compute the sqrt(2.0E+18) to
48640 significant digits. You then want the result
487output with 40 decimal places. But the output
488from above is :
489
490[1.4142135623730950488016887242096980785697E+9]
491[1414213562.3730950488016887242096980785697000000000]
492
493
494Why are there 9 '0' in the fixed point formatted string??
495
496The sqrt calculation computed 40 significant digits relative
497to the number in EXPONENTIAL format. When the number is
498output in exponential format, the 40 digits are as expected
499with an exponent of 'E+9'.
500
501The same number formatted as fixed point appears to be an
502error.  Remember, we computed 40 significant digits.  However,
503the result has an exponent of '+9'. So, 9 of the digits are
504needed *before* the decimal point. In our calculation, only
50531 digits of precision remain from our original 40. We then
506asked the fixed point formatting function to format 40 digits.
507Only 31 are left so 9 zeros are used as pad at the end to
508fulfill the 40 places asked for.
509
510Keep this in mind if you truly desire more accurate results
511in fixed point formatting and your result contains a large
512positive exponent.
513
514**************************************************************************
515
516MAPM C++ WRAPPER CLASS:
517
518Orion Sky Lawlor (olawlor@acm.org) has added a very nice C++ wrapper
519class to m_apm.h. This C++ class will have no effect if you just use
520a normal C compiler. The library will operate as before with no user
521impacts.
522
523For now, I recommend compiling the library as 'C'. The library will
524compile as a C++ library, but then it is likely that you would not
525be able to use the library in a straight C application.  Since the C++
526wrapper class works very nicely as is, there is no pressing need to compile
527the library as C++.
528
529See the file 'cpp_function.ref' to see a description of how to use
530the MAPM class.
531
532To compile and link the C++ demo program: (assuming the library is
533already built)
534
535UNIX:
536
537g++ cpp_demo.cpp lib_mapm.a -s -o cpp_demo -lm
538
539GCC for DOS: (gxx (or g++) is the C++ compiler)
540
541gxx cpp_demo.cpp lib_mapm.a -s -o cpp_demo.exe -lm
542
543
544Using the C++ wrapper allows you to do things like:
545
546// Compute the factorial of the integer n
547
548MAPM factorial(MAPM n)
549{
550	MAPM i;
551	MAPM product = 1;
552	for (i=2; i <= n; i++)
553		product *= i;
554	return product;
555}
556
557
558The syntax is the same as if you were just writing normal code, but all
559the computations will be performed with the high precision math library,
560using the new 'datatype' MAPM.
561
562The default precision of the computations is as follows:
563
564Addition, subtraction, and multiplication will maintain ALL significant
565digits.
566
567All other operations (divide, sin, etc) will use the following rules:
568
5691) if the operation uses only one input value [y = sin(x)], the result 'y'
570   will be the same precision as 'x', with a minimum of 30 digits if 'x' is
571   less than 30 digits.
572
5732) if the operation uses two input values [z = atan2(y,x)], the result 'z'
574   will be the max digits of 'x' or 'y' with a minimum of 30.
575
576The default precision is 30 digits. You can change the precision at
577any time with the function 'm_apm_cpp_precision'. (See function.ref)
578
579---->        m_apm_cpp_precision(80);
580
581will result in all operations being accurate to a minimum of 80 significant
582digits. If any operand contains more than the minimum number of digits, then
583the result will contain the max number of digits of the operands.
584
585
586NOTE!: Some real life use with the C++ wrapper has revealed a certain
587       tendency for a program to become quite slow after many iterations
588       (like in a for/while loop).  After a little debug, the reason
589       became clear. Remember that multiplication will maintain ALL
590       significant digits :
591
592       20 digit number x 20 digit number =  40 digits
593       40 digit number x 40 digit number =  80 digits
594       80 digit number x 80 digit number = 160 digits
595       etc.
596
597       So after numerous iterations, the number of significant digits
598       was growing without bound. The easy way to fix the problem is
599       to simply *round* your result after a multiply or some other
600       complex operation. For example:
601
602       #define MAX_DIGITS 256
603
604       p1 = (p0 * sin(b1) * exp(1 + u1)) / sqrt(1 + b1);
605       p1 = p1.round(MAX_DIGITS);
606
607       If you 'round' as shown here, your program will likely be
608       nearly as fast as a straight 'C' version.
609
610
611NOTE #2!
612
613       Reference the following code snippet:
614
615       ...
616
617       MAPM  pi1, pi2;
618       char  obuf[256];
619
620       m_apm_cpp_precision(62);
621
622       pi1 = 2 * asin("1");
623       pi2 = 2 * asin(1.0);
624
625       pi1.toString(obuf, 60);   printf("PI1 = [%s] \n",obuf);
626       pi2.toString(obuf, 60);   printf("PI2 = [%s] \n",obuf);
627
628       ...
629
630       On my system, the output is :
631
632PI1 = [3.141592653589793238462643383279502884197169399375105820974945E+0]
633PI2 = [3.141592653589790000000000000000000000000000000000000000000000E+0]
634
635
636       PI2 only has 15 significant digits! This is due to how the second
637       asin is called. It is called with a 'double' as the argument, hence
638       the compiler will use the default double asin function from the
639       standard library. This is likely not the intent but this would be
640       easy to miss if this was a complex calculation and we didn't know
641       the 'right' answer.
642
643       In order to force the use of the overloaded MAPM functions, call the
644       MAPM functions with a quoted string as the argument (if the argument
645       is a constant and not a variable).
646
647       This would also work (though it seems less elegant ...) :
648
649       MAPM t = 1;
650       pi2 = 2 * asin(t);
651
652-----------
653
654If you have any questions or problems with the C++ wrapper, please let
655me know. I am not very C++ proficient, but I'd still like to know about any
656problems.  Orion Sky Lawlor (olawlor@acm.org) is the one who implemented
657the MAPM class, so he'll have to resolve any real hardcore problems, if you
658have any.
659
660**************************************************************************
661
662TESTING :
663
664Testing the library was interesting.  How do I know it's right?  Since I
665test the library against the standard C runtime math library (see below)
666I have high confidence that the MAPM library is giving correct results.
667The logic here is the basic algorithms are independent of the number of
668digits calculated, more digits just takes longer.
669
670The MAPM library has been tested in the following environments :
671
672Linux i486 / gcc 2.7.2.3, gcc 2.95.2
673Linux i686 / gcc 2.91.66, gcc 2.95.2, gcc 2.95.3, gcc 3.0.4, gcc 3.2.3
674Linux i686 / Intel Linux C/C++ Compiler Verison 7.0 / 8.1
675FreeBSD 4.1 / gcc 2.95.1
676FreeBSD 4.8 / gcc 2.95.4
677FreeBSD 5.x / gcc 2.95.2, gcc 2.95.3
678Redhat Linux 8.2 / gcc 3.2
679HP-UX 9.x /gcc 2.5.8
680HP-UX 10.x / gcc 2.95.2
681Sun 5.5.1  (output from uname), gcc 3.1.1
682Sun Solaris 2.6 / gcc 2.95.1, gcc 2.95.3, gcc 3.2.3
683MAC OSX / gcc ?
684DOS 5.0 using GCC 2.8.1 for DOS
685DOS 5.0 using GCC 2.95.2 for DOS
686DOS ??? using Borland Turbo C++ 3.0
687WIN NT+SP5 using Borland C++ 5.02 IDE, 5.2 & 5.5 command line.
688WIN 2000 using National Instruments LabWindows CVI 6.0
689WIN98 using Borland C++ 5.5 command line.
690WIN98 & NT & 2000 & XP using LCC-WIN32 Ver 3.2, 3.3
691WIN98 & NT using Watcom C 11.x
692WIN95 & NT using Open Watcom 1.0
693WIN95 & WIN98 using MINGW-32 version mingw-1.0.1-20010726
694WIN95 & WIN2000 using DEV-CPP 5.0 Beta 8, 4.9.8.0
695MINGW-32 with gcc 3.2   (mingw special 20020817-1)
696MINGW-32 with gcc 3.2.3 (mingw special 20030504-1)
697WINXP & MINGW-32 with gcc 3.4.5
698WINXP & Digital Mars Compiler 8.49
699WIN?? using Metrowerks CodeWarrior Pro 7.0 for Windows
700DOS 5.0 using Microsoft C 5.1 and 8.00c  (16 bit EXEs)
701WIN98 & NT using Microsoft Visual C++ 6.0
702WIN98 & NT using Microsoft Visual C++ 7.x (VS.NET 2003 except for
703                                           known compiler bug C2676)
704
705
706
707As a general rule, when calculating a quantity to a given number of decimal
708places, I calculated 4-6 extra digits and then rounded the result to what was
709asked for. I decided to be conservative and give a correct answer rather than
710to be faster and possibly have the last 2-3 digits in error. Also, some of
711the functions call other functions (calculating arc-cos will call cos, log
712will call exp, etc.) so I had to calculate a few extra digits in each iteration
713to guarantee the loops terminated correctly.
714
7151)  I debugged the 4 basic math operations. I threw numerous test cases at
716    each of the operations until I knew they were correct.
717
718    Also note that the math.h type functions all call the 4 basic operations
719    numerous times. So if all the math.h functions work, it is highly
720    probable the 4 basic math operations work also.
721
7222)  'math.h' type functions.
723
724     SQRT:     Not real hard to check. Single stepping through the iterative
725	       loop showed it was always converging to the sqrt.
726
727     CBRT:     Similar to sqrt, single stepping through the iterative loop
728	       showed it was always converging to the cube root.
729
730     EXP:      I wrote a separate algorithm which expanded the Taylor series
731	       manually and compared the results against the library.
732
733     POW:      Straightforward since this just calls 'EXP'.
734
735     LOG:      I wrote a separate algorithm which expanded the Taylor series
736	       manually and compared the results against the library. This
737	       took a long time to execute since the normal series converges
738	       VERY slowly for the log function. This is why the LOG function
739	       uses an iterative algorithm.
740
741     LOG10:    Straightforward since this just calls 'LOG'.
742
743     SIN/COS:  I wrote a separate algorithm which expanded the Taylor series
744	       manually and compared the results against the library.
745
746     TAN:      Straightforward since this just calls 'SIN' and 'COS'.
747
748     ARC-x:    Single stepping through the iterative loop showed the arc
749	       family of functions were always converging. Also used these
750	       to compute PI. The value of PI is now known to many, many
751	       digits. I computed PI to 1000+ digits by computing:
752
753	       6 * arcsin(0.5)  and  4 * arctan(1)  and  3 * arccos(0.5)
754
755	       and compared the output to the published 'real' values of PI.
756
757	       The arc family of functions exercise considerable portions
758	       of the library.
759
760  HYPERBOLIC:  The hyperbolic functions just use exp, log, and the 4 basic
761  	       math operations. All of these functions simply use other
762	       existing functions in the library.
763
764     FINALLY:  Run the program 'validate'. This program compares the first
765	       13-14 significant digits of the standard C library against
766	       the MAPM library. If this program passes, you can feel
767	       confident that the MAPM library is giving correct results.
768
769**************************************************************************
770