As a special service "Fossies" has tried to format the requested text file into HTML format (style: standard) with prefixed line numbers.
Alternatively you can here view or download the uninterpreted source code file.

1 ************************************************************************** 2 3 MAPM 4 5 Version 4.9.5a 6 7 February 21, 2010 8 9 Michael C. Ring 10 11 ringx004@umn.edu 12 13 Latest release will be available at 14 http://tc.umn.edu/~ringx004 15 16 *************************************************************************** 17 * * 18 * Copyright (C) 1999 - 2010 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 44 Mike's Arbitrary Precision Math Library is a set of functions that 45 allow the user to perform math to any level of accuracy that is 46 desired. The inspiration for this library was Lloyd Zusman's similar 47 APM package that was released in ~1988. I borrowed some of his ideas 48 in my implementation, creating a new data type (MAPM) and how the data 49 type was used by the user. However, there were a few things I wanted 50 my library to do that the original library did not : 51 52 1) 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 57 2) 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 64 3) 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 74 NOTE: 75 76 There is a COMPILER BUG in Microsoft's Visual C++ 7.x (VS.NET 2003) which 77 prevents a C++ MAPM application from compiling. 78 79 This only affects C++ applications. C applications are OK. 80 81 The compiler bug creates an error C2676 similar to this: 82 83 <path>...\mapm-49\M_APM.H(###) : error C2676: binary operator '-': 'MAPM' 84 doesn't define this operator or a conversion to a suitable type for the 85 predefined operator. 86 87 To work around this bug, go to http://www.microsoft.com 88 89 In the upper right corner of web page, search for "814455". 90 91 The results of the search will point you to an article on how to work 92 around the problem. 93 94 ************************************************************************** 95 96 NOTE: MAPM Library History can now be found in 'history.txt' 97 98 ************************************************************************** 99 100 ANOTHER NOTE: For the Windows/DOS distribution, the filename convention 101 will always be in 8.3 format. This is because there are some users who 102 still use 16 bit DOS .... 103 104 (I really wasn't sure what to call this library. 'Arbitrary Precision Math' 105 is a defacto standard for what this does, but that name was already taken, 106 so I just put an 'M' in front of it ...) 107 108 ************************************************************************** 109 110 MAPM LIBRARY NUMERICAL LIMITATIONS: 111 112 A general floating point number is of the form: 113 114 Sn.nnnnnnnnE+yyyy ex: -4.318384739357974916E+6215 115 Sn.nnnnnnnnE-yyyy ex: +8.208237913789131096523645193E-12873 116 117 'S' is the sign, + or -. 118 119 In MAPM, a number (n.nnnn...) may contain up to ( INT_MAX - 1 ) digits. 120 121 For example, an MAPM number with a 16 bit integer may contain 2^15 or 32,767 122 digits. An MAPM number with a 32 bit integer may contain 2^31 or 2,147,483,647 123 digits. All MAPM numbers are stored in RAM, there is no "data on disk" option. 124 So to represent the maximum number of digits on a 32 bit CPU will require 125 greater than 2 Gig of RAM. 126 127 If you have a CPU with 64 bit ints, then the limitation is 2^63 or 128 9,223,372,036,854,775,807. It should be a very long time before computers 129 have this much RAM in them. 130 131 For the exponent (yyyy), the limitations are also INT_MAX and INT_MIN. 132 133 For a 16 bit CPU, the largest number you can represent is approx 134 0.9999999....E+32767. (H) 135 136 For a 16 bit CPU, the smallest number you can represent (other than 0) 137 is approx 0.1E-32767. (L) 138 139 For a 32 bit CPU, the largest number you can represent is approx 140 0.9999999....E+2147483647. (H) 141 142 For a 32 bit CPU, the smallest number you can represent (other than 0) 143 is approx 0.1E-2147483647. (L) 144 145 The 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 156 MAPM can represent real numbers from -H to -L, 0.0, +L to +H. 157 158 The number of significant digits is typically only limited by available RAM. 159 160 In 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 162 detect integer overflow and underflow. Per K&R C, the result of integer 163 overflow/underflow is implementation dependent. In assembly language, when 164 you add two numbers, you have access to a "carry flag" to see if an overflow 165 occurred. C has no analogous operator to a carry flag. 166 167 It is up to the user to decide if the results are valid for a given operation. 168 In a 32 bit environment, the limit is so large that this is likely not an 169 issue for most typical applications. However, it doesn't take much to overflow 170 a 16 bit int so care should taken in a 16 bit environment. 171 172 The 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 179 Since I don't have control over where this library could possibly run, 180 I chose to ignore the problem for now. If anyone has some suggestions 181 (that's portable), please let me know. 182 183 ************************************************************************** 184 185 KNOWN BUGS : (Other than integer overflow discussed above....) None 186 187 ************************************************************************** 188 189 IF YOU ARE IN A HURRY ... 190 191 UNIX: (assumes gcc compiler) 192 193 run make (build library + 4 executables) 194 195 run make -f makefile.osx (for MAC OSX) 196 197 --- OR --- 198 199 run: mklib (this will create the library, lib_mapm.a) 200 201 run: mkdemo (this will create 4 executables, 'calc', 'validate', 202 'primenum', and 'cpp_demo') 203 204 205 DOS / Win NT/9x (in a DOS window for NT/9x): 206 207 see the file 'README.DOS' for instructions. 208 209 210 ************************************************************************** 211 212 calc: 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 : 229 3.1415926535897932384626433832795028841971693993751058209749445923078164 230 231 calc -d70 -1 ac (arccos(-1) for fastest possible way) 232 233 validate : 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 245 primenum: 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 257 To use the library, simply include 'm_apm.h' and link your program 258 with the library, libmapm.a (unix) or mapm.lib (dos). 259 260 Note: If your system creates libraries with a '.a' extension, then the 261 library will be named libmapm.a. (This conforms to typical unix naming 262 conventions). 263 264 Note: If your system creates libraries with a '.lib' extension, then the 265 library will be named mapm.lib. 266 267 For unix builds, you also may need to specify the math library (-lm) when 268 linking. The reason is some of the MAPM functions use an iterative algorithm. 269 When you use an iterative solution, you have to supply an initial guess. I 270 use the standard math library to generate this initial guess. I debated 271 whether this library should be stand-alone, i.e. generate it's own initial 272 guesses with some algorithm. In the end, I decided using the standard math 273 library would not be a big inconvienence and also it was just too tempting 274 having an immediate 15 digits of precision. When you prime the iterative 275 routines with 15 accurate digits, the MAPM functions converge faster. 276 277 See the file 'algorithms.used' to see a description of the algorithms I 278 used in the library. Some I derived on my own, others I borrowed from people 279 smarter than me. Version 2 of the library supports a 'fast' multiplication 280 algorithm. The algorithm used is described in the algorithms.used file. A 281 considerable amount of time went into finding fast algorithms for the 282 library. However, some could possibly be even better. If anyone has a 283 more efficient algorithm for any of these functions, I would like to here 284 from you. 285 286 See the file 'function.ref' to see a description of all the functions in 287 the library and the calling conventions, prototypes, etc. 288 289 See the file 'struct.ref' which documents how I store the numbers internally 290 in the MAPM data structure. This will not be needed for normal use, but it 291 will be very useful if you need to change/add to the existing library. 292 293 ************************************************************************** 294 295 USING MAPM IN A MULTI-THREADED APPLICATION : 296 297 Note that the default MAPM library is NOT thread safe. MAPM internal data 298 structures could get corrupted if multiple MAPM functions are active at the 299 same time. The user should guarantee that only one thread is performing 300 MAPM functions. This can usually be achieved by a call to the operating 301 system to obtain a 'semaphore', 'mutex', or 'critical code section' so the 302 operating system will guarantee that only one MAPM thread will be active 303 at a time. 304 305 The necessary function wrappers for thread safe operation can be found in 306 the sub-directory 'multi_thread' (unix) or 'multithd' (Win/Dos). For now, 307 only Microsoft Visual C++ 6.0 is supported. 308 309 ************************************************************************** 310 311 QUICK TEMPLATE FOR NORMAL USE : 312 313 The MAPM math is done on a new data type called "M_APM". This is 314 actually a pointer to a structure, but the contents of the structure 315 should never be manipulated: all operations on MAPM entities are done 316 through a functional interface. 317 318 The MAPM routines will automatically allocate enough space in their 319 results to hold the proper number of digits. 320 321 The caller must initialize all MAPM values before the routines can 322 operate on them (including the values intended to contain results of 323 calculations). Once this initialization is done, the user never needs 324 to worry about resizing the MAPM values, as this is handled inside the 325 MAPM routines and is totally invisible to the user. 326 327 The result of a MAPM operation cannot be one of the other MAPM operands. 328 If you want this to be the case, you must put the result into a 329 temporary variable and then assign (m_apm_copy) it to the appropriate operand. 330 331 All of the MAPM math functions begin with m_apm_*. 332 333 There are some predefined constants for your use. In case it's not obvious, 334 these should never appear as a 'result' parameter in a function call. 335 336 The following constants are available : (declared in m_apm.h) 337 338 MM_Zero MM_One MM_Two MM_Three 339 MM_Four MM_Five MM_Ten 340 MM_PI MM_HALF_PI MM_2_PI MM_E 341 MM_LOG_E_BASE_10 MM_LOG_10_BASE_E MM_LOG_2_BASE_E MM_LOG_3_BASE_E 342 343 344 The non-integer constants above (PI, log(2), etc) are accurate to 128 345 decimal places. The file mapmcnst.c contains these constants. I've 346 included 512 digit constants in this file also that are commented out. 347 If you need more than 512 digits, you can simply use the 'calc' program 348 to generate more precise constants (or create more precise constants at 349 run-time in your app). The number of significant digits in the constants 350 should be 6-8 more than the value specified in the #define. 351 352 353 Basic 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 364 M_APM area_mapm; /* declare your variables */ 365 M_APM tmp_mapm; 366 M_APM array_mapm[10]; /* can use normal array notation */ 367 M_APM array2d_mapm[10][10]; 368 369 area_mapm = m_apm_init() /* init your variables */ 370 tmp_mapm = m_apm_init(); 371 372 for (i=0; i < M; i++) /* must init every element of the array */ 373 array_mapm[i] = m_apm_init(); 374 375 for (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 388 m_apm_set_string(tmp_mapm, "5.3286E-7"); 389 m_apm_set_long(array_mapm[6], -872253L); 390 m_apm_set_double(array2d_mapm[3][7], -529.4486711); 391 392 393 -------- (2) -------- 394 395 do your math ... 396 397 m_apm_add(cc_mapm, aa_mapm, MM_PI); 398 m_apm_divide(bb_mapm, DECIMAL_PLACES, aa_mapm, MM_LOG_2_BASE_E); 399 m_apm_sin(bb_mapm, DECIMAL_PLACES, aa_mapm); 400 whatever ... 401 402 403 -------- (3) -------- 404 405 There are 5 total functions for converting an M_APM number into something 406 useful. (See the file 'function.ref' for full function descriptions) 407 408 For these 5 functions, M_APM number -> string is the conversion 409 410 =================== 411 ==== METHOD 1 ==== : floating point values (m_apm_to_string) 412 =================== 413 414 the format will be in scientific (exponential) notation 415 416 output string = [-]n.nnnnnE+x or ...E-x 417 418 where 'n' are digits and the exponent will be always be present, 419 including E+0 420 421 it's easy to convert this to a double: 422 423 double 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) 430 the format will be in fixed point notation 431 432 output string = [-]mmm.nnnnnn 433 434 where 'm' & 'n' are digits. 435 436 437 =================== 438 ==== METHOD 3 ==== : integer values (m_apm_to_integer_string) 439 =================== 440 441 the format will simply be digits with a possible leading '-' sign. 442 443 output string = [-]nnnnnn 444 445 where 'n' are digits. 446 447 it's easy to convert this to a long : 448 long mtmp = atol(out_buffer); 449 450 ... or an int : 451 int itmp = atoi(out_buffer); 452 453 Note that if the M_APM number has a fractional portion, the fraction 454 will be truncated and only the integer portion will be output. 455 456 457 char out_buffer[1024]; 458 459 m_apm_to_string(out_buffer, DECIMAL_PLACES, mapm_number); 460 461 m_apm_to_fixpt_string(out_buffer, DECIMAL_PLACES, mapm_number); 462 463 m_apm_to_integer_string(out_buffer, mapm_number); 464 465 466 ************************************************************************** 467 468 ********************************************************* 469 **** NOTES on the fixed point formatting functions **** 470 ********************************************************* 471 472 Assume 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 485 It is desired to compute the sqrt(2.0E+18) to 486 40 significant digits. You then want the result 487 output with 40 decimal places. But the output 488 from above is : 489 490 [1.4142135623730950488016887242096980785697E+9] 491 [1414213562.3730950488016887242096980785697000000000] 492 493 494 Why are there 9 '0' in the fixed point formatted string?? 495 496 The sqrt calculation computed 40 significant digits relative 497 to the number in EXPONENTIAL format. When the number is 498 output in exponential format, the 40 digits are as expected 499 with an exponent of 'E+9'. 500 501 The same number formatted as fixed point appears to be an 502 error. Remember, we computed 40 significant digits. However, 503 the result has an exponent of '+9'. So, 9 of the digits are 504 needed *before* the decimal point. In our calculation, only 505 31 digits of precision remain from our original 40. We then 506 asked the fixed point formatting function to format 40 digits. 507 Only 31 are left so 9 zeros are used as pad at the end to 508 fulfill the 40 places asked for. 509 510 Keep this in mind if you truly desire more accurate results 511 in fixed point formatting and your result contains a large 512 positive exponent. 513 514 ************************************************************************** 515 516 MAPM C++ WRAPPER CLASS: 517 518 Orion Sky Lawlor (olawlor@acm.org) has added a very nice C++ wrapper 519 class to m_apm.h. This C++ class will have no effect if you just use 520 a normal C compiler. The library will operate as before with no user 521 impacts. 522 523 For now, I recommend compiling the library as 'C'. The library will 524 compile as a C++ library, but then it is likely that you would not 525 be able to use the library in a straight C application. Since the C++ 526 wrapper class works very nicely as is, there is no pressing need to compile 527 the library as C++. 528 529 See the file 'cpp_function.ref' to see a description of how to use 530 the MAPM class. 531 532 To compile and link the C++ demo program: (assuming the library is 533 already built) 534 535 UNIX: 536 537 g++ cpp_demo.cpp lib_mapm.a -s -o cpp_demo -lm 538 539 GCC for DOS: (gxx (or g++) is the C++ compiler) 540 541 gxx cpp_demo.cpp lib_mapm.a -s -o cpp_demo.exe -lm 542 543 544 Using the C++ wrapper allows you to do things like: 545 546 // Compute the factorial of the integer n 547 548 MAPM 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 558 The syntax is the same as if you were just writing normal code, but all 559 the computations will be performed with the high precision math library, 560 using the new 'datatype' MAPM. 561 562 The default precision of the computations is as follows: 563 564 Addition, subtraction, and multiplication will maintain ALL significant 565 digits. 566 567 All other operations (divide, sin, etc) will use the following rules: 568 569 1) 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 573 2) 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 576 The default precision is 30 digits. You can change the precision at 577 any time with the function 'm_apm_cpp_precision'. (See function.ref) 578 579 ----> m_apm_cpp_precision(80); 580 581 will result in all operations being accurate to a minimum of 80 significant 582 digits. If any operand contains more than the minimum number of digits, then 583 the result will contain the max number of digits of the operands. 584 585 586 NOTE!: 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 611 NOTE #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 632 PI1 = [3.141592653589793238462643383279502884197169399375105820974945E+0] 633 PI2 = [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 654 If you have any questions or problems with the C++ wrapper, please let 655 me know. I am not very C++ proficient, but I'd still like to know about any 656 problems. Orion Sky Lawlor (olawlor@acm.org) is the one who implemented 657 the MAPM class, so he'll have to resolve any real hardcore problems, if you 658 have any. 659 660 ************************************************************************** 661 662 TESTING : 663 664 Testing the library was interesting. How do I know it's right? Since I 665 test the library against the standard C runtime math library (see below) 666 I have high confidence that the MAPM library is giving correct results. 667 The logic here is the basic algorithms are independent of the number of 668 digits calculated, more digits just takes longer. 669 670 The MAPM library has been tested in the following environments : 671 672 Linux i486 / gcc 2.7.2.3, gcc 2.95.2 673 Linux i586 / gcc 2.95.3, gcc 3.3.6 674 Linux i686 / gcc 2.91.66, gcc 2.95.2, gcc 2.95.3, gcc 3.0.4, gcc 3.2.3 675 Linux i686 / gcc 3.3.6, gcc 3.4.6, gcc 4.2.4 676 Linux i686 / Intel Linux C/C++ Compiler Verison 7.0 / 8.1 677 FreeBSD 4.1 / gcc 2.95.1 678 FreeBSD 4.8 / gcc 2.95.4 679 FreeBSD 5.x / gcc 2.95.2, gcc 2.95.3 680 Redhat Linux 8.2 / gcc 3.2 681 RHEL4 / gcc 3.4.3 682 HP-UX 9.x /gcc 2.5.8 683 HP-UX 10.x / gcc 2.95.2 684 Sun 5.5.1 (output from uname), gcc 3.1.1 685 Sun Solaris 2.6 / gcc 2.95.1, gcc 2.95.3, gcc 3.2.3 686 MAC OSX / gcc ? 687 DOS 5.0 using GCC 2.8.1 for DOS 688 DOS 5.0 using GCC 2.95.2 for DOS 689 DOS ??? using Borland Turbo C++ 3.0 690 WIN NT+SP5 using Borland C++ 5.02 IDE, 5.2 & 5.5 command line. 691 WIN 2000 using National Instruments LabWindows CVI 6.0 692 WIN98 using Borland C++ 5.5 command line. 693 WIN98 & NT & 2000 & XP using LCC-WIN32 Ver 3.2, 3.3 694 WIN98 & NT using Watcom C 11.x 695 WIN95 & NT using Open Watcom 1.0 696 WIN95 & WIN98 using MINGW-32 version mingw-1.0.1-20010726 697 WIN95 & WIN2000 using DEV-CPP 5.0 Beta 8, 4.9.8.0 698 MINGW-32 with gcc 3.2 (mingw special 20020817-1) 699 MINGW-32 with gcc 3.2.3 (mingw special 20030504-1) 700 WINXP & MINGW-32 with gcc 3.4.5 701 WINXP & Digital Mars Compiler 8.49 702 WIN?? using Metrowerks CodeWarrior Pro 7.0 for Windows 703 DOS 5.0 using Microsoft C 5.1 and 8.00c (16 bit EXEs) 704 WIN98 & NT using Microsoft Visual C++ 6.0 705 WIN98 & NT using Microsoft Visual C++ 7.x (VS.NET 2003 except for 706 known compiler bug C2676) 707 708 709 710 As a general rule, when calculating a quantity to a given number of decimal 711 places, I calculated 4-6 extra digits and then rounded the result to what was 712 asked for. I decided to be conservative and give a correct answer rather than 713 to be faster and possibly have the last 2-3 digits in error. Also, some of 714 the functions call other functions (calculating arc-cos will call cos, log 715 will call exp, etc.) so I had to calculate a few extra digits in each iteration 716 to guarantee the loops terminated correctly. 717 718 1) I debugged the 4 basic math operations. I threw numerous test cases at 719 each of the operations until I knew they were correct. 720 721 Also note that the math.h type functions all call the 4 basic operations 722 numerous times. So if all the math.h functions work, it is highly 723 probable the 4 basic math operations work also. 724 725 2) 'math.h' type functions. 726 727 SQRT: Not real hard to check. Single stepping through the iterative 728 loop showed it was always converging to the sqrt. 729 730 CBRT: Similar to sqrt, single stepping through the iterative loop 731 showed it was always converging to the cube root. 732 733 EXP: I wrote a separate algorithm which expanded the Taylor series 734 manually and compared the results against the library. 735 736 POW: Straightforward since this just calls 'EXP'. 737 738 LOG: I wrote a separate algorithm which expanded the Taylor series 739 manually and compared the results against the library. This 740 took a long time to execute since the normal series converges 741 VERY slowly for the log function. This is why the LOG function 742 uses an iterative algorithm. 743 744 LOG10: Straightforward since this just calls 'LOG'. 745 746 SIN/COS: I wrote a separate algorithm which expanded the Taylor series 747 manually and compared the results against the library. 748 749 TAN: Straightforward since this just calls 'SIN' and 'COS'. 750 751 ARC-x: Single stepping through the iterative loop showed the arc 752 family of functions were always converging. Also used these 753 to compute PI. The value of PI is now known to many, many 754 digits. I computed PI to 1000+ digits by computing: 755 756 6 * arcsin(0.5) and 4 * arctan(1) and 3 * arccos(0.5) 757 758 and compared the output to the published 'real' values of PI. 759 760 The arc family of functions exercise considerable portions 761 of the library. 762 763 HYPERBOLIC: The hyperbolic functions just use exp, log, and the 4 basic 764 math operations. All of these functions simply use other 765 existing functions in the library. 766 767 FINALLY: Run the program 'validate'. This program compares the first 768 13-14 significant digits of the standard C library against 769 the MAPM library. If this program passes, you can feel 770 confident that the MAPM library is giving correct results. 771 772 **************************************************************************