mp.c
1 //! @file mp.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 INT, REAL and COMPLEX routines.
25
26 // Multiprecision calculations are useful in these cases:
27 //
28 // Ill-conditioned linear systems
29 // Summation of large series
30 // Long-time or large-scale simulations
31 // Small-scale phenomena
32 // 'Experimental mathematics'
33 //
34 // The routines in this library follow algorithms as described in the
35 // literature, notably
36 //
37 // D.M. Smith, "Efficient Multiple-Precision Evaluation of Elementary Functions"
38 // Mathematics of Computation 52 (1989) 131-134
39 //
40 // D.M. Smith, "A Multiple-Precision Division Algorithm"
41 // Mathematics of Computation 66 (1996) 157-163
42 //
43 // The GNU MPFR library documentation
44 //
45 // Multiprecision libraries are (freely) available, but this one is particularly
46 // designed to work with Algol68G. It implements following modes:
47 //
48 // LONG INT, LONG REAL, LONG COMPLEX, LONG BITS
49 // LONG LONG INT, LONG LONG REAL, LONG LONG COMPLEX, LONG LONG BITS
50 //
51 // Note that recent implementations of GCC make available 64-bit LONG INT and
52 // 128-bit LONG REAL. This suits many multiprecision needs already.
53 // On such platforms, below code is used for LONG LONG modes only.
54 // Now that 64-bit integers are commonplace, this library has been adapted to
55 // exploit them. Having some more digits per word gives performance gain.
56 // This is implemented through macros to keep the library compatible with
57 // old platforms with 32-bit integers and 64-bit doubles.
58 //
59 // Currently, LONG modes have a fixed precision, and LONG LONG modes have
60 // user-definable precision. Precisions span about 30 decimal digits for
61 // LONG modes up to (default) about 60 decimal digits for LONG LONG modes, a
62 // range that is thought to be adequate for most multiprecision applications.
63 //
64 // Although the maximum length of a number is in principle unbound, this
65 // implementation is not designed for more than a few hundred decimal places.
66 // At higher precisions, expect a performance penalty with respect to
67 // state of the art implementations that may for instance use convolution for
68 // multiplication.
69 //
70 // This library takes a sloppy approach towards LONG INT and LONG BITS which are
71 // implemented as LONG REAL and truncated where appropriate. This keeps the code
72 // short at the penalty of some performance loss.
73 //
74 // As is common practice, mp numbers are represented by a row of digits
75 // in a large base. Layout of a mp number "z" is:
76 //
77 // MP_T *z;
78 //
79 // MP_STATUS (z) Status word
80 // MP_EXPONENT (z) Exponent with base MP_RADIX
81 // MP_DIGIT (z, 1 .. N) Digits 1 .. N
82 //
83 // Note that this library assumes IEEE 754 compatible implementation of
84 // type "double". It also assumes a 32- (or 64-) bit type "int".
85 //
86 // Most legacy multiple precision libraries stored numbers as [] int*4.
87 // However, since division and multiplication are O(N ** 2) operations, it is
88 // advantageous to keep the base as high as possible. Modern computers handle
89 // doubles at similar or better speed as integers, therefore this library
90 // opts for storing numbers as [] words were a word is real*8 (legacy) or
91 // int*8 (on f.i. ix86 processors that have real*10), trading space for speed.
92 //
93 // Set a base such that "base^2" can be exactly represented by a word.
94 // To facilitate transput, we require a base that is a power of 10.
95 //
96 // If we choose the base right then in multiplication and division we do not need
97 // to normalise intermediate results at each step since a number of additions
98 // can be made before overflow occurs. That is why we specify "MAX_REPR_INT".
99 //
100 // Mind that the precision of a mp number is at worst just
101 // (LONG_MP_DIGITS - 1) * LOG_MP_RADIX + 1, since the most significant mp digit
102 // is also in range [0 .. MP_RADIX>. Do not specify less than 2 digits.
103 //
104 // Since this software is distributed without any warranty, it is your
105 // responsibility to validate the behaviour of the routines and their accuracy
106 // using the source code provided. See the GNU General Public License for details.
107
108 #include "a68g.h"
109 #include "a68g-genie.h"
110 #include "a68g-prelude.h"
111 #include "a68g-double.h"
112 #include "a68g-mp.h"
113
114 // Internal mp constants.
115
116 //! @brief Set number of digits for long long numbers.
117
118 void set_long_mp_digits (int n)
119 {
120 A68_MP (varying_mp_digits) = n;
121 }
122
123 //! @brief Convert precision to digits for long long number.
124
125 int width_to_mp_digits (int n)
126 {
127 return 1 + n / LOG_MP_RADIX;
128 }
129
130 //! @brief Unformatted write of z to stdout; debugging routine.
131
132 #if !defined (BUILD_WIN32)
133
134 void raw_write_mp (char *str, MP_T * z, int digs)
135 {
136 fprintf (stdout, "\n(%d digits)%s", digs, str);
137 for (unt i = 1; i <= digs; i++) {
138 #if (A68_LEVEL >= 3)
139 fprintf (stdout, " %09lld", (MP_INT_T) MP_DIGIT (z, i));
140 #else
141 fprintf (stdout, " %07d", (MP_INT_T) MP_DIGIT (z, i));
142 #endif
143 }
144 fprintf (stdout, " E" A68_LD, (MP_INT_T) MP_EXPONENT (z));
145 fprintf (stdout, " S" A68_LD, (MP_INT_T) MP_STATUS (z));
146 fprintf (stdout, "\n");
147 ASSERT (fflush (stdout) == 0);
148 }
149
150 #endif
151
152 //! @brief Whether z is a valid representation for its mode.
153
154 BOOL_T check_mp_int (MP_T * z, MOID_T * m)
155 {
156 if (m == M_LONG_INT || m == M_LONG_BITS) {
157 return (BOOL_T) ((MP_EXPONENT (z) >= (MP_T) 0) && (MP_EXPONENT (z) < (MP_T) LONG_MP_DIGITS));
158 } else if (m == M_LONG_LONG_INT || m == M_LONG_LONG_BITS) {
159 return (BOOL_T) ((MP_EXPONENT (z) >= (MP_T) 0) && (MP_EXPONENT (z) < (MP_T) A68_MP (varying_mp_digits)));
160 }
161 return A68_FALSE;
162 }
163
164 //! @brief |x|
165
166 MP_T *abs_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
167 {
168 (void) p;
169 if (x != z) {
170 (void) move_mp (z, x, digs);
171 }
172 MP_DIGIT (z, 1) = ABS (MP_DIGIT (z, 1));
173 MP_STATUS (z) = (MP_T) INIT_MASK;
174 return z;
175 }
176
177 //! @brief -x
178
179 MP_T *minus_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
180 {
181 (void) p;
182 if (x != z) {
183 (void) move_mp (z, x, digs);
184 }
185 MP_DIGIT (z, 1) = -(MP_DIGIT (z, 1));
186 MP_STATUS (z) = (MP_T) INIT_MASK;
187 return z;
188 }
189
190 //! @brief 1 - x
191
192 MP_T *one_minus_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
193 {
194 ADDR_T pop_sp = A68_SP;
195 (void) sub_mp (p, z, mp_one (digs), x, digs);
196 A68_SP = pop_sp;
197 return z;
198 }
199
200 //! @brief x - 1
201
202 MP_T *minus_one_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
203 {
204 ADDR_T pop_sp = A68_SP;
205 (void) sub_mp (p, z, x, mp_one (digs), digs);
206 A68_SP = pop_sp;
207 return z;
208 }
209
210 //! @brief x + 1
211
212 MP_T *plus_one_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
213 {
214 ADDR_T pop_sp = A68_SP;
215 (void) add_mp (p, z, x, mp_one (digs), digs);
216 A68_SP = pop_sp;
217 return z;
218 }
219
220 //! @brief Test whether x = y.
221
222 BOOL_T same_mp (NODE_T * p, MP_T * x, MP_T * y, int digs)
223 {
224 (void) p;
225 if ((MP_STATUS (x) == MP_STATUS (y)) && (MP_EXPONENT (x) == MP_EXPONENT (y))) {
226 for (unt k = digs; k >= 1; k--) {
227 if (MP_DIGIT (x, k) != MP_DIGIT (y, k)) {
228 return A68_FALSE;
229 }
230 }
231 return A68_TRUE;
232 } else {
233 return A68_FALSE;
234 }
235 }
236
237 //! @brief Align 10-base z in a MP_RADIX mantissa.
238
239 MP_T *align_mp (MP_T * z, INT_T * expo, int digs)
240 {
241 INT_T shift;
242 if (*expo >= 0) {
243 shift = LOG_MP_RADIX - (*expo) % LOG_MP_RADIX - 1;
244 (*expo) /= LOG_MP_RADIX;
245 } else {
246 shift = (-(*expo) - 1) % LOG_MP_RADIX;
247 (*expo) = ((*expo) + 1) / LOG_MP_RADIX;
248 (*expo)--;
249 }
250 // Optimising below code does not make the library noticeably faster.
251 for (unt i = 1; i <= shift; i++) {
252 INT_T carry = 0;
253 for (INT_T j = 1; j <= digs; j++) {
254 MP_INT_T k = ((MP_INT_T) MP_DIGIT (z, j)) % 10;
255 MP_DIGIT (z, j) = (MP_T) ((MP_INT_T) (MP_DIGIT (z, j) / 10) + carry * (MP_RADIX / 10));
256 carry = k;
257 }
258 }
259 return z;
260 }
261
262 //! @brief Transform string into multi-precision number.
263
264 MP_T *strtomp (NODE_T * p, MP_T * z, char *s, int digs)
265 {
266 BOOL_T ok = A68_TRUE;
267 errno = 0;
268 SET_MP_ZERO (z, digs);
269 while (IS_SPACE (s[0])) {
270 s++;
271 }
272 // Get the sign.
273 int sign = (s[0] == '-' ? -1 : 1);
274 if (s[0] == '+' || s[0] == '-') {
275 s++;
276 }
277 // Scan mantissa digs and put them into "z".
278 while (s[0] == '0') {
279 s++;
280 }
281 int i = 0, dig = 1;
282 INT_T sum = 0, dot = -1, one = -1, pow = 0, W = MP_RADIX / 10;
283 while (s[i] != NULL_CHAR && dig <= digs && (IS_DIGIT (s[i]) || s[i] == POINT_CHAR)) {
284 if (s[i] == POINT_CHAR) {
285 dot = i;
286 } else {
287 int value = (int) s[i] - (int) '0';
288 if (one < 0 && value > 0) {
289 one = pow;
290 }
291 sum += W * value;
292 if (one >= 0) {
293 W /= 10;
294 }
295 pow++;
296 if (W < 1) {
297 MP_DIGIT (z, dig++) = (MP_T) sum;
298 sum = 0;
299 W = MP_RADIX / 10;
300 }
301 }
302 i++;
303 }
304 // Store the last digs.
305 if (dig <= digs) {
306 MP_DIGIT (z, dig++) = (MP_T) sum;
307 }
308 // See if there is an exponent.
309 INT_T expo;
310 if (s[i] != NULL_CHAR && TO_UPPER (s[i]) == TO_UPPER (EXPONENT_CHAR)) {
311 char *end;
312 expo = (int) strtol (&(s[++i]), &end, 10);
313 ok = (BOOL_T) (end[0] == NULL_CHAR);
314 } else {
315 expo = 0;
316 ok = (BOOL_T) (s[i] == NULL_CHAR);
317 }
318 // Calculate effective exponent.
319 if (dot >= 0) {
320 if (one > dot) {
321 expo -= one - dot + 1;
322 } else {
323 expo += dot - 1;
324 }
325 } else {
326 expo += pow - 1;
327 }
328 (void) align_mp (z, &expo, digs);
329 MP_EXPONENT (z) = (MP_DIGIT (z, 1) == 0 ? 0 : (MP_T) expo);
330 MP_DIGIT (z, 1) *= sign;
331 check_mp_exp (p, z);
332 if (errno == 0 && ok) {
333 return z;
334 } else {
335 return NaN_MP;
336 }
337 }
338
339 //! @brief Convert integer to multi-precison number.
340
341 MP_T *int_to_mp (NODE_T * p, MP_T * z, INT_T k, int digs)
342 {
343 int sign_k = 1;
344 if (k < 0) {
345 k = -k;
346 sign_k = -1;
347 }
348 int m = k, n = 0;
349 while ((m /= MP_RADIX) != 0) {
350 n++;
351 }
352 set_mp (z, 0, n, digs);
353 for (unt j = 1 + n; j >= 1; j--) {
354 MP_DIGIT (z, j) = (MP_T) (k % MP_RADIX);
355 k /= MP_RADIX;
356 }
357 MP_DIGIT (z, 1) = sign_k * MP_DIGIT (z, 1);
358 check_mp_exp (p, z);
359 return z;
360 }
361
362 //! @brief Convert unt to multi-precison number.
363
364 MP_T *unt_to_mp (NODE_T * p, MP_T * z, UNSIGNED_T k, int digs)
365 {
366 int m = k, n = 0;
367 while ((m /= MP_RADIX) != 0) {
368 n++;
369 }
370 set_mp (z, 0, n, digs);
371 for (unt j = 1 + n; j >= 1; j--) {
372 MP_DIGIT (z, j) = (MP_T) (k % MP_RADIX);
373 k /= MP_RADIX;
374 }
375 check_mp_exp (p, z);
376 return z;
377 }
378
379 //! @brief Convert multi-precision number to integer.
380
381 INT_T mp_to_int (NODE_T * p, MP_T * z, int digs)
382 {
383 // This routines looks a lot like "strtol".
384 INT_T expo = (int) MP_EXPONENT (z), sum = 0, weight = 1;
385 if (expo >= digs) {
386 diagnostic (A68_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, MOID (p));
387 exit_genie (p, A68_RUNTIME_ERROR);
388 }
389 BOOL_T negative = (BOOL_T) (MP_DIGIT (z, 1) < 0);
390 if (negative) {
391 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
392 }
393 for (int j = 1 + expo; j >= 1; j--) {
394 if ((MP_INT_T) MP_DIGIT (z, j) > A68_MAX_INT / weight) {
395 diagnostic (A68_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, M_INT);
396 exit_genie (p, A68_RUNTIME_ERROR);
397 }
398 INT_T term = (MP_INT_T) MP_DIGIT (z, j) * weight;
399 if (sum > A68_MAX_INT - term) {
400 diagnostic (A68_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, M_INT);
401 exit_genie (p, A68_RUNTIME_ERROR);
402 }
403 sum += term;
404 weight *= MP_RADIX;
405 }
406 return negative ? -sum : sum;
407 }
408
409 //! @brief Convert REAL_T to multi-precison number.
410
411 MP_T *real_to_mp (NODE_T * p, MP_T * z, REAL_T x, int digs)
412 {
413 SET_MP_ZERO (z, digs);
414 if (x == 0.0) {
415 return z;
416 }
417 // Small integers can be done better by int_to_mp.
418 if (ABS (x) < MP_RADIX && trunc (x) == x) {
419 return int_to_mp (p, z, (INT_T) trunc (x), digs);
420 }
421 int sign_x = SIGN (x);
422 // Scale to [0, 0.1>.
423 REAL_T a = ABS (x);
424 INT_T expo = (int) log10 (a);
425 a /= ten_up (expo);
426 expo--;
427 if (a >= 1) {
428 a /= 10;
429 expo++;
430 }
431 // Transport digs of x to the mantissa of z.
432 INT_T sum = 0, weight = (MP_RADIX / 10);
433 int j = 1;
434 for (unt k = 0; a != 0.0 && j <= digs && k < REAL_DIGITS; k++) {
435 REAL_T u = a * 10;
436 REAL_T v = floor (u);
437 a = u - v;
438 sum += weight * (INT_T) v;
439 weight /= 10;
440 if (weight < 1) {
441 MP_DIGIT (z, j++) = (MP_T) sum;
442 sum = 0;
443 weight = (MP_RADIX / 10);
444 }
445 }
446 // Store the last digs.
447 if (j <= digs) {
448 MP_DIGIT (z, j) = (MP_T) sum;
449 }
450 (void) align_mp (z, &expo, digs);
451 MP_EXPONENT (z) = (MP_T) expo;
452 MP_DIGIT (z, 1) *= sign_x;
453 check_mp_exp (p, z);
454 return z;
455 }
456
457 //! @brief Convert multi-precision number to real.
458
459 REAL_T mp_to_real (NODE_T * p, MP_T * z, int digs)
460 {
461 // This routine looks a lot like "strtod".
462 (void) p;
463 if (MP_EXPONENT (z) * (MP_T) LOG_MP_RADIX <= (MP_T) REAL_MIN_10_EXP) {
464 return 0;
465 } else {
466 REAL_T terms[1 + MP_MAX_DIGITS];
467 REAL_T weight = ten_up ((int) (MP_EXPONENT (z) * LOG_MP_RADIX));
468 unt lim = MIN (digs, MP_MAX_DIGITS);
469 for (unt k = 1; k <= lim; k++) {
470 terms[k] = ABS (MP_DIGIT (z, k)) * weight;
471 weight /= MP_RADIX;
472 }
473 // Sum terms from small to large.
474 REAL_T sum = 0;
475 for (unt k = lim; k >= 1; k--) {
476 sum += terms[k];
477 }
478 CHECK_REAL (p, sum);
479 return MP_DIGIT (z, 1) >= 0 ? sum : -sum;
480 }
481 }
482
483 //! @brief Normalise positive intermediate, fast.
484
485 static inline void norm_mp_light (MP_T * w, int k, int digs)
486 {
487 // Bring every digit back to [0 .. MP_RADIX>.
488 MP_T *z = &MP_DIGIT (w, digs);
489 for (unt j = digs; j >= k; j--, z--) {
490 if (z[0] >= MP_RADIX) {
491 z[0] -= (MP_T) MP_RADIX;
492 z[-1] += 1;
493 } else if (z[0] < 0) {
494 z[0] += (MP_T) MP_RADIX;
495 z[-1] -= 1;
496 }
497 }
498 }
499
500 //! @brief Normalise positive intermediate.
501
502 static inline void norm_mp (MP_T * w, int k, int digs)
503 {
504 // Bring every digit back to [0 .. MP_RADIX>.
505 unt j; MP_T *z;
506 for (j = digs, z = &MP_DIGIT (w, digs); j >= k; j--, z--) {
507 if (z[0] >= (MP_T) MP_RADIX) {
508 MP_T carry = (MP_T) ((MP_INT_T) (z[0] / (MP_T) MP_RADIX));
509 z[0] -= carry * (MP_T) MP_RADIX;
510 z[-1] += carry;
511 } else if (z[0] < (MP_T) 0) {
512 MP_T carry = (MP_T) 1 + (MP_T) ((MP_INT_T) ((-z[0] - 1) / (MP_T) MP_RADIX));
513 z[0] += carry * (MP_T) MP_RADIX;
514 z[-1] -= carry;
515 }
516 }
517 }
518
519 //! @brief Round multi-precision number.
520
521 static inline void round_internal_mp (MP_T * z, MP_T * w, int digs)
522 {
523 // Assume that w has precision of at least 2 + digs.
524 int last = (MP_DIGIT (w, 1) == 0 ? 2 + digs : 1 + digs);
525 if (MP_DIGIT (w, last) >= MP_RADIX / 2) {
526 MP_DIGIT (w, last - 1) += 1;
527 }
528 if (MP_DIGIT (w, last - 1) >= MP_RADIX) {
529 norm_mp (w, 2, last); // Hardly ever happens - no need to optimise
530 }
531 if (MP_DIGIT (w, 1) == 0) {
532 (void) move_mp_part (&MP_DIGIT (z, 1), &MP_DIGIT (w, 2), digs);
533 MP_EXPONENT (z) = MP_EXPONENT (w) - 1;
534 } else if (z != w) {
535 (void) move_mp_part (&MP_EXPONENT (z), &MP_EXPONENT (w), (1 + digs));
536 }
537 // Zero is zero is zero.
538 if (MP_DIGIT (z, 1) == 0) {
539 MP_EXPONENT (z) = (MP_T) 0;
540 }
541 }
542
543 //! @brief Truncate at decimal point.
544
545 MP_T *trunc_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
546 {
547 if (MP_EXPONENT (x) < 0) {
548 SET_MP_ZERO (z, digs);
549 } else if (MP_EXPONENT (x) >= (MP_T) digs) {
550 errno = EDOM;
551 diagnostic (A68_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, (IS (MOID (p), PROC_SYMBOL) ? SUB_MOID (p) : MOID (p)));
552 exit_genie (p, A68_RUNTIME_ERROR);
553 } else {
554 (void) move_mp (z, x, digs);
555 for (int k = (int) (MP_EXPONENT (x) + 2); k <= digs; k++) {
556 MP_DIGIT (z, k) = (MP_T) 0;
557 }
558 }
559 return z;
560 }
561
562 //! @brief Floor - largest integer smaller than x.
563
564 MP_T *floor_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
565 {
566 if (MP_EXPONENT (x) < 0) {
567 SET_MP_ZERO (z, digs);
568 } else if (MP_EXPONENT (x) >= (MP_T) digs) {
569 errno = EDOM;
570 diagnostic (A68_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, (IS (MOID (p), PROC_SYMBOL) ? SUB_MOID (p) : MOID (p)));
571 exit_genie (p, A68_RUNTIME_ERROR);
572 } else {
573 (void) move_mp (z, x, digs);
574 for (int k = (int) (MP_EXPONENT (x) + 2); k <= digs; k++) {
575 MP_DIGIT (z, k) = (MP_T) 0;
576 }
577 }
578 if (MP_DIGIT (x, 1) < 0 && ! same_mp (p, z, x, digs)) {
579 (void) minus_one_mp (p, z, z, digs);
580 }
581 return z;
582 }
583
584 BOOL_T is_int_mp (NODE_T *p, MP_T *z, int digs)
585 {
586 ADDR_T pop_sp = A68_SP;
587 MP_T *y = nil_mp (p, digs);
588 trunc_mp (p, y, z, digs);
589 BOOL_T tst = same_mp (p, y, z, digs);
590 A68_SP = pop_sp;
591 return tst;
592 }
593
594 //! @brief Shorten and round.
595
596 MP_T *shorten_mp (NODE_T * p, MP_T * z, int digs, MP_T * x, int digs_x)
597 {
598 if (digs > digs_x) {
599 return lengthen_mp (p, z, digs, x, digs_x);
600 } else if (digs == digs_x) {
601 return move_mp (z, x, digs);
602 } else {
603 // Reserve extra digs for proper rounding.
604 ADDR_T pop_sp = A68_SP;
605 int digs_h = digs + 2;
606 BOOL_T negative = (BOOL_T) (MP_DIGIT (x, 1) < 0);
607 MP_T *w = nil_mp (p, digs_h);
608 if (negative) {
609 MP_DIGIT (x, 1) = -MP_DIGIT (x, 1);
610 }
611 MP_STATUS (w) = (MP_T) 0;
612 MP_EXPONENT (w) = MP_EXPONENT (x) + 1;
613 MP_DIGIT (w, 1) = (MP_T) 0;
614 (void) move_mp_part (&MP_DIGIT (w, 2), &MP_DIGIT (x, 1), digs + 1);
615 round_internal_mp (z, w, digs);
616 if (negative) {
617 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
618 }
619 A68_SP = pop_sp;
620 return z;
621 }
622 }
623
624 //! @brief Lengthen x and assign to z.
625
626 MP_T *lengthen_mp (NODE_T * p, MP_T * z, int digs_z, MP_T * x, int digs_x)
627 {
628 if (digs_z < digs_x) {
629 return shorten_mp (p, z, digs_z, x, digs_x);
630 } else if (digs_z == digs_x) {
631 return move_mp (z, x, digs_z);
632 } else {
633 if (z != x) {
634 (void) move_mp_part (&MP_DIGIT (z, 1), &MP_DIGIT (x, 1), digs_x);
635 MP_EXPONENT (z) = MP_EXPONENT (x);
636 MP_STATUS (z) = MP_STATUS (x);
637 }
638 for (unt j = 1 + digs_x; j <= digs_z; j++) {
639 MP_DIGIT (z, j) = (MP_T) 0;
640 }
641 }
642 return z;
643 }
644
645 //! @brief Set "z" to the sum of positive "x" and positive "y".
646
647 MP_T *add_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
648 {
649 MP_STATUS (z) = (MP_T) INIT_MASK;
650 // Trivial cases.
651 if (MP_DIGIT (x, 1) == (MP_T) 0) {
652 (void) move_mp (z, y, digs);
653 return z;
654 } else if (MP_DIGIT (y, 1) == 0) {
655 (void) move_mp (z, x, digs);
656 return z;
657 }
658 // We want positive arguments.
659 ADDR_T pop_sp = A68_SP;
660 MP_T x_1 = MP_DIGIT (x, 1), y_1 = MP_DIGIT (y, 1);
661 MP_DIGIT (x, 1) = ABS (x_1);
662 MP_DIGIT (y, 1) = ABS (y_1);
663 if (x_1 >= 0 && y_1 < 0) {
664 (void) sub_mp (p, z, x, y, digs);
665 } else if (x_1 < 0 && y_1 >= 0) {
666 (void) sub_mp (p, z, y, x, digs);
667 } else if (x_1 < 0 && y_1 < 0) {
668 (void) add_mp (p, z, x, y, digs);
669 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
670 } else {
671 // Add.
672 int digs_h = 2 + digs;
673 MP_T *w = nil_mp (p, digs_h);
674 if (MP_EXPONENT (x) == MP_EXPONENT (y)) {
675 MP_EXPONENT (w) = (MP_T) 1 + MP_EXPONENT (x);
676 for (unt j = 1; j <= digs; j++) {
677 MP_DIGIT (w, j + 1) = MP_DIGIT (x, j) + MP_DIGIT (y, j);
678 }
679 MP_DIGIT (w, digs_h) = (MP_T) 0;
680 } else if (MP_EXPONENT (x) > MP_EXPONENT (y)) {
681 int shl_y = (int) MP_EXPONENT (x) - (int) MP_EXPONENT (y);
682 MP_EXPONENT (w) = (MP_T) 1 + MP_EXPONENT (x);
683 for (unt j = 1; j < digs_h; j++) {
684 int i_y = j - shl_y;
685 MP_T x_j = (j > digs ? 0 : MP_DIGIT (x, j));
686 MP_T y_j = (i_y <= 0 || i_y > digs ? 0 : MP_DIGIT (y, i_y));
687 MP_DIGIT (w, j + 1) = x_j + y_j;
688 }
689 } else {
690 int shl_x = (int) MP_EXPONENT (y) - (int) MP_EXPONENT (x);
691 MP_EXPONENT (w) = (MP_T) 1 + MP_EXPONENT (y);
692 for (unt j = 1; j < digs_h; j++) {
693 int i_x = j - shl_x;
694 MP_T x_j = (i_x <= 0 || i_x > digs ? 0 : MP_DIGIT (x, i_x));
695 MP_T y_j = (j > digs ? 0 : MP_DIGIT (y, j));
696 MP_DIGIT (w, j + 1) = x_j + y_j;
697 }
698 }
699 norm_mp_light (w, 2, digs_h);
700 round_internal_mp (z, w, digs);
701 check_mp_exp (p, z);
702 }
703 // Restore and exit.
704 A68_SP = pop_sp;
705 MP_T z_1 = MP_DIGIT (z, 1);
706 MP_DIGIT (x, 1) = x_1;
707 MP_DIGIT (y, 1) = y_1;
708 MP_DIGIT (z, 1) = z_1; // In case z IS x OR z IS y
709 return z;
710 }
711
712 //! @brief Set "z" to the difference of positive "x" and positive "y".
713
714 MP_T *sub_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
715 {
716 MP_STATUS (z) = (MP_T) INIT_MASK;
717 // Trivial cases.
718 if (MP_DIGIT (x, 1) == (MP_T) 0) {
719 (void) move_mp (z, y, digs);
720 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
721 return z;
722 } else if (MP_DIGIT (y, 1) == (MP_T) 0) {
723 (void) move_mp (z, x, digs);
724 return z;
725 }
726 // We want positive arguments.
727 ADDR_T pop_sp = A68_SP;
728 MP_T x_1 = MP_DIGIT (x, 1), y_1 = MP_DIGIT (y, 1);
729 MP_DIGIT (x, 1) = ABS (x_1);
730 MP_DIGIT (y, 1) = ABS (y_1);
731 if (x_1 >= 0 && y_1 < 0) {
732 (void) add_mp (p, z, x, y, digs);
733 } else if (x_1 < 0 && y_1 >= 0) {
734 (void) add_mp (p, z, y, x, digs);
735 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
736 } else if (x_1 < 0 && y_1 < 0) {
737 (void) sub_mp (p, z, y, x, digs);
738 } else {
739 // Subtract.
740 BOOL_T negative = A68_FALSE;
741 int fnz, digs_h = 2 + digs;
742 MP_T *w = nil_mp (p, digs_h);
743 if (MP_EXPONENT (x) == MP_EXPONENT (y)) {
744 MP_EXPONENT (w) = (MP_T) 1 + MP_EXPONENT (x);
745 for (unt j = 1; j <= digs; j++) {
746 MP_DIGIT (w, j + 1) = MP_DIGIT (x, j) - MP_DIGIT (y, j);
747 }
748 MP_DIGIT (w, digs_h) = (MP_T) 0;
749 } else if (MP_EXPONENT (x) > MP_EXPONENT (y)) {
750 int shl_y = (int) MP_EXPONENT (x) - (int) MP_EXPONENT (y);
751 MP_EXPONENT (w) = (MP_T) 1 + MP_EXPONENT (x);
752 for (unt j = 1; j < digs_h; j++) {
753 int i_y = j - shl_y;
754 MP_T x_j = (j > digs ? 0 : MP_DIGIT (x, j));
755 MP_T y_j = (i_y <= 0 || i_y > digs ? 0 : MP_DIGIT (y, i_y));
756 MP_DIGIT (w, j + 1) = x_j - y_j;
757 }
758 } else {
759 int shl_x = (int) MP_EXPONENT (y) - (int) MP_EXPONENT (x);
760 MP_EXPONENT (w) = (MP_T) 1 + MP_EXPONENT (y);
761 for (unt j = 1; j < digs_h; j++) {
762 int i_x = j - shl_x;
763 MP_T x_j = (i_x <= 0 || i_x > digs ? 0 : MP_DIGIT (x, i_x));
764 MP_T y_j = (j > digs ? 0 : MP_DIGIT (y, j));
765 MP_DIGIT (w, j + 1) = x_j - y_j;
766 }
767 }
768 // Correct if we subtract large from small.
769 if (MP_DIGIT (w, 2) <= 0) {
770 fnz = -1;
771 for (unt j = 2; j <= digs_h && fnz < 0; j++) {
772 if (MP_DIGIT (w, j) != 0) {
773 fnz = j;
774 }
775 }
776 negative = (BOOL_T) (MP_DIGIT (w, fnz) < 0);
777 if (negative) {
778 for (unt j = fnz; j <= digs_h; j++) {
779 MP_DIGIT (w, j) = -MP_DIGIT (w, j);
780 }
781 }
782 }
783 // Normalise.
784 norm_mp_light (w, 2, digs_h);
785 fnz = -1;
786 for (unt j = 1; j <= digs_h && fnz < 0; j++) {
787 if (MP_DIGIT (w, j) != 0) {
788 fnz = j;
789 }
790 }
791 if (fnz > 1) {
792 int j2 = fnz - 1;
793 for (int j = 1; j <= digs_h - j2; j++) {
794 MP_DIGIT (w, j) = MP_DIGIT (w, j + j2);
795 MP_DIGIT (w, j + j2) = (MP_T) 0;
796 }
797 MP_EXPONENT (w) -= j2;
798 }
799 // Round.
800 round_internal_mp (z, w, digs);
801 if (negative) {
802 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
803 }
804 check_mp_exp (p, z);
805 }
806 // Restore and exit.
807 A68_SP = pop_sp;
808 MP_T z_1 = MP_DIGIT (z, 1);
809 MP_DIGIT (x, 1) = x_1;
810 MP_DIGIT (y, 1) = y_1;
811 MP_DIGIT (z, 1) = z_1; // In case z IS x OR z IS y
812 return z;
813 }
814
815 //! @brief Set "z" to the product of "x" and "y".
816
817 MP_T *mul_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
818 {
819 if (IS_ZERO_MP (x) || IS_ZERO_MP (y)) {
820 SET_MP_ZERO (z, digs);
821 return z;
822 }
823 // Grammar school algorithm with intermittent normalisation.
824 ADDR_T pop_sp = A68_SP;
825 int digs_h = 2 + digs;
826 MP_T x_1 = MP_DIGIT (x, 1), y_1 = MP_DIGIT (y, 1);
827 MP_DIGIT (x, 1) = ABS (x_1);
828 MP_DIGIT (y, 1) = ABS (y_1);
829 MP_STATUS (z) = (MP_T) INIT_MASK;
830 MP_T *w = lit_mp (p, 0, MP_EXPONENT (x) + MP_EXPONENT (y) + 1, digs_h);
831 int oflow = (int) FLOOR_MP ((MP_REAL_T) MAX_REPR_INT / (2 * MP_REAL_RADIX * MP_REAL_RADIX)) - 1;
832 for (unt i = digs; i >= 1; i--) {
833 MP_T yi = MP_DIGIT (y, i);
834 if (yi != 0) {
835 int k = digs_h - i;
836 int j = (k > digs ? digs : k);
837 MP_T *u = &MP_DIGIT (w, i + j), *v = &MP_DIGIT (x, j);
838 if ((digs - i + 1) % oflow == 0) {
839 norm_mp (w, 2, digs_h);
840 }
841 while (j-- >= 1) {
842 (u--)[0] += yi * (v--)[0];
843 }
844 }
845 }
846 norm_mp (w, 2, digs_h);
847 round_internal_mp (z, w, digs);
848 // Restore and exit.
849 A68_SP = pop_sp;
850 MP_T z_1 = MP_DIGIT (z, 1);
851 MP_DIGIT (x, 1) = x_1;
852 MP_DIGIT (y, 1) = y_1;
853 MP_DIGIT (z, 1) = ((x_1 * y_1) >= 0 ? z_1 : -z_1);
854 check_mp_exp (p, z);
855 return z;
856 }
857
858 //! @brief Set "z" to the quotient of "x" and "y".
859
860 MP_T *div_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
861 {
862 // This routine is based on
863 //
864 // D. M. Smith, "A Multiple-Precision Division Algorithm"
865 // Mathematics of Computation 66 (1996) 157-163.
866 //
867 // This is O(N^2) but runs faster than straightforward methods by skipping
868 // most of the intermediate normalisation and recovering from wrong
869 // guesses without separate correction steps.
870 //
871 // Depending on application, div_mp cost is circa 3 times that of mul_mp.
872 // Therefore Newton-Raphson division makes no sense here.
873 //
874 if (IS_ZERO_MP (y)) {
875 errno = ERANGE;
876 return NaN_MP;
877 }
878 // Determine normalisation interval assuming that q < 2b in each step.
879 #if (A68_LEVEL <= 2)
880 int oflow = (int) FLOOR_MP ((MP_REAL_T) MAX_REPR_INT / (3 * MP_REAL_RADIX * MP_REAL_RADIX)) - 1;
881 #else
882 int oflow = (int) FLOOR_MP ((MP_REAL_T) MAX_REPR_INT / (2 * MP_REAL_RADIX * MP_REAL_RADIX)) - 1;
883 #endif
884 //
885 MP_T x_1 = MP_DIGIT (x, 1), y_1 = MP_DIGIT (y, 1);
886 MP_DIGIT (x, 1) = ABS (x_1);
887 MP_DIGIT (y, 1) = ABS (y_1);
888 MP_STATUS (z) = (MP_T) INIT_MASK;
889 // Slight optimisation when the denominator has few digits.
890 int nzdigs = digs;
891 while (MP_DIGIT (y, nzdigs) == 0 && nzdigs > 1) {
892 nzdigs--;
893 }
894 if (nzdigs == 1 && MP_EXPONENT (y) == 0) {
895 (void) div_mp_digit (p, z, x, MP_DIGIT (y, 1), digs);
896 MP_T z_1 = MP_DIGIT (z, 1);
897 MP_DIGIT (x, 1) = x_1;
898 MP_DIGIT (y, 1) = y_1;
899 MP_DIGIT (z, 1) = ((x_1 * y_1) >= 0 ? z_1 : -z_1);
900 check_mp_exp (p, z);
901 return z;
902 }
903 // Working nominator in which the quotient develops.
904 ADDR_T pop_sp = A68_SP;
905 int wdigs = 4 + digs;
906 MP_T *w = lit_mp (p, 0, MP_EXPONENT (x) - MP_EXPONENT (y), wdigs);
907 (void) move_mp_part (&MP_DIGIT (w, 2), &MP_DIGIT (x, 1), digs);
908 // Estimate the denominator. For small MP_RADIX add: MP_DIGIT (y, 4) / MP_REAL_RADIX.
909 MP_REAL_T den = (MP_DIGIT (y, 1) * MP_REAL_RADIX + MP_DIGIT (y, 2)) * MP_REAL_RADIX + MP_DIGIT (y, 3);
910 MP_T *t = &MP_DIGIT (w, 2);
911 for (unt k = 1, len = digs + 2, first = 3; k <= digs + 2; k++, len++, first++, t++) {
912 // Estimate quotient digit.
913 MP_REAL_T q, nom = ((t[-1] * MP_REAL_RADIX + t[0]) * MP_REAL_RADIX + t[1]) * MP_REAL_RADIX + (wdigs >= (first + 2) ? t[2] : 0);
914 if (nom == 0) {
915 q = 0;
916 } else {
917 // Correct the nominator.
918 q = (MP_T) (MP_INT_T) (nom / den);
919 int lim = MINIMUM (len, wdigs);
920 if (nzdigs <= lim - first + 1) {
921 lim = first + nzdigs - 1;
922 }
923 MP_T *u = t, *v = &MP_DIGIT (y, 1);
924 for (unt j = first; j <= lim; j++) {
925 (u++)[0] -= q * (v++)[0];
926 }
927 }
928 t[0] += t[-1] * MP_RADIX;
929 t[-1] = q;
930 if (k % oflow == 0 || k == digs + 2) {
931 norm_mp (w, first, wdigs);
932 }
933 }
934 norm_mp (w, 2, digs);
935 round_internal_mp (z, w, digs);
936 // Restore and exit.
937 A68_SP = pop_sp;
938 MP_T z_1 = MP_DIGIT (z, 1);
939 MP_DIGIT (x, 1) = x_1;
940 MP_DIGIT (y, 1) = y_1;
941 MP_DIGIT (z, 1) = ((x_1 * y_1) >= 0 ? z_1 : -z_1);
942 check_mp_exp (p, z);
943 return z;
944 }
945
946 //! @brief Set "z" to the integer quotient of "x" and "y".
947
948 MP_T *over_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
949 {
950 if (MP_DIGIT (y, 1) == 0) {
951 errno = ERANGE;
952 return NaN_MP;
953 }
954 int digs_g = FUN_DIGITS (digs);
955 ADDR_T pop_sp = A68_SP;
956 MP_T *x_g = len_mp (p, x, digs, digs_g);
957 MP_T *y_g = len_mp (p, y, digs, digs_g);
958 MP_T *z_g = nil_mp (p, digs_g);
959 (void) div_mp (p, z_g, x_g, y_g, digs_g);
960 trunc_mp (p, z_g, z_g, digs_g);
961 (void) shorten_mp (p, z, digs, z_g, digs_g);
962 MP_STATUS (z) = (MP_T) INIT_MASK;
963 // Restore and exit.
964 A68_SP = pop_sp;
965 return z;
966 }
967
968 //! @brief Set "z" to x mod y.
969
970 MP_T *mod_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
971 {
972 if (MP_DIGIT (y, 1) == 0) {
973 errno = EDOM;
974 return NaN_MP;
975 }
976 int digs_g = FUN_DIGITS (digs);
977 ADDR_T pop_sp = A68_SP;
978 MP_T *x_g = len_mp (p, x, digs, digs_g);
979 MP_T *y_g = len_mp (p, y, digs, digs_g);
980 MP_T *z_g = nil_mp (p, digs_g);
981 // x mod y = x - y * trunc (x / y).
982 (void) over_mp (p, z_g, x_g, y_g, digs_g);
983 (void) mul_mp (p, z_g, y_g, z_g, digs_g);
984 (void) sub_mp (p, z_g, x_g, z_g, digs_g);
985 (void) shorten_mp (p, z, digs, z_g, digs_g);
986 // Restore and exit.
987 A68_SP = pop_sp;
988 return z;
989 }
990
991 //! @brief Set "z" to the product of x and digit y.
992
993 MP_T *mul_mp_digit (NODE_T * p, MP_T * z, MP_T * x, MP_T y, int digs)
994 {
995 // This is an O(N) routine for multiplication by a short value.
996 MP_T x_1 = MP_DIGIT (x, 1);
997 int digs_h = 2 + digs;
998 ADDR_T pop_sp = A68_SP;
999 MP_DIGIT (x, 1) = ABS (x_1);
1000 MP_STATUS (z) = (MP_T) INIT_MASK;
1001 MP_T y_1 = y;
1002 y = ABS (y_1);
1003 if (y == 2) {
1004 (void) add_mp (p, z, x, x, digs);
1005 } else {
1006 MP_T *w = lit_mp (p, 0, MP_EXPONENT (x) + 1, digs_h);
1007 MP_T *u = &MP_DIGIT (w, 1 + digs), *v = &MP_DIGIT (x, digs);
1008 int j = digs;
1009 while (j-- >= 1) {
1010 (u--)[0] += y * (v--)[0];
1011 }
1012 norm_mp (w, 2, digs_h);
1013 round_internal_mp (z, w, digs);
1014 }
1015 // Restore and exit.
1016 A68_SP = pop_sp;
1017 MP_T z_1 = MP_DIGIT (z, 1);
1018 MP_DIGIT (x, 1) = x_1;
1019 MP_DIGIT (z, 1) = ((x_1 * y_1) >= 0 ? z_1 : -z_1);
1020 check_mp_exp (p, z);
1021 return z;
1022 }
1023
1024 //! @brief Set "z" to x/2.
1025
1026 MP_T *half_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1027 {
1028 MP_T *w, z_1, x_1 = MP_DIGIT (x, 1), *u, *v;
1029 int j, digs_h = 2 + digs;
1030 ADDR_T pop_sp = A68_SP;
1031 MP_DIGIT (x, 1) = ABS (x_1);
1032 MP_STATUS (z) = (MP_T) INIT_MASK;
1033 // Calculate x * 0.5.
1034 w = lit_mp (p, 0, MP_EXPONENT (x), digs_h);
1035 j = digs;
1036 u = &MP_DIGIT (w, 1 + digs);
1037 v = &MP_DIGIT (x, digs);
1038 while (j-- >= 1) {
1039 (u--)[0] += (MP_RADIX / 2) * (v--)[0];
1040 }
1041 norm_mp (w, 2, digs_h);
1042 round_internal_mp (z, w, digs);
1043 // Restore and exit.
1044 A68_SP = pop_sp;
1045 z_1 = MP_DIGIT (z, 1);
1046 MP_DIGIT (x, 1) = x_1;
1047 MP_DIGIT (z, 1) = (x_1 >= 0 ? z_1 : -z_1);
1048 check_mp_exp (p, z);
1049 return z;
1050 }
1051
1052 //! @brief Set "z" to x/10.
1053
1054 MP_T *tenth_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1055 {
1056 MP_T *w, z_1, x_1 = MP_DIGIT (x, 1), *u, *v;
1057 int j, digs_h = 2 + digs;
1058 ADDR_T pop_sp = A68_SP;
1059 MP_DIGIT (x, 1) = ABS (x_1);
1060 MP_STATUS (z) = (MP_T) INIT_MASK;
1061 // Calculate x * 0.1.
1062 w = lit_mp (p, 0, MP_EXPONENT (x), digs_h);
1063 j = digs;
1064 u = &MP_DIGIT (w, 1 + digs);
1065 v = &MP_DIGIT (x, digs);
1066 while (j-- >= 1) {
1067 (u--)[0] += (MP_RADIX / 10) * (v--)[0];
1068 }
1069 norm_mp (w, 2, digs_h);
1070 round_internal_mp (z, w, digs);
1071 // Restore and exit.
1072 A68_SP = pop_sp;
1073 z_1 = MP_DIGIT (z, 1);
1074 MP_DIGIT (x, 1) = x_1;
1075 MP_DIGIT (z, 1) = (x_1 >= 0 ? z_1 : -z_1);
1076 check_mp_exp (p, z);
1077 return z;
1078 }
1079
1080 //! @brief Set "z" to the quotient of x and digit y.
1081
1082 MP_T *div_mp_digit (NODE_T * p, MP_T * z, MP_T * x, MP_T y, int digs)
1083 {
1084 if (y == 0) {
1085 errno = ERANGE;
1086 return NaN_MP;
1087 }
1088 // Determine normalisation interval assuming that q < 2b in each step.
1089 #if (A68_LEVEL <= 2)
1090 int oflow = (int) FLOOR_MP ((MP_REAL_T) MAX_REPR_INT / (3 * MP_REAL_RADIX * MP_REAL_RADIX)) - 1;
1091 #else
1092 int oflow = (int) FLOOR_MP ((MP_REAL_T) MAX_REPR_INT / (2 * MP_REAL_RADIX * MP_REAL_RADIX)) - 1;
1093 #endif
1094 // Work with positive operands.
1095 ADDR_T pop_sp = A68_SP;
1096 MP_T x_1 = MP_DIGIT (x, 1), y_1 = y;
1097 MP_DIGIT (x, 1) = ABS (x_1);
1098 MP_STATUS (z) = (MP_T) INIT_MASK;
1099 y = ABS (y_1);
1100 //
1101 if (y == 2) {
1102 (void) half_mp (p, z, x, digs);
1103 } else if (y == 10) {
1104 (void) tenth_mp (p, z, x, digs);
1105 } else {
1106 int wdigs = 4 + digs;
1107 MP_T *w = lit_mp (p, 0, MP_EXPONENT (x), wdigs);
1108 (void) move_mp_part (&MP_DIGIT (w, 2), &MP_DIGIT (x, 1), digs);
1109 // Estimate the denominator.
1110 MP_REAL_T den = (MP_REAL_T) y * MP_REAL_RADIX * MP_REAL_RADIX;
1111 MP_T *t = &MP_DIGIT (w, 2);
1112 for (unt k = 1, first = 3; k <= digs + 2; k++, first++, t++) {
1113 // Estimate quotient digit and correct.
1114 MP_REAL_T nom = ((t[-1] * MP_REAL_RADIX + t[0]) * MP_REAL_RADIX + t[1]) * MP_REAL_RADIX + (wdigs >= (first + 2) ? t[2] : 0);
1115 MP_REAL_T q = (MP_T) (MP_INT_T) (nom / den);
1116 t[0] += t[-1] * MP_RADIX - q * y;
1117 t[-1] = q;
1118 if (k % oflow == 0 || k == digs + 2) {
1119 norm_mp (w, first, wdigs);
1120 }
1121 }
1122 norm_mp (w, 2, digs);
1123 round_internal_mp (z, w, digs);
1124 }
1125 // Restore and exit.
1126 A68_SP = pop_sp;
1127 MP_T z_1 = MP_DIGIT (z, 1);
1128 MP_DIGIT (x, 1) = x_1;
1129 MP_DIGIT (z, 1) = ((x_1 * y_1) >= 0 ? z_1 : -z_1);
1130 check_mp_exp (p, z);
1131 return z;
1132 }
1133
1134 //! @brief Set "z" to the integer quotient of "x" and "y".
1135
1136 MP_T *over_mp_digit (NODE_T * p, MP_T * z, MP_T * x, MP_T y, int digs)
1137 {
1138 if (y == 0) {
1139 errno = ERANGE;
1140 return NaN_MP;
1141 }
1142 int digs_g = FUN_DIGITS (digs);
1143 ADDR_T pop_sp = A68_SP;
1144 MP_T *x_g = len_mp (p, x, digs, digs_g);
1145 MP_T *z_g = nil_mp (p, digs_g);
1146 (void) div_mp_digit (p, z_g, x_g, y, digs_g);
1147 trunc_mp (p, z_g, z_g, digs_g);
1148 (void) shorten_mp (p, z, digs, z_g, digs_g);
1149 // Restore and exit.
1150 A68_SP = pop_sp;
1151 return z;
1152 }
1153
1154 //! @brief Set "z" to the reciprocal of "x".
1155
1156 MP_T *rec_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1157 {
1158 if (IS_ZERO_MP (x)) {
1159 errno = ERANGE;
1160 return NaN_MP;
1161 }
1162 ADDR_T pop_sp = A68_SP;
1163 (void) div_mp (p, z, mp_one (digs), x, digs);
1164 A68_SP = pop_sp;
1165 return z;
1166 }
1167
1168 //! @brief LONG REAL long pi
1169
1170 void genie_pi_mp (NODE_T * p)
1171 {
1172 int digs = DIGITS (MOID (p));
1173 MP_T *z = nil_mp (p, digs);
1174 (void) mp_pi (p, z, MP_PI, digs);
1175 MP_STATUS (z) = (MP_T) INIT_MASK;
1176 }
1177
1178 //! @brief Set "z" to "x" ** "n".
1179
1180 MP_T *pow_mp_int (NODE_T * p, MP_T * z, MP_T * x, INT_T n, int digs)
1181 {
1182 ADDR_T pop_sp = A68_SP;
1183 int bit, digs_g = FUN_DIGITS (digs);
1184 BOOL_T negative;
1185 MP_T *x_g = len_mp (p, x, digs, digs_g);
1186 MP_T *z_g = lit_mp (p, 1, 0, digs_g);
1187 negative = (BOOL_T) (n < 0);
1188 if (negative) {
1189 n = -n;
1190 }
1191 bit = 1;
1192 while ((int) bit <= (int) n) {
1193 if (n & bit) {
1194 (void) mul_mp (p, z_g, z_g, x_g, digs_g);
1195 }
1196 (void) mul_mp (p, x_g, x_g, x_g, digs_g);
1197 bit <<= 1;
1198 }
1199 (void) shorten_mp (p, z, digs, z_g, digs_g);
1200 A68_SP = pop_sp;
1201 if (negative) {
1202 (void) rec_mp (p, z, z, digs);
1203 }
1204 check_mp_exp (p, z);
1205 return z;
1206 }
1207
1208 //! @brief Set "z" to "x" ** "y".
1209
1210 MP_T *pow_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
1211 {
1212 PRELUDE_ERROR (ln_mp (p, z, x, digs) == NaN_MP, p, ERROR_INVALID_ARGUMENT, MOID (p));
1213 (void) mul_mp (p, z, y, z, digs);
1214 (void) exp_mp (p, z, z, digs);
1215 return z;
1216 }
1217
1218 //! @brief Set "z" to 10 ** "n".
1219
1220 MP_T *ten_up_mp (NODE_T * p, MP_T * z, int n, int digs)
1221 {
1222 #if (A68_LEVEL >= 3)
1223 static MP_T y[LOG_MP_RADIX] = { 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000 };
1224 #else
1225 static MP_T y[LOG_MP_RADIX] = { 1, 10, 100, 1000, 10000, 100000, 1000000 };
1226 #endif
1227 if (n >= 0) {
1228 set_mp (z, y[n % LOG_MP_RADIX], n / LOG_MP_RADIX, digs);
1229 } else {
1230 set_mp (z, y[(LOG_MP_RADIX + n % LOG_MP_RADIX) % LOG_MP_RADIX], (n + 1) / LOG_MP_RADIX - 1, digs);
1231 }
1232 check_mp_exp (p, z);
1233 return z;
1234 }
1235
1236 //! @brief Comparison of "x" and "y".
1237
1238 void eq_mp (NODE_T * p, A68_BOOL * z, MP_T * x, MP_T * y, int digs)
1239 {
1240 ADDR_T pop_sp = A68_SP;
1241 MP_T *v = nil_mp (p, digs);
1242 (void) sub_mp (p, v, x, y, digs);
1243 STATUS (z) = INIT_MASK;
1244 VALUE (z) = (MP_DIGIT (v, 1) == 0 ? A68_TRUE : A68_FALSE);
1245 A68_SP = pop_sp;
1246 }
1247
1248 //! @brief Comparison of "x" and "y".
1249
1250 void ne_mp (NODE_T * p, A68_BOOL * z, MP_T * x, MP_T * y, int digs)
1251 {
1252 ADDR_T pop_sp = A68_SP;
1253 MP_T *v = nil_mp (p, digs);
1254 (void) sub_mp (p, v, x, y, digs);
1255 STATUS (z) = INIT_MASK;
1256 VALUE (z) = (MP_DIGIT (v, 1) != 0 ? A68_TRUE : A68_FALSE);
1257 A68_SP = pop_sp;
1258 }
1259
1260 //! @brief Comparison of "x" and "y".
1261
1262 void lt_mp (NODE_T * p, A68_BOOL * z, MP_T * x, MP_T * y, int digs)
1263 {
1264 ADDR_T pop_sp = A68_SP;
1265 MP_T *v = nil_mp (p, digs);
1266 (void) sub_mp (p, v, x, y, digs);
1267 STATUS (z) = INIT_MASK;
1268 VALUE (z) = (MP_DIGIT (v, 1) < 0 ? A68_TRUE : A68_FALSE);
1269 A68_SP = pop_sp;
1270 }
1271
1272 //! @brief Comparison of "x" and "y".
1273
1274 void le_mp (NODE_T * p, A68_BOOL * z, MP_T * x, MP_T * y, int digs)
1275 {
1276 ADDR_T pop_sp = A68_SP;
1277 MP_T *v = nil_mp (p, digs);
1278 (void) sub_mp (p, v, x, y, digs);
1279 STATUS (z) = INIT_MASK;
1280 VALUE (z) = (MP_DIGIT (v, 1) <= 0 ? A68_TRUE : A68_FALSE);
1281 A68_SP = pop_sp;
1282 }
1283
1284 //! @brief Comparison of "x" and "y".
1285
1286 void gt_mp (NODE_T * p, A68_BOOL * z, MP_T * x, MP_T * y, int digs)
1287 {
1288 ADDR_T pop_sp = A68_SP;
1289 MP_T *v = nil_mp (p, digs);
1290 (void) sub_mp (p, v, x, y, digs);
1291 STATUS (z) = INIT_MASK;
1292 VALUE (z) = (MP_DIGIT (v, 1) > 0 ? A68_TRUE : A68_FALSE);
1293 A68_SP = pop_sp;
1294 }
1295
1296 //! @brief Comparison of "x" and "y".
1297
1298 void ge_mp (NODE_T * p, A68_BOOL * z, MP_T * x, MP_T * y, int digs)
1299 {
1300 ADDR_T pop_sp = A68_SP;
1301 MP_T *v = nil_mp (p, digs);
1302 (void) sub_mp (p, v, x, y, digs);
1303 STATUS (z) = INIT_MASK;
1304 VALUE (z) = (MP_DIGIT (v, 1) >= 0 ? A68_TRUE : A68_FALSE);
1305 A68_SP = pop_sp;
1306 }
1307
1308 //! @brief round (x).
1309
1310 MP_T *round_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1311 {
1312 MP_T *y = nil_mp (p, digs);
1313 SET_MP_HALF (y, digs);
1314 if (MP_DIGIT (x, 1) >= 0) {
1315 (void) add_mp (p, z, x, y, digs);
1316 (void) trunc_mp (p, z, z, digs);
1317 } else {
1318 (void) sub_mp (p, z, x, y, digs);
1319 (void) trunc_mp (p, z, z, digs);
1320 }
1321 MP_STATUS (z) = (MP_T) INIT_MASK;
1322 return z;
1323 }
1324
1325 //! @brief Entier (x).
1326
1327 MP_T *entier_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1328 {
1329 if (MP_DIGIT (x, 1) >= 0) {
1330 (void) trunc_mp (p, z, x, digs);
1331 } else {
1332 MP_T *y = nil_mp (p, digs);
1333 (void) move_mp (y, z, digs);
1334 (void) trunc_mp (p, z, x, digs);
1335 (void) sub_mp (p, y, y, z, digs);
1336 if (MP_DIGIT (y, 1) != 0) {
1337 SET_MP_ONE (y, digs);
1338 (void) sub_mp (p, z, z, y, digs);
1339 }
1340 }
1341 MP_STATUS (z) = (MP_T) INIT_MASK;
1342 return z;
1343 }
