1 2 /* 3 * M_APM - mapm_cpi.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: mapm_cpi.c,v 1.4 2007/12/03 01:34:29 mike Exp $ 23 * 24 * This file contains the PI related functions. 25 * 26 * $Log: mapm_cpi.c,v $ 27 * Revision 1.4 2007/12/03 01:34:29 mike 28 * Update license 29 * 30 * Revision 1.3 2002/11/05 23:10:14 mike 31 * streamline the PI AGM algorithm 32 * 33 * Revision 1.2 2002/11/03 21:56:21 mike 34 * Updated function parameters to use the modern style 35 * 36 * Revision 1.1 2001/03/25 21:01:53 mike 37 * Initial revision 38 */ 39 40 #include "m_apm_lc.h" 41 42 /****************************************************************************/ 43 /* 44 * check if our local copy of PI is precise enough 45 * for our purpose. if not, calculate PI so it's 46 * as precise as desired, accurate to 'places' decimal 47 * places. 48 */ 49 void M_check_PI_places(int places) 50 { 51 int dplaces; 52 53 dplaces = places + 2; 54 55 if (dplaces > MM_lc_PI_digits) 56 { 57 MM_lc_PI_digits = dplaces + 2; 58 59 /* compute PI using the AGM (see right below) */ 60 61 M_calculate_PI_AGM(MM_lc_PI, (dplaces + 5)); 62 63 m_apm_multiply(MM_lc_HALF_PI, MM_0_5, MM_lc_PI); 64 m_apm_multiply(MM_lc_2_PI, MM_Two, MM_lc_PI); 65 } 66 } 67 /****************************************************************************/ 68 /* 69 * Calculate PI using the AGM (Arithmetic-Geometric Mean) 70 * 71 * Init : A0 = 1 72 * B0 = 1 / sqrt(2) 73 * Sum = 1 74 * 75 * Iterate: n = 1... 76 * 77 * 78 * A = 0.5 * [ A + B ] 79 * n n-1 n-1 80 * 81 * 82 * B = sqrt [ A * B ] 83 * n n-1 n-1 84 * 85 * 86 * 87 * C = 0.5 * [ A - B ] 88 * n n-1 n-1 89 * 90 * 91 * 2 n+1 92 * Sum = Sum - C * 2 93 * n 94 * 95 * 96 * At the end when C is 'small enough' : 97 * n 98 * 99 * 2 100 * PI = 4 * A / Sum 101 * n+1 102 * 103 * -OR- 104 * 105 * 2 106 * PI = ( A + B ) / Sum 107 * n n 108 * 109 */ 110 void M_calculate_PI_AGM(M_APM outv, int places) 111 { 112 M_APM tmp1, tmp2, a0, b0, c0, a1, b1, sum, pow_2; 113 int dplaces, nn; 114 115 tmp1 = M_get_stack_var(); 116 tmp2 = M_get_stack_var(); 117 a0 = M_get_stack_var(); 118 b0 = M_get_stack_var(); 119 c0 = M_get_stack_var(); 120 a1 = M_get_stack_var(); 121 b1 = M_get_stack_var(); 122 sum = M_get_stack_var(); 123 pow_2 = M_get_stack_var(); 124 125 dplaces = places + 16; 126 127 m_apm_copy(a0, MM_One); 128 m_apm_copy(sum, MM_One); 129 m_apm_copy(pow_2, MM_Four); 130 m_apm_sqrt(b0, dplaces, MM_0_5); /* sqrt(0.5) */ 131 132 while (TRUE) 133 { 134 m_apm_add(tmp1, a0, b0); 135 m_apm_multiply(a1, MM_0_5, tmp1); 136 137 m_apm_multiply(tmp1, a0, b0); 138 m_apm_sqrt(b1, dplaces, tmp1); 139 140 m_apm_subtract(tmp1, a0, b0); 141 m_apm_multiply(c0, MM_0_5, tmp1); 142 143 /* 144 * the net 'PI' calculated from this iteration will 145 * be accurate to ~4 X the value of (c0)'s exponent. 146 * this was determined experimentally. 147 */ 148 149 nn = -4 * c0->m_apm_exponent; 150 151 m_apm_multiply(tmp1, c0, c0); 152 m_apm_multiply(tmp2, tmp1, pow_2); 153 m_apm_subtract(tmp1, sum, tmp2); 154 m_apm_round(sum, dplaces, tmp1); 155 156 if (nn >= dplaces) 157 break; 158 159 m_apm_copy(a0, a1); 160 m_apm_copy(b0, b1); 161 162 m_apm_multiply(tmp1, pow_2, MM_Two); 163 m_apm_copy(pow_2, tmp1); 164 } 165 166 m_apm_add(tmp1, a1, b1); 167 m_apm_multiply(tmp2, tmp1, tmp1); 168 m_apm_divide(tmp1, dplaces, tmp2, sum); 169 m_apm_round(outv, places, tmp1); 170 171 M_restore_stack(9); 172 } 173 /****************************************************************************/ 174