"Fossies" - the Fresh Open Source Software Archive 
Member "mapm_4.9.5a/validate.c" (21 Feb 2010, 63752 Bytes) of package /linux/misc/old/mapm-4.9.5a.tar.gz:
As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style:
standard) with prefixed line numbers and
code folding option.
Alternatively you can here
view or
download the uninterpreted source code file.
For more information about "validate.c" see the
Fossies "Dox" file reference documentation.
1
2 /*
3 * M_APM - validate.c
4 *
5 * Copyright (C) 1999 - 2007 Michael C. Ring
6 *
7 * Permission to use, copy, and distribute this software and its
8 * documentation for any purpose with or without fee is hereby granted,
9 * provided that the above copyright notice appear in all copies and
10 * that both that copyright notice and this permission notice appear
11 * in supporting documentation.
12 *
13 * Permission to modify the software is granted. Permission to distribute
14 * the modified code is granted. Modifications are to be distributed by
15 * using the file 'license.txt' as a template to modify the file header.
16 * 'license.txt' is available in the official MAPM distribution.
17 *
18 * This software is provided "as is" without express or implied warranty.
19 */
20
21 /*
22 * $Id: validate.c,v 2.6 2007/12/03 02:01:56 mike Exp $
23 *
24 * This file contains the validation test program. It compares the
25 * M_APM library to the standard C lbrary math functions. It also
26 * does a number of high precision calculations and compares the
27 * results to known quantities, like PI, log(2), etc.
28 *
29 * $Log: validate.c,v $
30 * Revision 2.6 2007/12/03 02:01:56 mike
31 * Update license
32 *
33 * Revision 2.5 2003/12/29 05:16:35 mike
34 * create exp_2 to exercise new code in _exp
35 *
36 * Revision 2.4 2003/07/28 19:37:47 mike
37 * change 16 bit constant
38 *
39 * Revision 2.3 2003/07/26 19:30:21 mike
40 * end one test early in 16 bit environments.
41 *
42 * Revision 2.2 2003/07/25 18:08:53 mike
43 * fix pow test when x==0
44 *
45 * Revision 2.1 2003/07/24 21:11:48 mike
46 * update comments
47 *
48 * Revision 2.0 2003/07/24 21:00:49 mike
49 * initial version
50 */
51
52 #include <stdio.h>
53 #include <stdlib.h>
54 #include <string.h>
55 #include <math.h>
56 #include "m_apm.h"
57
58 extern void factorial_local(M_APM, M_APM);
59 extern double cbrt_local(double);
60 extern double asinh_local(double);
61 extern double acosh_local(double);
62 extern double atanh_local(double);
63 extern M_APM M_get_stack_var(void);
64 extern void M_restore_stack(int);
65 extern void M_log_AGM_R_func(M_APM, int, M_APM, M_APM);
66 extern int M_get_sizeof_int(void);
67
68 #define DECIMAL_PLACES 22
69
70 /****************************************************************************/
71
72 void library_version_check()
73 {
74 char version_info[80];
75
76 m_apm_lib_short_version(version_info);
77 m_apm_lib_version(version_info);
78
79 fprintf(stdout,"\n%s\n\n",version_info);
80
81 /*
82 * check MAPM version info from m_apm.h to what was compiled into
83 * the library. This check could be used in your application to
84 * confirm that you are using the correct include file (m_apm.h).
85 */
86
87 if (strcmp(version_info, MAPM_LIB_VERSION) != 0)
88 {
89 fprintf(stdout,
90 "ERROR: MAPM Library compiled with different include file\n");
91
92 fprintf(stdout, "Include file version : %s \n", MAPM_LIB_SHORT_VERSION);
93 fprintf(stdout, "Library file version : %s \n",
94 m_apm_lib_short_version(version_info));
95 exit(2);
96 }
97 }
98
99 /****************************************************************************/
100
101 /*
102 * in case some compilers don't have a native cube root or
103 * inverse hyperbolic functions, we'll just make our own.
104 *
105 * note that we are not doing any decent error checking on
106 * the inputs, these functions are here just to support this
107 * program.
108 */
109
110 double cbrt_local(double x)
111 {
112 if (x == 0.0)
113 return(0.0);
114 else
115 {
116 if (x < 0.0)
117 return(-pow(-x, 0.333333333333333333));
118 else
119 return(pow(x, 0.333333333333333333));
120 }
121 }
122
123 /****************************************************************************/
124
125 double asinh_local(double x)
126 {
127 return(log(x + sqrt(x * x + 1.0)));
128 }
129
130 /****************************************************************************/
131
132 double acosh_local(double x)
133 {
134 if (x >= 1.0)
135 return(log(x + sqrt(x * x - 1.0)));
136 else
137 return(0.0);
138 }
139
140 /****************************************************************************/
141
142 double atanh_local(double x)
143 {
144 if (fabs(x) < 1.0)
145 return(0.5 * log((1.0 + x) / (1.0 - x)));
146 else
147 return(0.0);
148 }
149
150 /****************************************************************************/
151 /*
152 *
153 * slow but obvious version for testing
154 *
155 */
156 void factorial_local(M_APM result, M_APM ainput)
157 {
158 M_APM ftmp1, ftmp2, ftmp3;
159
160 if (m_apm_compare(ainput, MM_One) <= 0)
161 {
162 m_apm_copy(result, MM_One);
163 return;
164 }
165
166 ftmp1 = m_apm_init();
167 ftmp2 = m_apm_init();
168 ftmp3 = m_apm_init();
169
170 m_apm_copy(result, ainput);
171 m_apm_copy(ftmp1, ainput);
172
173 while (1)
174 {
175 m_apm_subtract(ftmp2, ftmp1, MM_One);
176 m_apm_multiply(ftmp3, ftmp2, result);
177
178 if (m_apm_compare(ftmp2, MM_Two) <= 0)
179 {
180 m_apm_copy(result, ftmp3);
181 break;
182 }
183
184 m_apm_subtract(ftmp1, ftmp2, MM_One);
185 m_apm_multiply(result, ftmp3, ftmp1);
186
187 if (m_apm_compare(ftmp1, MM_Two) <= 0)
188 break;
189 }
190
191 m_apm_free(ftmp3);
192 m_apm_free(ftmp2);
193 m_apm_free(ftmp1);
194 }
195
196 /****************************************************************************/
197
198 /*
199 * Traditional GCD from Knuth, The Art of Computer Programming:
200 *
201 * To compute GCD(u,v)
202 *
203 * A1:
204 * if (v == 0) return (u)
205 * A2:
206 * t = u mod v
207 * u = v
208 * v = t
209 * goto A1
210 */
211
212 void gcd_local(M_APM r, M_APM u, M_APM v)
213 {
214 M_APM tmpD, tmpN, tmpU, tmpV;
215
216 tmpD = M_get_stack_var();
217 tmpN = M_get_stack_var();
218 tmpU = M_get_stack_var();
219 tmpV = M_get_stack_var();
220
221 m_apm_absolute_value(tmpU, u);
222 m_apm_absolute_value(tmpV, v);
223
224 while (1)
225 {
226 if (tmpV->m_apm_sign == 0)
227 break;
228
229 m_apm_integer_div_rem(tmpD, tmpN, tmpU, tmpV);
230 m_apm_copy(tmpU, tmpV);
231 m_apm_copy(tmpV, tmpN);
232 }
233
234 m_apm_copy(r, tmpU);
235 M_restore_stack(4);
236 }
237
238 /****************************************************************************/
239
240 /*
241 * calculate log using an AGM algorithm.
242 */
243
244 void log_local(M_APM outval, int places, M_APM inval)
245 {
246 int dplaces;
247 M_APM tmp6, tmp7, tmp8, tmp9;
248
249 dplaces = places + 8;
250
251 tmp6 = M_get_stack_var();
252 tmp7 = M_get_stack_var();
253 tmp8 = M_get_stack_var();
254 tmp9 = M_get_stack_var();
255
256 m_apm_copy(tmp7, MM_One);
257 tmp7->m_apm_exponent = -places;
258
259 M_log_AGM_R_func(tmp8, dplaces, MM_One, tmp7);
260
261 m_apm_divide(tmp6, dplaces, tmp7, inval);
262
263 M_log_AGM_R_func(tmp9, dplaces, MM_One, tmp6);
264
265 m_apm_subtract(tmp6, tmp9, tmp8);
266 m_apm_round(outval, places, tmp6);
267
268 M_restore_stack(4);
269 }
270
271 /****************************************************************************/
272
273 void mapm_show_random()
274 {
275 int i, j;
276 char buffer2[32];
277 M_APM mapm1;
278
279 mapm1 = m_apm_init();
280
281 fprintf(stdout,"Output the first 10 random numbers ...\n");
282
283 for (i=1; i <= 10; i++)
284 {
285 m_apm_get_random(mapm1);
286 m_apm_to_string(buffer2, 14, mapm1);
287
288 fprintf(stdout,"Random Num %2d : %s \n",i,buffer2);
289 }
290
291 for (j=0; j < 2; j++)
292 {
293 fprintf(stdout,
294 "Set random seed to 2003 and output next 5 random numbers\n");
295
296 m_apm_set_random_seed("2003");
297
298 for (i=1; i <= 5; i++)
299 {
300 m_apm_get_random(mapm1);
301 m_apm_to_string(buffer2, 14, mapm1);
302
303 fprintf(stdout,"Random Num %2d : %s \n",i,buffer2);
304 }
305 }
306
307 m_apm_free(mapm1);
308 }
309
310 /****************************************************************************/
311
312 int mapm_test_sqrt()
313 {
314 int pass;
315 double xx, xx1, xdelta, yy, yy_mapm_lib, ydiff, tolerance;
316 char buffer[32], cbuf[32];
317 M_APM mapm1, mapm2;
318
319 fprintf(stdout,"Validating the SQRT function ... \n");
320
321 mapm1 = m_apm_init();
322 mapm2 = m_apm_init();
323
324 xx = 3.184E-7;
325 xdelta = 2.83576;
326 tolerance = 1.0E-13;
327 pass = 0;
328
329 while (1)
330 {
331 sprintf(cbuf, "%.6E", xx);
332 xx1 = atof(cbuf);
333
334 yy = sqrt(xx1);
335
336 m_apm_set_string(mapm1, cbuf);
337 m_apm_sqrt(mapm2, DECIMAL_PLACES, mapm1);
338 m_apm_to_string(buffer, DECIMAL_PLACES, mapm2);
339 yy_mapm_lib = atof(buffer);
340
341 if (fabs(yy) > 1.0)
342 ydiff = (yy_mapm_lib - yy) / yy;
343 else
344 ydiff = yy_mapm_lib - yy;
345
346 if (fabs(ydiff) > tolerance)
347 {
348 pass = 1;
349 fprintf(stdout, "x = %17.10E C-math-lib = %.15E \n", xx1, yy);
350 fprintf(stdout, "DIFF = %17.10E M_APM_LIB = %s \n", ydiff, buffer);
351 }
352
353 if (xx > 3.4E6)
354 break;
355
356 xx *= xdelta;
357 }
358
359 if (pass == 0)
360 fprintf(stdout,"... SQRT function passes\n");
361
362 m_apm_free(mapm1);
363 m_apm_free(mapm2);
364
365 return(pass);
366 }
367
368 /****************************************************************************/
369
370 int mapm_test_cbrt()
371 {
372 int sflag, pass;
373 double xx, xx1, xdelta, yy, yy_mapm_lib, ydiff, tolerance;
374 char buffer[32], cbuf[32];
375 M_APM mapm1, mapm2;
376
377 fprintf(stdout,"Validating the CBRT function ... \n");
378
379 mapm1 = m_apm_init();
380 mapm2 = m_apm_init();
381
382 xx = 3.9521E-8;
383 xdelta = 2.71329;
384 tolerance = 1.0E-13;
385 pass = 0;
386 sflag = 1;
387
388 while (1)
389 {
390 sprintf(cbuf, "%.6E", (sflag * xx));
391 xx1 = atof(cbuf);
392
393 yy = cbrt_local(xx1);
394
395 m_apm_set_string(mapm1, cbuf);
396 m_apm_cbrt(mapm2, DECIMAL_PLACES, mapm1);
397 m_apm_to_string(buffer, DECIMAL_PLACES, mapm2);
398 yy_mapm_lib = atof(buffer);
399
400 if (fabs(yy) > 1.0)
401 ydiff = (yy_mapm_lib - yy) / yy;
402 else
403 ydiff = yy_mapm_lib - yy;
404
405 if (fabs(ydiff) > tolerance)
406 {
407 pass = 1;
408 fprintf(stdout, "x = %17.10E C-math-lib = %.15E \n", xx1, yy);
409 fprintf(stdout, "DIFF = %17.10E M_APM_LIB = %s \n", ydiff, buffer);
410 }
411
412 if (xx > 4.1E7)
413 break;
414
415 xx *= xdelta;
416 sflag = -sflag;
417 }
418
419 if (pass == 0)
420 fprintf(stdout,"... CBRT function passes\n");
421
422 m_apm_free(mapm1);
423 m_apm_free(mapm2);
424
425 return(pass);
426 }
427
428 /****************************************************************************/
429
430 int mapm_test_sin()
431 {
432 int pass;
433 double xx, xdelta, yy, yy_mapm_lib, ydiff, tolerance;
434 char buffer[32];
435 M_APM mapm1, mapm2;
436
437 fprintf(stdout, "Validating the SIN function ... \n");
438
439 mapm1 = m_apm_init();
440 mapm2 = m_apm_init();
441
442 xx = -7.214;
443 xdelta = 0.817;
444 tolerance = 1.0E-13;
445 pass = 0;
446
447 while (1)
448 {
449 yy = sin(xx);
450
451 m_apm_set_double(mapm1, xx);
452 m_apm_sin(mapm2, DECIMAL_PLACES, mapm1);
453 m_apm_to_string(buffer, DECIMAL_PLACES, mapm2);
454 yy_mapm_lib = atof(buffer);
455
456 ydiff = yy_mapm_lib - yy;
457
458 if (fabs(ydiff) > tolerance)
459 {
460 pass = 1;
461 fprintf(stdout,"x = %17.10E C-math-lib = %.15E \n", xx, yy);
462 fprintf(stdout,"DIFF = %17.10E M_APM_LIB = %s \n", ydiff, buffer);
463 }
464
465 xx += xdelta;
466 if (xx > 7.21)
467 break;
468 }
469
470 if (pass == 0)
471 fprintf(stdout,"... SIN function passes\n");
472
473 m_apm_free(mapm1);
474 m_apm_free(mapm2);
475
476 return(pass);
477 }
478
479 /****************************************************************************/
480
481 int mapm_test_cos()
482 {
483 int pass;
484 double xx, xdelta, yy, yy_mapm_lib, ydiff, tolerance;
485 char buffer[32];
486 M_APM mapm1, mapm2;
487
488 fprintf(stdout, "Validating the COS function ... \n");
489
490 mapm1 = m_apm_init();
491 mapm2 = m_apm_init();
492
493 xx = -7.315;
494 xdelta = 0.813;
495 tolerance = 1.0E-13;
496 pass = 0;
497
498 while (1)
499 {
500 yy = cos(xx);
501
502 m_apm_set_double(mapm1, xx);
503 m_apm_cos(mapm2, DECIMAL_PLACES, mapm1);
504 m_apm_to_string(buffer, DECIMAL_PLACES, mapm2);
505 yy_mapm_lib = atof(buffer);
506
507 ydiff = yy_mapm_lib - yy;
508
509 if (fabs(ydiff) > tolerance)
510 {
511 pass = 1;
512 fprintf(stdout,"x = %17.10E C-math-lib = %.15E \n", xx, yy);
513 fprintf(stdout,"DIFF = %17.10E M_APM_LIB = %s \n", ydiff, buffer);
514 }
515
516 xx += xdelta;
517 if (xx > 7.21)
518 break;
519 }
520
521 if (pass == 0)
522 fprintf(stdout,"... COS function passes\n");
523
524 m_apm_free(mapm1);
525 m_apm_free(mapm2);
526
527 return(pass);
528 }
529
530 /****************************************************************************/
531
532 int mapm_test_sin_cos()
533 {
534 int pass;
535 double xx, xdelta, yy_sin, yy_cos, yy_mapm_lib, ydiff, tolerance;
536 char buffer[32];
537 M_APM mapm1, sinval, cosval;
538
539 fprintf(stdout, "Validating the SIN_COS function ... \n");
540
541 mapm1 = m_apm_init();
542 sinval = m_apm_init();
543 cosval = m_apm_init();
544
545 xx = -7.715;
546 xdelta = 0.737;
547 tolerance = 1.0E-13;
548 pass = 0;
549
550 while (1)
551 {
552 yy_sin = sin(xx);
553 yy_cos = cos(xx);
554
555 m_apm_set_double(mapm1, xx);
556 m_apm_sin_cos(sinval, cosval, DECIMAL_PLACES, mapm1);
557
558 m_apm_to_string(buffer, DECIMAL_PLACES, sinval);
559 yy_mapm_lib = atof(buffer);
560
561 ydiff = yy_mapm_lib - yy_sin;
562
563 if (fabs(ydiff) > tolerance)
564 {
565 pass = 1;
566 fprintf(stdout,"x = %17.10E C-math-lib = %.15E \n",xx,yy_sin);
567 fprintf(stdout,"DIFF = %17.10E M_APM_LIB = %s \n",ydiff,buffer);
568 }
569
570 m_apm_to_string(buffer, DECIMAL_PLACES, cosval);
571 yy_mapm_lib = atof(buffer);
572
573 ydiff = yy_mapm_lib - yy_cos;
574
575 if (fabs(ydiff) > tolerance)
576 {
577 pass = 1;
578 fprintf(stdout,"x = %17.10E C-math-lib = %.15E \n",xx,yy_cos);
579 fprintf(stdout,"DIFF = %17.10E M_APM_LIB = %s \n",ydiff,buffer);
580 }
581
582 xx += xdelta;
583 if (xx > 7.61)
584 break;
585 }
586
587 if (pass == 0)
588 fprintf(stdout,"... SIN_COS function passes\n");
589
590 m_apm_free(mapm1);
591 m_apm_free(sinval);
592 m_apm_free(cosval);
593
594 return(pass);
595 }
596
597 /****************************************************************************/
598
599 int mapm_test_tan()
600 {
601 int pass;
602 double xx, xdelta, yy, yy_mapm_lib, ydiff, tolerance;
603 char buffer[32];
604 M_APM mapm1, mapm2;
605
606 fprintf(stdout, "Validating the TAN function ... \n");
607
608 mapm1 = m_apm_init();
609 mapm2 = m_apm_init();
610
611 xx = -7.315;
612 xdelta = 0.673;
613 tolerance = 1.0E-13;
614 pass = 0;
615
616 while (1)
617 {
618 yy = tan(xx);
619
620 m_apm_set_double(mapm1, xx);
621 m_apm_tan(mapm2, DECIMAL_PLACES, mapm1);
622 m_apm_to_string(buffer, DECIMAL_PLACES, mapm2);
623 yy_mapm_lib = atof(buffer);
624
625 if (fabs(yy) > 1.0)
626 ydiff = (yy_mapm_lib - yy) / yy;
627 else
628 ydiff = yy_mapm_lib - yy;
629
630 if (fabs(ydiff) > tolerance)
631 {
632 pass = 1;
633 fprintf(stdout,"x = %17.10E C-math-lib = %.15E \n", xx, yy);
634 fprintf(stdout,"DIFF = %17.10E M_APM_LIB = %s \n", ydiff, buffer);
635 }
636
637 xx += xdelta;
638 if (xx > 7.21)
639 break;
640 }
641
642 if (pass == 0)
643 fprintf(stdout,"... TAN function passes\n");
644
645 m_apm_free(mapm1);
646 m_apm_free(mapm2);
647
648 return(pass);
649 }
650
651 /****************************************************************************/
652
653 int mapm_test_asin()
654 {
655 int pass;
656 double xx, xdelta, yy, yy_mapm_lib, ydiff, tolerance;
657 char buffer[32];
658 M_APM mapm1, mapm2;
659
660 fprintf(stdout, "Validating the ASIN function ... \n");
661
662 mapm1 = m_apm_init();
663 mapm2 = m_apm_init();
664
665 xx = -0.997;
666 xdelta = 0.1073;
667 tolerance = 1.0E-13;
668 pass = 0;
669
670 while (1)
671 {
672 yy = asin(xx);
673
674 m_apm_set_double(mapm1, xx);
675 m_apm_asin(mapm2, DECIMAL_PLACES, mapm1);
676 m_apm_to_string(buffer, DECIMAL_PLACES, mapm2);
677 yy_mapm_lib = atof(buffer);
678
679 ydiff = yy_mapm_lib - yy;
680
681 if (fabs(ydiff) > tolerance)
682 {
683 pass = 1;
684 fprintf(stdout,"x = %17.10E C-math-lib = %.15E \n", xx, yy);
685 fprintf(stdout,"DIFF = %17.10E M_APM_LIB = %s \n", ydiff, buffer);
686 }
687
688 xx += xdelta;
689 if (xx > 0.99)
690 break;
691 }
692
693 if (pass == 0)
694 fprintf(stdout,"... ASIN function passes\n");
695
696 m_apm_free(mapm1);
697 m_apm_free(mapm2);
698
699 return(pass);
700 }
701
702 /****************************************************************************/
703
704 int mapm_test_acos()
705 {
706 int pass;
707 double xx, xdelta, yy, yy_mapm_lib, ydiff, tolerance;
708 char buffer[32];
709 M_APM mapm1, mapm2;
710
711 fprintf(stdout, "Validating the ACOS function ... \n");
712
713 mapm1 = m_apm_init();
714 mapm2 = m_apm_init();
715
716 xx = -0.998;
717 xdelta = 0.09373;
718 tolerance = 1.0E-13;
719 pass = 0;
720
721 while (1)
722 {
723 yy = acos(xx);
724
725 m_apm_set_double(mapm1, xx);
726 m_apm_acos(mapm2, DECIMAL_PLACES, mapm1);
727 m_apm_to_string(buffer, DECIMAL_PLACES, mapm2);
728 yy_mapm_lib = atof(buffer);
729
730 ydiff = yy_mapm_lib - yy;
731
732 if (fabs(ydiff) > tolerance)
733 {
734 pass = 1;
735 fprintf(stdout,"x = %17.10E C-math-lib = %.15E \n", xx, yy);
736 fprintf(stdout,"DIFF = %17.10E M_APM_LIB = %s \n", ydiff, buffer);
737 }
738
739 xx += xdelta;
740 if (xx > 0.99)
741 break;
742 }
743
744 if (pass == 0)
745 fprintf(stdout,"... ACOS function passes\n");
746
747 m_apm_free(mapm1);
748 m_apm_free(mapm2);
749
750 return(pass);
751 }
752
753 /****************************************************************************/
754
755 int mapm_test_atan()
756 {
757 int pass;
758 double xx, xdelta, yy, yy_mapm_lib, ydiff, tolerance;
759 char buffer[32];
760 M_APM mapm1, mapm2;
761
762 fprintf(stdout, "Validating the ATAN function ... \n");
763
764 mapm1 = m_apm_init();
765 mapm2 = m_apm_init();
766
767 xx = -6.8215;
768 xdelta = 0.3967;
769 tolerance = 1.0E-13;
770 pass = 0;
771
772 while (1)
773 {
774 yy = atan(xx);
775
776 m_apm_set_double(mapm1, xx);
777 m_apm_atan(mapm2, DECIMAL_PLACES, mapm1);
778 m_apm_to_string(buffer, DECIMAL_PLACES, mapm2);
779 yy_mapm_lib = atof(buffer);
780
781 ydiff = yy_mapm_lib - yy;
782
783 if (fabs(ydiff) > tolerance)
784 {
785 pass = 1;
786 fprintf(stdout,"x = %17.10E C-math-lib = %.15E \n", xx, yy);
787 fprintf(stdout,"DIFF = %17.10E M_APM_LIB = %s \n", ydiff, buffer);
788 }
789
790 if (xx > 6.5)
791 break;
792
793 xx += xdelta;
794 }
795
796 if (pass == 0)
797 fprintf(stdout,"... ATAN function passes\n");
798
799 m_apm_free(mapm1);
800 m_apm_free(mapm2);
801
802 return(pass);
803 }
804
805 /****************************************************************************/
806
807 int mapm_test_atan2()
808 {
809 int i, j, pass;
810 double xx1, xx2, yy, yy_mapm_lib, ydiff, tolerance;
811 double xarray[10], yarray[10];
812 char buffer[32];
813 M_APM mapmX, mapmY, mapm0;
814
815 fprintf(stdout, "Validating the ATAN2 function ... \n");
816
817 mapm0 = m_apm_init();
818 mapmX = m_apm_init();
819 mapmY = m_apm_init();
820
821 pass = 0;
822 tolerance = 1.0E-13;
823
824 xarray[1] = -1.638;
825 xarray[2] = -2.179;
826 xarray[3] = 0.0;
827 xarray[4] = 0.842;
828 xarray[5] = 1.976;
829 xarray[6] = 9.237;
830 xarray[7] = -8.324;
831 xarray[8] = -1.287E-9;
832 xarray[9] = 74.861;
833
834 yarray[1] = -0.374;
835 yarray[2] = -4.127;
836 yarray[3] = 0.0;
837 yarray[4] = 0.732;
838 yarray[5] = 3.261;
839 yarray[6] = -6.305;
840 yarray[7] = 11.427;
841 yarray[8] = -0.8827;
842 yarray[9] = 7.385E10;
843
844 for (i=1; i <= 9; i++)
845 {
846 for (j=1; j <= 9; j++)
847 {
848 if (i != 3 || j != 3) /* domain error if both == 0.0 */
849 {
850 xx1 = xarray[i];
851 xx2 = yarray[j];
852
853 yy = atan2(xx2, xx1);
854
855 m_apm_set_double(mapmX, xx1);
856 m_apm_set_double(mapmY, xx2);
857 m_apm_atan2(mapm0, DECIMAL_PLACES, mapmY, mapmX);
858 m_apm_to_string(buffer, DECIMAL_PLACES, mapm0);
859 yy_mapm_lib = atof(buffer);
860
861 ydiff = yy_mapm_lib - yy;
862
863 if (fabs(ydiff) > tolerance)
864 {
865 pass = 1;
866
867 fprintf(stdout,
868 "x2 = %10.3E x1 = %10.3E C-math-lib = %.15E \n",xx2,xx1,yy);
869 fprintf(stdout,
870 "DIFF = %17.10E M_APM_LIB = %s \n",ydiff,buffer);
871 }
872 }
873 }
874 }
875
876 if (pass == 0)
877 fprintf(stdout,"... ATAN2 function passes\n");
878
879 m_apm_free(mapm0);
880 m_apm_free(mapmX);
881 m_apm_free(mapmY);
882
883 return(pass);
884 }
885
886 /****************************************************************************/
887
888 int mapm_test_sinh()
889 {
890 int pass;
891 double xx, xdelta, yy, yy_mapm_lib, ydiff, tolerance;
892 char buffer[32];
893 M_APM mapm1, mapm2;
894
895 fprintf(stdout, "Validating the SINH function ... \n");
896
897 mapm1 = m_apm_init();
898 mapm2 = m_apm_init();
899
900 xx = -4.2215;
901 xdelta = 0.3967;
902 tolerance = 1.0E-13;
903 pass = 0;
904
905 while (1)
906 {
907 yy = sinh(xx);
908
909 m_apm_set_double(mapm1, xx);
910 m_apm_sinh(mapm2, DECIMAL_PLACES, mapm1);
911 m_apm_to_string(buffer, DECIMAL_PLACES, mapm2);
912 yy_mapm_lib = atof(buffer);
913
914 if (fabs(yy) > 1.0)
915 ydiff = (yy_mapm_lib - yy) / yy;
916 else
917 ydiff = yy_mapm_lib - yy;
918
919 if (fabs(ydiff) > tolerance)
920 {
921 pass = 1;
922 fprintf(stdout,"x = %17.10E C-math-lib = %.15E \n", xx, yy);
923 fprintf(stdout,"DIFF = %17.10E M_APM_LIB = %s \n", ydiff, buffer);
924 }
925
926 if (xx > 4.2)
927 break;
928
929 xx += xdelta;
930 }
931
932 if (pass == 0)
933 fprintf(stdout,"... SINH function passes\n");
934
935 m_apm_free(mapm1);
936 m_apm_free(mapm2);
937
938 return(pass);
939 }
940
941 /****************************************************************************/
942
943 int mapm_test_cosh()
944 {
945 int pass;
946 double xx, xdelta, yy, yy_mapm_lib, ydiff, tolerance;
947 char buffer[32];
948 M_APM mapm1, mapm2;
949
950 fprintf(stdout, "Validating the COSH function ... \n");
951
952 mapm1 = m_apm_init();
953 mapm2 = m_apm_init();
954
955 xx = -4.2215;
956 xdelta = 0.3816;
957 tolerance = 1.0E-13;
958 pass = 0;
959
960 while (1)
961 {
962 yy = cosh(xx);
963
964 m_apm_set_double(mapm1, xx);
965 m_apm_cosh(mapm2, DECIMAL_PLACES, mapm1);
966 m_apm_to_string(buffer, DECIMAL_PLACES, mapm2);
967 yy_mapm_lib = atof(buffer);
968
969 if (fabs(yy) > 1.0)
970 ydiff = (yy_mapm_lib - yy) / yy;
971 else
972 ydiff = yy_mapm_lib - yy;
973
974 if (fabs(ydiff) > tolerance)
975 {
976 pass = 1;
977 fprintf(stdout,"x = %17.10E C-math-lib = %.15E \n", xx, yy);
978 fprintf(stdout,"DIFF = %17.10E M_APM_LIB = %s \n", ydiff, buffer);
979 }
980
981 if (xx > 4.2)
982 break;
983
984 xx += xdelta;
985 }
986
987 if (pass == 0)
988 fprintf(stdout,"... COSH function passes\n");
989
990 m_apm_free(mapm1);
991 m_apm_free(mapm2);
992
993 return(pass);
994 }
995
996 /****************************************************************************/
997
998 int mapm_test_tanh()
999 {
1000 int pass;
1001 double xx, xdelta, yy, yy_mapm_lib, ydiff, tolerance;
1002 char buffer[32];
1003 M_APM mapm1, mapm2;
1004
1005 fprintf(stdout, "Validating the TANH function ... \n");
1006
1007 mapm1 = m_apm_init();
1008 mapm2 = m_apm_init();
1009
1010 xx = -2.1083;
1011 xdelta = 0.10769;
1012 tolerance = 1.0E-13;
1013 pass = 0;
1014
1015 while (1)
1016 {
1017 yy = tanh(xx);
1018
1019 m_apm_set_double(mapm1, xx);
1020 m_apm_tanh(mapm2, DECIMAL_PLACES, mapm1);
1021 m_apm_to_string(buffer, DECIMAL_PLACES, mapm2);
1022 yy_mapm_lib = atof(buffer);
1023
1024 ydiff = yy_mapm_lib - yy;
1025
1026 if (fabs(ydiff) > tolerance)
1027 {
1028 pass = 1;
1029 fprintf(stdout,"x = %17.10E C-math-lib = %.15E \n", xx, yy);
1030 fprintf(stdout,"DIFF = %17.10E M_APM_LIB = %s \n", ydiff, buffer);
1031 }
1032
1033 if (xx > 2.1)
1034 break;
1035
1036 xx += xdelta;
1037 }
1038
1039 if (pass == 0)
1040 fprintf(stdout,"... TANH function passes\n");
1041
1042 m_apm_free(mapm1);
1043 m_apm_free(mapm2);
1044
1045 return(pass);
1046 }
1047
1048 /****************************************************************************/
1049
1050 int mapm_test_asinh()
1051 {
1052 int pass;
1053 double xx, xdelta, yy, yy_mapm_lib, ydiff, tolerance;
1054 char buffer[32];
1055 M_APM mapm1, mapm2;
1056
1057 fprintf(stdout, "Validating the ASINH function ... \n");
1058
1059 mapm1 = m_apm_init();
1060 mapm2 = m_apm_init();
1061
1062 xx = -27.2215;
1063 xdelta = 3.1817;
1064 tolerance = 1.0E-13;
1065 pass = 0;
1066
1067 while (1)
1068 {
1069 yy = asinh_local(xx);
1070
1071 m_apm_set_double(mapm1, xx);
1072 m_apm_asinh(mapm2, DECIMAL_PLACES, mapm1);
1073 m_apm_to_string(buffer, DECIMAL_PLACES, mapm2);
1074 yy_mapm_lib = atof(buffer);
1075
1076 if (fabs(yy) > 1.0)
1077 ydiff = (yy_mapm_lib - yy) / yy;
1078 else
1079 ydiff = yy_mapm_lib - yy;
1080
1081 if (fabs(ydiff) > tolerance)
1082 {
1083 pass = 1;
1084 fprintf(stdout,"x = %17.10E C-math-lib = %.15E \n", xx, yy);
1085 fprintf(stdout,"DIFF = %17.10E M_APM_LIB = %s \n", ydiff, buffer);
1086 }
1087
1088 if (xx > 27.2)
1089 break;
1090
1091 xx += xdelta;
1092 }
1093
1094 if (pass == 0)
1095 fprintf(stdout,"... ASINH function passes\n");
1096
1097 m_apm_free(mapm1);
1098 m_apm_free(mapm2);
1099
1100 return(pass);
1101 }
1102
1103 /****************************************************************************/
1104
1105 int mapm_test_acosh()
1106 {
1107 int pass;
1108 double xx, xdelta, yy, yy_mapm_lib, ydiff, tolerance;
1109 char buffer[32];
1110 M_APM mapm1, mapm2;
1111
1112 fprintf(stdout, "Validating the ACOSH function ... \n");
1113
1114 mapm1 = m_apm_init();
1115 mapm2 = m_apm_init();
1116
1117 xx = 1.02;
1118 xdelta = 0.7816;
1119 tolerance = 1.0E-13;
1120 pass = 0;
1121
1122 while (1)
1123 {
1124 yy = acosh_local(xx);
1125
1126 m_apm_set_double(mapm1, xx);
1127 m_apm_acosh(mapm2, DECIMAL_PLACES, mapm1);
1128 m_apm_to_string(buffer, DECIMAL_PLACES, mapm2);
1129 yy_mapm_lib = atof(buffer);
1130
1131 if (fabs(yy) > 1.0)
1132 ydiff = (yy_mapm_lib - yy) / yy;
1133 else
1134 ydiff = yy_mapm_lib - yy;
1135
1136 if (fabs(ydiff) > tolerance)
1137 {
1138 pass = 1;
1139 fprintf(stdout,"x = %17.10E C-math-lib = %.15E \n", xx, yy);
1140 fprintf(stdout,"DIFF = %17.10E M_APM_LIB = %s \n", ydiff, buffer);
1141 }
1142
1143 if (xx > 30.0)
1144 break;
1145
1146 xx += xdelta;
1147 }
1148
1149 if (pass == 0)
1150 fprintf(stdout,"... acosh function passes\n");
1151
1152 m_apm_free(mapm1);
1153 m_apm_free(mapm2);
1154
1155 return(pass);
1156 }
1157
1158 /****************************************************************************/
1159
1160 int mapm_test_atanh()
1161 {
1162 int pass;
1163 double xx, xdelta, yy, yy_mapm_lib, ydiff, tolerance;
1164 char buffer[32];
1165 M_APM mapm1, mapm2;
1166
1167 fprintf(stdout, "Validating the ATANH function ... \n");
1168
1169 mapm1 = m_apm_init();
1170 mapm2 = m_apm_init();
1171
1172 xx = -0.998;
1173 xdelta = 0.1167;
1174 tolerance = 1.0E-13;
1175 pass = 0;
1176
1177 while (1)
1178 {
1179 yy = atanh_local(xx);
1180
1181 m_apm_set_double(mapm1, xx);
1182 m_apm_atanh(mapm2, DECIMAL_PLACES, mapm1);
1183 m_apm_to_string(buffer, DECIMAL_PLACES, mapm2);
1184 yy_mapm_lib = atof(buffer);
1185
1186 if (fabs(yy) > 1.0)
1187 ydiff = (yy_mapm_lib - yy) / yy;
1188 else
1189 ydiff = yy_mapm_lib - yy;
1190
1191 if (fabs(ydiff) > tolerance)
1192 {
1193 pass = 1;
1194 fprintf(stdout,"x = %17.10E C-math-lib = %.15E \n", xx, yy);
1195 fprintf(stdout,"DIFF = %17.10E M_APM_LIB = %s \n", ydiff, buffer);
1196 }
1197
1198 xx += xdelta;
1199
1200 if (xx > 0.99)
1201 break;
1202 }
1203
1204 if (pass == 0)
1205 fprintf(stdout,"... ATANH function passes\n");
1206
1207 m_apm_free(mapm1);
1208 m_apm_free(mapm2);
1209
1210 return(pass);
1211 }
1212
1213 /****************************************************************************/
1214
1215 int mapm_test_exp()
1216 {
1217 int pass;
1218 double xx, xdelta, yy, yy_mapm_lib, ydiff, tolerance;
1219 char buffer[32];
1220 M_APM mapm1, mapm2;
1221
1222 fprintf(stdout, "Validating the EXP function ... \n");
1223
1224 mapm1 = m_apm_init();
1225 mapm2 = m_apm_init();
1226
1227 xx = -42.215;
1228 xdelta = 3.387;
1229 tolerance = 1.0E-13;
1230 pass = 0;
1231
1232 while (1)
1233 {
1234 yy = exp(xx);
1235
1236 m_apm_set_double(mapm1, xx);
1237 m_apm_exp(mapm2, DECIMAL_PLACES, mapm1);
1238 m_apm_to_string(buffer, DECIMAL_PLACES, mapm2);
1239 yy_mapm_lib = atof(buffer);
1240
1241 ydiff = (yy_mapm_lib - yy) / yy;
1242
1243 if (fabs(ydiff) > tolerance)
1244 {
1245 pass = 1;
1246 fprintf(stdout,"x = %17.10E C-math-lib = %.15E \n", xx, yy);
1247 fprintf(stdout,"DIFF = %17.10E M_APM_LIB = %s \n", ydiff, buffer);
1248 }
1249
1250 if (xx > 44.2)
1251 break;
1252
1253 xx += xdelta;
1254 }
1255
1256 if (pass == 0)
1257 fprintf(stdout,"... EXP function passes\n");
1258
1259 m_apm_free(mapm1);
1260 m_apm_free(mapm2);
1261
1262 return(pass);
1263 }
1264
1265 /****************************************************************************/
1266
1267 int mapm_test_exp_2()
1268 {
1269 int kk, dplaces, pass;
1270 M_APM mapm1, mapm2, mapm3, mapm4;
1271 char *cp;
1272 static char *str1[4] = { "5.2107E-4", "-7.138927E-5", "8.21735E-6", "0.0" };
1273 static char *str2[4] = { "128.00047", "511.99874", ".03125012", ".01562483" };
1274
1275 fprintf(stdout,
1276 "Validating calculations involving EXP & LOG ... \n");
1277
1278 mapm1 = m_apm_init();
1279 mapm2 = m_apm_init();
1280 mapm3 = m_apm_init();
1281 mapm4 = m_apm_init();
1282 pass = 0;
1283
1284 dplaces = 250;
1285
1286 for (kk=0; kk < 4; kk++)
1287 {
1288 cp = str1[kk];
1289
1290 m_apm_set_string(mapm4, cp);
1291
1292 fprintf(stdout,
1293 "Calculate log(exp(%s)) to %d decimal places ... \n", cp, dplaces);
1294
1295 m_apm_exp(mapm3, (dplaces + 8), mapm4);
1296 m_apm_add(mapm1, mapm3, MM_Zero);
1297 m_apm_log(mapm2, dplaces, mapm1);
1298 m_apm_add(mapm1, MM_Zero, mapm2);
1299
1300 if (m_apm_compare(mapm1, mapm4) == 0)
1301 {
1302 fprintf(stdout, "Verified %s == log(exp(%s)) \n", cp, cp);
1303 }
1304 else
1305 {
1306 pass = 1;
1307 fprintf(stdout,
1308 "***** FAILED to verify %s == log(exp(%s)) \n", cp, cp);
1309 }
1310 }
1311
1312
1313 for (kk=0; kk < 4; kk++)
1314 {
1315 cp = str2[kk];
1316
1317 m_apm_set_string(mapm1, cp); /* close to exact power of 2 */
1318
1319 fprintf(stdout,
1320 "Calculate exp(log(%s)) to %d decimal places ... \n", cp, dplaces);
1321
1322 m_apm_log(mapm2, (dplaces + 8), mapm1);
1323 m_apm_exp(mapm3, dplaces, mapm2);
1324
1325 if (m_apm_compare(mapm1, mapm3) == 0)
1326 {
1327 fprintf(stdout, "Verified %s == exp(log(%s)) \n", cp, cp);
1328 }
1329 else
1330 {
1331 pass = 1;
1332 fprintf(stdout,
1333 "***** FAILED to verify %s == exp(log(%s)) \n", cp, cp);
1334 }
1335 }
1336
1337
1338 if (pass == 0)
1339 fprintf(stdout,"... EXP & LOG Calculations pass\n");
1340
1341 m_apm_free(mapm1);
1342 m_apm_free(mapm2);
1343 m_apm_free(mapm3);
1344 m_apm_free(mapm4);
1345
1346 return(pass);
1347 }
1348
1349 /****************************************************************************/
1350
1351 int mapm_test_log()
1352 {
1353 int pass;
1354 double xx, xdelta, yy, yy_mapm_lib, ydiff, tolerance;
1355 char buffer[32];
1356 M_APM mapm1, mapm2;
1357
1358 fprintf(stdout, "Validating the LOG function ... \n");
1359
1360 mapm1 = m_apm_init();
1361 mapm2 = m_apm_init();
1362
1363 xx = 1.7235E-20;
1364 xdelta = 8.3816;
1365 tolerance = 1.0E-13;
1366 pass = 0;
1367
1368 while (1)
1369 {
1370 yy = log(xx);
1371
1372 m_apm_set_double(mapm1, xx);
1373 m_apm_log(mapm2, DECIMAL_PLACES, mapm1);
1374 m_apm_to_string(buffer, DECIMAL_PLACES, mapm2);
1375 yy_mapm_lib = atof(buffer);
1376
1377 if (fabs(yy) > 1.0)
1378 ydiff = (yy_mapm_lib - yy) / yy;
1379 else
1380 ydiff = yy_mapm_lib - yy;
1381
1382 if (fabs(ydiff) > tolerance)
1383 {
1384 pass = 1;
1385 fprintf(stdout,"x = %17.10E C-math-lib = %.15E \n", xx, yy);
1386 fprintf(stdout,"DIFF = %17.10E M_APM_LIB = %s \n", ydiff, buffer);
1387 }
1388
1389 if (xx > 2.5E20)
1390 break;
1391
1392 xx *= xdelta;
1393 }
1394
1395 if (pass == 0)
1396 fprintf(stdout,"... LOG function passes\n");
1397
1398 m_apm_free(mapm1);
1399 m_apm_free(mapm2);
1400
1401 return(pass);
1402 }
1403
1404 /****************************************************************************/
1405
1406 int mapm_test_log10()
1407 {
1408 int pass;
1409 double xx, xdelta, yy, yy_mapm_lib, ydiff, tolerance;
1410 char buffer[32];
1411 M_APM mapm1, mapm2;
1412
1413 fprintf(stdout, "Validating the LOG10 function ... \n");
1414
1415 mapm1 = m_apm_init();
1416 mapm2 = m_apm_init();
1417
1418 xx = 3.4235E-22;
1419 xdelta = 12.7816;
1420 tolerance = 1.0E-13;
1421 pass = 0;
1422
1423 while (1)
1424 {
1425 yy = log10(xx);
1426
1427 m_apm_set_double(mapm1, xx);
1428 m_apm_log10(mapm2, DECIMAL_PLACES, mapm1);
1429 m_apm_to_string(buffer, DECIMAL_PLACES, mapm2);
1430 yy_mapm_lib = atof(buffer);
1431
1432 if (fabs(yy) > 1.0)
1433 ydiff = (yy_mapm_lib - yy) / yy;
1434 else
1435 ydiff = yy_mapm_lib - yy;
1436
1437 if (fabs(ydiff) > tolerance)
1438 {
1439 pass = 1;
1440 fprintf(stdout,"x = %17.10E C-math-lib = %.15E \n", xx, yy);
1441 fprintf(stdout,"DIFF = %17.10E M_APM_LIB = %s \n", ydiff, buffer);
1442 }
1443
1444 if (xx > 2.5E22)
1445 break;
1446
1447 xx *= xdelta;
1448 }
1449
1450 if (pass == 0)
1451 fprintf(stdout,"... LOG10 function passes\n");
1452
1453 m_apm_free(mapm1);
1454 m_apm_free(mapm2);
1455
1456 return(pass);
1457 }
1458
1459 /****************************************************************************/
1460
1461 int mapm_test_pow()
1462 {
1463 int i, pass;
1464 double xx1, xx2, yy, yy_mapm_lib, ydiff, tolerance;
1465 double xarray[22], yarray[22];
1466 char buffer[32];
1467 M_APM mapm0, mapmX, mapmY;
1468
1469 fprintf(stdout, "Validating the POW function ... \n");
1470
1471 mapm0 = m_apm_init();
1472 mapmX = m_apm_init();
1473 mapmY = m_apm_init();
1474
1475 pass = 0;
1476 tolerance = 1.0E-13;
1477
1478 xarray[1] = 1.638; yarray[1] = 7.14;
1479 xarray[2] = 8.321e-12; yarray[2] = -5.2;
1480 xarray[3] = 1.0; yarray[3] = 4.19;
1481 xarray[4] = 7.3231e17; yarray[4] = 2.11;
1482 xarray[5] = 2.0; yarray[5] = 32.0;
1483 xarray[6] = 9.997e4; yarray[6] = -1.0;
1484 xarray[7] = 1.0528e3; yarray[7] = -3.83e-2;
1485 xarray[8] = 6.319e-6; yarray[8] = 3.17;
1486 xarray[9] = 2.108e9; yarray[9] = -2.0;
1487 xarray[10] = 8.417e-3; yarray[10] = 0.5;
1488 xarray[11] = 8.417e-3; yarray[11] = -6.2109;
1489 xarray[12] = 5.347e7; yarray[12] = 1.0;
1490 xarray[13] = 10.0; yarray[13] = -17.0;
1491 xarray[14] = 10.0; yarray[14] = 23.37;
1492 xarray[15] = 9.0e2; yarray[15] = 0.0;
1493 xarray[16] = 2.0; yarray[16] = -3.62;
1494 xarray[17] = 2.0; yarray[17] = 23.017;
1495 xarray[18] = 8.27693; yarray[18] = 5.19832;
1496 xarray[19] = 0.0; yarray[19] = 32.02384;
1497 xarray[20] = 78.23; yarray[20] = 9.2104;
1498
1499 for (i=1; i <= 20; i++)
1500 {
1501 xx1 = xarray[i];
1502 xx2 = yarray[i];
1503 yy = pow(xx1, xx2);
1504
1505 m_apm_set_double(mapmX, xx1);
1506 m_apm_set_double(mapmY, xx2);
1507 m_apm_pow(mapm0, DECIMAL_PLACES, mapmX, mapmY);
1508 m_apm_to_string(buffer, DECIMAL_PLACES, mapm0);
1509 yy_mapm_lib = atof(buffer);
1510
1511 if (fabs(yy) > 1.0)
1512 ydiff = (yy_mapm_lib - yy) / yy;
1513 else
1514 ydiff = yy_mapm_lib - yy;
1515
1516 if (fabs(ydiff) > tolerance)
1517 {
1518 pass = 1;
1519 fprintf(stdout,
1520 "x1 = %10.3E x2 = %10.3E C-math-lib = %.15E \n",xx1,xx2,yy);
1521 fprintf(stdout,
1522 "DIFF = %17.10E M_APM_LIB = %s \n",ydiff,buffer);
1523 }
1524 }
1525
1526 if (pass == 0)
1527 fprintf(stdout,"... POW function passes\n");
1528
1529 m_apm_free(mapm0);
1530 m_apm_free(mapmX);
1531 m_apm_free(mapmY);
1532
1533 return(pass);
1534 }
1535
1536 /****************************************************************************/
1537
1538 int mapm_test_factorial()
1539 {
1540 int i, pass;
1541 M_APM mapm1, mapm2, mapmi;
1542
1543 fprintf(stdout, "Validating the FACTORIAL function ... \n");
1544
1545 pass = 0;
1546
1547 mapmi = m_apm_init();
1548 mapm1 = m_apm_init();
1549 mapm2 = m_apm_init();
1550
1551 for (i=0; i <= 15; i++)
1552 {
1553 m_apm_set_long(mapmi, (long)i);
1554
1555 factorial_local(mapm2, mapmi);
1556 m_apm_factorial(mapm1, mapmi);
1557
1558 if (m_apm_compare(mapm1, mapm2) != 0)
1559 {
1560 pass = 1;
1561 fprintf(stdout, "FAIL: miscompare of factorial %d \n", i);
1562 }
1563 }
1564
1565
1566 m_apm_set_string(mapmi, "37");
1567
1568 factorial_local(mapm2, mapmi);
1569 m_apm_factorial(mapm1, mapmi);
1570
1571 if (m_apm_compare(mapm1, mapm2) != 0)
1572 {
1573 pass = 1;
1574 fprintf(stdout, "FAIL: miscompare of factorial 37 \n");
1575 }
1576
1577
1578 m_apm_set_string(mapmi, "738");
1579
1580 factorial_local(mapm2, mapmi);
1581 m_apm_factorial(mapm1, mapmi);
1582
1583 if (m_apm_compare(mapm1, mapm2) != 0)
1584 {
1585 pass = 1;
1586 fprintf(stdout, "FAIL: miscompare of factorial 738 \n");
1587 }
1588
1589
1590 m_apm_set_string(mapmi, "143");
1591
1592 factorial_local(mapm2, mapmi);
1593 m_apm_factorial(mapm1, mapmi);
1594
1595 if (m_apm_compare(mapm1, mapm2) != 0)
1596 {
1597 pass = 1;
1598 fprintf(stdout, "FAIL: miscompare of factorial 143 \n");
1599 }
1600
1601
1602 if (pass == 0)
1603 fprintf(stdout,"... FACTORIAL function passes\n");
1604
1605 m_apm_free(mapmi);
1606 m_apm_free(mapm1);
1607 m_apm_free(mapm2);
1608
1609 return(pass);
1610 }
1611
1612 /****************************************************************************/
1613
1614 int mapm_test_gcd()
1615 {
1616 int i, pass;
1617 M_APM mapmt, mapm1, mapm2, mapmu, mapmv;
1618
1619 static char *str_u[15] = { "0", "108", "821632",
1620 "15688", "2149236", "67368481619705856",
1621 "66424636", "6058695",
1622 "24791601236051755008",
1623 "436138461283596140227642238435328",
1624 "982147792896",
1625 "33184916509275",
1626 "8947823498749749279247894",
1627 "5823614", "912759" };
1628
1629 static char *str_v[15] = { "151222", "0", "912374",
1630 "3388608", "93863556", "3595508517961728",
1631 "9996388", "17793890499",
1632 "7938882807659495424",
1633 "8728893958572592934385063297024",
1634 "3004167094272",
1635 "3497310791015625",
1636 "783487423789",
1637 "1135823614", "2759" };
1638
1639 fprintf(stdout, "Validating the GCD/LCM functions ... \n");
1640
1641 mapm1 = m_apm_init();
1642 mapm2 = m_apm_init();
1643 mapmt = m_apm_init();
1644 mapmu = m_apm_init();
1645 mapmv = m_apm_init();
1646
1647 pass = 0;
1648
1649 for (i=0; i < 15; i++)
1650 {
1651 m_apm_set_string(mapmu, str_u[i]);
1652 m_apm_set_string(mapmv, str_v[i]);
1653
1654 gcd_local(mapm1, mapmu, mapmv);
1655 m_apm_gcd(mapm2, mapmu, mapmv);
1656
1657 if (m_apm_compare(mapm1, mapm2) != 0)
1658 {
1659 pass = 1;
1660 fprintf(stdout, "FAIL: miscompare of GCD %d \n", i);
1661 }
1662
1663 /*
1664 * LCM(u, v) = (u * v) / GCD(u, v)
1665 */
1666
1667 m_apm_multiply(mapmt, mapmu, mapmv);
1668 m_apm_divide(mapm2, 200, mapmt, mapm1);
1669
1670 m_apm_lcm(mapm1, mapmu, mapmv);
1671
1672 if (m_apm_compare(mapm1, mapm2) != 0)
1673 {
1674 pass = 1;
1675 fprintf(stdout, "FAIL: miscompare of LCM %d \n", i);
1676 }
1677 }
1678
1679 if (pass == 0)
1680 fprintf(stdout,"... GCD/LCM functions pass\n");
1681
1682 m_apm_free(mapm1);
1683 m_apm_free(mapm2);
1684 m_apm_free(mapmt);
1685 m_apm_free(mapmu);
1686 m_apm_free(mapmv);
1687
1688 return(pass);
1689 }
1690
1691 /****************************************************************************/
1692
1693 int mapm_test_PI()
1694 {
1695 int i, dplaces, pass;
1696 M_APM mapmt, mapm1, mapm2, mapm3, mapm_pi;
1697
1698 fprintf(stdout, "Validating calculations involving PI ... \n");
1699
1700 mapm1 = m_apm_init();
1701 mapm2 = m_apm_init();
1702 mapm3 = m_apm_init();
1703 mapmt = m_apm_init();
1704 mapm_pi = m_apm_init();
1705
1706 pass = 0;
1707
1708 dplaces = 300;
1709
1710 for (i=0; i < 2; i++)
1711 {
1712 fprintf(stdout, "Calculate PI to %d decimal places ... \n", dplaces);
1713
1714 /* PI = acos(-1) */
1715
1716 m_apm_negate(mapmt, MM_One);
1717 m_apm_arccos(mapm_pi, dplaces, mapmt);
1718
1719 /*
1720 *
1721 */
1722
1723 fprintf(stdout,
1724 "Calculate 3 * acos(0.5) to %d decimal places ... \n", dplaces);
1725
1726 /* PI = 3 * acos(0.5) */
1727
1728 m_apm_set_string(mapmt, "0.5");
1729 m_apm_arccos(mapm1, (dplaces + 8), mapmt);
1730
1731 m_apm_multiply(mapm2, MM_Three, mapm1);
1732 m_apm_round(mapm1, dplaces, mapm2);
1733
1734 if (m_apm_compare(mapm1, mapm_pi) == 0)
1735 {
1736 fprintf(stdout, "Verified PI == 3 * acos(0.5) \n");
1737 }
1738 else
1739 {
1740 pass = 1;
1741 fprintf(stdout, "***** FAILED to verify PI == 3 * acos(0.5) \n");
1742 }
1743
1744 /*
1745 *
1746 */
1747
1748 fprintf(stdout,
1749 "Calculate 4 * atan(1) to %d decimal places ... \n", dplaces);
1750
1751 /* PI = 4 * atan(1) */
1752
1753 m_apm_arctan(mapm1, (dplaces + 8), MM_One);
1754
1755 m_apm_multiply(mapm2, mapm1, MM_Four);
1756 m_apm_round(mapm1, dplaces, mapm2);
1757
1758 if (m_apm_compare(mapm1, mapm_pi) == 0)
1759 {
1760 fprintf(stdout, "Verified PI == 4 * atan(1) \n");
1761 }
1762 else
1763 {
1764 pass = 1;
1765 fprintf(stdout, "***** FAILED to verify PI == 4 * atan(1) \n");
1766 }
1767
1768 /*
1769 *
1770 */
1771
1772 fprintf(stdout,
1773 "Calculate 1.2 * acos(-sqrt(3) / 2) to %d decimal places ... \n", dplaces);
1774
1775 /* PI = 1.2 * acos(-sqrt(3) / 2) */
1776
1777 m_apm_sqrt(mapm1, (dplaces + 8), MM_Three);
1778 m_apm_divide(mapm2, (dplaces + 8), mapm1, MM_Two);
1779 m_apm_negate(mapmt, mapm2);
1780 m_apm_arccos(mapm1, (dplaces + 8), mapmt);
1781
1782 m_apm_set_string(mapm2, "1.2");
1783 m_apm_multiply(mapmt, mapm1, mapm2);
1784 m_apm_round(mapm1, dplaces, mapmt);
1785
1786 if (m_apm_compare(mapm1, mapm_pi) == 0)
1787 {
1788 fprintf(stdout, "Verified PI == 1.2 * acos(-sqrt(3) / 2) \n");
1789 }
1790 else
1791 {
1792 pass = 1;
1793 fprintf(stdout,
1794 "***** FAILED to verify PI == 1.2 * acos(-sqrt(3) / 2) \n");
1795 }
1796
1797
1798 /*
1799 *
1800 */
1801
1802 fprintf(stdout,
1803 "Calculate 12 * atan(2 - sqrt(3)) to %d decimal places ... \n", dplaces);
1804
1805 /* PI = 12 * atan(2 - sqrt(3)) */
1806 /* sqrt(3) = 10 ^ (log10(3) / 2) */
1807
1808 m_apm_log10(mapm1, (dplaces + 8), MM_Three);
1809 m_apm_divide(mapm2, (dplaces + 8), mapm1, MM_Two);
1810 m_apm_pow(mapm1, (dplaces + 8), MM_Ten, mapm2);
1811
1812 m_apm_subtract(mapm2, MM_Two, mapm1);
1813
1814 m_apm_atan(mapm1, (dplaces + 8), mapm2);
1815
1816 m_apm_set_long(mapm2, 12L);
1817 m_apm_multiply(mapmt, mapm1, mapm2);
1818 m_apm_round(mapm1, dplaces, mapmt);
1819
1820 if (m_apm_compare(mapm1, mapm_pi) == 0)
1821 {
1822 fprintf(stdout, "Verified PI == 12 * atan(2 - sqrt(3)) \n");
1823 }
1824 else
1825 {
1826 pass = 1;
1827 fprintf(stdout,
1828 "***** FAILED to verify PI == 12 * atan(2 - sqrt(3)) \n");
1829 }
1830
1831 /*
1832 *
1833 */
1834
1835 fprintf(stdout,
1836 "Calculate cos(PI / 6) to %d decimal places ... \n", dplaces);
1837
1838 m_apm_set_string(mapm1, "6");
1839 m_apm_divide(mapmt, dplaces, mapm_pi, mapm1);
1840 m_apm_cos(mapm2, (dplaces - 8), mapmt);
1841
1842 fprintf(stdout,
1843 "Calculate sqrt(0.75) to %d decimal places ... \n", dplaces);
1844
1845 m_apm_set_string(mapmt, "0.75");
1846 m_apm_sqrt(mapm1, (dplaces - 8), mapmt);
1847
1848 if (m_apm_compare(mapm1, mapm2) == 0)
1849 {
1850 fprintf(stdout, "Verified cos(PI / 6) == sqrt(0.75) \n");
1851 }
1852 else
1853 {
1854 pass = 1;
1855 fprintf(stdout, "***** FAILED to verify cos(PI / 6) == sqrt(0.75) \n");
1856 }
1857
1858 /*
1859 *
1860 */
1861
1862 fprintf(stdout,
1863 "Calculate tan(5 * PI / 12) to %d decimal places ... \n", dplaces);
1864
1865 m_apm_multiply(mapm1, MM_Five, mapm_pi);
1866 m_apm_set_string(mapmt, "12");
1867 m_apm_divide(mapm2, dplaces, mapm1, mapmt);
1868 m_apm_tan(mapm1, (dplaces - 8), mapm2);
1869
1870 fprintf(stdout,
1871 "Calculate 2 + sqrt(3) to %d decimal places ... \n", dplaces);
1872
1873 m_apm_sqrt(mapmt, (dplaces - 8), MM_Three);
1874 m_apm_add(mapm2, MM_Two, mapmt);
1875
1876 if (m_apm_compare(mapm1, mapm2) == 0)
1877 {
1878 fprintf(stdout, "Verified tan(5 * PI / 12) == 2 + sqrt(3) \n");
1879 }
1880 else
1881 {
1882 pass = 1;
1883 fprintf(stdout,
1884 "***** FAILED to verify tan(5 * PI / 12) == 2 + sqrt(3) \n");
1885 }
1886
1887 /*
1888 *
1889 */
1890
1891 fprintf(stdout,
1892 "Calculate 2 * (asin(x) + acos(x)) to %d decimal places ... \n", dplaces);
1893
1894 m_apm_set_string(mapm2, "12.730348");
1895 m_apm_cbrt(mapm1, 96, mapm2);
1896 m_apm_divide(mapm3, 92, mapm1, MM_Four);
1897
1898 m_apm_asin(mapm1, (dplaces + 8), mapm3);
1899 m_apm_acos(mapm2, (dplaces + 8), mapm3);
1900
1901 m_apm_add(mapmt, mapm1, mapm2);
1902 m_apm_multiply(mapm1, MM_Two, mapmt);
1903 m_apm_round(mapm3, dplaces, mapm1);
1904
1905 if (m_apm_compare(mapm3, mapm_pi) == 0)
1906 {
1907 fprintf(stdout, "Verified PI = 2 * (asin(x) + acos(x)) \n");
1908 }
1909 else
1910 {
1911 pass = 1;
1912 fprintf(stdout,
1913 "***** FAILED to verify PI = 2 * (asin(x) + acos(x)) \n");
1914 }
1915
1916 dplaces += 500;
1917 }
1918
1919 if (pass == 0)
1920 fprintf(stdout,"... PI Calculations pass\n");
1921
1922 m_apm_free(mapm1);
1923 m_apm_free(mapm2);
1924 m_apm_free(mapm3);
1925 m_apm_free(mapmt);
1926 m_apm_free(mapm_pi);
1927
1928 return(pass);
1929 }
1930
1931 /****************************************************************************/
1932
1933 int mapm_test_exp_pwr2()
1934 {
1935 int i, dplaces, pass;
1936 M_APM mapm_E, mapm1, mapm2, mapm3, mapm4;
1937
1938 fprintf(stdout,
1939 "Validating calculations involving E & Integer POW_NR ... \n");
1940
1941 mapm_E = m_apm_init();
1942 mapm1 = m_apm_init();
1943 mapm2 = m_apm_init();
1944 mapm3 = m_apm_init();
1945 mapm4 = m_apm_init();
1946 pass = 0;
1947
1948 dplaces = 250;
1949
1950 for (i=0; i < 2; i++)
1951 {
1952 fprintf(stdout,
1953 "Calculate E (2.71828...) to %d decimal places ... \n", dplaces);
1954
1955 /* E = exp(1) */
1956
1957 m_apm_exp(mapm_E, (dplaces + 8), MM_One);
1958
1959 /*
1960 *
1961 */
1962
1963 fprintf(stdout, "Calculate exp(0.5) to %d decimal places ... \n", dplaces);
1964
1965 m_apm_set_string(mapm3, "0.5");
1966 m_apm_exp(mapm1, dplaces, mapm3);
1967
1968 fprintf(stdout, "Calculate sqrt(E) to %d decimal places ... \n", dplaces);
1969 m_apm_sqrt(mapm2, dplaces, mapm_E);
1970
1971 if (m_apm_compare(mapm1, mapm2) == 0)
1972 {
1973 fprintf(stdout, "Verified exp(0.5) == sqrt(E) \n");
1974 }
1975 else
1976 {
1977 pass = 1;
1978 fprintf(stdout, "***** FAILED to verify exp(0.5) == sqrt(E) \n");
1979 }
1980
1981 /*
1982 *
1983 */
1984
1985 fprintf(stdout, "Calculate exp(15) to %d decimal places ... \n", dplaces);
1986
1987 m_apm_set_string(mapm3, "15");
1988 m_apm_exp(mapm1, dplaces, mapm3);
1989
1990 m_apm_integer_pow_nr(mapm3, mapm_E, 15);
1991 m_apm_round(mapm2, dplaces, mapm3);
1992
1993 if (m_apm_compare(mapm1, mapm2) == 0)
1994 {
1995 fprintf(stdout, "Verified exp(15) == E ^ 15 \n");
1996 }
1997 else
1998 {
1999 pass = 1;
2000 fprintf(stdout, "***** FAILED to verify exp(15) == E ^ 15 \n");
2001 }
2002
2003 /*
2004 *
2005 */
2006
2007 fprintf(stdout, "Calculate exp(2) to %d decimal places ... \n", dplaces);
2008
2009 m_apm_exp(mapm1, dplaces, MM_Two);
2010
2011 m_apm_integer_pow_nr(mapm3, mapm_E, 2);
2012 m_apm_round(mapm2, dplaces, mapm3);
2013
2014 if (m_apm_compare(mapm1, mapm2) == 0)
2015 {
2016 fprintf(stdout, "Verified exp(2) == E ^ 2 \n");
2017 }
2018 else
2019 {
2020 pass = 1;
2021 fprintf(stdout, "***** FAILED to verify exp(2) == E ^ 2 \n");
2022 }
2023
2024 /*
2025 *
2026 */
2027
2028 fprintf(stdout, "Calculate exp(-6.5) to %d decimal places ... \n", dplaces);
2029
2030 m_apm_set_string(mapm3, "-6.5");
2031 m_apm_exp(mapm1, dplaces, mapm3);
2032
2033 m_apm_integer_pow_nr(mapm3, mapm_E, 13);
2034 m_apm_sqrt(mapm4, (dplaces + 8), mapm3);
2035 m_apm_reciprocal(mapm2, dplaces, mapm4);
2036
2037 if (m_apm_compare(mapm1, mapm2) == 0)
2038 {
2039 fprintf(stdout, "Verified exp(-6.5) == 1 / sqrt(E ^ 13) \n");
2040 }
2041 else
2042 {
2043 pass = 1;
2044 fprintf(stdout,
2045 "***** FAILED to verify exp(-6.5) == 1 / sqrt(E ^ 13) \n");
2046 }
2047
2048 /*
2049 *
2050 */
2051
2052 fprintf(stdout, "Calculate exp(8) to %d decimal places ... \n", dplaces);
2053
2054 m_apm_set_string(mapm3, "8");
2055 m_apm_exp(mapm1, dplaces, mapm3);
2056
2057 m_apm_integer_pow_nr(mapm3, mapm_E, 24);
2058 m_apm_cbrt(mapm2, dplaces, mapm3);
2059
2060 if (m_apm_compare(mapm1, mapm2) == 0)
2061 {
2062 fprintf(stdout, "Verified exp(8) == cbrt(E ^ 24) \n");
2063 }
2064 else
2065 {
2066 pass = 1;
2067 fprintf(stdout,
2068 "***** FAILED to verify exp(8) == cbrt(E ^ 24) \n");
2069 }
2070
2071 dplaces += 350;
2072
2073 if (M_get_sizeof_int() <= 2) /* end now for 16 bit ints */
2074 break;
2075 }
2076
2077 if (pass == 0)
2078 fprintf(stdout,"... E & Integer POW_NR Calculations pass\n");
2079
2080 m_apm_free(mapm1);
2081 m_apm_free(mapm2);
2082 m_apm_free(mapm3);
2083 m_apm_free(mapm4);
2084 m_apm_free(mapm_E);
2085
2086 return(pass);
2087 }
2088
2089 /****************************************************************************/
2090
2091 int mapm_test_log_2()
2092 {
2093 int i, dplaces, pass;
2094 M_APM mapm1, mapm2, mapm3, mapm4;
2095
2096 fprintf(stdout,
2097 "Validating log calculations with more digits ... \n");
2098
2099 mapm1 = m_apm_init();
2100 mapm2 = m_apm_init();
2101 mapm3 = m_apm_init();
2102 mapm4 = m_apm_init();
2103
2104 pass = 0;
2105
2106 dplaces = 220;
2107
2108 for (i=0; i < 2; i++)
2109 {
2110 fprintf(stdout,
2111 "Calculate log(3.6107) to %d decimal places ... \n", dplaces);
2112
2113 m_apm_set_string(mapm4, "3.6107");
2114 m_apm_log(mapm1, dplaces, mapm4);
2115
2116 fprintf(stdout,
2117 "Calculate log(3.6107) with an AGM algorithm to %d decimal places ... \n",
2118 dplaces);
2119
2120 log_local(mapm2, dplaces, mapm4);
2121
2122 if (m_apm_compare(mapm1, mapm2) == 0)
2123 {
2124 fprintf(stdout, "Verified log(3.6107) with an AGM algorithm \n");
2125 }
2126 else
2127 {
2128 pass = 1;
2129 fprintf(stdout,
2130 "***** FAILED to verify log(3.6107) with an AGM algorithm \n");
2131 }
2132
2133 /*
2134 *
2135 */
2136
2137 fprintf(stdout,
2138 "Calculate log(63.1874) to %d decimal places ... \n", dplaces);
2139
2140 m_apm_set_string(mapm4, "63.1874");
2141 m_apm_log(mapm1, dplaces, mapm4);
2142
2143 fprintf(stdout,
2144 "Calculate log(63.1874) with an AGM algorithm to %d decimal places ... \n",
2145 dplaces);
2146
2147 log_local(mapm2, dplaces, mapm4);
2148
2149 if (m_apm_compare(mapm1, mapm2) == 0)
2150 {
2151 fprintf(stdout, "Verified log(63.1874) with an AGM algorithm \n");
2152 }
2153 else
2154 {
2155 pass = 1;
2156 fprintf(stdout,
2157 "***** FAILED to verify log(63.1874) with an AGM algorithm \n");
2158 }
2159
2160 /*
2161 *
2162 */
2163
2164 fprintf(stdout,
2165 "Calculate log(184.9536) to %d decimal places ... \n", dplaces);
2166
2167 m_apm_set_string(mapm3, "184.9536");
2168 m_apm_log(mapm1, dplaces, mapm3);
2169
2170 fprintf(stdout,
2171 "Calculate log(184.9536) with an AGM algorithm to %d decimal places ... \n",
2172 dplaces);
2173
2174 log_local(mapm2, dplaces, mapm3);
2175
2176 if (m_apm_compare(mapm1, mapm2) == 0)
2177 {
2178 fprintf(stdout, "Verified log(184.9536) with an AGM algorithm \n");
2179 }
2180 else
2181 {
2182 pass = 1;
2183 fprintf(stdout,
2184 "***** FAILED to verify log(184.9536) with an AGM algorithm \n");
2185 }
2186
2187 dplaces += 480;
2188 }
2189
2190 if (pass == 0)
2191 fprintf(stdout,"... log calculations pass\n");
2192
2193 m_apm_free(mapm1);
2194 m_apm_free(mapm2);
2195 m_apm_free(mapm3);
2196 m_apm_free(mapm4);
2197
2198 return(pass);
2199 }
2200
2201 /****************************************************************************/
2202
2203 int mapm_test_floor_ceil()
2204 {
2205 int ict, pass;
2206 M_APM mapm0, mapm1, mapm2, mapmr;
2207
2208 fprintf(stdout, "Validating the FLOOR/CEIL functions ... \n");
2209
2210 mapm0 = m_apm_init();
2211 mapm1 = m_apm_init();
2212 mapm2 = m_apm_init();
2213 mapmr = m_apm_init();
2214 pass = 0;
2215
2216 m_apm_set_string(mapm0, "26.90253936868423");
2217 m_apm_set_long(mapmr, 26L);
2218
2219 ict = 50;
2220
2221 while (1)
2222 {
2223 m_apm_floor(mapm1, mapm0);
2224
2225 if (m_apm_compare(mapm1, mapmr) != 0)
2226 {
2227 pass = 1;
2228 fprintf(stdout,
2229 "***** FAILED \'floor\' iteration %d \n", ict);
2230 }
2231
2232 m_apm_add(mapm2, MM_One, mapmr);
2233
2234 m_apm_ceil(mapm1, mapm0);
2235
2236 if (m_apm_compare(mapm1, mapm2) != 0)
2237 {
2238 pass = 1;
2239 fprintf(stdout,
2240 "***** FAILED \'ceil\' iteration %d \n", ict);
2241 }
2242
2243 m_apm_subtract(mapm2, mapm0, MM_One);
2244 m_apm_copy(mapm0, mapm2);
2245
2246 m_apm_subtract(mapm2, mapmr, MM_One);
2247 m_apm_copy(mapmr, mapm2);
2248
2249 if (--ict == 0)
2250 break;
2251 }
2252
2253 /*
2254 * the following verifies that an exact int input
2255 * yields the same int as a result.
2256 */
2257
2258 m_apm_set_string(mapm0, "7892379239879273");
2259 m_apm_floor(mapm1, mapm0);
2260
2261 if (m_apm_compare(mapm1, mapm0) != 0)
2262 {
2263 pass = 1;
2264 fprintf(stdout,
2265 "***** FAILED \'floor\' int constant \n");
2266 }
2267
2268 m_apm_ceil(mapm1, mapm0);
2269
2270 if (m_apm_compare(mapm1, mapm0) != 0)
2271 {
2272 pass = 1;
2273 fprintf(stdout,
2274 "***** FAILED \'ceil\' int constant \n");
2275 }
2276
2277
2278 m_apm_set_string(mapm0, "-834734734781238480");
2279 m_apm_floor(mapm1, mapm0);
2280
2281 if (m_apm_compare(mapm1, mapm0) != 0)
2282 {
2283 pass = 1;
2284 fprintf(stdout,
2285 "***** FAILED \'floor\' int constant \n");
2286 }
2287
2288 m_apm_ceil(mapm1, mapm0);
2289
2290 if (m_apm_compare(mapm1, mapm0) != 0)
2291 {
2292 pass = 1;
2293 fprintf(stdout,
2294 "***** FAILED \'ceil\' int constant \n");
2295 }
2296
2297
2298 if (pass == 0)
2299 fprintf(stdout,"... FLOOR/CEIL functions pass\n");
2300
2301 m_apm_free(mapm0);
2302 m_apm_free(mapm1);
2303 m_apm_free(mapm2);
2304 m_apm_free(mapmr);
2305
2306 return(pass);
2307 }
2308
2309 /****************************************************************************/
2310
2311 int mapm_test_log_near_1()
2312 {
2313 int pass;
2314 char cbuf[32];
2315 M_APM mapm1, mapm2, mapmr;
2316
2317 fprintf(stdout, "Validating the LOG_NEAR_1 function ... \n");
2318
2319 mapm1 = m_apm_init();
2320 mapm2 = m_apm_init();
2321 mapmr = m_apm_init();
2322 pass = 0;
2323
2324 strcpy(cbuf, "1.000086347");
2325 m_apm_set_string(mapm1, cbuf);
2326
2327 m_apm_set_string(mapmr,
2328 "8.6343272312377051473047957048191159556789410110004420720659541743787375E-5");
2329
2330 m_apm_log(mapm2, 70, mapm1);
2331
2332 if (m_apm_compare(mapm2, mapmr) != 0)
2333 {
2334 pass = 1;
2335 fprintf(stdout,
2336 "***** FAILED \'log_near_1\' #1\n");
2337 }
2338
2339
2340 strcpy(cbuf, "0.99997206");
2341 m_apm_set_string(mapm1, cbuf);
2342
2343 m_apm_set_string(mapmr,
2344 "-2.7940390329070546415846319819251647850034442803891202616388034131872278E-5");
2345
2346 m_apm_log(mapm2, 70, mapm1);
2347
2348 if (m_apm_compare(mapm2, mapmr) != 0)
2349 {
2350 pass = 1;
2351 fprintf(stdout,
2352 "***** FAILED \'log_near_1\' #2\n");
2353 }
2354
2355 if (pass == 0)
2356 fprintf(stdout,"... LOG_NEAR_1 function passes\n");
2357
2358 m_apm_free(mapm1);
2359 m_apm_free(mapm2);
2360 m_apm_free(mapmr);
2361
2362 return(pass);
2363 }
2364
2365 /****************************************************************************/
2366
2367 int mapm_test_asin_near_0()
2368 {
2369 int dplaces, pass;
2370 char cbuf[32];
2371 M_APM mapm1, mapm2, mapm3, mapm4, mapm_pi2;
2372
2373 fprintf(stdout, "Validating the INV TRIG NEAR 0 functions ... \n");
2374
2375 mapm1 = m_apm_init();
2376 mapm2 = m_apm_init();
2377 mapm3 = m_apm_init();
2378 mapm4 = m_apm_init();
2379 mapm_pi2 = m_apm_init();
2380 pass = 0;
2381 dplaces = 320;
2382
2383 /*
2384 * exercise the special cases in the inverse sin/cos/tan functions:
2385 *
2386 * atan at 0, ~0
2387 * asin at +1, -1, 0, ~0
2388 * acos at +1, -1, 0, ~0
2389 */
2390
2391 fprintf(stdout, "Calculate PI/2 to %d decimal places ... \n", dplaces);
2392
2393 /* PI/2 = acos(-1) / 2 */
2394
2395 m_apm_negate(mapm1, MM_One);
2396 m_apm_acos(mapm2, (dplaces + 8), mapm1);
2397 m_apm_divide(mapm_pi2, dplaces, mapm2, MM_Two);
2398
2399 /* verify atan(0) == 0 */
2400
2401 m_apm_atan(mapm2, dplaces, MM_Zero);
2402
2403 if (m_apm_compare(mapm2, MM_Zero) != 0)
2404 {
2405 pass = 1;
2406 fprintf(stdout, "***** FAILED atan(0) \n");
2407 }
2408
2409 /* check 'arctan_near_0' function */
2410
2411 strcpy(cbuf, "-6.1423379E-7");
2412 m_apm_set_string(mapm1, cbuf);
2413 m_apm_set_string(mapm3,
2414 "-6.1423378999992275331387438418270080480455433775102970001007888571614664E-7");
2415
2416 m_apm_atan(mapm2, 70, mapm1);
2417
2418 if (m_apm_compare(mapm2, mapm3) != 0)
2419 {
2420 pass = 1;
2421 fprintf(stdout, "***** FAILED atan_near_0 \n");
2422 }
2423
2424 /* verify asin(0) == 0 */
2425
2426 m_apm_asin(mapm2, dplaces, MM_Zero);
2427
2428 if (m_apm_compare(mapm2, MM_Zero) != 0)
2429 {
2430 pass = 1;
2431 fprintf(stdout, "***** FAILED asin(0) \n");
2432 }
2433
2434 /* verify asin(1) == PI/2 */
2435
2436 m_apm_asin(mapm2, dplaces, MM_One);
2437
2438 if (m_apm_compare(mapm2, mapm_pi2) != 0)
2439 {
2440 pass = 1;
2441 fprintf(stdout, "***** FAILED asin(1) = PI / 2 \n");
2442 }
2443
2444 /* verify asin(-1) == -PI/2 */
2445
2446 m_apm_negate(mapm3, MM_One);
2447 m_apm_asin(mapm2, dplaces, mapm3);
2448 m_apm_negate(mapm1, mapm_pi2);
2449
2450 if (m_apm_compare(mapm2, mapm1) != 0)
2451 {
2452 pass = 1;
2453 fprintf(stdout, "***** FAILED asin(-1) = -PI / 2 \n");
2454 }
2455
2456 /* check 'arcsin_near_0' function */
2457
2458 strcpy(cbuf, "2.9402782364E-6");
2459 m_apm_set_string(mapm1, cbuf);
2460 m_apm_set_string(mapm3,
2461 "2.9402782364042365665958943683202914596951723300870666725716858547635946E-6");
2462
2463 m_apm_asin(mapm2, 70, mapm1);
2464
2465 if (m_apm_compare(mapm2, mapm3) != 0)
2466 {
2467 pass = 1;
2468 fprintf(stdout, "***** FAILED asin_near_0 \n");
2469 }
2470
2471 /* verify acos(1) == 0 */
2472
2473 m_apm_acos(mapm2, dplaces, MM_One);
2474
2475 if (m_apm_compare(mapm2, MM_Zero) != 0)
2476 {
2477 pass = 1;
2478 fprintf(stdout, "***** FAILED acos(1) = 0 \n");
2479 }
2480
2481 /* verify acos(0) == PI / 2 */
2482
2483 m_apm_acos(mapm2, dplaces, MM_Zero);
2484
2485 if (m_apm_compare(mapm2, mapm_pi2) != 0)
2486 {
2487 pass = 1;
2488 fprintf(stdout, "***** FAILED acos(0) = PI / 2 \n");
2489 }
2490
2491
2492 /* check 'arccos_near_0' function */
2493
2494 strcpy(cbuf, "-7.2804311E-5");
2495 m_apm_set_string(mapm1, cbuf);
2496 m_apm_set_string(mapm3,
2497 "1.5708691311059609353812900497247032950115203749494458830119113661688305");
2498
2499 m_apm_acos(mapm2, 70, mapm1);
2500
2501 if (m_apm_compare(mapm2, mapm3) != 0)
2502 {
2503 pass = 1;
2504 fprintf(stdout, "***** FAILED acos_near_0 \n");
2505 }
2506
2507
2508 /* verify (acos(x) + acos(-x)) / 2 = PI / 2 where x = 8.635254948248924E-7 */
2509
2510 strcpy(cbuf, "8.635254948248924E-7");
2511 m_apm_set_string(mapm1, cbuf);
2512 m_apm_negate(mapm2, mapm1);
2513
2514 m_apm_acos(mapm3, (dplaces + 8), mapm1);
2515 m_apm_acos(mapm4, (dplaces + 8), mapm2);
2516 m_apm_add(mapm1, mapm3, mapm4);
2517 m_apm_divide(mapm2, dplaces, mapm1, MM_Two);
2518
2519 if (m_apm_compare(mapm2, mapm_pi2) != 0)
2520 {
2521 pass = 1;
2522 fprintf(stdout, "***** FAILED (acos(x) + acos(-x)) / 2 == PI / 2 \n");
2523 }
2524
2525 if (pass == 0)
2526 fprintf(stdout,"... INV TRIG NEAR 0 functions pass\n");
2527
2528 m_apm_free(mapm1);
2529 m_apm_free(mapm2);
2530 m_apm_free(mapm3);
2531 m_apm_free(mapm4);
2532 m_apm_free(mapm_pi2);
2533
2534 return(pass);
2535 }
2536
2537 /****************************************************************************/
2538
2539 /*
2540 * verify sin^2 + cos^2 = 1
2541 */
2542
2543 int mapm_test_sin2_cos2()
2544 {
2545 int dplaces, pass;
2546 M_APM mapm1, mapm2, mapm3, mapm4, mapm5;
2547
2548 fprintf(stdout, "Validating SIN^2 + COS^2 = 1 ... \n");
2549
2550 mapm1 = m_apm_init();
2551 mapm2 = m_apm_init();
2552 mapm3 = m_apm_init();
2553 mapm4 = m_apm_init();
2554 mapm5 = m_apm_init();
2555
2556 pass = 0;
2557 dplaces = 150;
2558
2559 while (1)
2560 {
2561 m_apm_set_string(mapm3, "0.6216566124128346");
2562
2563 fprintf(stdout,
2564 "Calculate sin(0.6216566124128346) to %d decimal places ... \n", dplaces);
2565
2566 m_apm_sin(mapm1, (dplaces + 6), mapm3);
2567
2568 fprintf(stdout,
2569 "Calculate cos(0.6216566124128346) to %d decimal places ... \n", dplaces);
2570
2571 m_apm_cos(mapm2, (dplaces + 6), mapm3);
2572
2573 m_apm_multiply(mapm4, mapm1, mapm1);
2574 m_apm_multiply(mapm5, mapm2, mapm2);
2575
2576 m_apm_add(mapm3, mapm4, mapm5);
2577 m_apm_round(mapm4, dplaces, mapm3);
2578
2579 if (m_apm_compare(mapm4, MM_One) != 0)
2580 {
2581 pass = 1;
2582 fprintf(stdout, "***** FAIL: SIN^2 + COS^2 != 1 \n");
2583 }
2584
2585 m_apm_set_string(mapm3, "2.1073765187265308");
2586
2587 fprintf(stdout,
2588 "Calculate sin(2.1073765187265308) to %d decimal places ... \n", dplaces);
2589 m_apm_sin(mapm1, (dplaces + 6), mapm3);
2590
2591 fprintf(stdout,
2592 "Calculate cos(2.1073765187265308) to %d decimal places ... \n", dplaces);
2593 m_apm_cos(mapm2, (dplaces + 6), mapm3);
2594
2595 m_apm_multiply(mapm4, mapm1, mapm1);
2596 m_apm_multiply(mapm5, mapm2, mapm2);
2597
2598 m_apm_add(mapm3, mapm4, mapm5);
2599 m_apm_round(mapm4, dplaces, mapm3);
2600
2601 if (m_apm_compare(mapm4, MM_One) != 0)
2602 {
2603 pass = 1;
2604 fprintf(stdout, "***** FAIL: SIN^2 + COS^2 != 1 \n");
2605 }
2606
2607 if (dplaces > 480)
2608 break;
2609
2610 dplaces += 175;
2611 }
2612
2613 if (pass == 0)
2614 fprintf(stdout,"... SIN^2 + COS^2 = 1 pass\n");
2615
2616 m_apm_free(mapm1);
2617 m_apm_free(mapm2);
2618 m_apm_free(mapm3);
2619 m_apm_free(mapm4);
2620 m_apm_free(mapm5);
2621
2622 return(pass);
2623 }
2624
2625 /****************************************************************************/
2626
2627 /*
2628 * this test forces the FFT multiply to malloc additional working
2629 * arrays for temp data
2630 */
2631
2632 int mapm_test_fft_digits()
2633 {
2634 int digits, pass;
2635 M_APM mapm1, mapm2, mapm3, mapm4;
2636
2637 if (M_get_sizeof_int() >= 4)
2638 digits = 12000;
2639 else
2640 digits = 900;
2641
2642 fprintf(stdout, "Validating FFT multiply with %d digits ... \n", digits);
2643
2644 mapm1 = m_apm_init();
2645 mapm2 = m_apm_init();
2646 mapm3 = m_apm_init();
2647 mapm4 = m_apm_init();
2648 pass = 0;
2649
2650 /*
2651 * let a = sqrt(27.3)
2652 * let b = cbrt(-19.6)
2653 *
2654 * then verify ((a * b) ^ 2) * b = (-19.6 * 27.3)
2655 */
2656
2657 m_apm_set_string(mapm1, "27.3");
2658 m_apm_set_string(mapm2, "-19.6");
2659
2660 fprintf(stdout, "Calculate sqrt(27.3) to %d decimal places ... \n", digits);
2661
2662 m_apm_sqrt(mapm3, (digits + 6), mapm1);
2663
2664 fprintf(stdout, "Calculate cbrt(-19.6) to %d decimal places ... \n", digits);
2665
2666 m_apm_cbrt(mapm4, (digits + 6), mapm2);
2667
2668 m_apm_multiply(mapm1, mapm3, mapm4);
2669 m_apm_multiply(mapm2, mapm1, mapm1);
2670 m_apm_multiply(mapm1, mapm2, mapm4);
2671
2672 m_apm_round(mapm2, digits, mapm1);
2673
2674 m_apm_set_string(mapm1, "-535.08"); /* 27.3 * -19.6 */
2675
2676 if (m_apm_compare(mapm2, mapm1) != 0)
2677 {
2678 pass = 1;
2679 fprintf(stdout, "***** FAILED: FFT multiply with %d digits\n", digits);
2680 }
2681
2682 if (pass == 0)
2683 fprintf(stdout,"... %d digit FFT multiply passes \n", digits);
2684
2685 m_apm_free(mapm1);
2686 m_apm_free(mapm2);
2687 m_apm_free(mapm3);
2688 m_apm_free(mapm4);
2689
2690 return(pass);
2691 }
2692
2693 /****************************************************************************/
2694
2695 /*
2696 * this code forces the _add function to init the
2697 * static variables in that module. no pass/fail needed here.
2698 */
2699
2700 int mapm_test_add_init()
2701 {
2702 M_APM mapm0;
2703
2704 mapm0 = m_apm_init();
2705
2706 m_apm_add(mapm0, MM_Ten, MM_Three);
2707
2708 m_apm_free(mapm0);
2709
2710 return(0);
2711 }
2712
2713 /****************************************************************************/
2714
2715 /*
2716 * this test forces relatively rare corner conditions in the
2717 * Knuth division algorithm.
2718 */
2719
2720 int mapm_test_special_div()
2721 {
2722 int digits, pass;
2723 M_APM mapm1, mapm2, mapm3, mapm4, mapm5, mapm6;
2724
2725 digits = 170;
2726
2727 mapm1 = m_apm_init();
2728 mapm2 = m_apm_init();
2729 mapm3 = m_apm_init();
2730 mapm4 = m_apm_init();
2731 mapm5 = m_apm_init();
2732 mapm6 = m_apm_init();
2733 pass = 0;
2734
2735 fprintf(stdout, "Calculate sqrt(3) to %d decimal places ... \n", digits);
2736
2737 m_apm_sqrt(mapm4, (digits + 6), MM_Three);
2738
2739 /*
2740 * let N = sqrt(3)
2741 *
2742 * compute M = N / (N * 1.00000001)
2743 *
2744 * then M * 1.0000001 should == 1
2745 */
2746
2747 m_apm_set_string(mapm1, "1.00000001");
2748 m_apm_multiply(mapm3, mapm1, mapm4);
2749
2750 m_apm_divide(mapm2, (digits + 6), mapm4, mapm4); /* n / n == 1 */
2751
2752 if (m_apm_compare(mapm2, MM_One) != 0)
2753 {
2754 pass = 1;
2755 fprintf(stdout,
2756 "***** FAILED: Special Division Check #1 %d digits\n", digits);
2757 }
2758
2759 m_apm_negate(mapm6, mapm3);
2760 m_apm_set_long(mapm5, -1L);
2761 m_apm_divide(mapm2, (digits + 6), mapm6, mapm3); /* n / n == 1 */
2762
2763 if (m_apm_compare(mapm2, mapm5) != 0)
2764 {
2765 pass = 1;
2766 fprintf(stdout,
2767 "***** FAILED: Special Division Check #2 %d digits\n", digits);
2768 }
2769
2770 m_apm_divide(mapm2, (digits + 6), mapm4, mapm3);
2771
2772 m_apm_multiply(mapm5, mapm2, mapm1);
2773 m_apm_round(mapm2, digits, mapm5);
2774
2775 if (m_apm_compare(mapm2, MM_One) != 0)
2776 {
2777 pass = 1;
2778 fprintf(stdout,
2779 "***** FAILED: Special Division Check #3 %d digits\n", digits);
2780 }
2781
2782 if (pass == 0)
2783 fprintf(stdout,"... %d digit Special Division Check passes \n", digits);
2784
2785 m_apm_free(mapm1);
2786 m_apm_free(mapm2);
2787 m_apm_free(mapm3);
2788 m_apm_free(mapm4);
2789 m_apm_free(mapm5);
2790 m_apm_free(mapm6);
2791
2792 return(pass);
2793 }
2794
2795 /****************************************************************************/
2796
2797 int mapm_test_special_fpf()
2798 {
2799 int digits, kk, pass;
2800 char *out_string, buf[64];
2801 M_APM mapm1, mapm2, mapm3, mapm4;
2802
2803 digits = 156;
2804
2805 mapm1 = m_apm_init();
2806 mapm2 = m_apm_init();
2807 mapm3 = m_apm_init();
2808 mapm4 = m_apm_init();
2809 pass = 0;
2810
2811 fprintf(stdout, "Validating Fixed Point Formatting ... \n");
2812
2813 strcpy(buf, " +12.57564854193421");
2814 buf[2] = '\t'; /* put some tabs in the string */
2815 buf[7] = '\t';
2816 m_apm_set_string(mapm1, buf);
2817 m_apm_sqrt(mapm3, digits, mapm1);
2818
2819 digits = 150;
2820
2821 while (1)
2822 {
2823 out_string = m_apm_to_fixpt_stringexp(digits, mapm3, '.', 0, 0);
2824
2825 if (out_string == NULL)
2826 {
2827 fprintf(stderr, "VALIDATE: Out of Memory \n");
2828 exit(100);
2829 }
2830
2831 m_apm_set_string(mapm2, out_string);
2832
2833 m_apm_subtract(mapm1, mapm3, mapm2);
2834
2835 kk = m_apm_exponent(mapm1);
2836
2837 if (digits <= 50)
2838 fprintf(stdout, "[%s] \n", out_string);
2839
2840 if ((-kk) <= digits)
2841 {
2842 pass = 1;
2843 fprintf(stdout,
2844 "***** FAILED: Fixed Point Formatting: digits = %d \n", digits);
2845 }
2846
2847 free(out_string);
2848
2849 if (digits == -1)
2850 break;
2851
2852 digits--;
2853 }
2854
2855 if ((out_string = (char *)malloc(1024)) == NULL)
2856 {
2857 fprintf(stderr, "VALIDATE: Out of Memory \n");
2858 exit(100);
2859 }
2860
2861 m_apm_set_string(mapm1, "7.42812659128912E+1469");
2862 m_apm_sqrt(mapm2, 750, mapm1);
2863 m_apm_floor(mapm1, mapm2);
2864 m_apm_set_string(mapm3, "1.0E+23");
2865 m_apm_multiply(mapm4, mapm3, mapm1);
2866 m_apm_to_integer_string(out_string, mapm4);
2867 m_apm_set_string(mapm3, out_string);
2868
2869 kk = m_apm_compare(mapm3, mapm4);
2870
2871 if (kk != 0)
2872 {
2873 pass = 1;
2874 fprintf(stdout, "***** FAILED: Fixed Point Formatting \n");
2875 }
2876
2877 free(out_string);
2878
2879 if (pass == 0)
2880 fprintf(stdout,"... Fixed Point Formatting passes \n");
2881
2882 m_apm_free(mapm1);
2883 m_apm_free(mapm2);
2884 m_apm_free(mapm3);
2885 m_apm_free(mapm4);
2886
2887 return(pass);
2888 }
2889
2890 /****************************************************************************/
2891
2892 int mapm_test_special_set()
2893 {
2894 int digits, kk, pass, alldigits;
2895 char *out_string, buf[64];
2896 M_APM mapm1, mapm2, mapm3;
2897
2898 digits = 350;
2899 alldigits = -1;
2900
2901 mapm1 = m_apm_init();
2902 mapm2 = m_apm_init();
2903 mapm3 = m_apm_init();
2904 pass = 0;
2905
2906 fprintf(stdout, "Validating Set String & Fixed Point Formatting ... \n");
2907
2908 strcpy(buf, " -7.52164854093421E+12");
2909 m_apm_set_string(mapm1, buf);
2910
2911 while (1)
2912 {
2913 m_apm_cbrt(mapm3, digits, mapm1);
2914
2915 if (digits > 1400)
2916 break;
2917
2918 out_string = m_apm_to_fixpt_stringexp(alldigits, mapm3, '.', 0, 0);
2919
2920 if (out_string == NULL)
2921 {
2922 fprintf(stderr, "VALIDATE: Out of Memory \n");
2923 exit(100);
2924 }
2925
2926 m_apm_set_string(mapm2, out_string);
2927
2928 kk = m_apm_compare(mapm2, mapm3);
2929
2930 if (kk != 0)
2931 {
2932 pass = 1;
2933 fprintf(stdout,
2934 "***** FAILED: Set String & Fixed Point Formatting \n");
2935 }
2936
2937 free(out_string);
2938
2939 digits += 265;
2940 }
2941
2942 if (pass == 0)
2943 fprintf(stdout,"... Set String & Fixed Point Formatting passes \n");
2944
2945 m_apm_free(mapm1);
2946 m_apm_free(mapm2);
2947 m_apm_free(mapm3);
2948
2949 return(pass);
2950 }
2951
2952 /****************************************************************************/
2953
2954 /*
2955 *
2956 * start of main test procedure
2957 *
2958 */
2959
2960 int main(int argc, char *argv[])
2961 {
2962 int total_pass;
2963 M_APM tmp;
2964
2965 total_pass = 0;
2966
2967 tmp = m_apm_init();
2968
2969 library_version_check();
2970
2971 mapm_show_random();
2972 m_apm_trim_mem_usage();
2973 m_apm_set_random_seed("649378126582104");
2974 mapm_show_random();
2975
2976 total_pass |= mapm_test_sqrt();
2977 total_pass |= mapm_test_cbrt();
2978
2979 total_pass |= mapm_test_sin();
2980 total_pass |= mapm_test_cos();
2981 total_pass |= mapm_test_sin_cos();
2982 total_pass |= mapm_test_sin2_cos2();
2983 total_pass |= mapm_test_tan();
2984
2985 total_pass |= mapm_test_asin();
2986 total_pass |= mapm_test_acos();
2987 total_pass |= mapm_test_atan();
2988 total_pass |= mapm_test_atan2();
2989 total_pass |= mapm_test_asin_near_0();
2990
2991 m_apm_trim_mem_usage();
2992
2993 total_pass |= mapm_test_sinh();
2994 total_pass |= mapm_test_cosh();
2995 total_pass |= mapm_test_tanh();
2996
2997 total_pass |= mapm_test_asinh();
2998 total_pass |= mapm_test_acosh();
2999 total_pass |= mapm_test_atanh();
3000
3001 m_apm_free(tmp);
3002 m_apm_free_all_mem();
3003 tmp = m_apm_init();
3004
3005 total_pass |= mapm_test_exp();
3006 total_pass |= mapm_test_exp_2();
3007 total_pass |= mapm_test_log();
3008 total_pass |= mapm_test_log10();
3009 total_pass |= mapm_test_log_near_1();
3010 total_pass |= mapm_test_pow();
3011
3012 total_pass |= mapm_test_factorial();
3013 total_pass |= mapm_test_gcd();
3014 total_pass |= mapm_test_floor_ceil();
3015
3016 m_apm_trim_mem_usage();
3017
3018 total_pass |= mapm_test_PI();
3019 total_pass |= mapm_test_exp_pwr2();
3020 total_pass |= mapm_test_log_2();
3021 total_pass |= mapm_test_fft_digits();
3022
3023 m_apm_free(tmp);
3024 m_apm_free_all_mem();
3025
3026 total_pass |= mapm_test_add_init();
3027
3028 m_apm_trim_mem_usage();
3029 tmp = m_apm_init();
3030
3031 total_pass |= mapm_test_special_div();
3032 total_pass |= mapm_test_special_fpf();
3033 total_pass |= mapm_test_special_set();
3034
3035 if (total_pass == 0)
3036 fprintf(stdout, "\nvalidate : PASS \n");
3037 else
3038 fprintf(stdout, "\nvalidate : FAIL \n");
3039
3040 m_apm_free(tmp);
3041 m_apm_free_all_mem();
3042
3043 exit(total_pass);
3044 }
3045
3046 /****************************************************************************/
3047