mp-pi.c
1 //! @file mp-pi.c
2 //! @author J. Marcel van der Veer
3
4 //! @section Copyright
5 //!
6 //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
7 //! Copyright 2001-2026 J. Marcel van der Veer [algol68g@xs4all.nl].
8
9 //! @section License
10 //!
11 //! This program is free software; you can redistribute it and/or modify it
12 //! under the terms of the GNU General Public License as published by the
13 //! Free Software Foundation; either version 3 of the License, or
14 //! (at your option) any later version.
15 //!
16 //! This program is distributed in the hope that it will be useful, but
17 //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
19 //! more details. You should have received a copy of the GNU General Public
20 //! License along with this program. If not, see [http://www.gnu.org/licenses/].
21
22 //! @section Synopsis
23 //!
24 //! [LONG] LONG REAL value of pi by AGM.
25
26 #include "a68g.h"
27 #include "a68g-mp.h"
28
29 //! @brief Return "pi" with "digs" precision, using a Borwein & Borwein AGM.
30
31 MP_T *mp_pi(NODE_T *p, MP_T *api, int mod, int digs)
32 {
33 int gdigs = FUN_DIGITS(digs);
34 //
35 if (gdigs > A68G_MP(mp_pi_size))
36 {
37 //
38 // No luck with the cached value, so we compute a longer "pi".
39 // Calculate "pi" using a Borwein & Borwein AGM algorithm that doubles
40 // the numbers of digits on each iteration.
41 //
42 // J. M. Borwein, P. B. Borwein. "Pi and the AGM – A Study in Analytic Number
43 // Theory and Computational Complexity", Wiley, New York, 1987.
44 //
45 a68g_free(A68G_MP(mp_pi));
46 a68g_free(A68G_MP(mp_half_pi));
47 a68g_free(A68G_MP(mp_two_pi));
48 a68g_free(A68G_MP(mp_sqrt_two_pi));
49 a68g_free(A68G_MP(mp_sqrt_pi));
50 a68g_free(A68G_MP(mp_ln_pi));
51 a68g_free(A68G_MP(mp_180_over_pi));
52 a68g_free(A68G_MP(mp_pi_over_180));
53 ADDR_T pop_sp = A68G_SP;
54 MP_T *pi_g = nil_mp(p, gdigs);
55 MP_T *two = lit_mp(p, 2, 0, gdigs);
56 MP_T *x_g = lit_mp(p, 2, 0, gdigs);
57 MP_T *y_g = nil_mp(p, gdigs);
58 MP_T *u_g = nil_mp(p, gdigs);
59 MP_T *v_g = nil_mp(p, gdigs);
60 (void)sqrt_mp(p, x_g, x_g, gdigs);
61 (void)add_mp(p, pi_g, x_g, two, gdigs);
62 (void)sqrt_mp(p, y_g, x_g, gdigs);
63 BOOL_T iterate = A68G_TRUE;
64 while (iterate)
65 {
66 // New x.
67 (void)sqrt_mp(p, u_g, x_g, gdigs);
68 (void)rec_mp(p, v_g, u_g, gdigs);
69 (void)add_mp(p, u_g, u_g, v_g, gdigs);
70 (void)half_mp(p, x_g, u_g, gdigs);
71 // New pi.
72 (void)plus_one_mp(p, u_g, x_g, gdigs);
73 (void)plus_one_mp(p, v_g, y_g, gdigs);
74 (void)div_mp(p, u_g, u_g, v_g, gdigs);
75 (void)mul_mp(p, v_g, pi_g, u_g, gdigs);
76 // Done yet?
77 if (same_mp(p, v_g, pi_g, gdigs))
78 {
79 iterate = A68G_FALSE;
80 }
81 else
82 {
83 (void)move_mp(pi_g, v_g, gdigs);
84 // New y.
85 (void)sqrt_mp(p, u_g, x_g, gdigs);
86 (void)rec_mp(p, v_g, u_g, gdigs);
87 (void)mul_mp(p, u_g, y_g, u_g, gdigs);
88 (void)add_mp(p, u_g, u_g, v_g, gdigs);
89 (void)plus_one_mp(p, v_g, y_g, gdigs);
90 (void)div_mp(p, y_g, u_g, v_g, gdigs);
91 }
92 }
93 // Keep the result for future restore.
94 (void)shorten_mp(p, api, digs, pi_g, gdigs);
95 A68G_MP(mp_pi) = (MP_T *)get_heap_space(SIZE_MP(digs));
96 (void)move_mp(A68G_MP(mp_pi), api, digs);
97 A68G_MP(mp_half_pi) = (MP_T *)get_heap_space(SIZE_MP(digs));
98 (void)half_mp(p, A68G_MP(mp_half_pi), api, digs);
99 A68G_MP(mp_sqrt_pi) = (MP_T *)get_heap_space(SIZE_MP(digs));
100 (void)sqrt_mp(p, A68G_MP(mp_sqrt_pi), api, digs);
101 A68G_MP(mp_ln_pi) = (MP_T *)get_heap_space(SIZE_MP(digs));
102 (void)ln_mp(p, A68G_MP(mp_ln_pi), api, digs);
103 A68G_MP(mp_two_pi) = (MP_T *)get_heap_space(SIZE_MP(digs));
104 (void)mul_mp_digit(p, A68G_MP(mp_two_pi), api, (MP_T)2, digs);
105 A68G_MP(mp_sqrt_two_pi) = (MP_T *)get_heap_space(SIZE_MP(digs));
106 (void)sqrt_mp(p, A68G_MP(mp_sqrt_two_pi), A68G_MP(mp_two_pi), digs);
107 A68G_MP(mp_pi_over_180) = (MP_T *)get_heap_space(SIZE_MP(digs));
108 (void)div_mp_digit(p, A68G_MP(mp_pi_over_180), api, 180, digs);
109 A68G_MP(mp_180_over_pi) = (MP_T *)get_heap_space(SIZE_MP(digs));
110 (void)rec_mp(p, A68G_MP(mp_180_over_pi), A68G_MP(mp_pi_over_180), digs);
111 A68G_MP(mp_pi_size) = gdigs;
112 A68G_SP = pop_sp;
113 }
114 switch (mod)
115 {
116 case MP_PI:
117 {
118 return move_mp(api, A68G_MP(mp_pi), digs);
119 }
120 case MP_HALF_PI:
121 {
122 return move_mp(api, A68G_MP(mp_half_pi), digs);
123 }
124 case MP_TWO_PI:
125 {
126 return move_mp(api, A68G_MP(mp_two_pi), digs);
127 }
128 case MP_SQRT_TWO_PI:
129 {
130 return move_mp(api, A68G_MP(mp_sqrt_two_pi), digs);
131 }
132 case MP_SQRT_PI:
133 {
134 return move_mp(api, A68G_MP(mp_sqrt_pi), digs);
135 }
136 case MP_LN_PI:
137 {
138 return move_mp(api, A68G_MP(mp_ln_pi), digs);
139 }
140 case MP_180_OVER_PI:
141 {
142 return move_mp(api, A68G_MP(mp_180_over_pi), digs);
143 }
144 case MP_PI_OVER_180:
145 {
146 return move_mp(api, A68G_MP(mp_pi_over_180), digs);
147 }
148 default:
149 {
150 return NaN_MP; // Should not be here.
151 }
152 }
153 }
|
© 2002-2026 J.M. van der Veer (jmvdveer@xs4all.nl)
|
