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
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 = A68_SP;
41 // Note we need double precision!
42 BOOL_T siga = A68_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 A68_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 = A68_SP;
102 int gdigs;
103 BOOL_T siga = A68_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 (A68_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, gdigs);
119 #else
120 REAL_T y = ABS (mp_to_real (p, x, digs));
121 if (y < erf (4.0)) {
122 y = a68_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 A68_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 = A68_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 A68_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 (A68_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, __func__);
203 // error = 1 / (sqrt (b) * a68_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 A68_MP (mp_gamma_size) = b;
208 A68_MP (mp_gam_ck) = (MP_T **) get_heap_space ((size_t) ((b + 1) * sizeof (MP_T *)));
209 A68_MP (mp_gam_ck)[0] = (MP_T *) get_heap_space (SIZE_MP (gdigs));
210 (void) mp_pi (p, (A68_MP (mp_gam_ck)[0]), MP_SQRT_TWO_PI, gdigs);
211 ADDR_T pop_sp = A68_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 A68_MP (mp_gam_ck)[k] = (MP_T *) get_heap_space (SIZE_MP (gdigs));
229 (void) move_mp ((A68_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 }
233 A68_SP = pop_sp;
234 }
235 }
236
237 MP_T *mp_spouge_sum (NODE_T *p, MP_T *sum, MP_T *x_g, int gdigs)
238 {
239 ADDR_T pop_sp = A68_SP;
240 int a = A68_MP (mp_gamma_size);
241 MP_T *da = nil_mp (p, gdigs);
242 MP_T *dz = nil_mp (p, gdigs);
243 (void) move_mp (sum, A68_MP (mp_gam_ck)[0], gdigs);
244 // Sum small to large to preserve precision.
245 for (int k = a - 1; k > 0; k--) {
246 set_mp (da, k, 0, gdigs);
247 (void) add_mp (p, dz, x_g, da, gdigs);
248 (void) div_mp (p, dz, A68_MP (mp_gam_ck)[k], dz, gdigs);
249 (void) add_mp (p, sum, sum, dz, gdigs);
250 }
251 A68_SP = pop_sp;
252 return sum;
253 }
254
255 //! @brief PROC (LONG REAL) LONG REAL gamma
256
257 MP_T *gamma_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
258 {
259 mp_gamma_table (p, digs);
260 int gdigs = GAMMA_PRECISION (digs);
261 // Set up coefficient table.
262 ADDR_T pop_sp = A68_SP;
263 if (MP_DIGIT (x, 1) < 0) {
264 // G(1-x)G(x) = pi / sin (pi*x)
265 MP_T *pi = nil_mp (p, digs);
266 MP_T *sz = nil_mp (p, digs);
267 MP_T *xm = nil_mp (p, digs);
268 (void) mp_pi (p, pi, MP_PI, digs);
269 (void) one_minus_mp (p, xm, x, digs);
270 (void) gamma_mp (p, xm, xm, digs);
271 (void) sinpi_mp (p, sz, x, digs);
272 (void) mul_mp (p, sz, sz, xm, digs);
273 (void) div_mp (p, z, pi, sz, digs);
274 A68_SP = pop_sp;
275 return z;
276 }
277 int a = A68_MP (mp_gamma_size);
278 // Compute Spouge's Gamma
279 MP_T *sum = nil_mp (p, gdigs);
280 MP_T *x_g = len_mp (p, x, digs, gdigs);
281 (void) minus_one_mp (p, x_g, x_g, gdigs);
282 (void) mp_spouge_sum (p, sum, x_g, gdigs);
283 // (z+a)^(z+0.5)*exp(-(z+a)) * Sum
284 MP_T *fac = nil_mp (p, gdigs);
285 MP_T *dz = nil_mp (p, gdigs);
286 MP_T *az = nil_mp (p, gdigs);
287 MP_T *da = nil_mp (p, gdigs);
288 MP_T *hlf = nil_mp (p, gdigs);
289 SET_MP_HALF (hlf, gdigs);
290 set_mp (da, a, 0, gdigs);
291 (void) add_mp (p, az, x_g, da, gdigs);
292 (void) add_mp (p, dz, x_g, hlf, gdigs);
293 (void) pow_mp (p, fac, az, dz, gdigs);
294 (void) minus_mp (p, az, az, gdigs);
295 (void) exp_mp (p, dz, az, gdigs);
296 (void) mul_mp (p, fac, fac, dz, gdigs);
297 (void) mul_mp (p, fac, sum, fac, gdigs);
298 (void) shorten_mp (p, z, digs, fac, gdigs);
299 A68_SP = pop_sp;
300 return z;
301 }
302
303 //! @brief PROC (LONG REAL) LONG REAL ln gamma
304
305 MP_T *lngamma_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
306 {
307 mp_gamma_table (p, digs);
308 int gdigs = GAMMA_PRECISION (digs);
309 // Set up coefficient table.
310 ADDR_T pop_sp = A68_SP;
311 if (MP_DIGIT (x, 1) < 0) {
312 // G(1-x)G(x) = pi / sin (pi*x)
313 MP_T *sz = nil_mp (p, digs);
314 MP_T *dz = nil_mp (p, digs);
315 MP_T *xm = nil_mp (p, digs);
316 (void) mp_pi (p, dz, MP_LN_PI, digs);
317 (void) sinpi_mp (p, sz, x, digs);
318 (void) ln_mp (p, sz, sz, digs);
319 (void) sub_mp (p, dz, dz, sz, digs);
320 (void) one_minus_mp (p, xm, x, digs);
321 (void) lngamma_mp (p, xm, xm, digs);
322 (void) sub_mp (p, z, dz, xm, digs);
323 A68_SP = pop_sp;
324 return z;
325 }
326 int a = A68_MP (mp_gamma_size);
327 // Compute Spouge's ln Gamma
328 MP_T *sum = nil_mp (p, gdigs);
329 MP_T *x_g = len_mp (p, x, digs, gdigs);
330 (void) minus_one_mp (p, x_g, x_g, gdigs);
331 (void) mp_spouge_sum (p, sum, x_g, gdigs);
332 // (x+0.5) * ln (x+a) - (x+a) + ln Sum
333 MP_T *da = nil_mp (p, gdigs);
334 MP_T *hlf = nil_mp (p, gdigs);
335 SET_MP_HALF (hlf, gdigs);
336 MP_T *fac = nil_mp (p, gdigs);
337 MP_T *dz = nil_mp (p, gdigs);
338 MP_T *az = nil_mp (p, gdigs);
339 set_mp (da, a, 0, gdigs);
340 (void) add_mp (p, az, x_g, da, gdigs);
341 (void) ln_mp (p, fac, az, gdigs);
342 (void) add_mp (p, dz, x_g, hlf, gdigs);
343 (void) mul_mp (p, fac, fac, dz, gdigs);
344 (void) sub_mp (p, fac, fac, az, gdigs);
345 (void) ln_mp (p, dz, sum, gdigs);
346 (void) add_mp (p, fac, fac, dz, gdigs);
347 (void) shorten_mp (p, z, digs, fac, gdigs);
348 A68_SP = pop_sp;
349 return z;
350 }
351
352 //! @brief PROC (LONG REAL, LONG REAL) LONG REAL ln beta
353
354 MP_T *lnbeta_mp (NODE_T * p, MP_T * z, MP_T * a, MP_T *b, int digs)
355 {
356 ADDR_T pop_sp = A68_SP;
357 MP_T *aa = nil_mp (p, digs);
358 MP_T *bb = nil_mp (p, digs);
359 MP_T *ab = nil_mp (p, digs);
360 (void) lngamma_mp (p, aa, a, digs);
361 (void) lngamma_mp (p, bb, b, digs);
362 (void) add_mp (p, ab, a, b, digs);
363 (void) lngamma_mp (p, ab, ab, digs);
364 (void) add_mp (p, z, aa, bb, digs);
365 (void) sub_mp (p, z, z, ab, digs);
366 A68_SP = pop_sp;
367 return z;
368 }
369
370 //! @brief PROC (LONG REAL, LONG REAL) LONG REAL beta
371
372 MP_T *beta_mp (NODE_T * p, MP_T * z, MP_T * a, MP_T *b, int digs)
373 {
374 ADDR_T pop_sp = A68_SP;
375 MP_T *u = nil_mp (p, digs);
376 lnbeta_mp (p, u, a, b, digs);
377 exp_mp (p, z, u, digs);
378 A68_SP = pop_sp;
379 return z;
380 }
381
382 //! @brief PROC (LONG REAL, LONG REAL, LONG REAL) LONG REAL beta inc
383
384 MP_T *beta_inc_mp (NODE_T * p, MP_T * z, MP_T * s, MP_T *t, MP_T *x, int digs)
385 {
386 // Incomplete beta function I{x}(s, t).
387 // Continued fraction, see dlmf.nist.gov/8.17; Lentz's algorithm.
388 ADDR_T pop_sp = A68_SP;
389 A68_BOOL gt;
390 MP_T *one = lit_mp (p, 1, 0, digs);
391 gt_mp (p, >, x, one, digs);
392 if (MP_DIGIT (x, 1) < 0 || VALUE (>)) {
393 errno = EDOM;
394 A68_SP = pop_sp;
395 return NaN_MP;
396 }
397 if (same_mp (p, x, one, digs)) {
398 SET_MP_ONE (z, digs);
399 A68_SP = pop_sp;
400 return z;
401 }
402 // Rapid convergence when x <= (s+1)/((s+1)+(t+1)) or else recursion.
403 {
404 MP_T *u = nil_mp (p, digs), *v = nil_mp (p, digs), *w = nil_mp (p, digs);
405 (void) plus_one_mp (p, u, s, digs);
406 (void) plus_one_mp (p, v, t, digs);
407 (void) add_mp (p, w, u, v, digs);
408 (void) div_mp (p, u, u, w, digs);
409 gt_mp (p, >, x, u, digs);
410 if (VALUE (>)) {
411 // B{x}(s, t) = 1 - B{1-x}(t, s)
412 (void) one_minus_mp (p, x, x, digs);
413 PRELUDE_ERROR (beta_inc_mp (p, z, s, t, x, digs) == NaN_MP, p, ERROR_INVALID_ARGUMENT, MOID (p));
414 (void) one_minus_mp (p, z, z, digs);
415 A68_SP = pop_sp;
416 return z;
417 }
418 }
419 // Lentz's algorithm for continued fraction.
420 A68_SP = pop_sp;
421 int gdigs = FUN_DIGITS (digs);
422 const INT_T lim = gdigs * LOG_MP_RADIX;
423 BOOL_T cont = A68_TRUE;
424 MP_T *F = lit_mp (p, 1, 0, gdigs);
425 MP_T *T = lit_mp (p, 1, 0, gdigs);
426 MP_T *W = lit_mp (p, 1, 0, gdigs);
427 MP_T *c = lit_mp (p, 1, 0, gdigs);
428 MP_T *d = nil_mp (p, gdigs);
429 MP_T *m = nil_mp (p, gdigs);
430 MP_T *s_g = len_mp (p, s, digs, gdigs);
431 MP_T *t_g = len_mp (p, t, digs, gdigs);
432 MP_T *x_g = len_mp (p, x, digs, gdigs);
433 MP_T *u = lit_mp (p, 1, 0, gdigs);
434 MP_T *v = nil_mp (p, gdigs);
435 MP_T *w = nil_mp (p, gdigs);
436 for (int N = 0; cont && N < lim; N++) {
437 if (N == 0) {
438 SET_MP_ONE (T, gdigs);
439 } else if (N % 2 == 0) {
440 // d{2m} := x m(t-m)/((s+2m-1)(s+2m))
441 // T = x * m * (t - m) / (s + 2.0q * m - 1.0q) / (s + 2.0q * m);
442 (void) sub_mp (p, u, t_g, m, gdigs);
443 (void) mul_mp (p, u, u, m, gdigs);
444 (void) mul_mp (p, u, u, x_g, gdigs);
445 (void) add_mp (p, v, m, m, gdigs);
446 (void) add_mp (p, v, s_g, v, gdigs);
447 (void) minus_one_mp (p, v, v, gdigs);
448 (void) add_mp (p, w, m, m, gdigs);
449 (void) add_mp (p, w, s_g, w, gdigs);
450 (void) div_mp (p, T, u, v, gdigs);
451 (void) div_mp (p, T, T, w, gdigs);
452 } else {
453 // d{2m+1} := -x (s+m)(s+t+m)/((s+2m+1)(s+2m))
454 // T = -x * (s + m) * (s + t + m) / (s + 2.0q * m + 1.0q) / (s + 2.0q * m);
455 (void) add_mp (p, u, s_g, m, gdigs);
456 (void) add_mp (p, v, u, t_g, gdigs);
457 (void) mul_mp (p, u, u, v, gdigs);
458 (void) mul_mp (p, u, u, x_g, gdigs);
459 (void) minus_mp (p, u, u, gdigs);
460 (void) add_mp (p, v, m, m, gdigs);
461 (void) add_mp (p, v, s_g, v, gdigs);
462 (void) plus_one_mp (p, v, v, gdigs);
463 (void) add_mp (p, w, m, m, gdigs);
464 (void) add_mp (p, w, s_g, w, gdigs);
465 (void) div_mp (p, T, u, v, gdigs);
466 (void) div_mp (p, T, T, w, gdigs);
467 (void) plus_one_mp (p, m, m, gdigs);
468 }
469 // d = 1.0q / (T * d + 1.0q);
470 (void) mul_mp (p, d, T, d, gdigs);
471 (void) plus_one_mp (p, d, d, gdigs);
472 (void) rec_mp (p, d, d, gdigs);
473 // c = T / c + 1.0q;
474 (void) div_mp (p, c, T, c, gdigs);
475 (void) plus_one_mp (p, c, c, gdigs);
476 // F *= c * d;
477 (void) mul_mp (p, F, F, c, gdigs);
478 (void) mul_mp (p, F, F, d, gdigs);
479 if (same_mp (p, F, W, gdigs)) {
480 cont = A68_FALSE;
481 } else {
482 (void) move_mp (W, F, gdigs);
483 }
484 }
485 minus_one_mp (p, F, F, gdigs);
486 // I{x}(s,t)=x^s(1-x)^t / s / B(s,t) F
487 (void) pow_mp (p, u, x_g, s_g, gdigs);
488 (void) one_minus_mp (p, v, x_g, gdigs);
489 (void) pow_mp (p, v, v, t_g, gdigs);
490 (void) beta_mp (p, w, s_g, t_g, gdigs);
491 (void) mul_mp (p, m, u, v, gdigs);
492 (void) div_mp (p, m, m, s_g, gdigs);
493 (void) div_mp (p, m, m, w, gdigs);
494 (void) mul_mp (p, m, m, F, gdigs);
495 (void) shorten_mp (p, z, digs, m, gdigs);
496 A68_SP = pop_sp;
497 return z;
498 }
