mp-gamma.c
1 //! @file mp-gamma.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@algol68genie.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 error, gamma and beta functions.
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 PROC (LONG REAL) LONG REAL erf
33
34 MP_T *erf_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
35 {
36 if (IS_ZERO_MP (x)) {
37 SET_MP_ZERO (z, digs);
38 return z;
39 } else {
40 ADDR_T pop_sp = A68G_SP;
41 // Note we need double precision!
42 BOOL_T siga = A68G_TRUE;
43 int sign = MP_SIGN (x);
44 int gdigs = FUN_DIGITS (2 * digs);
45 MP_T *y_g = nil_mp (p, gdigs);
46 MP_T *z_g = len_mp (p, x, digs, gdigs);
47 (void) abs_mp (p, z_g, z_g, gdigs);
48 (void) set_mp (y_g, gdigs * LOG_MP_RADIX, 0, gdigs);
49 (void) sqrt_mp (p, y_g, y_g, gdigs);
50 (void) sub_mp (p, y_g, z_g, y_g, gdigs);
51 if (MP_SIGN (y_g) >= 0) {
52 SET_MP_ONE (z, digs);
53 } else {
54 // Taylor expansion.
55 MP_T *p_g = nil_mp (p, gdigs);
56 MP_T *r_g = nil_mp (p, gdigs);
57 MP_T *s_g = nil_mp (p, gdigs);
58 MP_T *t_g = nil_mp (p, gdigs);
59 MP_T *u_g = nil_mp (p, gdigs);
60 (void) mul_mp (p, y_g, z_g, z_g, gdigs);
61 SET_MP_ONE (s_g, gdigs);
62 SET_MP_ONE (t_g, gdigs);
63 for (int k = 1; siga; k++) {
64 (void) mul_mp (p, t_g, y_g, t_g, gdigs);
65 (void) div_mp_digit (p, t_g, t_g, (MP_T) k, gdigs);
66 (void) div_mp_digit (p, u_g, t_g, (MP_T) (2 * k + 1), gdigs);
67 if (EVEN (k)) {
68 (void) add_mp (p, s_g, s_g, u_g, gdigs);
69 } else {
70 (void) sub_mp (p, s_g, s_g, u_g, gdigs);
71 }
72 siga = (MP_EXPONENT (s_g) - MP_EXPONENT (u_g)) < gdigs;
73 }
74 (void) mul_mp (p, r_g, z_g, s_g, gdigs);
75 (void) mul_mp_digit (p, r_g, r_g, 2, gdigs);
76 (void) mp_pi (p, p_g, MP_SQRT_PI, gdigs);
77 (void) div_mp (p, r_g, r_g, p_g, gdigs);
78 (void) shorten_mp (p, z, digs, r_g, gdigs);
79 }
80 if (sign < 0) {
81 (void) minus_mp (p, z, z, digs);
82 }
83 A68G_SP = pop_sp;
84 return z;
85 }
86 }
87
88 //! @brief PROC (LONG REAL) LONG REAL erfc
89
90 MP_T *erfc_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
91 {
92 (void) erf_mp (p, z, x, digs);
93 (void) one_minus_mp (p, z, z, digs);
94 return z;
95 }
96
97 //! @brief PROC (LONG REAL) LONG REAL inverf
98
99 MP_T *inverf_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
100 {
101 ADDR_T pop_sp = A68G_SP;
102 int gdigs;
103 BOOL_T siga = A68G_TRUE;
104 // Precision adapts to the argument, but not too much.
105 // If this is not precise enough, you need more digs
106 // in your entire calculation, not just in this routine.
107 // Calculate an initial Newton-Raphson estimate while at it.
108 #if (A68G_LEVEL >= 3)
109 DOUBLE_T y = ABS (mp_to_double (p, x, digs));
110 if (y < erf_double (5.0q)) {
111 y = inverf_double (y);
112 gdigs = FUN_DIGITS (digs);
113 } else {
114 y = 5.0q;
115 gdigs = FUN_DIGITS (2 * digs);
116 }
117 MP_T *z_g = nil_mp (p, gdigs);
118 (void) double_to_mp (p, z_g, y, A68G_DOUBLE_DIG, A68G_FALSE, gdigs);
119 #else
120 REAL_T y = ABS (mp_to_real (p, x, digs));
121 if (y < erf (4.0)) {
122 y = a68g_inverf_real (y);
123 gdigs = FUN_DIGITS (digs);
124 } else {
125 y = 4.0;
126 gdigs = FUN_DIGITS (2 * digs);
127 }
128 MP_T *z_g = nil_mp (p, gdigs);
129 (void) real_to_mp (p, z_g, y, gdigs);
130 #endif
131 MP_T *x_g = len_mp (p, x, digs, gdigs);
132 MP_T *y_g = nil_mp (p, gdigs);
133 int sign = MP_SIGN (x);
134 (void) abs_mp (p, x_g, x_g, gdigs);
135 SET_MP_ONE (y_g, gdigs);
136 (void) sub_mp (p, y_g, x_g, y_g, gdigs);
137 if (MP_SIGN (y_g) >= 0) {
138 errno = EDOM;
139 A68G_SP = pop_sp;
140 return NaN_MP;
141 } else {
142 // Newton-Raphson.
143 MP_T *d_g = nil_mp (p, gdigs);
144 MP_T *f_g = nil_mp (p, gdigs);
145 MP_T *p_g = nil_mp (p, gdigs);
146 // sqrt (pi) / 2
147 (void) mp_pi (p, p_g, MP_SQRT_PI, gdigs);
148 (void) half_mp (p, p_g, p_g, gdigs);
149 // nrdigs prevents end-less iteration
150 int nrdigs;
151 for (nrdigs = 0; nrdigs < digs && siga; nrdigs++) {
152 (void) move_mp (y_g, z_g, gdigs);
153 (void) erf_mp (p, f_g, z_g, gdigs);
154 (void) sub_mp (p, f_g, f_g, x_g, gdigs);
155 (void) mul_mp (p, d_g, z_g, z_g, gdigs);
156 (void) minus_mp (p, d_g, d_g, gdigs);
157 (void) exp_mp (p, d_g, d_g, gdigs);
158 (void) div_mp (p, f_g, f_g, d_g, gdigs);
159 (void) mul_mp (p, f_g, f_g, p_g, gdigs);
160 (void) sub_mp (p, z_g, z_g, f_g, gdigs);
161 (void) sub_mp (p, y_g, z_g, y_g, gdigs);
162 if (IS_ZERO_MP (y_g)) {
163 siga = A68G_FALSE;
164 } else {
165 siga = ABS (MP_EXPONENT (y_g)) < digs;
166 }
167 }
168 (void) shorten_mp (p, z, digs, z_g, gdigs);
169 }
170 if (sign < 0) {
171 (void) minus_mp (p, z, z, digs);
172 }
173 A68G_SP = pop_sp;
174 return z;
175 }
176
177 //! @brief PROC (LONG REAL) LONG REAL inverfc
178
179 MP_T *inverfc_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
180 {
181 (void) one_minus_mp (p, z, x, digs);
182 (void) inverf_mp (p, z, z, digs);
183 return z;
184 }
185
186 // Reference:
187 // John L. Spouge. "Computation of the Gamma, Digamma, and Trigamma Functions".
188 // SIAM Journal on Numerical Analysis. 31 (3) [1994]
189 // Spouge's algorithm sums terms of greatly varying magnitude.
190
191 #define GAMMA_PRECISION(z) (2 * (z))
192
193 //! brief Set up gamma coefficient table
194
195 void mp_gamma_table (NODE_T *p, int digs)
196 {
197 if (A68G_MP (mp_gamma_size) <= 0) {
198 int b = 1;
199 int gdigs = GAMMA_PRECISION (digs);
200 REAL_T log_lim = -digs * LOG_MP_RADIX, log_error;
201 do {
202 ABEND (b >= MP_RADIX, ERROR_HIGH_PRECISION, NO_TEXT);
203 // error = 1 / (sqrt (b) * a68g_x_up_y (2 * M_PI, b + 0.5));
204 log_error = -(log10 (b) / 2 + (b + 0.5) * log10 (2 *M_PI));
205 b += 1;
206 } while (log_error > log_lim);
207 A68G_MP (mp_gamma_size) = b;
208 A68G_MP (mp_gam_ck) = (MP_T **) get_heap_space ((b + 1) * sizeof (MP_T *));
209 A68G_MP (mp_gam_ck)[0] = (MP_T *) get_heap_space (SIZE_MP (gdigs));
210 (void) mp_pi (p, (A68G_MP (mp_gam_ck)[0]), MP_SQRT_TWO_PI, gdigs);
211 ADDR_T pop_sp = A68G_SP;
212 MP_T *ak = nil_mp (p, gdigs);
213 MP_T *db = lit_mp (p, b, 0, gdigs);
214 MP_T *ck = nil_mp (p, gdigs);
215 MP_T *dk = nil_mp (p, gdigs);
216 MP_T *dz = nil_mp (p, gdigs);
217 MP_T *hlf = nil_mp (p, gdigs);
218 MP_T *fac = lit_mp (p, 1, 0, gdigs);
219 SET_MP_HALF (hlf, gdigs);
220 for (int k = 1; k < b; k++) {
221 set_mp (dk, k, 0, gdigs);
222 (void) sub_mp (p, ak, db, dk, gdigs);
223 (void) sub_mp (p, dz, dk, hlf, gdigs);
224 (void) pow_mp (p, ck, ak, dz, gdigs);
225 (void) exp_mp (p, dz, ak, gdigs);
226 (void) mul_mp (p, ck, ck, dz, gdigs);
227 (void) div_mp (p, ck, ck, fac, gdigs);
228 A68G_MP (mp_gam_ck)[k] = (MP_T *) get_heap_space (SIZE_MP (gdigs));
229 (void) move_mp ((A68G_MP (mp_gam_ck)[k]), ck, gdigs);
230 (void) mul_mp (p, fac, fac, dk, gdigs);
231 (void) minus_mp (p, fac, fac, gdigs);
232 errno = 0;
233 }
234 A68G_SP = pop_sp;
235 }
236 }
237
238 MP_T *mp_spouge_sum (NODE_T *p, MP_T *sum, MP_T *x_g, int gdigs)
239 {
240 ADDR_T pop_sp = A68G_SP;
241 int a = A68G_MP (mp_gamma_size);
242 MP_T *da = nil_mp (p, gdigs);
243 MP_T *dz = nil_mp (p, gdigs);
244 (void) move_mp (sum, A68G_MP (mp_gam_ck)[0], gdigs);
245 // Sum small to large to preserve precision.
246 for (int k = a - 1; k > 0; k--) {
247 set_mp (da, k, 0, gdigs);
248 (void) add_mp (p, dz, x_g, da, gdigs);
249 (void) div_mp (p, dz, A68G_MP (mp_gam_ck)[k], dz, gdigs);
250 (void) add_mp (p, sum, sum, dz, gdigs);
251 }
252 A68G_SP = pop_sp;
253 return sum;
254 }
255
256 //! @brief PROC (LONG REAL) LONG REAL gamma
257
258 MP_T *gamma_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
259 {
260 mp_gamma_table (p, digs);
261 int gdigs = GAMMA_PRECISION (digs);
262 // Set up coefficient table.
263 ADDR_T pop_sp = A68G_SP;
264 if (MP_DIGIT (x, 1) < 0) {
265 // G(1-x)G(x) = pi / sin (pi*x)
266 MP_T *pi = nil_mp (p, digs);
267 MP_T *sz = nil_mp (p, digs);
268 MP_T *xm = nil_mp (p, digs);
269 (void) mp_pi (p, pi, MP_PI, digs);
270 (void) one_minus_mp (p, xm, x, digs);
271 (void) gamma_mp (p, xm, xm, digs);
272 (void) sinpi_mp (p, sz, x, digs);
273 (void) mul_mp (p, sz, sz, xm, digs);
274 (void) div_mp (p, z, pi, sz, digs);
275 A68G_SP = pop_sp;
276 return z;
277 }
278 int a = A68G_MP (mp_gamma_size);
279 // Compute Spouge's Gamma
280 MP_T *sum = nil_mp (p, gdigs);
281 MP_T *x_g = len_mp (p, x, digs, gdigs);
282 (void) minus_one_mp (p, x_g, x_g, gdigs);
283 (void) mp_spouge_sum (p, sum, x_g, gdigs);
284 // (z+a)^(z+0.5)*exp(-(z+a)) * Sum
285 MP_T *fac = nil_mp (p, gdigs);
286 MP_T *dz = nil_mp (p, gdigs);
287 MP_T *az = nil_mp (p, gdigs);
288 MP_T *da = nil_mp (p, gdigs);
289 MP_T *hlf = nil_mp (p, gdigs);
290 SET_MP_HALF (hlf, gdigs);
291 set_mp (da, a, 0, gdigs);
292 (void) add_mp (p, az, x_g, da, gdigs);
293 (void) add_mp (p, dz, x_g, hlf, gdigs);
294 (void) pow_mp (p, fac, az, dz, gdigs);
295 (void) minus_mp (p, az, az, gdigs);
296 (void) exp_mp (p, dz, az, gdigs);
297 (void) mul_mp (p, fac, fac, dz, gdigs);
298 (void) mul_mp (p, fac, sum, fac, gdigs);
299 (void) shorten_mp (p, z, digs, fac, gdigs);
300 A68G_SP = pop_sp;
301 return z;
302 }
303
304 //! @brief PROC (LONG REAL) LONG REAL ln gamma
305
306 MP_T *lngamma_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
307 {
308 mp_gamma_table (p, digs);
309 errno = 0;
310 int gdigs = GAMMA_PRECISION (digs);
311 // Set up coefficient table.
312 ADDR_T pop_sp = A68G_SP;
313 if (MP_DIGIT (x, 1) < 0) {
314 // G(1-x)G(x) = pi / sin (pi*x)
315 MP_T *sz = nil_mp (p, digs);
316 MP_T *dz = nil_mp (p, digs);
317 MP_T *xm = nil_mp (p, digs);
318 (void) mp_pi (p, dz, MP_LN_PI, digs);
319 (void) sinpi_mp (p, sz, x, digs);
320 (void) ln_mp (p, sz, sz, digs);
321 (void) sub_mp (p, dz, dz, sz, digs);
322 (void) one_minus_mp (p, xm, x, digs);
323 (void) lngamma_mp (p, xm, xm, digs);
324 (void) sub_mp (p, z, dz, xm, digs);
325 A68G_SP = pop_sp;
326 return z;
327 }
328 int a = A68G_MP (mp_gamma_size);
329 // Compute Spouge's ln Gamma
330 MP_T *sum = nil_mp (p, gdigs);
331 MP_T *x_g = len_mp (p, x, digs, gdigs);
332 (void) minus_one_mp (p, x_g, x_g, gdigs);
333 (void) mp_spouge_sum (p, sum, x_g, gdigs);
334 // (x+0.5) * ln (x+a) - (x+a) + ln Sum
335 MP_T *da = nil_mp (p, gdigs);
336 MP_T *hlf = nil_mp (p, gdigs);
337 SET_MP_HALF (hlf, gdigs);
338 MP_T *fac = nil_mp (p, gdigs);
339 MP_T *dz = nil_mp (p, gdigs);
340 MP_T *az = nil_mp (p, gdigs);
341 set_mp (da, a, 0, gdigs);
342 (void) add_mp (p, az, x_g, da, gdigs);
343 (void) ln_mp (p, fac, az, gdigs);
344 (void) add_mp (p, dz, x_g, hlf, gdigs);
345 (void) mul_mp (p, fac, fac, dz, gdigs);
346 (void) sub_mp (p, fac, fac, az, gdigs);
347 (void) ln_mp (p, dz, sum, gdigs);
348 (void) add_mp (p, fac, fac, dz, gdigs);
349 (void) shorten_mp (p, z, digs, fac, gdigs);
350 A68G_SP = pop_sp;
351 return z;
352 }
353
354 //! @brief PROC (LONG REAL, LONG REAL) LONG REAL ln beta
355
356 MP_T *lnbeta_mp (NODE_T * p, MP_T * z, MP_T * a, MP_T *b, int digs)
357 {
358 ADDR_T pop_sp = A68G_SP;
359 MP_T *aa = nil_mp (p, digs);
360 MP_T *bb = nil_mp (p, digs);
361 MP_T *ab = nil_mp (p, digs);
362 (void) lngamma_mp (p, aa, a, digs);
363 (void) lngamma_mp (p, bb, b, digs);
364 (void) add_mp (p, ab, a, b, digs);
365 (void) lngamma_mp (p, ab, ab, digs);
366 (void) add_mp (p, z, aa, bb, digs);
367 (void) sub_mp (p, z, z, ab, digs);
368 A68G_SP = pop_sp;
369 return z;
370 }
371
372 //! @brief PROC (LONG REAL, LONG REAL) LONG REAL beta
373
374 MP_T *beta_mp (NODE_T * p, MP_T * z, MP_T * a, MP_T *b, int digs)
375 {
376 ADDR_T pop_sp = A68G_SP;
377 MP_T *u = nil_mp (p, digs);
378 lnbeta_mp (p, u, a, b, digs);
379 exp_mp (p, z, u, digs);
380 A68G_SP = pop_sp;
381 return z;
382 }
383
384 //! @brief PROC (LONG REAL, LONG REAL, LONG REAL) LONG REAL beta inc
385
386 MP_T *beta_inc_mp (NODE_T * p, MP_T * z, MP_T * s, MP_T *t, MP_T *x, int digs)
387 {
388 // Incomplete beta function I{x}(s, t).
389 // Continued fraction, see dlmf.nist.gov/8.17; Lentz's algorithm.
390 ADDR_T pop_sp = A68G_SP;
391 A68G_BOOL gt;
392 MP_T *one = lit_mp (p, 1, 0, digs);
393 gt_mp (p, >, x, one, digs);
394 if (MP_DIGIT (x, 1) < 0 || VALUE (>)) {
395 errno = EDOM;
396 A68G_SP = pop_sp;
397 return NaN_MP;
398 }
399 if (same_mp (p, x, one, digs)) {
400 SET_MP_ONE (z, digs);
401 A68G_SP = pop_sp;
402 return z;
403 }
404 // Rapid convergence when x <= (s+1)/((s+1)+(t+1)) or else recursion.
405 {
406 MP_T *u = nil_mp (p, digs), *v = nil_mp (p, digs), *w = nil_mp (p, digs);
407 (void) plus_one_mp (p, u, s, digs);
408 (void) plus_one_mp (p, v, t, digs);
409 (void) add_mp (p, w, u, v, digs);
410 (void) div_mp (p, u, u, w, digs);
411 gt_mp (p, >, x, u, digs);
412 if (VALUE (>)) {
413 // B{x}(s, t) = 1 - B{1-x}(t, s)
414 (void) one_minus_mp (p, x, x, digs);
415 PRELUDE_ERROR (beta_inc_mp (p, z, s, t, x, digs) == NaN_MP, p, ERROR_INVALID_ARGUMENT, MOID (p));
416 (void) one_minus_mp (p, z, z, digs);
417 A68G_SP = pop_sp;
418 return z;
419 }
420 }
421 // Lentz's algorithm for continued fraction.
422 A68G_SP = pop_sp;
423 int gdigs = FUN_DIGITS (digs);
424 const INT_T lim = gdigs * LOG_MP_RADIX;
425 BOOL_T cont = A68G_TRUE;
426 MP_T *F = lit_mp (p, 1, 0, gdigs);
427 MP_T *T = lit_mp (p, 1, 0, gdigs);
428 MP_T *W = lit_mp (p, 1, 0, gdigs);
429 MP_T *c = lit_mp (p, 1, 0, gdigs);
430 MP_T *d = nil_mp (p, gdigs);
431 MP_T *m = nil_mp (p, gdigs);
432 MP_T *s_g = len_mp (p, s, digs, gdigs);
433 MP_T *t_g = len_mp (p, t, digs, gdigs);
434 MP_T *x_g = len_mp (p, x, digs, gdigs);
435 MP_T *u = lit_mp (p, 1, 0, gdigs);
436 MP_T *v = nil_mp (p, gdigs);
437 MP_T *w = nil_mp (p, gdigs);
438 for (int N = 0; cont && N < lim; N++) {
439 if (N == 0) {
440 SET_MP_ONE (T, gdigs);
441 } else if (N % 2 == 0) {
442 // d{2m} := x m(t-m)/((s+2m-1)(s+2m))
443 // T = x * m * (t - m) / (s + 2.0q * m - 1.0q) / (s + 2.0q * m);
444 (void) sub_mp (p, u, t_g, m, gdigs);
445 (void) mul_mp (p, u, u, m, gdigs);
446 (void) mul_mp (p, u, u, x_g, gdigs);
447 (void) add_mp (p, v, m, m, gdigs);
448 (void) add_mp (p, v, s_g, v, gdigs);
449 (void) minus_one_mp (p, v, v, gdigs);
450 (void) add_mp (p, w, m, m, gdigs);
451 (void) add_mp (p, w, s_g, w, gdigs);
452 (void) div_mp (p, T, u, v, gdigs);
453 (void) div_mp (p, T, T, w, gdigs);
454 } else {
455 // d{2m+1} := -x (s+m)(s+t+m)/((s+2m+1)(s+2m))
456 // T = -x * (s + m) * (s + t + m) / (s + 2.0q * m + 1.0q) / (s + 2.0q * m);
457 (void) add_mp (p, u, s_g, m, gdigs);
458 (void) add_mp (p, v, u, t_g, gdigs);
459 (void) mul_mp (p, u, u, v, gdigs);
460 (void) mul_mp (p, u, u, x_g, gdigs);
461 (void) minus_mp (p, u, u, gdigs);
462 (void) add_mp (p, v, m, m, gdigs);
463 (void) add_mp (p, v, s_g, v, gdigs);
464 (void) plus_one_mp (p, v, v, gdigs);
465 (void) add_mp (p, w, m, m, gdigs);
466 (void) add_mp (p, w, s_g, w, gdigs);
467 (void) div_mp (p, T, u, v, gdigs);
468 (void) div_mp (p, T, T, w, gdigs);
469 (void) plus_one_mp (p, m, m, gdigs);
470 }
471 // d = 1.0q / (T * d + 1.0q);
472 (void) mul_mp (p, d, T, d, gdigs);
473 (void) plus_one_mp (p, d, d, gdigs);
474 (void) rec_mp (p, d, d, gdigs);
475 // c = T / c + 1.0q;
476 (void) div_mp (p, c, T, c, gdigs);
477 (void) plus_one_mp (p, c, c, gdigs);
478 // F *= c * d;
479 (void) mul_mp (p, F, F, c, gdigs);
480 (void) mul_mp (p, F, F, d, gdigs);
481 if (same_mp (p, F, W, gdigs)) {
482 cont = A68G_FALSE;
483 } else {
484 (void) move_mp (W, F, gdigs);
485 }
486 }
487 minus_one_mp (p, F, F, gdigs);
488 // I{x}(s,t)=x^s(1-x)^t / s / B(s,t) F
489 (void) pow_mp (p, u, x_g, s_g, gdigs);
490 (void) one_minus_mp (p, v, x_g, gdigs);
491 (void) pow_mp (p, v, v, t_g, gdigs);
492 (void) beta_mp (p, w, s_g, t_g, gdigs);
493 (void) mul_mp (p, m, u, v, gdigs);
494 (void) div_mp (p, m, m, s_g, gdigs);
495 (void) div_mp (p, m, m, w, gdigs);
496 (void) mul_mp (p, m, m, F, gdigs);
497 (void) shorten_mp (p, z, digs, m, gdigs);
498 A68G_SP = pop_sp;
499 return z;
500 }
© 2001-2026 J.M. van der Veer
jmvdveer@algol68genie.nl
