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