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-genie.h"
28 #include "a68g-prelude.h"
29 #include "a68g-double.h"
30 #include "a68g-mp.h"
31
32 //! @brief Return "pi" with "digs" precision, using Borwein & Borwein AGM.
33
34 MP_T *mp_pi (NODE_T * p, MP_T * api, int mod, int digs)
35 {
36 int gdigs = FUN_DIGITS (digs);
37 if (gdigs > A68_MP (mp_pi_size)) {
38 //
39 // No luck with the kept value, hence we generate a longer "pi".
40 // Calculate "pi" using the Borwein & Borwein AGM algorithm.
41 // This AGM doubles the numbers of digits every iteration.
42 //
43 a68_free (A68_MP (mp_pi));
44 a68_free (A68_MP (mp_half_pi));
45 a68_free (A68_MP (mp_two_pi));
46 a68_free (A68_MP (mp_sqrt_two_pi));
47 a68_free (A68_MP (mp_sqrt_pi));
48 a68_free (A68_MP (mp_ln_pi));
49 a68_free (A68_MP (mp_180_over_pi));
50 a68_free (A68_MP (mp_pi_over_180));
51 //
52 ADDR_T pop_sp = A68_SP;
53 MP_T *pi_g = nil_mp (p, gdigs);
54 MP_T *two = lit_mp (p, 2, 0, gdigs);
55 MP_T *x_g = lit_mp (p, 2, 0, gdigs);
56 MP_T *y_g = nil_mp (p, gdigs);
57 MP_T *u_g = nil_mp (p, gdigs);
58 MP_T *v_g = nil_mp (p, gdigs);
59 (void) sqrt_mp (p, x_g, x_g, gdigs);
60 (void) add_mp (p, pi_g, x_g, two, gdigs);
61 (void) sqrt_mp (p, y_g, x_g, gdigs);
62 BOOL_T iterate = A68_TRUE;
63 while (iterate) {
64 // New x.
65 (void) sqrt_mp (p, u_g, x_g, gdigs);
66 (void) rec_mp (p, v_g, u_g, gdigs);
67 (void) add_mp (p, u_g, u_g, v_g, gdigs);
68 (void) half_mp (p, x_g, u_g, gdigs);
69 // New pi.
70 (void) plus_one_mp (p, u_g, x_g, gdigs);
71 (void) plus_one_mp (p, v_g, y_g, gdigs);
72 (void) div_mp (p, u_g, u_g, v_g, gdigs);
73 (void) mul_mp (p, v_g, pi_g, u_g, gdigs);
74 // Done yet?.
75 if (same_mp (p, v_g, pi_g, gdigs)) {
76 iterate = A68_FALSE;
77 } else {
78 (void) move_mp (pi_g, v_g, gdigs);
79 // New y.
80 (void) sqrt_mp (p, u_g, x_g, gdigs);
81 (void) rec_mp (p, v_g, u_g, gdigs);
82 (void) mul_mp (p, u_g, y_g, u_g, gdigs);
83 (void) add_mp (p, u_g, u_g, v_g, gdigs);
84 (void) plus_one_mp (p, v_g, y_g, gdigs);
85 (void) div_mp (p, y_g, u_g, v_g, gdigs);
86 }
87 }
88 // Keep the result for future restore.
89 (void) shorten_mp (p, api, digs, pi_g, gdigs);
90 A68_MP (mp_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
91 (void) move_mp (A68_MP (mp_pi), api, digs);
92 A68_MP (mp_half_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
93 (void) half_mp (p, A68_MP (mp_half_pi), api, digs);
94 A68_MP (mp_sqrt_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
95 (void) sqrt_mp (p, A68_MP (mp_sqrt_pi), api, digs);
96 A68_MP (mp_ln_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
97 (void) ln_mp (p, A68_MP (mp_ln_pi), api, digs);
98 A68_MP (mp_two_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
99 (void) mul_mp_digit (p, A68_MP (mp_two_pi), api, (MP_T) 2, digs);
100 A68_MP (mp_sqrt_two_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
101 (void) sqrt_mp (p, A68_MP (mp_sqrt_two_pi), A68_MP (mp_two_pi), digs);
102 A68_MP (mp_pi_over_180) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
103 (void) div_mp_digit (p, A68_MP (mp_pi_over_180), api, 180, digs);
104 A68_MP (mp_180_over_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
105 (void) rec_mp (p, A68_MP (mp_180_over_pi), A68_MP (mp_pi_over_180), digs);
106 A68_MP (mp_pi_size) = gdigs;
107 A68_SP = pop_sp;
108 }
109 switch (mod) {
110 case MP_PI: return move_mp (api, A68_MP (mp_pi), digs);
111 case MP_HALF_PI: return move_mp (api, A68_MP (mp_half_pi), digs);
112 case MP_TWO_PI: return move_mp (api, A68_MP (mp_two_pi), digs);
113 case MP_SQRT_TWO_PI: return move_mp (api, A68_MP (mp_sqrt_two_pi), digs);
114 case MP_SQRT_PI: return move_mp (api, A68_MP (mp_sqrt_pi), digs);
115 case MP_LN_PI: return move_mp (api, A68_MP (mp_ln_pi), digs);
116 case MP_180_OVER_PI: return move_mp (api, A68_MP (mp_180_over_pi), digs);
117 case MP_PI_OVER_180: return move_mp (api, A68_MP (mp_pi_over_180), digs);
118 default: return NaN_MP; // Should not be here.
119 }
120 }
121
