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-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 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), 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 (unt 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_double_real (p, x, digs));
111 if (y < erfq (5.0q)) {
112 y = inverf_double_real (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) double_real_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 for (unt k = 1; k < b; k++) {
223 set_mp (dk, k, 0, gdigs);
224 (void) sub_mp (p, ak, db, dk, gdigs);
225 (void) sub_mp (p, dz, dk, hlf, gdigs);
226 (void) pow_mp (p, ck, ak, dz, gdigs);
227 (void) exp_mp (p, dz, ak, gdigs);
228 (void) mul_mp (p, ck, ck, dz, gdigs);
229 (void) div_mp (p, ck, ck, fac, gdigs);
230 A68_MP (mp_gam_ck)[k] = (MP_T *) get_heap_space (SIZE_MP (gdigs));
231 (void) move_mp ((A68_MP(mp_gam_ck)[k]), ck, gdigs);
232 (void) mul_mp (p, fac, fac, dk, gdigs);
233 (void) minus_mp (p, fac, fac, gdigs);
234 }
235 A68_SP = pop_sp;
236 }
237 }
238
239 MP_T *mp_spouge_sum (NODE_T *p, MP_T *sum, MP_T *x_g, int gdigs)
240 {
241 ADDR_T pop_sp = A68_SP;
242 int a = A68_MP (mp_gamma_size);
243 MP_T *da = nil_mp (p, gdigs);
244 MP_T *dz = nil_mp (p, gdigs);
245 (void) move_mp (sum, A68_MP (mp_gam_ck)[0], gdigs);
246 // Sum small to large to preserve precision.
247 for (int k = a - 1; k > 0; k--) {
248 set_mp (da, k, 0, gdigs);
249 (void) add_mp (p, dz, x_g, da, gdigs);
250 (void) div_mp (p, dz, A68_MP (mp_gam_ck)[k], dz, gdigs);
251 (void) add_mp (p, sum, sum, dz, gdigs);
252 }
253 A68_SP = pop_sp;
254 return sum;
255 }
256
257 //! @brief PROC (LONG REAL) LONG REAL gamma
258
259 MP_T *gamma_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
260 {
261 mp_gamma_table (p, digs);
262 int gdigs = GAMMA_PRECISION (digs);
263 // Set up coefficient table.
264 ADDR_T pop_sp = A68_SP;
265 if (MP_DIGIT (x, 1) < 0) {
266 // G(1-x)G(x) = pi / sin (pi*x)
267 MP_T *pi = nil_mp (p, digs);
268 MP_T *sz = nil_mp (p, digs);
269 MP_T *xm = nil_mp (p, digs);
270 (void) mp_pi (p, pi, MP_PI, digs);
271 (void) one_minus_mp (p, xm, x, digs);
272 (void) gamma_mp (p, xm, xm, digs);
273 (void) sinpi_mp (p, sz, x, digs);
274 (void) mul_mp (p, sz, sz, xm, digs);
275 (void) div_mp (p, z, pi, sz, digs);
276 A68_SP = pop_sp;
277 return z;
278 }
279 int a = A68_MP (mp_gamma_size);
280 // Compute Spouge's Gamma
281 MP_T *sum = nil_mp (p, gdigs);
282 MP_T *x_g = len_mp (p, x, digs, gdigs);
283 (void) minus_one_mp (p, x_g, x_g, gdigs);
284 (void) mp_spouge_sum (p, sum, x_g, gdigs);
285 // (z+a)^(z+0.5)*exp(-(z+a)) * Sum
286 MP_T *fac = nil_mp (p, gdigs);
287 MP_T *dz = nil_mp (p, gdigs);
288 MP_T *az = nil_mp (p, gdigs);
289 MP_T *da = nil_mp (p, gdigs);
290 MP_T *hlf = nil_mp (p, gdigs);
291 SET_MP_HALF (hlf, gdigs);
292 set_mp (da, a, 0, gdigs);
293 (void) add_mp (p, az, x_g, da, gdigs);
294 (void) add_mp (p, dz, x_g, hlf, gdigs);
295 (void) pow_mp (p, fac, az, dz, gdigs);
296 (void) minus_mp (p, az, az, gdigs);
297 (void) exp_mp (p, dz, az, gdigs);
298 (void) mul_mp (p, fac, fac, dz, gdigs);
299 (void) mul_mp (p, fac, sum, fac, gdigs);
300 (void) shorten_mp (p, z, digs, fac, gdigs);
301 A68_SP = pop_sp;
302 return z;
303 }
304
305 //! @brief PROC (LONG REAL) LONG REAL ln gamma
306
307 MP_T *lngamma_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
308 {
309 mp_gamma_table (p, digs);
310 int gdigs = GAMMA_PRECISION (digs);
311 // Set up coefficient table.
312 ADDR_T pop_sp = A68_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 A68_SP = pop_sp;
326 return z;
327 }
328 int a = A68_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 A68_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 = A68_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 A68_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 = A68_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 A68_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 = A68_SP;
391 A68_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 A68_SP = pop_sp;
397 return NaN_MP;
398 }
399 if (same_mp (p, x, one, digs)) {
400 SET_MP_ONE (z, digs);
401 A68_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 A68_SP = pop_sp;
418 return z;
419 }
420 }
421 // Lentz's algorithm for continued fraction.
422 A68_SP = pop_sp;
423 int gdigs = FUN_DIGITS (digs);
424 const INT_T lim = gdigs * LOG_MP_RADIX;
425 BOOL_T cont = A68_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 (unt 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 = A68_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 A68_SP = pop_sp;
499 return z;
500 }
