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-2026 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 A68G_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_LINUX)
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 (A68G_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" A68G_LD, (MP_INT_T) MP_EXPONENT (z));
139 fprintf (stdout, " S" A68G_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) A68G_MP (varying_mp_digits)));
154 }
155 return A68G_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 = A68G_SP;
189 (void) sub_mp (p, z, mp_one (digs), x, digs);
190 A68G_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 = A68G_SP;
199 (void) sub_mp (p, z, x, mp_one (digs), digs);
200 A68G_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 = A68G_SP;
209 (void) add_mp (p, z, x, mp_one (digs), digs);
210 A68G_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 A68G_FALSE;
223 }
224 }
225 return A68G_TRUE;
226 } else {
227 return A68G_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 = A68G_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_T 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 (A68G_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, MOID (p));
381 exit_genie (p, A68G_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) > A68G_MAX_INT / weight) {
389 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, M_INT);
390 exit_genie (p, A68G_RUNTIME_ERROR);
391 }
392 INT_T term = (MP_INT_T) MP_DIGIT (z, j) * weight;
393 if (sum > A68G_MAX_INT - term) {
394 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, M_INT);
395 exit_genie (p, A68G_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 (a68g_isinf_real (x)) {
409 if (x == A68G_POSINF_REAL) {
410 MP_STATUS (z) = (PLUS_INF_MASK | INIT_MASK);
411 } else {
412 MP_STATUS (z) = (MINUS_INF_MASK | INIT_MASK);
413 }
414 return z;
415 }
416 if (x == 0.0) {
417 return z;
418 }
419 // Small integers can be done better by int_to_mp.
420 if (ABS (x) < MP_RADIX && trunc (x) == x) {
421 return int_to_mp (p, z, (INT_T) trunc (x), digs);
422 }
423 int sign_x = SIGN (x);
424 // Scale to [0, 0.1>.
425 REAL_T a = ABS (x);
426 INT_T expo = (int) log10 (a);
427 a /= ten_up (expo);
428 expo--;
429 if (a >= 1) {
430 a /= 10;
431 expo++;
432 }
433 // Transport digs of x to the mantissa of z.
434 INT_T sum = 0, weight = (MP_RADIX / 10);
435 int j = 1;
436 for (int k = 0; a != 0.0 && j <= digs && k < A68G_REAL_DIG; k++) {
437 REAL_T u = a * 10;
438 REAL_T v = floor (u);
439 a = u - v;
440 sum += weight * (INT_T) v;
441 weight /= 10;
442 if (weight < 1) {
443 MP_DIGIT (z, j++) = (MP_T) sum;
444 sum = 0;
445 weight = (MP_RADIX / 10);
446 }
447 }
448 // Store the last digs.
449 if (j <= digs) {
450 MP_DIGIT (z, j) = (MP_T) sum;
451 }
452 (void) align_mp (z, &expo, digs);
453 MP_EXPONENT (z) = (MP_T) expo;
454 MP_DIGIT (z, 1) *= sign_x;
455 check_mp_exp (p, z);
456 return z;
457 }
458
459 //! @brief Convert multi-precision number to real.
460
461 REAL_T mp_to_real (NODE_T * p, MP_T * z, int digs)
462 {
463 // This routine looks a lot like "strtod".
464 (void) p;
465 if (PLUS_INF_MP (z)) {
466 return A68G_POSINF_REAL;
467 }
468 if (MINUS_INF_MP (z)) {
469 return A68G_MININF_REAL;
470 }
471 if (MP_EXPONENT (z) * (MP_T) LOG_MP_RADIX <= (MP_T) A68G_REAL_MIN_EXP) {
472 return 0;
473 } else {
474 REAL_T terms[1 + MP_MAX_DIGITS];
475 REAL_T weight = ten_up ((int) (MP_EXPONENT (z) * LOG_MP_RADIX));
476 int lim = MIN (digs, MP_MAX_DIGITS);
477 for (int k = 1; k <= lim; k++) {
478 terms[k] = ABS (MP_DIGIT (z, k)) * weight;
479 weight /= MP_RADIX;
480 }
481 // Sum terms from small to large.
482 REAL_T sum = 0;
483 for (int k = lim; k >= 1; k--) {
484 sum += terms[k];
485 }
486 CHECK_REAL (p, sum);
487 return MP_DIGIT (z, 1) >= 0 ? sum : -sum;
488 }
489 }
490
491 //! @brief Normalise positive intermediate, fast.
492
493 static inline void norm_mp_light (MP_T * w, int k, int digs)
494 {
495 // Bring every digit back to [0 .. MP_RADIX>.
496 MP_T *z = &MP_DIGIT (w, digs);
497 for (int j = digs; j >= k; j--, z--) {
498 if (z[0] >= MP_RADIX) {
499 z[0] -= (MP_T) MP_RADIX;
500 z[-1] += 1;
501 } else if (z[0] < 0) {
502 z[0] += (MP_T) MP_RADIX;
503 z[-1] -= 1;
504 }
505 }
506 }
507
508 //! @brief Normalise positive intermediate.
509
510 static inline void norm_mp (MP_T * w, int k, int digs)
511 {
512 // Bring every digit back to [0 .. MP_RADIX>.
513 MP_T *z = &MP_DIGIT (w, digs);
514 for (int j = digs; j >= k; j--, z--) {
515 if (z[0] >= (MP_T) MP_RADIX) {
516 MP_T carry = (MP_T) ((MP_INT_T) (z[0] / (MP_T) MP_RADIX));
517 z[0] -= carry * (MP_T) MP_RADIX;
518 z[-1] += carry;
519 } else if (z[0] < (MP_T) 0) {
520 MP_T carry = (MP_T) 1 + (MP_T) ((MP_INT_T) ((-z[0] - 1) / (MP_T) MP_RADIX));
521 z[0] += carry * (MP_T) MP_RADIX;
522 z[-1] -= carry;
523 }
524 }
525 }
526
527 //! @brief Round multi-precision number.
528
529 static inline void round_internal_mp (MP_T * z, MP_T * w, int digs)
530 {
531 // Assume that w has precision of at least 2 + digs.
532 int last = (MP_DIGIT (w, 1) == 0 ? 2 + digs : 1 + digs);
533 if (MP_DIGIT (w, last) >= MP_RADIX / 2) {
534 MP_DIGIT (w, last - 1) += 1;
535 }
536 if (MP_DIGIT (w, last - 1) >= MP_RADIX) {
537 norm_mp (w, 2, last); // Hardly ever happens - no need to optimise
538 }
539 if (MP_DIGIT (w, 1) == 0) {
540 (void) move_mp_part (&MP_DIGIT (z, 1), &MP_DIGIT (w, 2), digs);
541 MP_EXPONENT (z) = MP_EXPONENT (w) - 1;
542 } else if (z != w) {
543 (void) move_mp_part (&MP_EXPONENT (z), &MP_EXPONENT (w), (1 + digs));
544 }
545 // Zero is zero is zero.
546 if (MP_DIGIT (z, 1) == 0) {
547 MP_EXPONENT (z) = (MP_T) 0;
548 }
549 }
550
551 //! @brief Truncate at decimal point.
552
553 MP_T *trunc_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
554 {
555 if (MP_EXPONENT (x) < 0) {
556 SET_MP_ZERO (z, digs);
557 } else if (MP_EXPONENT (x) >= (MP_T) digs) {
558 errno = EDOM;
559 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, (IS (MOID (p), PROC_SYMBOL) ? SUB_MOID (p) : MOID (p)));
560 exit_genie (p, A68G_RUNTIME_ERROR);
561 } else {
562 (void) move_mp (z, x, digs);
563 for (int k = (int) (MP_EXPONENT (x) + 2); k <= digs; k++) {
564 MP_DIGIT (z, k) = (MP_T) 0;
565 }
566 }
567 return z;
568 }
569
570 //! @brief Floor - largest integer smaller than x.
571
572 MP_T *floor_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
573 {
574 if (MP_EXPONENT (x) < 0) {
575 SET_MP_ZERO (z, digs);
576 } else if (MP_EXPONENT (x) >= (MP_T) digs) {
577 errno = EDOM;
578 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, (IS (MOID (p), PROC_SYMBOL) ? SUB_MOID (p) : MOID (p)));
579 exit_genie (p, A68G_RUNTIME_ERROR);
580 } else {
581 (void) move_mp (z, x, digs);
582 for (int k = (int) (MP_EXPONENT (x) + 2); k <= digs; k++) {
583 MP_DIGIT (z, k) = (MP_T) 0;
584 }
585 }
586 if (MP_DIGIT (x, 1) < 0 && ! same_mp (p, z, x, digs)) {
587 (void) minus_one_mp (p, z, z, digs);
588 }
589 return z;
590 }
591
592 //! @brief Shorten and round.
593
594 MP_T *shorten_mp (NODE_T * p, MP_T * z, int digs, MP_T * x, int digs_x)
595 {
596 if (digs > digs_x) {
597 return lengthen_mp (p, z, digs, x, digs_x);
598 } else if (digs == digs_x) {
599 return move_mp (z, x, digs);
600 } else {
601 // Reserve extra digs for proper rounding.
602 ADDR_T pop_sp = A68G_SP;
603 int digs_h = digs + 2;
604 BOOL_T negative = (BOOL_T) (MP_DIGIT (x, 1) < 0);
605 MP_T *w = nil_mp (p, digs_h);
606 if (negative) {
607 MP_DIGIT (x, 1) = -MP_DIGIT (x, 1);
608 }
609 MP_STATUS (w) = (MP_T) 0;
610 MP_EXPONENT (w) = MP_EXPONENT (x) + 1;
611 MP_DIGIT (w, 1) = (MP_T) 0;
612 (void) move_mp_part (&MP_DIGIT (w, 2), &MP_DIGIT (x, 1), digs + 1);
613 round_internal_mp (z, w, digs);
614 if (negative) {
615 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
616 }
617 MP_STATUS (z) = MP_STATUS (x);
618 A68G_SP = pop_sp;
619 return z;
620 }
621 }
622
623 //! @brief Lengthen x and assign to z.
624
625 MP_T *lengthen_mp (NODE_T * p, MP_T * z, int digs_z, MP_T * x, int digs_x)
626 {
627 if (digs_z < digs_x) {
628 return shorten_mp (p, z, digs_z, x, digs_x);
629 } else if (digs_z == digs_x) {
630 return move_mp (z, x, digs_z);
631 } else {
632 if (z != x) {
633 (void) move_mp_part (&MP_DIGIT (z, 1), &MP_DIGIT (x, 1), digs_x);
634 MP_EXPONENT (z) = MP_EXPONENT (x);
635 MP_STATUS (z) = MP_STATUS (x);
636 }
637 for (int j = 1 + digs_x; j <= digs_z; j++) {
638 MP_DIGIT (z, j) = (MP_T) 0;
639 }
640 }
641 return z;
642 }
643
644 //! @brief Set "z" to the sum of positive "x" and positive "y".
645
646 MP_T *add_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
647 {
648 MP_STATUS (z) = (MP_T) INIT_MASK;
649 // Trivial cases.
650 if (MP_DIGIT (x, 1) == (MP_T) 0) {
651 (void) move_mp (z, y, digs);
652 return z;
653 } else if (MP_DIGIT (y, 1) == 0) {
654 (void) move_mp (z, x, digs);
655 return z;
656 }
657 // We want positive arguments.
658 ADDR_T pop_sp = A68G_SP;
659 MP_T x_1 = MP_DIGIT (x, 1), y_1 = MP_DIGIT (y, 1);
660 MP_DIGIT (x, 1) = ABS (x_1);
661 MP_DIGIT (y, 1) = ABS (y_1);
662 if (x_1 >= 0 && y_1 < 0) {
663 (void) sub_mp (p, z, x, y, digs);
664 } else if (x_1 < 0 && y_1 >= 0) {
665 (void) sub_mp (p, z, y, x, digs);
666 } else if (x_1 < 0 && y_1 < 0) {
667 (void) add_mp (p, z, x, y, digs);
668 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
669 } else {
670 // Add.
671 int digs_h = 2 + digs;
672 MP_T *w = nil_mp (p, digs_h);
673 if (MP_EXPONENT (x) == MP_EXPONENT (y)) {
674 MP_EXPONENT (w) = (MP_T) 1 + MP_EXPONENT (x);
675 for (int j = 1; j <= digs; j++) {
676 MP_DIGIT (w, j + 1) = MP_DIGIT (x, j) + MP_DIGIT (y, j);
677 }
678 MP_DIGIT (w, digs_h) = (MP_T) 0;
679 } else if (MP_EXPONENT (x) > MP_EXPONENT (y)) {
680 int shl_y = (int) MP_EXPONENT (x) - (int) MP_EXPONENT (y);
681 MP_EXPONENT (w) = (MP_T) 1 + MP_EXPONENT (x);
682 for (int j = 1; j < digs_h; j++) {
683 int i_y = j - shl_y;
684 MP_T x_j = (j > digs ? 0 : MP_DIGIT (x, j));
685 MP_T y_j = (i_y <= 0 || i_y > digs ? 0 : MP_DIGIT (y, i_y));
686 MP_DIGIT (w, j + 1) = x_j + y_j;
687 }
688 } else {
689 int shl_x = (int) MP_EXPONENT (y) - (int) MP_EXPONENT (x);
690 MP_EXPONENT (w) = (MP_T) 1 + MP_EXPONENT (y);
691 for (int j = 1; j < digs_h; j++) {
692 int i_x = j - shl_x;
693 MP_T x_j = (i_x <= 0 || i_x > digs ? 0 : MP_DIGIT (x, i_x));
694 MP_T y_j = (j > digs ? 0 : MP_DIGIT (y, j));
695 MP_DIGIT (w, j + 1) = x_j + y_j;
696 }
697 }
698 norm_mp_light (w, 2, digs_h);
699 round_internal_mp (z, w, digs);
700 check_mp_exp (p, z);
701 }
702 // Restore and exit.
703 A68G_SP = pop_sp;
704 MP_T z_1 = MP_DIGIT (z, 1);
705 MP_DIGIT (x, 1) = x_1;
706 MP_DIGIT (y, 1) = y_1;
707 MP_DIGIT (z, 1) = z_1; // In case z IS x OR z IS y
708 return z;
709 }
710
711 //! @brief Set "z" to the difference of positive "x" and positive "y".
712
713 MP_T *sub_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
714 {
715 MP_STATUS (z) = (MP_T) INIT_MASK;
716 // Trivial cases.
717 if (MP_DIGIT (x, 1) == (MP_T) 0) {
718 (void) move_mp (z, y, digs);
719 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
720 return z;
721 } else if (MP_DIGIT (y, 1) == (MP_T) 0) {
722 (void) move_mp (z, x, digs);
723 return z;
724 }
725 // We want positive arguments.
726 ADDR_T pop_sp = A68G_SP;
727 MP_T x_1 = MP_DIGIT (x, 1), y_1 = MP_DIGIT (y, 1);
728 MP_DIGIT (x, 1) = ABS (x_1);
729 MP_DIGIT (y, 1) = ABS (y_1);
730 if (x_1 >= 0 && y_1 < 0) {
731 (void) add_mp (p, z, x, y, digs);
732 } else if (x_1 < 0 && y_1 >= 0) {
733 (void) add_mp (p, z, y, x, digs);
734 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
735 } else if (x_1 < 0 && y_1 < 0) {
736 (void) sub_mp (p, z, y, x, digs);
737 } else {
738 // Subtract.
739 BOOL_T negative = A68G_FALSE;
740 int fnz, digs_h = 2 + digs;
741 MP_T *w = nil_mp (p, digs_h);
742 if (MP_EXPONENT (x) == MP_EXPONENT (y)) {
743 MP_EXPONENT (w) = (MP_T) 1 + MP_EXPONENT (x);
744 for (int j = 1; j <= digs; j++) {
745 MP_DIGIT (w, j + 1) = MP_DIGIT (x, j) - MP_DIGIT (y, j);
746 }
747 MP_DIGIT (w, digs_h) = (MP_T) 0;
748 } else if (MP_EXPONENT (x) > MP_EXPONENT (y)) {
749 int shl_y = (int) MP_EXPONENT (x) - (int) MP_EXPONENT (y);
750 MP_EXPONENT (w) = (MP_T) 1 + MP_EXPONENT (x);
751 for (int j = 1; j < digs_h; j++) {
752 int i_y = j - shl_y;
753 MP_T x_j = (j > digs ? 0 : MP_DIGIT (x, j));
754 MP_T y_j = (i_y <= 0 || i_y > digs ? 0 : MP_DIGIT (y, i_y));
755 MP_DIGIT (w, j + 1) = x_j - y_j;
756 }
757 } else {
758 int shl_x = (int) MP_EXPONENT (y) - (int) MP_EXPONENT (x);
759 MP_EXPONENT (w) = (MP_T) 1 + MP_EXPONENT (y);
760 for (int j = 1; j < digs_h; j++) {
761 int i_x = j - shl_x;
762 MP_T x_j = (i_x <= 0 || i_x > digs ? 0 : MP_DIGIT (x, i_x));
763 MP_T y_j = (j > digs ? 0 : MP_DIGIT (y, j));
764 MP_DIGIT (w, j + 1) = x_j - y_j;
765 }
766 }
767 // Correct if we subtract large from small.
768 if (MP_DIGIT (w, 2) <= 0) {
769 fnz = -1;
770 for (int j = 2; j <= digs_h && fnz < 0; j++) {
771 if (MP_DIGIT (w, j) != 0) {
772 fnz = j;
773 }
774 }
775 negative = (BOOL_T) (MP_DIGIT (w, fnz) < 0);
776 if (negative) {
777 for (int j = fnz; j <= digs_h; j++) {
778 MP_DIGIT (w, j) = -MP_DIGIT (w, j);
779 }
780 }
781 }
782 // Normalise.
783 norm_mp_light (w, 2, digs_h);
784 fnz = -1;
785 for (int j = 1; j <= digs_h && fnz < 0; j++) {
786 if (MP_DIGIT (w, j) != 0) {
787 fnz = j;
788 }
789 }
790 if (fnz > 1) {
791 int j2 = fnz - 1;
792 for (int j = 1; j <= digs_h - j2; j++) {
793 MP_DIGIT (w, j) = MP_DIGIT (w, j + j2);
794 MP_DIGIT (w, j + j2) = (MP_T) 0;
795 }
796 MP_EXPONENT (w) -= j2;
797 }
798 // Round.
799 round_internal_mp (z, w, digs);
800 if (negative) {
801 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
802 }
803 check_mp_exp (p, z);
804 }
805 // Restore and exit.
806 A68G_SP = pop_sp;
807 MP_T z_1 = MP_DIGIT (z, 1);
808 MP_DIGIT (x, 1) = x_1;
809 MP_DIGIT (y, 1) = y_1;
810 MP_DIGIT (z, 1) = z_1; // In case z IS x OR z IS y
811 return z;
812 }
813
814 //! @brief Set "z" to the product of "x" and "y".
815
816 MP_T *mul_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
817 {
818 if (IS_ZERO_MP (x) || IS_ZERO_MP (y)) {
819 SET_MP_ZERO (z, digs);
820 return z;
821 }
822 // Grammar school algorithm with intermittent normalisation.
823 ADDR_T pop_sp = A68G_SP;
824 int digs_h = 2 + digs;
825 MP_T x_1 = MP_DIGIT (x, 1), y_1 = MP_DIGIT (y, 1);
826 MP_DIGIT (x, 1) = ABS (x_1);
827 MP_DIGIT (y, 1) = ABS (y_1);
828 MP_STATUS (z) = (MP_T) INIT_MASK;
829 MP_T *w = lit_mp (p, 0, MP_EXPONENT (x) + MP_EXPONENT (y) + 1, digs_h);
830 int oflow = (int) FLOOR_MP ((MP_REAL_T) MAX_REPR_INT / (2 * MP_REAL_RADIX * MP_REAL_RADIX)) - 1;
831 for (int i = digs; i >= 1; i--) {
832 MP_T yi = MP_DIGIT (y, i);
833 if (yi != 0) {
834 int k = digs_h - i;
835 int j = (k > digs ? digs : k);
836 MP_T *u = &MP_DIGIT (w, i + j), *v = &MP_DIGIT (x, j);
837 if ((digs - i + 1) % oflow == 0) {
838 norm_mp (w, 2, digs_h);
839 }
840 while (j-- >= 1) {
841 (u--)[0] += yi * (v--)[0];
842 }
843 }
844 }
845 norm_mp (w, 2, digs_h);
846 round_internal_mp (z, w, digs);
847 // Restore and exit.
848 A68G_SP = pop_sp;
849 MP_T z_1 = MP_DIGIT (z, 1);
850 MP_DIGIT (x, 1) = x_1;
851 MP_DIGIT (y, 1) = y_1;
852 MP_DIGIT (z, 1) = ((x_1 * y_1) >= 0 ? z_1 : -z_1);
853 check_mp_exp (p, z);
854 return z;
855 }
856
857 //! @brief Set "z" to the quotient of "x" and "y".
858
859 MP_T *div_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
860 {
861 // This routine is based on
862 //
863 // D. M. Smith, "A Multiple-Precision Division Algorithm"
864 // Mathematics of Computation 66 (1996) 157-163.
865 //
866 // This is O(N^2) but runs faster than straightforward methods by skipping
867 // most of the intermediate normalisation and recovering from wrong
868 // guesses without separate correction steps.
869 // Depending on application, div_mp cost is circa 3 times that of mul_mp.
870 // Therefore Newton-Raphson division makes no sense here.
871 if (IS_ZERO_MP (y)) {
872 errno = ERANGE;
873 return NaN_MP;
874 }
875 // Determine normalisation interval assuming that q < 2b in each step.
876 #if (A68G_LEVEL <= 2)
877 int oflow = (int) FLOOR_MP ((MP_REAL_T) MAX_REPR_INT / (3 * MP_REAL_RADIX * MP_REAL_RADIX)) - 1;
878 #else
879 int oflow = (int) FLOOR_MP ((MP_REAL_T) MAX_REPR_INT / (2 * MP_REAL_RADIX * MP_REAL_RADIX)) - 1;
880 #endif
881 MP_T x_1 = MP_DIGIT (x, 1), y_1 = MP_DIGIT (y, 1);
882 MP_DIGIT (x, 1) = ABS (x_1);
883 MP_DIGIT (y, 1) = ABS (y_1);
884 MP_STATUS (z) = (MP_T) INIT_MASK;
885 // Slight optimisation when the denominator has few digits.
886 int nzdigs = digs;
887 while (MP_DIGIT (y, nzdigs) == 0 && nzdigs > 1) {
888 nzdigs--;
889 }
890 if (nzdigs == 1 && MP_EXPONENT (y) == 0) {
891 (void) div_mp_digit (p, z, x, MP_DIGIT (y, 1), digs);
892 MP_T z_1 = MP_DIGIT (z, 1);
893 MP_DIGIT (x, 1) = x_1;
894 MP_DIGIT (y, 1) = y_1;
895 MP_DIGIT (z, 1) = ((x_1 * y_1) >= 0 ? z_1 : -z_1);
896 check_mp_exp (p, z);
897 return z;
898 }
899 // Working nominator in which the quotient develops.
900 ADDR_T pop_sp = A68G_SP;
901 int wdigs = 4 + digs;
902 MP_T *w = lit_mp (p, 0, MP_EXPONENT (x) - MP_EXPONENT (y), wdigs);
903 (void) move_mp_part (&MP_DIGIT (w, 2), &MP_DIGIT (x, 1), digs);
904 // Estimate the denominator. For small MP_RADIX add: MP_DIGIT (y, 4) / MP_REAL_RADIX.
905 MP_REAL_T den = (MP_DIGIT (y, 1) * MP_REAL_RADIX + MP_DIGIT (y, 2)) * MP_REAL_RADIX + MP_DIGIT (y, 3);
906 MP_T *t = &MP_DIGIT (w, 2);
907 for (int k = 1, len = digs + 2, first = 3; k <= digs + 2; k++, len++, first++, t++) {
908 // Estimate quotient digit.
909 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);
910 if (nom == 0) {
911 q = 0;
912 } else {
913 // Correct the nominator.
914 q = (MP_T) (MP_INT_T) (nom / den);
915 int lim = MINIMUM (len, wdigs);
916 if (nzdigs <= lim - first + 1) {
917 lim = first + nzdigs - 1;
918 }
919 MP_T *u = t, *v = &MP_DIGIT (y, 1);
920 for (int j = first; j <= lim; j++) {
921 (u++)[0] -= q * (v++)[0];
922 }
923 }
924 t[0] += t[-1] * MP_RADIX;
925 t[-1] = q;
926 if (k % oflow == 0 || k == digs + 2) {
927 norm_mp (w, first, wdigs);
928 }
929 }
930 norm_mp (w, 2, digs);
931 round_internal_mp (z, w, digs);
932 // Restore and exit.
933 A68G_SP = pop_sp;
934 MP_T z_1 = MP_DIGIT (z, 1);
935 MP_DIGIT (x, 1) = x_1;
936 MP_DIGIT (y, 1) = y_1;
937 MP_DIGIT (z, 1) = ((x_1 * y_1) >= 0 ? z_1 : -z_1);
938 check_mp_exp (p, z);
939 return z;
940 }
941
942 //! @brief Set "z" to the integer quotient of "x" and "y".
943
944 MP_T *over_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
945 {
946 if (MP_DIGIT (y, 1) == 0) {
947 errno = ERANGE;
948 return NaN_MP;
949 }
950 ADDR_T pop_sp = A68G_SP;
951 int digs_g = FUN_DIGITS (digs);
952 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);
953 (void) div_mp (p, z_g, x_g, y_g, digs_g);
954 trunc_mp (p, z_g, z_g, digs_g);
955 (void) shorten_mp (p, z, digs, z_g, digs_g);
956 MP_STATUS (z) = (MP_T) INIT_MASK;
957 A68G_SP = pop_sp;
958 return z;
959 }
960
961 //! @brief Set "z" to x mod y.
962
963 MP_T *mod_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
964 {
965 if (MP_DIGIT (y, 1) == 0) {
966 errno = EDOM;
967 return NaN_MP;
968 }
969 int digs_g = FUN_DIGITS (digs);
970 ADDR_T pop_sp = A68G_SP;
971 MP_T *x_g = len_mp (p, x, digs, digs_g);
972 MP_T *y_g = len_mp (p, y, digs, digs_g);
973 MP_T *z_g = nil_mp (p, digs_g);
974 // x mod y = x - y * trunc (x / y).
975 (void) over_mp (p, z_g, x_g, y_g, digs_g);
976 (void) mul_mp (p, z_g, y_g, z_g, digs_g);
977 (void) sub_mp (p, z_g, x_g, z_g, digs_g);
978 (void) shorten_mp (p, z, digs, z_g, digs_g);
979 // Restore and exit.
980 A68G_SP = pop_sp;
981 return z;
982 }
983
984 //! @brief Set "z" to the product of x and digit y.
985
986 MP_T *mul_mp_digit (NODE_T * p, MP_T * z, MP_T * x, MP_T y, int digs)
987 {
988 // This is an O(N) routine for multiplication by a short value.
989 MP_T x_1 = MP_DIGIT (x, 1);
990 int digs_h = 2 + digs;
991 ADDR_T pop_sp = A68G_SP;
992 MP_DIGIT (x, 1) = ABS (x_1);
993 MP_STATUS (z) = (MP_T) INIT_MASK;
994 MP_T y_1 = y;
995 y = ABS (y_1);
996 if (y == 2) {
997 (void) add_mp (p, z, x, x, digs);
998 } else {
999 MP_T *w = lit_mp (p, 0, MP_EXPONENT (x) + 1, digs_h);
1000 MP_T *u = &MP_DIGIT (w, 1 + digs), *v = &MP_DIGIT (x, digs);
1001 int j = digs;
1002 while (j-- >= 1) {
1003 (u--)[0] += y * (v--)[0];
1004 }
1005 norm_mp (w, 2, digs_h);
1006 round_internal_mp (z, w, digs);
1007 }
1008 // Restore and exit.
1009 A68G_SP = pop_sp;
1010 MP_T z_1 = MP_DIGIT (z, 1);
1011 MP_DIGIT (x, 1) = x_1;
1012 MP_DIGIT (z, 1) = ((x_1 * y_1) >= 0 ? z_1 : -z_1);
1013 check_mp_exp (p, z);
1014 return z;
1015 }
1016
1017 //! @brief Set "z" to x/2.
1018
1019 MP_T *half_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1020 {
1021 ADDR_T pop_sp = A68G_SP;
1022 int digs_h = 2 + digs;
1023 MP_T x_1 = MP_DIGIT (x, 1);
1024 MP_DIGIT (x, 1) = ABS (x_1);
1025 MP_STATUS (z) = (MP_T) INIT_MASK;
1026 // Calculate x * 0.5.
1027 MP_T *w = lit_mp (p, 0, MP_EXPONENT (x), digs_h);
1028 MP_T *u = &MP_DIGIT (w, 1 + digs), *v = &MP_DIGIT (x, digs);
1029 int j = digs;
1030 while (j-- >= 1) {
1031 (u--)[0] += (MP_RADIX / 2) * (v--)[0];
1032 }
1033 norm_mp (w, 2, digs_h);
1034 round_internal_mp (z, w, digs);
1035 // Restore and exit.
1036 A68G_SP = pop_sp;
1037 MP_T z_1 = MP_DIGIT (z, 1);
1038 MP_DIGIT (x, 1) = x_1;
1039 MP_DIGIT (z, 1) = (x_1 >= 0 ? z_1 : -z_1);
1040 check_mp_exp (p, z);
1041 return z;
1042 }
1043
1044 //! @brief Set "z" to x/10.
1045
1046 MP_T *tenth_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1047 {
1048 ADDR_T pop_sp = A68G_SP;
1049 int digs_h = 2 + digs;
1050 MP_T x_1 = MP_DIGIT (x, 1);
1051 MP_DIGIT (x, 1) = ABS (x_1);
1052 MP_STATUS (z) = (MP_T) INIT_MASK;
1053 // Calculate x * 0.1.
1054 MP_T *w = lit_mp (p, 0, MP_EXPONENT (x), digs_h);
1055 MP_T *u = &MP_DIGIT (w, 1 + digs), *v = &MP_DIGIT (x, digs);
1056 int j = digs;
1057 while (j-- >= 1) {
1058 (u--)[0] += (MP_RADIX / 10) * (v--)[0];
1059 }
1060 norm_mp (w, 2, digs_h);
1061 round_internal_mp (z, w, digs);
1062 // Restore and exit.
1063 A68G_SP = pop_sp;
1064 MP_T z_1 = MP_DIGIT (z, 1);
1065 MP_DIGIT (x, 1) = x_1;
1066 MP_DIGIT (z, 1) = (x_1 >= 0 ? z_1 : -z_1);
1067 check_mp_exp (p, z);
1068 return z;
1069 }
1070
1071 //! @brief Set "z" to the quotient of x and digit y.
1072
1073 MP_T *div_mp_digit (NODE_T * p, MP_T * z, MP_T * x, MP_T y, int digs)
1074 {
1075 if (y == 0) {
1076 errno = ERANGE;
1077 return NaN_MP;
1078 }
1079 // Determine normalisation interval assuming that q < 2b in each step.
1080 #if (A68G_LEVEL <= 2)
1081 int oflow = (int) FLOOR_MP ((MP_REAL_T) MAX_REPR_INT / (3 * MP_REAL_RADIX * MP_REAL_RADIX)) - 1;
1082 #else
1083 int oflow = (int) FLOOR_MP ((MP_REAL_T) MAX_REPR_INT / (2 * MP_REAL_RADIX * MP_REAL_RADIX)) - 1;
1084 #endif
1085 // Work with positive operands.
1086 ADDR_T pop_sp = A68G_SP;
1087 MP_T x_1 = MP_DIGIT (x, 1), y_1 = y;
1088 MP_DIGIT (x, 1) = ABS (x_1);
1089 MP_STATUS (z) = (MP_T) INIT_MASK;
1090 y = ABS (y_1);
1091 if (y == 2) {
1092 (void) half_mp (p, z, x, digs);
1093 } else if (y == 10) {
1094 (void) tenth_mp (p, z, x, digs);
1095 } else {
1096 int wdigs = 4 + digs;
1097 MP_T *w = lit_mp (p, 0, MP_EXPONENT (x), wdigs);
1098 (void) move_mp_part (&MP_DIGIT (w, 2), &MP_DIGIT (x, 1), digs);
1099 // Estimate the denominator.
1100 MP_REAL_T den = (MP_REAL_T) y * MP_REAL_RADIX * MP_REAL_RADIX;
1101 MP_T *t = &MP_DIGIT (w, 2);
1102 for (int k = 1, first = 3; k <= digs + 2; k++, first++, t++) {
1103 // Estimate quotient digit and correct.
1104 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);
1105 MP_REAL_T q = (MP_T) (MP_INT_T) (nom / den);
1106 t[0] += t[-1] * MP_RADIX - q * y;
1107 t[-1] = q;
1108 if (k % oflow == 0 || k == digs + 2) {
1109 norm_mp (w, first, wdigs);
1110 }
1111 }
1112 norm_mp (w, 2, digs);
1113 round_internal_mp (z, w, digs);
1114 }
1115 // Restore and exit.
1116 A68G_SP = pop_sp;
1117 MP_T z_1 = MP_DIGIT (z, 1);
1118 MP_DIGIT (x, 1) = x_1;
1119 MP_DIGIT (z, 1) = ((x_1 * y_1) >= 0 ? z_1 : -z_1);
1120 check_mp_exp (p, z);
1121 return z;
1122 }
1123
1124 //! @brief Set "z" to the integer quotient of "x" and "y".
1125
1126 MP_T *over_mp_digit (NODE_T * p, MP_T * z, MP_T * x, MP_T y, int digs)
1127 {
1128 if (y == 0) {
1129 errno = ERANGE;
1130 return NaN_MP;
1131 }
1132 int digs_g = FUN_DIGITS (digs);
1133 ADDR_T pop_sp = A68G_SP;
1134 MP_T *x_g = len_mp (p, x, digs, digs_g), *z_g = nil_mp (p, digs_g);
1135 (void) div_mp_digit (p, z_g, x_g, y, digs_g);
1136 trunc_mp (p, z_g, z_g, digs_g);
1137 (void) shorten_mp (p, z, digs, z_g, digs_g);
1138 // Restore and exit.
1139 A68G_SP = pop_sp;
1140 return z;
1141 }
1142
1143 //! @brief Set "z" to the reciprocal of "x".
1144
1145 MP_T *rec_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1146 {
1147 if (IS_ZERO_MP (x)) {
1148 errno = ERANGE;
1149 return NaN_MP;
1150 }
1151 ADDR_T pop_sp = A68G_SP;
1152 (void) div_mp (p, z, mp_one (digs), x, digs);
1153 A68G_SP = pop_sp;
1154 return z;
1155 }
1156
1157 //! @brief LONG REAL long pi
1158
1159 void genie_pi_mp (NODE_T * p)
1160 {
1161 int digs = DIGITS (MOID (p));
1162 MP_T *z = nil_mp (p, digs);
1163 (void) mp_pi (p, z, MP_PI, digs);
1164 MP_STATUS (z) = (MP_T) INIT_MASK;
1165 }
1166
1167 //! @brief Set "z" to "x" ** "n".
1168
1169 MP_T *pow_mp_int (NODE_T * p, MP_T * z, MP_T * x, INT_T n, int digs)
1170 {
1171 ADDR_T pop_sp = A68G_SP;
1172 int bit, digs_g = FUN_DIGITS (digs);
1173 MP_T *x_g = len_mp (p, x, digs, digs_g), *z_g = lit_mp (p, 1, 0, digs_g);
1174 BOOL_T negative = (BOOL_T) (n < 0);
1175 if (negative) {
1176 n = -n;
1177 }
1178 bit = 1;
1179 while ((int) bit <= (int) n) {
1180 if (n & bit) {
1181 (void) mul_mp (p, z_g, z_g, x_g, digs_g);
1182 }
1183 (void) mul_mp (p, x_g, x_g, x_g, digs_g);
1184 bit <<= 1;
1185 }
1186 (void) shorten_mp (p, z, digs, z_g, digs_g);
1187 A68G_SP = pop_sp;
1188 if (negative) {
1189 (void) rec_mp (p, z, z, digs);
1190 }
1191 check_mp_exp (p, z);
1192 return z;
1193 }
1194
1195 //! @brief Set "z" to "x" ** "y".
1196
1197 MP_T *pow_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
1198 {
1199 PRELUDE_ERROR (ln_mp (p, z, x, digs) == NaN_MP, p, ERROR_INVALID_ARGUMENT, MOID (p));
1200 (void) mul_mp (p, z, y, z, digs);
1201 (void) exp_mp (p, z, z, digs);
1202 return z;
1203 }
1204
1205 //! @brief Set "z" to 10 ** "n".
1206
1207 MP_T *ten_up_mp (NODE_T * p, MP_T * z, int n, int digs)
1208 {
1209 #if (A68G_LEVEL >= 3)
1210 static MP_T y[LOG_MP_RADIX] = { 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000 };
1211 #else
1212 static MP_T y[LOG_MP_RADIX] = { 1, 10, 100, 1000, 10000, 100000, 1000000 };
1213 #endif
1214 if (n >= 0) {
1215 set_mp (z, y[n % LOG_MP_RADIX], n / LOG_MP_RADIX, digs);
1216 } else {
1217 set_mp (z, y[(LOG_MP_RADIX + n % LOG_MP_RADIX) % LOG_MP_RADIX], (n + 1) / LOG_MP_RADIX - 1, digs);
1218 }
1219 check_mp_exp (p, z);
1220 return z;
1221 }
1222
1223 //! @brief Comparison of "x" and "y".
1224
1225 void eq_mp (NODE_T * p, A68G_BOOL * z, MP_T * x, MP_T * y, int digs)
1226 {
1227 ADDR_T pop_sp = A68G_SP;
1228 MP_T *v = nil_mp (p, digs);
1229 (void) sub_mp (p, v, x, y, digs);
1230 STATUS (z) = INIT_MASK;
1231 VALUE (z) = (MP_DIGIT (v, 1) == 0 ? A68G_TRUE : A68G_FALSE);
1232 A68G_SP = pop_sp;
1233 }
1234
1235 //! @brief Comparison of "x" and "y".
1236
1237 void ne_mp (NODE_T * p, A68G_BOOL * z, MP_T * x, MP_T * y, int digs)
1238 {
1239 ADDR_T pop_sp = A68G_SP;
1240 MP_T *v = nil_mp (p, digs);
1241 (void) sub_mp (p, v, x, y, digs);
1242 STATUS (z) = INIT_MASK;
1243 VALUE (z) = (MP_DIGIT (v, 1) != 0 ? A68G_TRUE : A68G_FALSE);
1244 A68G_SP = pop_sp;
1245 }
1246
1247 //! @brief Comparison of "x" and "y".
1248
1249 void lt_mp (NODE_T * p, A68G_BOOL * z, MP_T * x, MP_T * y, int digs)
1250 {
1251 ADDR_T pop_sp = A68G_SP;
1252 MP_T *v = nil_mp (p, digs);
1253 (void) sub_mp (p, v, x, y, digs);
1254 STATUS (z) = INIT_MASK;
1255 VALUE (z) = (MP_DIGIT (v, 1) < 0 ? A68G_TRUE : A68G_FALSE);
1256 A68G_SP = pop_sp;
1257 }
1258
1259 //! @brief Comparison of "x" and "y".
1260
1261 void le_mp (NODE_T * p, A68G_BOOL * z, MP_T * x, MP_T * y, int digs)
1262 {
1263 ADDR_T pop_sp = A68G_SP;
1264 MP_T *v = nil_mp (p, digs);
1265 (void) sub_mp (p, v, x, y, digs);
1266 STATUS (z) = INIT_MASK;
1267 VALUE (z) = (MP_DIGIT (v, 1) <= 0 ? A68G_TRUE : A68G_FALSE);
1268 A68G_SP = pop_sp;
1269 }
1270
1271 //! @brief Comparison of "x" and "y".
1272
1273 void gt_mp (NODE_T * p, A68G_BOOL * z, MP_T * x, MP_T * y, int digs)
1274 {
1275 ADDR_T pop_sp = A68G_SP;
1276 MP_T *v = nil_mp (p, digs);
1277 (void) sub_mp (p, v, x, y, digs);
1278 STATUS (z) = INIT_MASK;
1279 VALUE (z) = (MP_DIGIT (v, 1) > 0 ? A68G_TRUE : A68G_FALSE);
1280 A68G_SP = pop_sp;
1281 }
1282
1283 //! @brief Comparison of "x" and "y".
1284
1285 void ge_mp (NODE_T * p, A68G_BOOL * z, MP_T * x, MP_T * y, int digs)
1286 {
1287 ADDR_T pop_sp = A68G_SP;
1288 MP_T *v = nil_mp (p, digs);
1289 (void) sub_mp (p, v, x, y, digs);
1290 STATUS (z) = INIT_MASK;
1291 VALUE (z) = (MP_DIGIT (v, 1) >= 0 ? A68G_TRUE : A68G_FALSE);
1292 A68G_SP = pop_sp;
1293 }
1294
1295 //! @brief round (x).
1296
1297 MP_T *round_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1298 {
1299 MP_T *y = nil_mp (p, digs);
1300 SET_MP_HALF (y, digs);
1301 if (MP_DIGIT (x, 1) >= 0) {
1302 (void) add_mp (p, z, x, y, digs);
1303 (void) trunc_mp (p, z, z, digs);
1304 } else {
1305 (void) sub_mp (p, z, x, y, digs);
1306 (void) trunc_mp (p, z, z, digs);
1307 }
1308 MP_STATUS (z) = (MP_T) INIT_MASK;
1309 return z;
1310 }
1311
1312 //! @brief Entier (x).
1313
1314 MP_T *entier_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1315 {
1316 if (MP_DIGIT (x, 1) >= 0) {
1317 (void) trunc_mp (p, z, x, digs);
1318 } else {
1319 MP_T *y = nil_mp (p, digs);
1320 (void) move_mp (y, z, digs);
1321 (void) trunc_mp (p, z, x, digs);
1322 (void) sub_mp (p, y, y, z, digs);
1323 if (MP_DIGIT (y, 1) != 0) {
1324 SET_MP_ONE (y, digs);
1325 (void) sub_mp (p, z, z, y, digs);
1326 }
1327 }
1328 MP_STATUS (z) = (MP_T) INIT_MASK;
1329 return z;
1330 }
|
© 2002-2026 J.M. van der Veer (jmvdveer@xs4all.nl)
|
