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