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