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 */
49void	M_check_PI_places(int places)
50{
51int     dplaces;
52
53dplaces = places + 2;
54
55if (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 */
110void	M_calculate_PI_AGM(M_APM outv, int places)
111{
112M_APM   tmp1, tmp2, a0, b0, c0, a1, b1, sum, pow_2;
113int     dplaces, nn;
114
115tmp1  = M_get_stack_var();
116tmp2  = M_get_stack_var();
117a0    = M_get_stack_var();
118b0    = M_get_stack_var();
119c0    = M_get_stack_var();
120a1    = M_get_stack_var();
121b1    = M_get_stack_var();
122sum   = M_get_stack_var();
123pow_2 = M_get_stack_var();
124
125dplaces = places + 16;
126
127m_apm_copy(a0, MM_One);
128m_apm_copy(sum, MM_One);
129m_apm_copy(pow_2, MM_Four);
130m_apm_sqrt(b0, dplaces, MM_0_5);        /* sqrt(0.5) */
131
132while (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
166m_apm_add(tmp1, a1, b1);
167m_apm_multiply(tmp2, tmp1, tmp1);
168m_apm_divide(tmp1, dplaces, tmp2, sum);
169m_apm_round(outv, places, tmp1);
170
171M_restore_stack(9);
172}
173/****************************************************************************/
174