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