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-2023 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 if (gdigs > A68_MP (mp_pi_size)) {
35 // No luck with the cached value, so we compute a longer "pi".
36 // Calculate "pi" using a Borwein & Borwein AGM algorithm that doubles
37 // the numbers of digits on each iteration.
38 //
39 // J. M. Borwein, P. B. Borwein. "Pi and the AGM – A Study in Analytic Number
40 // Theory and Computational Complexity", Wiley, New York, 1987.
41 a68_free (A68_MP (mp_pi));
42 a68_free (A68_MP (mp_half_pi));
43 a68_free (A68_MP (mp_two_pi));
44 a68_free (A68_MP (mp_sqrt_two_pi));
45 a68_free (A68_MP (mp_sqrt_pi));
46 a68_free (A68_MP (mp_ln_pi));
47 a68_free (A68_MP (mp_180_over_pi));
48 a68_free (A68_MP (mp_pi_over_180));
49 ADDR_T pop_sp = A68_SP;
50 MP_T *pi_g = nil_mp (p, gdigs);
51 MP_T *two = lit_mp (p, 2, 0, gdigs);
52 MP_T *x_g = lit_mp (p, 2, 0, gdigs);
53 MP_T *y_g = nil_mp (p, gdigs);
54 MP_T *u_g = nil_mp (p, gdigs);
55 MP_T *v_g = nil_mp (p, gdigs);
56 (void) sqrt_mp (p, x_g, x_g, gdigs);
57 (void) add_mp (p, pi_g, x_g, two, gdigs);
58 (void) sqrt_mp (p, y_g, x_g, gdigs);
59 BOOL_T iterate = A68_TRUE;
60 while (iterate) {
61 // New x.
62 (void) sqrt_mp (p, u_g, x_g, gdigs);
63 (void) rec_mp (p, v_g, u_g, gdigs);
64 (void) add_mp (p, u_g, u_g, v_g, gdigs);
65 (void) half_mp (p, x_g, u_g, gdigs);
66 // New pi.
67 (void) plus_one_mp (p, u_g, x_g, gdigs);
68 (void) plus_one_mp (p, v_g, y_g, gdigs);
69 (void) div_mp (p, u_g, u_g, v_g, gdigs);
70 (void) mul_mp (p, v_g, pi_g, u_g, gdigs);
71 // Done yet?.
72 if (same_mp (p, v_g, pi_g, gdigs)) {
73 iterate = A68_FALSE;
74 } else {
75 (void) move_mp (pi_g, v_g, gdigs);
76 // New y.
77 (void) sqrt_mp (p, u_g, x_g, gdigs);
78 (void) rec_mp (p, v_g, u_g, gdigs);
79 (void) mul_mp (p, u_g, y_g, u_g, gdigs);
80 (void) add_mp (p, u_g, u_g, v_g, gdigs);
81 (void) plus_one_mp (p, v_g, y_g, gdigs);
82 (void) div_mp (p, y_g, u_g, v_g, gdigs);
83 }
84 }
85 // Keep the result for future restore.
86 (void) shorten_mp (p, api, digs, pi_g, gdigs);
87 A68_MP (mp_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
88 (void) move_mp (A68_MP (mp_pi), api, digs);
89 A68_MP (mp_half_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
90 (void) half_mp (p, A68_MP (mp_half_pi), api, digs);
91 A68_MP (mp_sqrt_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
92 (void) sqrt_mp (p, A68_MP (mp_sqrt_pi), api, digs);
93 A68_MP (mp_ln_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
94 (void) ln_mp (p, A68_MP (mp_ln_pi), api, digs);
95 A68_MP (mp_two_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
96 (void) mul_mp_digit (p, A68_MP (mp_two_pi), api, (MP_T) 2, digs);
97 A68_MP (mp_sqrt_two_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
98 (void) sqrt_mp (p, A68_MP (mp_sqrt_two_pi), A68_MP (mp_two_pi), digs);
99 A68_MP (mp_pi_over_180) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
100 (void) div_mp_digit (p, A68_MP (mp_pi_over_180), api, 180, digs);
101 A68_MP (mp_180_over_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
102 (void) rec_mp (p, A68_MP (mp_180_over_pi), A68_MP (mp_pi_over_180), digs);
103 A68_MP (mp_pi_size) = gdigs;
104 A68_SP = pop_sp;
105 }
106 switch (mod) {
107 case MP_PI: {
108 return move_mp (api, A68_MP (mp_pi), digs);
109 }
110 case MP_HALF_PI: {
111 return move_mp (api, A68_MP (mp_half_pi), digs);
112 }
113 case MP_TWO_PI: {
114 return move_mp (api, A68_MP (mp_two_pi), digs);
115 }
116 case MP_SQRT_TWO_PI: {
117 return move_mp (api, A68_MP (mp_sqrt_two_pi), digs);
118 }
119 case MP_SQRT_PI: {
120 return move_mp (api, A68_MP (mp_sqrt_pi), digs);
121 }
122 case MP_LN_PI: {
123 return move_mp (api, A68_MP (mp_ln_pi), digs);
124 }
125 case MP_180_OVER_PI: {
126 return move_mp (api, A68_MP (mp_180_over_pi), digs);
127 }
128 case MP_PI_OVER_180: {
129 return move_mp (api, A68_MP (mp_pi_over_180), digs);
130 }
131 default: {
132 return NaN_MP; // Should not be here.
133 }
134 }
135 }
136
