Bug Summary

File:mp.c
Warning:line 505, column 30
Out of bound memory access (accessed memory precedes memory block)

Annotated Source Code

Press '?' to see keyboard shortcuts

clang -cc1 -cc1 -triple x86_64-pc-linux-gnu -analyze -disable-free -clear-ast-before-backend -disable-llvm-verifier -discard-value-names -main-file-name mp.c -analyzer-store=region -analyzer-opt-analyze-nested-blocks -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -mrelocation-model pic -pic-level 2 -pic-is-pie -mframe-pointer=all -fmath-errno -ffp-contract=on -fno-rounding-math -mconstructor-aliases -funwind-tables=2 -target-cpu x86-64 -tune-cpu generic -debugger-tuning=gdb -fcoverage-compilation-dir=/rootdir/src -resource-dir /usr/lib/llvm-14/lib/clang/14.0.6 -D HAVE_CONFIG_H -I . -I .. -D VERSION="1.25.0" -D LOCALE_DIR="/usr/local/share/locale" -D GETTEXT_PACKAGE="cafe-calc" -I /usr/include/ctk-3.0 -I /usr/include/pango-1.0 -I /usr/include/glib-2.0 -I /usr/lib/x86_64-linux-gnu/glib-2.0/include -I /usr/include/harfbuzz -I /usr/include/freetype2 -I /usr/include/libpng16 -I /usr/include/libmount -I /usr/include/blkid -I /usr/include/fribidi -I /usr/include/uuid -I /usr/include/cairo -I /usr/include/pixman-1 -I /usr/include/gdk-pixbuf-2.0 -I /usr/include/x86_64-linux-gnu -I /usr/include/gio-unix-2.0 -I /usr/include/atk-1.0 -I /usr/include/at-spi2-atk/2.0 -I /usr/include/at-spi-2.0 -I /usr/include/dbus-1.0 -I /usr/lib/x86_64-linux-gnu/dbus-1.0/include -I /usr/include/libxml2 -internal-isystem /usr/lib/llvm-14/lib/clang/14.0.6/include -internal-isystem /usr/local/include -internal-isystem /usr/lib/gcc/x86_64-linux-gnu/12/../../../../x86_64-linux-gnu/include -internal-externc-isystem /usr/include/x86_64-linux-gnu -internal-externc-isystem /include -internal-externc-isystem /usr/include -fdebug-compilation-dir=/rootdir/src -ferror-limit 19 -fgnuc-version=4.2.1 -analyzer-checker deadcode.DeadStores -analyzer-checker alpha.deadcode.UnreachableCode -analyzer-checker alpha.core.CastSize -analyzer-checker alpha.core.CastToStruct -analyzer-checker alpha.core.IdenticalExpr -analyzer-checker alpha.core.SizeofPtr -analyzer-checker alpha.security.ArrayBoundV2 -analyzer-checker alpha.security.MallocOverflow -analyzer-checker alpha.security.ReturnPtrRange -analyzer-checker alpha.unix.SimpleStream -analyzer-checker alpha.unix.cstring.BufferOverlap -analyzer-checker alpha.unix.cstring.NotNullTerminated -analyzer-checker alpha.unix.cstring.OutOfBounds -analyzer-checker alpha.core.FixedAddr -analyzer-output=html -faddrsig -D__GCC_HAVE_DWARF2_CFI_ASM=1 -o /rootdir/html-report/2022-12-29-234211-15618-1 -x c mp.c
1/*
2 * Copyright (C) 1987-2008 Sun Microsystems, Inc. All Rights Reserved.
3 * Copyright (C) 2008-2011 Robert Ancell
4 *
5 * This program is free software: you can redistribute it and/or modify it under
6 * the terms of the GNU General Public License as published by the Free Software
7 * Foundation, either version 2 of the License, or (at your option) any later
8 * version. See http://www.gnu.org/copyleft/gpl.html the full text of the
9 * license.
10 */
11
12#include <stdlib.h>
13#include <stdio.h>
14#include <math.h>
15#include <errno(*__errno_location ()).h>
16
17#include "mp.h"
18#include "mp-private.h"
19
20static MPNumber eulers_number;
21static gboolean have_eulers_number = FALSE(0);
22
23// FIXME: Re-add overflow and underflow detection
24
25char *mp_error = NULL((void*)0);
26
27/* THIS ROUTINE IS CALLED WHEN AN ERROR CONDITION IS ENCOUNTERED, AND
28 * AFTER A MESSAGE HAS BEEN WRITTEN TO STDERR.
29 */
30void
31mperr(const char *format, ...)
32{
33 char text[1024];
34 va_list args;
35
36 va_start(args, format)__builtin_va_start(args, format);
37 vsnprintf(text, 1024, format, args);
38 va_end(args)__builtin_va_end(args);
39
40 if (mp_error)
41 free(mp_error);
42 mp_error = strdup(text);
43}
44
45
46const char *
47mp_get_error()
48{
49 return mp_error;
50}
51
52
53void mp_clear_error()
54{
55 if (mp_error)
56 free(mp_error);
57 mp_error = NULL((void*)0);
58}
59
60
61/* ROUTINE CALLED BY MP_DIVIDE AND MP_SQRT TO ENSURE THAT
62 * RESULTS ARE REPRESENTED EXACTLY IN T-2 DIGITS IF THEY
63 * CAN BE. X IS AN MP NUMBER, I AND J ARE INTEGERS.
64 */
65static void
66mp_ext(int i, int j, MPNumber *x)
67{
68 int q, s;
69
70 if (mp_is_zero(x) || MP_T100 <= 2 || i == 0)
71 return;
72
73 /* COMPUTE MAXIMUM POSSIBLE ERROR IN THE LAST PLACE */
74 q = (j + 1) / i + 1;
75 s = MP_BASE10000 * x->fraction[MP_T100 - 2] + x->fraction[MP_T100 - 1];
76
77 /* SET LAST TWO DIGITS TO ZERO */
78 if (s <= q) {
79 x->fraction[MP_T100 - 2] = 0;
80 x->fraction[MP_T100 - 1] = 0;
81 return;
82 }
83
84 if (s + q < MP_BASE10000 * MP_BASE10000)
85 return;
86
87 /* ROUND UP HERE */
88 x->fraction[MP_T100 - 2] = MP_BASE10000 - 1;
89 x->fraction[MP_T100 - 1] = MP_BASE10000;
90
91 /* NORMALIZE X (LAST DIGIT B IS OK IN MP_MULTIPLY_INTEGER) */
92 mp_multiply_integer(x, 1, x);
93}
94
95
96void
97mp_get_eulers(MPNumber *z)
98{
99 if (!have_eulers_number) {
100 MPNumber t;
101 mp_set_from_integer(1, &t);
102 mp_epowy(&t, &eulers_number);
103 have_eulers_number = TRUE(!(0));
104 }
105 mp_set_from_mp(&eulers_number, z);
106}
107
108
109void
110mp_get_i(MPNumber *z)
111{
112 mp_set_from_integer(0, z);
113 z->im_sign = 1;
114 z->im_exponent = 1;
115 z->im_fraction[0] = 1;
116}
117
118
119void
120mp_abs(const MPNumber *x, MPNumber *z)
121{
122 if (mp_is_complex(x)){
123 MPNumber x_real, x_im;
124
125 mp_real_component(x, &x_real);
126 mp_imaginary_component(x, &x_im);
127
128 mp_multiply(&x_real, &x_real, &x_real);
129 mp_multiply(&x_im, &x_im, &x_im);
130 mp_add(&x_real, &x_im, z);
131 mp_root(z, 2, z);
132 }
133 else {
134 mp_set_from_mp(x, z);
135 if (z->sign < 0)
136 z->sign = -z->sign;
137 }
138}
139
140
141void
142mp_arg(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
143{
144 MPNumber x_real, x_im, pi;
145
146 if (mp_is_zero(x)) {
147 /* Translators: Error display when attempting to take argument of zero */
148 mperr(_("Argument not defined for zero")gettext ("Argument not defined for zero"));
149 mp_set_from_integer(0, z);
150 return;
151 }
152
153 mp_real_component(x, &x_real);
154 mp_imaginary_component(x, &x_im);
155 mp_get_pi(&pi);
156
157 if (mp_is_zero(&x_im)) {
158 if (mp_is_negative(&x_real))
159 convert_from_radians(&pi, MP_RADIANS, z);
160 else
161 mp_set_from_integer(0, z);
162 }
163 else if (mp_is_zero(&x_real)) {
164 mp_set_from_mp(&pi, z);
165 if (mp_is_negative(&x_im))
166 mp_divide_integer(z, -2, z);
167 else
168 mp_divide_integer(z, 2, z);
169 }
170 else if (mp_is_negative(&x_real)) {
171 mp_divide(&x_im, &x_real, z);
172 mp_atan(z, MP_RADIANS, z);
173 if (mp_is_negative(&x_im))
174 mp_subtract(z, &pi, z);
175 else
176 mp_add(z, &pi, z);
177 }
178 else {
179 mp_divide(&x_im, &x_real, z);
180 mp_atan(z, MP_RADIANS, z);
181 }
182
183 convert_from_radians(z, unit, z);
184}
185
186
187void
188mp_conjugate(const MPNumber *x, MPNumber *z)
189{
190 mp_set_from_mp(x, z);
191 z->im_sign = -z->im_sign;
192}
193
194
195void
196mp_real_component(const MPNumber *x, MPNumber *z)
197{
198 mp_set_from_mp(x, z);
199
200 /* Clear imaginary component */
201 z->im_sign = 0;
202 z->im_exponent = 0;
203 memset(z->im_fraction, 0, sizeof(int) * MP_SIZE1000);
204}
205
206
207void
208mp_imaginary_component(const MPNumber *x, MPNumber *z)
209{
210 /* Copy imaginary component to real component */
211 z->sign = x->im_sign;
212 z->exponent = x->im_exponent;
213 memcpy(z->fraction, x->im_fraction, sizeof(int) * MP_SIZE1000);
214
215 /* Clear (old) imaginary component */
216 z->im_sign = 0;
217 z->im_exponent = 0;
218 memset(z->im_fraction, 0, sizeof(int) * MP_SIZE1000);
219}
220
221
222static void
223mp_add_real(const MPNumber *x, int y_sign, const MPNumber *y, MPNumber *z)
224{
225 int sign_prod, i, c;
226 int exp_diff, med;
227 bool_Bool x_largest = false0;
228 const int *big_fraction, *small_fraction;
229 MPNumber x_copy, y_copy;
230
231 /* 0 + y = y */
232 if (mp_is_zero(x)) {
233 mp_set_from_mp(y, z);
234 z->sign = y_sign;
235 return;
236 }
237 /* x + 0 = x */
238 else if (mp_is_zero(y)) {
239 mp_set_from_mp(x, z);
240 return;
241 }
242
243 sign_prod = y_sign * x->sign;
244 exp_diff = x->exponent - y->exponent;
245 med = abs(exp_diff);
246 if (exp_diff < 0) {
247 x_largest = false0;
248 } else if (exp_diff > 0) {
249 x_largest = true1;
250 } else {
251 /* EXPONENTS EQUAL SO COMPARE SIGNS, THEN FRACTIONS IF NEC. */
252 if (sign_prod < 0) {
253 /* Signs are not equal. find out which mantissa is larger. */
254 int j;
255 for (j = 0; j < MP_T100; j++) {
256 int i = x->fraction[j] - y->fraction[j];
257 if (i == 0)
258 continue;
259 if (i < 0)
260 x_largest = false0;
261 else if (i > 0)
262 x_largest = true1;
263 break;
264 }
265
266 /* Both mantissas equal, so result is zero. */
267 if (j >= MP_T100) {
268 mp_set_from_integer(0, z);
269 return;
270 }
271 }
272 }
273
274 mp_set_from_mp(x, &x_copy);
275 mp_set_from_mp(y, &y_copy);
276 mp_set_from_integer(0, z);
277
278 if (x_largest) {
279 z->sign = x_copy.sign;
280 z->exponent = x_copy.exponent;
281 big_fraction = x_copy.fraction;
282 small_fraction = y_copy.fraction;
283 } else {
284 z->sign = y_sign;
285 z->exponent = y_copy.exponent;
286 big_fraction = y_copy.fraction;
287 small_fraction = x_copy.fraction;
288 }
289
290 /* CLEAR GUARD DIGITS TO RIGHT OF X DIGITS */
291 for(i = 3; i >= med; i--)
292 z->fraction[MP_T100 + i] = 0;
293
294 if (sign_prod >= 0) {
295 /* This is probably insufficient overflow detection, but it makes us
296 * not crash at least.
297 */
298 if (MP_T100 + 3 < med) {
299 mperr(_("Overflow: the result couldn't be calculated")gettext ("Overflow: the result couldn't be calculated"));
300 mp_set_from_integer(0, z);
301 return;
302 }
303
304 /* HERE DO ADDITION, EXPONENT(Y) >= EXPONENT(X) */
305 for (i = MP_T100 + 3; i >= MP_T100; i--)
306 z->fraction[i] = small_fraction[i - med];
307
308 c = 0;
309 for (; i >= med; i--) {
310 c = big_fraction[i] + small_fraction[i - med] + c;
311
312 if (c < MP_BASE10000) {
313 /* NO CARRY GENERATED HERE */
314 z->fraction[i] = c;
315 c = 0;
316 } else {
317 /* CARRY GENERATED HERE */
318 z->fraction[i] = c - MP_BASE10000;
319 c = 1;
320 }
321 }
322
323 for (; i >= 0; i--)
324 {
325 c = big_fraction[i] + c;
326 if (c < MP_BASE10000) {
327 z->fraction[i] = c;
328 i--;
329
330 /* NO CARRY POSSIBLE HERE */
331 for (; i >= 0; i--)
332 z->fraction[i] = big_fraction[i];
333
334 c = 0;
335 break;
336 }
337
338 z->fraction[i] = 0;
339 c = 1;
340 }
341
342 /* MUST SHIFT RIGHT HERE AS CARRY OFF END */
343 if (c != 0) {
344 for (i = MP_T100 + 3; i > 0; i--)
345 z->fraction[i] = z->fraction[i - 1];
346 z->fraction[0] = 1;
347 z->exponent++;
348 }
349 }
350 else {
351 c = 0;
352 for (i = MP_T100 + med - 1; i >= MP_T100; i--) {
353 /* HERE DO SUBTRACTION, ABS(Y) > ABS(X) */
354 z->fraction[i] = c - small_fraction[i - med];
355 c = 0;
356
357 /* BORROW GENERATED HERE */
358 if (z->fraction[i] < 0) {
359 c = -1;
360 z->fraction[i] += MP_BASE10000;
361 }
362 }
363
364 for(; i >= med; i--) {
365 c = big_fraction[i] + c - small_fraction[i - med];
366 if (c >= 0) {
367 /* NO BORROW GENERATED HERE */
368 z->fraction[i] = c;
369 c = 0;
370 } else {
371 /* BORROW GENERATED HERE */
372 z->fraction[i] = c + MP_BASE10000;
373 c = -1;
374 }
375 }
376
377 for (; i >= 0; i--) {
378 c = big_fraction[i] + c;
379
380 if (c >= 0) {
381 z->fraction[i] = c;
382 i--;
383
384 /* NO CARRY POSSIBLE HERE */
385 for (; i >= 0; i--)
386 z->fraction[i] = big_fraction[i];
387
388 break;
389 }
390
391 z->fraction[i] = c + MP_BASE10000;
392 c = -1;
393 }
394 }
395
396 mp_normalize(z);
397}
398
399
400static void
401mp_add_with_sign(const MPNumber *x, int y_sign, const MPNumber *y, MPNumber *z)
402{
403 if (mp_is_complex(x) || mp_is_complex(y)) {
404 MPNumber real_x, real_y, im_x, im_y, real_z, im_z;
405
406 mp_real_component(x, &real_x);
407 mp_imaginary_component(x, &im_x);
408 mp_real_component(y, &real_y);
409 mp_imaginary_component(y, &im_y);
410
411 mp_add_real(&real_x, y_sign * y->sign, &real_y, &real_z);
412 mp_add_real(&im_x, y_sign * y->im_sign, &im_y, &im_z);
413
414 mp_set_from_complex(&real_z, &im_z, z);
415 }
416 else
417 mp_add_real(x, y_sign * y->sign, y, z);
418}
419
420
421void
422mp_add(const MPNumber *x, const MPNumber *y, MPNumber *z)
423{
424 mp_add_with_sign(x, 1, y, z);
425}
426
427
428void
429mp_add_integer(const MPNumber *x, int64_t y, MPNumber *z)
430{
431 MPNumber t;
432 mp_set_from_integer(y, &t);
433 mp_add(x, &t, z);
434}
435
436
437void
438mp_add_fraction(const MPNumber *x, int64_t i, int64_t j, MPNumber *y)
439{
440 MPNumber t;
441 mp_set_from_fraction(i, j, &t);
442 mp_add(x, &t, y);
443}
444
445
446void
447mp_subtract(const MPNumber *x, const MPNumber *y, MPNumber *z)
448{
449 mp_add_with_sign(x, -1, y, z);
450}
451
452
453void
454mp_sgn(const MPNumber *x, MPNumber *z)
455{
456 if (mp_is_zero(x))
457 mp_set_from_integer(0, z);
458 else if (mp_is_negative(x))
459 mp_set_from_integer(-1, z);
460 else
461 mp_set_from_integer(1, z);
462}
463
464void
465mp_integer_component(const MPNumber *x, MPNumber *z)
466{
467 int i;
468
469 /* Clear fraction */
470 mp_set_from_mp(x, z);
471 for (i = z->exponent; i < MP_SIZE1000; i++)
472 z->fraction[i] = 0;
473 z->im_sign = 0;
474 z->im_exponent = 0;
475 memset(z->im_fraction, 0, sizeof(int) * MP_SIZE1000);
476}
477
478void
479mp_fractional_component(const MPNumber *x, MPNumber *z)
480{
481 int i, shift;
482
483 /* Fractional component of zero is 0 */
484 if (mp_is_zero(x)) {
14
←
Assuming the condition is false
→
15
←
Taking false branch
→
485 mp_set_from_integer(0, z);
486 return;
487 }
488
489 /* All fractional */
490 if (x->exponent <= 0) {
16
←
Assuming field 'exponent' is > 0
→
17
←
Taking false branch
→
491 mp_set_from_mp(x, z);
492 return;
493 }
494
495 /* Shift fractional component */
496 shift = x->exponent;
497 for (i = shift; i < MP_SIZE1000 && x->fraction[i] == 0; i++)
18
←
Assuming 'i' is >= MP_SIZE
→
498 shift++;
499 z->sign = x->sign;
500 z->exponent = x->exponent - shift;
501 for (i = 0; i < MP_SIZE1000; i++) {
19
←
Loop condition is true. Entering loop body
→
21
←
Loop condition is true. Entering loop body
→
502 if (i + shift >= MP_SIZE1000)
20
←
Taking true branch
→
22
←
Assuming the condition is false
→
23
←
Taking false branch
→
503 z->fraction[i] = 0;
504 else
505 z->fraction[i] = x->fraction[i + shift];
24
←
Out of bound memory access (accessed memory precedes memory block)
506 }
507 if (z->fraction[0] == 0)
508 z->sign = 0;
509
510 z->im_sign = 0;
511 z->im_exponent = 0;
512 memset(z->im_fraction, 0, sizeof(int) * MP_SIZE1000);
513}
514
515
516void
517mp_fractional_part(const MPNumber *x, MPNumber *z)
518{
519 MPNumber f;
520 mp_floor(x, &f);
521 mp_subtract(x, &f, z);
522}
523
524
525void
526mp_floor(const MPNumber *x, MPNumber *z)
527{
528 int i;
529 bool_Bool have_fraction = false0, is_negative;
530
531 /* Integer component of zero = 0 */
532 if (mp_is_zero(x)) {
533 mp_set_from_mp(x, z);
534 return;
535 }
536
537 /* If all fractional then no integer component */
538 if (x->exponent <= 0) {
539 mp_set_from_integer(0, z);
540 return;
541 }
542
543 is_negative = mp_is_negative(x);
544
545 /* Clear fraction */
546 mp_set_from_mp(x, z);
547 for (i = z->exponent; i < MP_SIZE1000; i++) {
548 if (z->fraction[i])
549 have_fraction = true1;
550 z->fraction[i] = 0;
551 }
552 z->im_sign = 0;
553 z->im_exponent = 0;
554 memset(z->im_fraction, 0, sizeof(int) * MP_SIZE1000);
555
556 if (have_fraction && is_negative)
557 mp_add_integer(z, -1, z);
558}
559
560
561void
562mp_ceiling(const MPNumber *x, MPNumber *z)
563{
564 MPNumber f;
565
566 mp_floor(x, z);
567 mp_fractional_component(x, &f);
568 if (mp_is_zero(&f))
569 return;
570 mp_add_integer(z, 1, z);
571}
572
573
574void
575mp_round(const MPNumber *x, MPNumber *z)
576{
577 MPNumber f, one;
578 bool_Bool do_floor;
579
580 do_floor = !mp_is_negative(x);
581
582 mp_fractional_component(x, &f);
583 mp_multiply_integer(&f, 2, &f);
584 mp_abs(&f, &f);
585 mp_set_from_integer(1, &one);
586 if (mp_is_greater_equal(&f, &one))
587 do_floor = !do_floor;
588
589 if (do_floor)
590 mp_floor(x, z);
591 else
592 mp_ceiling(x, z);
593}
594
595int
596mp_compare_mp_to_mp(const MPNumber *x, const MPNumber *y)
597{
598 int i;
599
600 if (x->sign != y->sign) {
601 if (x->sign > y->sign)
602 return 1;
603 else
604 return -1;
605 }
606
607 /* x = y = 0 */
608 if (mp_is_zero(x))
609 return 0;
610
611 /* See if numbers are of different magnitude */
612 if (x->exponent != y->exponent) {
613 if (x->exponent > y->exponent)
614 return x->sign;
615 else
616 return -x->sign;
617 }
618
619 /* Compare fractions */
620 for (i = 0; i < MP_SIZE1000; i++) {
621 if (x->fraction[i] == y->fraction[i])
622 continue;
623
624 if (x->fraction[i] > y->fraction[i])
625 return x->sign;
626 else
627 return -x->sign;
628 }
629
630 /* x = y */
631 return 0;
632}
633
634
635void
636mp_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
637{
638 int i, ie;
639 MPNumber t;
640
641 /* x/0 */
642 if (mp_is_zero(y)) {
643 /* Translators: Error displayed attempted to divide by zero */
644 mperr(_("Division by zero is undefined")gettext ("Division by zero is undefined"));
645 mp_set_from_integer(0, z);
646 return;
647 }
648
649 /* 0/y = 0 */
650 if (mp_is_zero(x)) {
651 mp_set_from_integer(0, z);
652 return;
653 }
654
655 /* z = x × y⁻¹ */
656 /* FIXME: Set exponent to zero to avoid overflow in mp_multiply??? */
657 mp_reciprocal(y, &t);
658 ie = t.exponent;
659 t.exponent = 0;
660 i = t.fraction[0];
661 mp_multiply(x, &t, z);
662 mp_ext(i, z->fraction[0], z);
663 z->exponent += ie;
664}
665
666
667static void
668mp_divide_integer_real(const MPNumber *x, int64_t y, MPNumber *z)
669{
670 int c, i, k, b2, c2, j1, j2;
671 MPNumber x_copy;
672
673 /* x/0 */
674 if (y == 0) {
675 /* Translators: Error displayed attempted to divide by zero */
676 mperr(_("Division by zero is undefined")gettext ("Division by zero is undefined"));
677 mp_set_from_integer(0, z);
678 return;
679 }
680
681 /* 0/y = 0 */
682 if (mp_is_zero(x)) {
683 mp_set_from_integer(0, z);
684 return;
685 }
686
687 /* Division by -1 or 1 just changes sign */
688 if (y == 1 || y == -1) {
689 if (y < 0)
690 mp_invert_sign(x, z);
691 else
692 mp_set_from_mp(x, z);
693 return;
694 }
695
696 /* Copy x as z may also refer to x */
697 mp_set_from_mp(x, &x_copy);
698 mp_set_from_integer(0, z);
699
700 if (y < 0) {
701 y = -y;
702 z->sign = -x_copy.sign;
703 }
704 else
705 z->sign = x_copy.sign;
706 z->exponent = x_copy.exponent;
707
708 c = 0;
709 i = 0;
710
711 /* IF y*B NOT REPRESENTABLE AS AN INTEGER HAVE TO SIMULATE
712 * LONG DIVISION. ASSUME AT LEAST 16-BIT WORD.
713 */
714
715 /* Computing MAX */
716 b2 = max(MP_BASE << 3, 32767 / MP_BASE)((10000 << 3) >= (32767 / 10000) ? (10000 << 3
) : (32767 / 10000));
717 if (y < b2) {
718 int kh, r1;
719
720 /* LOOK FOR FIRST NONZERO DIGIT IN QUOTIENT */
721 do {
722 c = MP_BASE10000 * c;
723 if (i < MP_T100)
724 c += x_copy.fraction[i];
725 i++;
726 r1 = c / y;
727 if (r1 < 0)
728 goto L210;
729 } while (r1 == 0);
730
731 /* ADJUST EXPONENT AND GET T+4 DIGITS IN QUOTIENT */
732 z->exponent += 1 - i;
733 z->fraction[0] = r1;
734 c = MP_BASE10000 * (c - y * r1);
735 kh = 1;
736 if (i < MP_T100) {
737 kh = MP_T100 + 1 - i;
738 for (k = 1; k < kh; k++) {
739 c += x_copy.fraction[i];
740 z->fraction[k] = c / y;
741 c = MP_BASE10000 * (c - y * z->fraction[k]);
742 i++;
743 }
744 if (c < 0)
745 goto L210;
746 }
747
748 for (k = kh; k < MP_T100 + 4; k++) {
749 z->fraction[k] = c / y;
750 c = MP_BASE10000 * (c - y * z->fraction[k]);
751 }
752 if (c < 0)
753 goto L210;
754
755 mp_normalize(z);
756 return;
757 }
758
759 /* HERE NEED SIMULATED DOUBLE-PRECISION DIVISION */
760 j1 = y / MP_BASE10000;
761 j2 = y - j1 * MP_BASE10000;
762
763 /* LOOK FOR FIRST NONZERO DIGIT */
764 c2 = 0;
765 do {
766 c = MP_BASE10000 * c + c2;
767 c2 = i < MP_T100 ? x_copy.fraction[i] : 0;
768 i++;
769 } while (c < j1 || (c == j1 && c2 < j2));
770
771 /* COMPUTE T+4 QUOTIENT DIGITS */
772 z->exponent += 1 - i;
773 i--;
774
775 /* MAIN LOOP FOR LARGE ABS(y) CASE */
776 for (k = 1; k <= MP_T100 + 4; k++) {
777 int ir, iq, iqj;
778
779 /* GET APPROXICAFE QUOTIENT FIRST */
780 ir = c / (j1 + 1);
781
782 /* NOW REDUCE SO OVERFLOW DOES NOT OCCUR */
783 iq = c - ir * j1;
784 if (iq >= b2) {
785 /* HERE IQ*B WOULD POSSIBLY OVERFLOW SO INCREASE IR */
786 ++ir;
787 iq -= j1;
788 }
789
790 iq = iq * MP_BASE10000 - ir * j2;
791 if (iq < 0) {
792 /* HERE IQ NEGATIVE SO IR WAS TOO LARGE */
793 ir--;
794 iq += y;
795 }
796
797 if (i < MP_T100)
798 iq += x_copy.fraction[i];
799 i++;
800 iqj = iq / y;
801
802 /* R(K) = QUOTIENT, C = REMAINDER */
803 z->fraction[k - 1] = iqj + ir;
804 c = iq - y * iqj;
805
806 if (c < 0)
807 goto L210;
808 }
809
810 mp_normalize(z);
811
812L210:
813 /* CARRY NEGATIVE SO OVERFLOW MUST HAVE OCCURRED */
814 mperr("*** INTEGER OVERFLOW IN MP_DIVIDE_INTEGER, B TOO LARGE ***");
815 mp_set_from_integer(0, z);
816}
817
818
819void
820mp_divide_integer(const MPNumber *x, int64_t y, MPNumber *z)
821{
822 if (mp_is_complex(x)) {
823 MPNumber re_z, im_z;
824
825 mp_real_component(x, &re_z);
826 mp_imaginary_component(x, &im_z);
827 mp_divide_integer_real(&re_z, y, &re_z);
828 mp_divide_integer_real(&im_z, y, &im_z);
829 mp_set_from_complex(&re_z, &im_z, z);
830 }
831 else
832 mp_divide_integer_real(x, y, z);
833}
834
835
836bool_Bool
837mp_is_integer(const MPNumber *x)
838{
839 MPNumber t1, t2, t3;
840
841 if (mp_is_complex(x))
842 return false0;
843
844 /* This fix is required for 1/3 repiprocal not being detected as an integer */
845 /* Multiplication and division by 10000 is used to get around a
846 * limitation to the "fix" for Sun bugtraq bug #4006391 in the
847 * mp_floor() routine in mp.c, when the exponent is less than 1.
848 */
849 mp_set_from_integer(10000, &t3);
850 mp_multiply(x, &t3, &t1);
851 mp_divide(&t1, &t3, &t1);
852 mp_floor(&t1, &t2);
853 return mp_is_equal(&t1, &t2);
854
855 /* Correct way to check for integer */
856 /*int i;
857
858 // Zero is an integer
859 if (mp_is_zero(x))
860 return true;
861
862 // Fractional
863 if (x->exponent <= 0)
864 return false;
865
866 // Look for fractional components
867 for (i = x->exponent; i < MP_SIZE; i++) {
868 if (x->fraction[i] != 0)
869 return false;
870 }
871
872 return true;*/
873}
874
875
876bool_Bool
877mp_is_positive_integer(const MPNumber *x)
878{
879 if (mp_is_complex(x))
880 return false0;
881 else
882 return x->sign >= 0 && mp_is_integer(x);
883}
884
885
886bool_Bool
887mp_is_natural(const MPNumber *x)
888{
889 if (mp_is_complex(x))
890 return false0;
891 else
892 return x->sign > 0 && mp_is_integer(x);
893}
894
895
896bool_Bool
897mp_is_complex(const MPNumber *x)
898{
899 return x->im_sign != 0;
900}
901
902
903bool_Bool
904mp_is_equal(const MPNumber *x, const MPNumber *y)
905{
906 return mp_compare_mp_to_mp(x, y) == 0;
907}
908
909
910/* Return e^x for |x| < 1 USING AN O(SQRT(T).M(T)) ALGORITHM
911 * DESCRIBED IN - R. P. BRENT, THE COMPLEXITY OF MULTIPLE-
912 * PRECISION ARITHMETIC (IN COMPLEXITY OF COMPUTATIONAL PROBLEM
913 * SOLVING, UNIV. OF QUEENSLAND PRESS, BRISBANE, 1976, 126-165).
914 * ASYMPTOTICALLY FASTER METHODS EXIST, BUT ARE NOT USEFUL
915 * UNLESS T IS VERY LARGE. SEE COMMENTS TO MP_ATAN AND MPPIGL.
916 */
917static void
918mp_exp(const MPNumber *x, MPNumber *z)
919{
920 int i, q;
921 float rlb;
922 MPNumber t1, t2;
923
924 /* e^0 = 1 */
925 if (mp_is_zero(x)) {
926 mp_set_from_integer(1, z);
927 return;
928 }
929
930 /* Only defined for |x| < 1 */
931 if (x->exponent > 0) {
932 mperr("*** ABS(X) NOT LESS THAN 1 IN CALL TO MP_EXP ***");
933 mp_set_from_integer(0, z);
934 return;
935 }
936
937 mp_set_from_mp(x, &t1);
938 rlb = log((float)MP_BASE10000);
939
940 /* Compute approxicafely optimal q (and divide x by 2^q) */
941 q = (int)(sqrt((float)MP_T100 * 0.48f * rlb) + (float) x->exponent * 1.44f * rlb);
942
943 /* HALVE Q TIMES */
944 if (q > 0) {
945 int ib, ic;
946
947 ib = MP_BASE10000 << 2;
948 ic = 1;
949 for (i = 1; i <= q; ++i) {
950 ic *= 2;
951 if (ic < ib && ic != MP_BASE10000 && i < q)
952 continue;
953 mp_divide_integer(&t1, ic, &t1);
954 ic = 1;
955 }
956 }
957
958 if (mp_is_zero(&t1)) {
959 mp_set_from_integer(0, z);
960 return;
961 }
962
963 /* Sum series, reducing t where possible */
964 mp_set_from_mp(&t1, z);
965 mp_set_from_mp(&t1, &t2);
966 for (i = 2; MP_T100 + t2.exponent - z->exponent > 0; i++) {
967 mp_multiply(&t1, &t2, &t2);
968 mp_divide_integer(&t2, i, &t2);
969 mp_add(&t2, z, z);
970 if (mp_is_zero(&t2))
971 break;
972 }
973
974 /* Apply (x+1)^2 - 1 = x(2 + x) for q iterations */
975 for (i = 1; i <= q; ++i) {
976 mp_add_integer(z, 2, &t1);
977 mp_multiply(&t1, z, z);
978 }
979
980 mp_add_integer(z, 1, z);
981}
982
983
984static void
985mp_epowy_real(const MPNumber *x, MPNumber *z)
986{
987 float r__1;
988 int i, ix, xs, tss;
989 float rx, rz;
990 MPNumber t1, t2;
991
992 /* e^0 = 1 */
993 if (mp_is_zero(x)) {
10
←
Taking false branch
→
994 mp_set_from_integer(1, z);
995 return;
996 }
997
998 /* If |x| < 1 use mp_exp */
999 if (x->exponent <= 0) {
11
←
Assuming field 'exponent' is > 0
→
12
←
Taking false branch
→
1000 mp_exp(x, z);
1001 return;
1002 }
1003
1004 /* NOW SAFE TO CONVERT X TO REAL */
1005 rx = mp_cast_to_float(x);
1006
1007 /* SAVE SIGN AND WORK WITH ABS(X) */
1008 xs = x->sign;
1009 mp_abs(x, &t2);
1010
1011 /* GET FRACTIONAL AND INTEGER PARTS OF ABS(X) */
1012 ix = mp_cast_to_int(&t2);
1013 mp_fractional_component(&t2, &t2);
13
←
Calling 'mp_fractional_component'
→
1014
1015 /* ATTACH SIGN TO FRACTIONAL PART AND COMPUTE EXP OF IT */
1016 t2.sign *= xs;
1017 mp_exp(&t2, z);
1018
1019 /* COMPUTE E-2 OR 1/E USING TWO EXTRA DIGITS IN CASE ABS(X) LARGE
1020 * (BUT ONLY ONE EXTRA DIGIT IF T < 4)
1021 */
1022 if (MP_T100 < 4)
1023 tss = MP_T100 + 1;
1024 else
1025 tss = MP_T100 + 2;
1026
1027 /* LOOP FOR E COMPUTATION. DECREASE T IF POSSIBLE. */
1028 /* Computing MIN */
1029 mp_set_from_integer(xs, &t1);
1030
1031 t2.sign = 0;
1032 for (i = 2 ; ; i++) {
1033 if (min(tss, tss + 2 + t1.exponent)((tss) <= (tss + 2 + t1.exponent) ? (tss) : (tss + 2 + t1.
exponent)) <= 2)
1034 break;
1035
1036 mp_divide_integer(&t1, i * xs, &t1);
1037 mp_add(&t2, &t1, &t2);
1038 if (mp_is_zero(&t1))
1039 break;
1040 }
1041
1042 /* RAISE E OR 1/E TO POWER IX */
1043 if (xs > 0)
1044 mp_add_integer(&t2, 2, &t2);
1045 mp_xpowy_integer(&t2, ix, &t2);
1046
1047 /* MULTIPLY EXPS OF INTEGER AND FRACTIONAL PARTS */
1048 mp_multiply(z, &t2, z);
1049
1050 /* CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS ABS(X) LARGE
1051 * (WHEN EXP MIGHT OVERFLOW OR UNDERFLOW)
1052 */
1053 if (fabs(rx) > 10.0f)
1054 return;
1055
1056 rz = mp_cast_to_float(z);
1057 r__1 = rz - exp(rx);
1058 if (fabs(r__1) < rz * 0.01f)
1059 return;
1060
1061 /* THE FOLLOWING MESSAGE MAY INDICATE THAT
1062 * B**(T-1) IS TOO SMALL, OR THAT M IS TOO SMALL SO THE
1063 * RESULT UNDERFLOWED.
1064 */
1065 mperr("*** ERROR OCCURRED IN MP_EPOWY, RESULT INCORRECT ***");
1066}
1067
1068
1069void
1070mp_epowy(const MPNumber *x, MPNumber *z)
1071{
1072 /* e^0 = 1 */
1073 if (mp_is_zero(x)) {
7
←
Taking false branch
→
1074 mp_set_from_integer(1, z);
1075 return;
1076 }
1077
1078 if (mp_is_complex(x)) {
8
←
Taking false branch
→
1079 MPNumber x_real, r, theta;
1080
1081 mp_real_component(x, &x_real);
1082 mp_imaginary_component(x, &theta);
1083
1084 mp_epowy_real(&x_real, &r);
1085 mp_set_from_polar(&r, MP_RADIANS, &theta, z);
1086 }
1087 else
1088 mp_epowy_real(x, z);
9
←
Calling 'mp_epowy_real'
→
1089}
1090
1091
1092/* RETURNS K = K/GCD AND L = L/GCD, WHERE GCD IS THE
1093 * GREATEST COMMON DIVISOR OF K AND L.
1094 * SAVE INPUT PARAMETERS IN LOCAL VARIABLES
1095 */
1096void
1097mp_gcd(int64_t *k, int64_t *l)
1098{
1099 int64_t i, j;
1100
1101 i = abs(*k);
1102 j = abs(*l);
1103 if (j == 0) {
1104 /* IF J = 0 RETURN (1, 0) UNLESS I = 0, THEN (0, 0) */
1105 *k = 1;
1106 *l = 0;
1107 if (i == 0)
1108 *k = 0;
1109 return;
1110 }
1111
1112 /* EUCLIDEAN ALGORITHM LOOP */
1113 do {
1114 i %= j;
1115 if (i == 0) {
1116 *k = *k / j;
1117 *l = *l / j;
1118 return;
1119 }
1120 j %= i;
1121 } while (j != 0);
1122
1123 /* HERE J IS THE GCD OF K AND L */
1124 *k = *k / i;
1125 *l = *l / i;
1126}
1127
1128
1129bool_Bool
1130mp_is_zero(const MPNumber *x)
1131{
1132 return x->sign == 0 && x->im_sign == 0;
1133}
1134
1135
1136bool_Bool
1137mp_is_negative(const MPNumber *x)
1138{
1139 return x->sign < 0;
1140}
1141
1142
1143bool_Bool
1144mp_is_greater_equal(const MPNumber *x, const MPNumber *y)
1145{
1146 return mp_compare_mp_to_mp(x, y) >= 0;
1147}
1148
1149
1150bool_Bool
1151mp_is_greater_than(const MPNumber *x, const MPNumber *y)
1152{
1153 return mp_compare_mp_to_mp(x, y) > 0;
1154}
1155
1156
1157bool_Bool
1158mp_is_less_equal(const MPNumber *x, const MPNumber *y)
1159{
1160 return mp_compare_mp_to_mp(x, y) <= 0;
1161}
1162
1163
1164/* RETURNS MP Y = LN(1+X) IF X IS AN MP NUMBER SATISFYING THE
1165 * CONDITION ABS(X) < 1/B, ERROR OTHERWISE.
1166 * USES NEWTONS METHOD TO SOLVE THE EQUATION
1167 * EXP1(-Y) = X, THEN REVERSES SIGN OF Y.
1168 */
1169static void
1170mp_lns(const MPNumber *x, MPNumber *z)
1171{
1172 int t, it0;
1173 MPNumber t1, t2, t3;
1174
1175 /* ln(1+0) = 0 */
1176 if (mp_is_zero(x)) {
1177 mp_set_from_integer(0, z);
1178 return;
1179 }
1180
1181 /* Get starting approximation -ln(1+x) ~= -x + x^2/2 - x^3/3 + x^4/4 */
1182 mp_set_from_mp(x, &t2);
1183 mp_divide_integer(x, 4, &t1);
1184 mp_add_fraction(&t1, -1, 3, &t1);
1185 mp_multiply(x, &t1, &t1);
1186 mp_add_fraction(&t1, 1, 2, &t1);
1187 mp_multiply(x, &t1, &t1);
1188 mp_add_integer(&t1, -1, &t1);
1189 mp_multiply(x, &t1, z);
1190
1191 /* Solve using Newtons method */
1192 it0 = t = 5;
1193 while(1)
1194 {
1195 int ts2, ts3;
1196
1197 /* t3 = (e^t3 - 1) */
1198 /* z = z - (t2 + t3 + (t2 * t3)) */
1199 mp_epowy(z, &t3);
1200 mp_add_integer(&t3, -1, &t3);
1201 mp_multiply(&t2, &t3, &t1);
1202 mp_add(&t3, &t1, &t3);
1203 mp_add(&t2, &t3, &t3);
1204 mp_subtract(z, &t3, z);
1205 if (t >= MP_T100)
1206 break;
1207
1208 /* FOLLOWING LOOP COMPUTES NEXT VALUE OF T TO USE.
1209 * BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE,
1210 * WE CAN ALMOST DOUBLE T EACH TIME.
1211 */
1212 ts3 = t;
1213 t = MP_T100;
1214 do {
1215 ts2 = t;
1216 t = (t + it0) / 2;
1217 } while (t > ts3);
1218 t = ts2;
1219 }
1220
1221 /* CHECK THAT NEWTON ITERATION WAS CONVERGING AS EXPECTED */
1222 if (t3.sign != 0 && t3.exponent << 1 > it0 - MP_T100) {
1223 mperr("*** ERROR OCCURRED IN MP_LNS, NEWTON ITERATION NOT CONVERGING PROPERLY ***");
1224 }
1225
1226 z->sign = -z->sign;
1227}
1228
1229
1230static void
1231mp_ln_real(const MPNumber *x, MPNumber *z)
1232{
1233 int e, k;
1234 double rx, rlx;
1235 MPNumber t1, t2;
1236
1237 /* LOOP TO GET APPROXICAFE LN(X) USING SINGLE-PRECISION */
1238 mp_set_from_mp(x, &t1);
1239 mp_set_from_integer(0, z);
1240 for(k = 0; k < 10; k++)
1241 {
1242 /* COMPUTE FINAL CORRECTION ACCURATELY USING MP_LNS */
1243 mp_add_integer(&t1, -1, &t2);
1244 if (mp_is_zero(&t2) || t2.exponent + 1 <= 0) {
1245 mp_lns(&t2, &t2);
1246 mp_add(z, &t2, z);
1247 return;
1248 }
1249
1250 /* REMOVE EXPONENT TO AVOID FLOATING-POINT OVERFLOW */
1251 e = t1.exponent;
1252 t1.exponent = 0;
1253 rx = mp_cast_to_double(&t1);
1254 t1.exponent = e;
1255 rlx = log(rx) + e * log(MP_BASE10000);
1256 mp_set_from_double(-(double)rlx, &t2);
1257
1258 /* UPDATE Z AND COMPUTE ACCURATE EXP OF APPROXICAFE LOG */
1259 mp_subtract(z, &t2, z);
1260 mp_epowy(&t2, &t2);
1261
1262 /* COMPUTE RESIDUAL WHOSE LOG IS STILL TO BE FOUND */
1263 mp_multiply(&t1, &t2, &t1);
1264 }
1265
1266 mperr("*** ERROR IN MP_LN, ITERATION NOT CONVERGING ***");
1267}
1268
1269
1270void
1271mp_ln(const MPNumber *x, MPNumber *z)
1272{
1273 /* ln(0) undefined */
1274 if (mp_is_zero(x)) {
1275 /* Translators: Error displayed when attempting to take logarithm of zero */
1276 mperr(_("Logarithm of zero is undefined")gettext ("Logarithm of zero is undefined"));
1277 mp_set_from_integer(0, z);
1278 return;
1279 }
1280
1281 /* ln(-x) complex */
1282 /* FIXME: Make complex numbers optional */
1283 /*if (mp_is_negative(x)) {
1284 // Translators: Error displayed attempted to take logarithm of negative value
1285 mperr(_("Logarithm of negative values is undefined"));
1286 mp_set_from_integer(0, z);
1287 return;
1288 }*/
1289
1290 if (mp_is_complex(x) || mp_is_negative(x)) {
1291 MPNumber r, theta, z_real;
1292
1293 /* ln(re^iθ) = e^(ln(r)+iθ) */
1294 mp_abs(x, &r);
1295 mp_arg(x, MP_RADIANS, &theta);
1296
1297 mp_ln_real(&r, &z_real);
1298 mp_set_from_complex(&z_real, &theta, z);
1299 }
1300 else
1301 mp_ln_real(x, z);
1302}
1303
1304
1305void
1306mp_logarithm(int64_t n, const MPNumber *x, MPNumber *z)
1307{
1308 MPNumber t1, t2;
1309
1310 /* log(0) undefined */
1311 if (mp_is_zero(x)) {
1312 /* Translators: Error displayed when attempting to take logarithm of zero */
1313 mperr(_("Logarithm of zero is undefined")gettext ("Logarithm of zero is undefined"));
1314 mp_set_from_integer(0, z);
1315 return;
1316 }
1317
1318 /* logn(x) = ln(x) / ln(n) */
1319 mp_set_from_integer(n, &t1);
1320 mp_ln(&t1, &t1);
1321 mp_ln(x, &t2);
1322 mp_divide(&t2, &t1, z);
1323}
1324
1325
1326bool_Bool
1327mp_is_less_than(const MPNumber *x, const MPNumber *y)
1328{
1329 return mp_compare_mp_to_mp(x, y) < 0;
1330}
1331
1332
1333static void
1334mp_multiply_real(const MPNumber *x, const MPNumber *y, MPNumber *z)
1335{
1336 int c, i, xi;
1337 MPNumber r;
1338
1339 /* x*0 = 0*y = 0 */
1340 if (x->sign == 0 || y->sign == 0) {
1341 mp_set_from_integer(0, z);
1342 return;
1343 }
1344
1345 z->sign = x->sign * y->sign;
1346 z->exponent = x->exponent + y->exponent;
1347 memset(&r, 0, sizeof(MPNumber));
1348
1349 /* PERFORM MULTIPLICATION */
1350 c = 8;
1351 for (i = 0; i < MP_T100; i++) {
1352 int j;
1353
1354 xi = x->fraction[i];
1355
1356 /* FOR SPEED, PUT THE NUMBER WITH MANY ZEROS FIRST */
1357 if (xi == 0)
1358 continue;
1359
1360 /* Computing MIN */
1361 for (j = 0; j < min(MP_T, MP_T + 3 - i)((100) <= (100 + 3 - i) ? (100) : (100 + 3 - i)); j++)
1362 r.fraction[i+j+1] += xi * y->fraction[j];
1363 c--;
1364 if (c > 0)
1365 continue;
1366
1367 /* CHECK FOR LEGAL BASE B DIGIT */
1368 if (xi < 0 || xi >= MP_BASE10000) {
1369 mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
1370 mp_set_from_integer(0, z);
1371 return;
1372 }
1373
1374 /* PROPAGATE CARRIES AT END AND EVERY EIGHTH TIME,
1375 * FASTER THAN DOING IT EVERY TIME.
1376 */
1377 for (j = MP_T100 + 3; j >= 0; j--) {
1378 int ri = r.fraction[j] + c;
1379 if (ri < 0) {
1380 mperr("*** INTEGER OVERFLOW IN MP_MULTIPLY, B TOO LARGE ***");
1381 mp_set_from_integer(0, z);
1382 return;
1383 }
1384 c = ri / MP_BASE10000;
1385 r.fraction[j] = ri - MP_BASE10000 * c;
1386 }
1387 if (c != 0) {
1388 mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
1389 mp_set_from_integer(0, z);
1390 return;
1391 }
1392 c = 8;
1393 }
1394
1395 if (c != 8) {
1396 if (xi < 0 || xi >= MP_BASE10000) {
1397 mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
1398 mp_set_from_integer(0, z);
1399 return;
1400 }
1401
1402 c = 0;
1403 for (i = MP_T100 + 3; i >= 0; i--) {
1404 int ri = r.fraction[i] + c;
1405 if (ri < 0) {
1406 mperr("*** INTEGER OVERFLOW IN MP_MULTIPLY, B TOO LARGE ***");
1407 mp_set_from_integer(0, z);
1408 return;
1409 }
1410 c = ri / MP_BASE10000;
1411 r.fraction[i] = ri - MP_BASE10000 * c;
1412 }
1413
1414 if (c != 0) {
1415 mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
1416 mp_set_from_integer(0, z);
1417 return;
1418 }
1419 }
1420
1421 /* Clear complex part */
1422 z->im_sign = 0;
1423 z->im_exponent = 0;
1424 memset(z->im_fraction, 0, sizeof(int) * MP_SIZE1000);
1425
1426 /* NORMALIZE AND ROUND RESULT */
1427 // FIXME: Use stack variable because of mp_normalize brokeness
1428 for (i = 0; i < MP_SIZE1000; i++)
1429 z->fraction[i] = r.fraction[i];
1430 mp_normalize(z);
1431}
1432
1433
1434void
1435mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
1436{
1437 /* x*0 = 0*y = 0 */
1438 if (mp_is_zero(x) || mp_is_zero(y)) {
1439 mp_set_from_integer(0, z);
1440 return;
1441 }
1442
1443 /* (a+bi)(c+di) = (ac-bd)+(ad+bc)i */
1444 if (mp_is_complex(x) || mp_is_complex(y)) {
1445 MPNumber real_x, real_y, im_x, im_y, t1, t2, real_z, im_z;
1446
1447 mp_real_component(x, &real_x);
1448 mp_imaginary_component(x, &im_x);
1449 mp_real_component(y, &real_y);
1450 mp_imaginary_component(y, &im_y);
1451
1452 mp_multiply_real(&real_x, &real_y, &t1);
1453 mp_multiply_real(&im_x, &im_y, &t2);
1454 mp_subtract(&t1, &t2, &real_z);
1455
1456 mp_multiply_real(&real_x, &im_y, &t1);
1457 mp_multiply_real(&im_x, &real_y, &t2);
1458 mp_add(&t1, &t2, &im_z);
1459
1460 mp_set_from_complex(&real_z, &im_z, z);
1461 }
1462 else {
1463 mp_multiply_real(x, y, z);
1464 }
1465}
1466
1467
1468static void
1469mp_multiply_integer_real(const MPNumber *x, int64_t y, MPNumber *z)
1470{
1471 int c, i;
1472 MPNumber x_copy;
1473
1474 /* x*0 = 0*y = 0 */
1475 if (mp_is_zero(x) || y == 0) {
1476 mp_set_from_integer(0, z);
1477 return;
1478 }
1479
1480 /* x*1 = x, x*-1 = -x */
1481 // FIXME: Why is this not working? mp_ext is using this function to do a normalization
1482 /*if (y == 1 || y == -1) {
1483 if (y < 0)
1484 mp_invert_sign(x, z);
1485 else
1486 mp_set_from_mp(x, z);
1487 return;
1488 }*/
1489
1490 /* Copy x as z may also refer to x */
1491 mp_set_from_mp(x, &x_copy);
1492 mp_set_from_integer(0, z);
1493
1494 if (y < 0) {
1495 y = -y;
1496 z->sign = -x_copy.sign;
1497 }
1498 else
1499 z->sign = x_copy.sign;
1500 z->exponent = x_copy.exponent + 4;
1501
1502 /* FORM PRODUCT IN ACCUMULATOR */
1503 c = 0;
1504
1505 /* IF y*B NOT REPRESENTABLE AS AN INTEGER WE HAVE TO SIMULATE
1506 * DOUBLE-PRECISION MULTIPLICATION.
1507 */
1508
1509 /* Computing MAX */
1510 if (y >= max(MP_BASE << 3, 32767 / MP_BASE)((10000 << 3) >= (32767 / 10000) ? (10000 << 3
) : (32767 / 10000))) {
1511 int64_t j1, j2;
1512
1513 /* HERE J IS TOO LARGE FOR SINGLE-PRECISION MULTIPLICATION */
1514 j1 = y / MP_BASE10000;
1515 j2 = y - j1 * MP_BASE10000;
1516
1517 /* FORM PRODUCT */
1518 for (i = MP_T100 + 3; i >= 0; i--) {
1519 int64_t c1, c2, is, ix, t;
1520
1521 c1 = c / MP_BASE10000;
1522 c2 = c - MP_BASE10000 * c1;
1523 ix = 0;
1524 if (i > 3)
1525 ix = x_copy.fraction[i - 4];
1526
1527 t = j2 * ix + c2;
1528 is = t / MP_BASE10000;
1529 c = j1 * ix + c1 + is;
1530 z->fraction[i] = t - MP_BASE10000 * is;
1531 }
1532 }
1533 else
1534 {
1535 int64_t ri = 0;
1536
1537 for (i = MP_T100 + 3; i >= 4; i--) {
1538 ri = y * x_copy.fraction[i - 4] + c;
1539 c = ri / MP_BASE10000;
1540 z->fraction[i] = ri - MP_BASE10000 * c;
1541 }
1542
1543 /* CHECK FOR INTEGER OVERFLOW */
1544 if (ri < 0) {
1545 mperr("*** INTEGER OVERFLOW IN mp_multiply_integer, B TOO LARGE ***");
1546 mp_set_from_integer(0, z);
1547 return;
1548 }
1549
1550 /* HAVE TO TREAT FIRST FOUR WORDS OF R SEPARATELY */
1551 for (i = 3; i >= 0; i--) {
1552 int t;
1553
1554 t = c;
1555 c = t / MP_BASE10000;
1556 z->fraction[i] = t - MP_BASE10000 * c;
1557 }
1558 }
1559
1560 /* HAVE TO SHIFT RIGHT HERE AS CARRY OFF END */
1561 while (c != 0) {
1562 int64_t t;
1563
1564 for (i = MP_T100 + 3; i >= 1; i--)
1565 z->fraction[i] = z->fraction[i - 1];
1566 t = c;
1567 c = t / MP_BASE10000;
1568 z->fraction[0] = t - MP_BASE10000 * c;
1569 z->exponent++;
1570 }
1571
1572 if (c < 0) {
1573 mperr("*** INTEGER OVERFLOW IN mp_multiply_integer, B TOO LARGE ***");
1574 mp_set_from_integer(0, z);
1575 return;
1576 }
1577
1578 z->im_sign = 0;
1579 z->im_exponent = 0;
1580 memset(z->im_fraction, 0, sizeof(int) * MP_SIZE1000);
1581 mp_normalize(z);
1582}
1583
1584
1585void
1586mp_multiply_integer(const MPNumber *x, int64_t y, MPNumber *z)
1587{
1588 if (mp_is_complex(x)) {
1589 MPNumber re_z, im_z;
1590 mp_real_component(x, &re_z);
1591 mp_imaginary_component(x, &im_z);
1592 mp_multiply_integer_real(&re_z, y, &re_z);
1593 mp_multiply_integer_real(&im_z, y, &im_z);
1594 mp_set_from_complex(&re_z, &im_z, z);
1595 }
1596 else
1597 mp_multiply_integer_real(x, y, z);
1598}
1599
1600
1601void
1602mp_multiply_fraction(const MPNumber *x, int64_t numerator, int64_t denominator, MPNumber *z)
1603{
1604 if (denominator == 0) {
1605 mperr(_("Division by zero is undefined")gettext ("Division by zero is undefined"));
1606 mp_set_from_integer(0, z);
1607 return;
1608 }
1609
1610 if (numerator == 0) {
1611 mp_set_from_integer(0, z);
1612 return;
1613 }
1614
1615 /* Reduce to lowest terms */
1616 mp_gcd(&numerator, &denominator);
1617 mp_divide_integer(x, denominator, z);
1618 mp_multiply_integer(z, numerator, z);
1619}
1620
1621
1622void
1623mp_invert_sign(const MPNumber *x, MPNumber *z)
1624{
1625 mp_set_from_mp(x, z);
1626 z->sign = -z->sign;
1627 z->im_sign = -z->im_sign;
1628}
1629
1630
1631// FIXME: Is r->fraction large enough? It seems to be in practise but it may be MP_T+4 instead of MP_T
1632// FIXME: There is some sort of stack corruption/use of unitialised variables here. Some functions are
1633// using stack variables as x otherwise there are corruption errors. e.g. "Cos(45) - 1/Sqrt(2) = -0"
1634// (try in scientific mode)
1635void
1636mp_normalize(MPNumber *x)
1637{
1638 int start_index;
1639
1640 /* Find first non-zero digit */
1641 for (start_index = 0; start_index < MP_SIZE1000 && x->fraction[start_index] == 0; start_index++);
1642
1643 /* Mark as zero */
1644 if (start_index >= MP_SIZE1000) {
1645 x->sign = 0;
1646 x->exponent = 0;
1647 return;
1648 }
1649
1650 /* Shift left so first digit is non-zero */
1651 if (start_index > 0) {
1652 int i;
1653
1654 x->exponent -= start_index;
1655 for (i = 0; (i + start_index) < MP_SIZE1000; i++)
1656 x->fraction[i] = x->fraction[i + start_index];
1657 for (; i < MP_SIZE1000; i++)
1658 x->fraction[i] = 0;
1659 }
1660}
1661
1662
1663static void
1664mp_pwr(const MPNumber *x, const MPNumber *y, MPNumber *z)
1665{
1666 MPNumber t;
1667
1668 /* (-x)^y imaginary */
1669 /* FIXME: Make complex numbers optional */
1670 /*if (x->sign < 0) {
1671 mperr(_("The power of negative numbers is only defined for integer exponents"));
1672 mp_set_from_integer(0, z);
1673 return;
1674 }*/
1675
1676 /* 0^y = 0, 0^-y undefined */
1677 if (mp_is_zero(x)) {
4
←
Taking false branch
→
1678 mp_set_from_integer(0, z);
1679 if (y->sign < 0)
1680 mperr(_("The power of zero is undefined for a negative exponent")gettext ("The power of zero is undefined for a negative exponent"
));
1681 return;
1682 }
1683
1684 /* x^0 = 1 */
1685 if (mp_is_zero(y)) {
5
←
Taking false branch
→
1686 mp_set_from_integer(1, z);
1687 return;
1688 }
1689
1690 mp_ln(x, &t);
1691 mp_multiply(y, &t, z);
1692 mp_epowy(z, z);
6
←
Calling 'mp_epowy'
→
1693}
1694
1695
1696static void
1697mp_reciprocal_real(const MPNumber *x, MPNumber *z)
1698{
1699 MPNumber t1, t2;
1700 int it0, t;
1701
1702 /* 1/0 invalid */
1703 if (mp_is_zero(x)) {
1704 mperr(_("Reciprocal of zero is undefined")gettext ("Reciprocal of zero is undefined"));
1705 mp_set_from_integer(0, z);
1706 return;
1707 }
1708
1709 /* Start by approximating value using floating point */
1710 mp_set_from_mp(x, &t1);
1711 t1.exponent = 0;
1712 mp_set_from_float(1.0f / mp_cast_to_float(&t1), &t1);
1713 t1.exponent -= x->exponent;
1714
1715 it0 = t = 3;
1716 while(1) {
1717 int ts2, ts3;
1718
1719 /* t1 = t1 - (t1 * ((x * t1) - 1)) (2*t1 - t1^2*x) */
1720 mp_multiply(x, &t1, &t2);
1721 mp_add_integer(&t2, -1, &t2);
1722 mp_multiply(&t1, &t2, &t2);
1723 mp_subtract(&t1, &t2, &t1);
1724 if (t >= MP_T100)
1725 break;
1726
1727 /* FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE
1728 * BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).
1729 */
1730 ts3 = t;
1731 t = MP_T100;
1732 do {
1733 ts2 = t;
1734 t = (t + it0) / 2;
1735 } while (t > ts3);
1736 t = min(ts2, MP_T)((ts2) <= (100) ? (ts2) : (100));
1737 }
1738
1739 /* RETURN IF NEWTON ITERATION WAS CONVERGING */
1740 if (t2.sign != 0 && (t1.exponent - t2.exponent) << 1 < MP_T100 - it0) {
1741 /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL,
1742 * OR THAT THE STARTING APPROXIMATION IS NOT ACCURATE ENOUGH.
1743 */
1744 mperr("*** ERROR OCCURRED IN MP_RECIPROCAL, NEWTON ITERATION NOT CONVERGING PROPERLY ***");
1745 }
1746
1747 mp_set_from_mp(&t1, z);
1748}
1749
1750
1751void
1752mp_reciprocal(const MPNumber *x, MPNumber *z)
1753{
1754 if (mp_is_complex(x)) {
1755 MPNumber t1, t2;
1756 MPNumber real_x, im_x;
1757
1758 mp_real_component(x, &real_x);
1759 mp_imaginary_component(x, &im_x);
1760
1761 /* 1/(a+bi) = (a-bi)/(a+bi)(a-bi) = (a-bi)/(a²+b²) */
1762 mp_multiply(&real_x, &real_x, &t1);
1763 mp_multiply(&im_x, &im_x, &t2);
1764 mp_add(&t1, &t2, &t1);
1765 mp_reciprocal_real(&t1, z);
1766 mp_conjugate(x, &t1);
1767 mp_multiply(&t1, z, z);
1768 }
1769 else
1770 mp_reciprocal_real(x, z);
1771}
1772
1773
1774static void
1775mp_root_real(const MPNumber *x, int64_t n, MPNumber *z)
1776{
1777 float approximation;
1778 int ex, np, it0, t;
1779 MPNumber t1, t2;
1780
1781 /* x^(1/1) = x */
1782 if (n == 1) {
1783 mp_set_from_mp(x, z);
1784 return;
1785 }
1786
1787 /* x^(1/0) invalid */
1788 if (n == 0) {
1789 mperr(_("Root must be non-zero")gettext ("Root must be non-zero"));
1790 mp_set_from_integer(0, z);
1791 return;
1792 }
1793
1794 np = abs(n);
1795
1796 /* LOSS OF ACCURACY IF NP LARGE, SO ONLY ALLOW NP <= MAX (B, 64) */
1797 if (np > max(MP_BASE, 64)((10000) >= (64) ? (10000) : (64))) {
1798 mperr("*** ABS(N) TOO LARGE IN CALL TO MP_ROOT ***");
1799 mp_set_from_integer(0, z);
1800 return;
1801 }
1802
1803 /* 0^(1/n) = 0 for positive n */
1804 if (mp_is_zero(x)) {
1805 mp_set_from_integer(0, z);
1806 if (n <= 0)
1807 mperr(_("Negative root of zero is undefined")gettext ("Negative root of zero is undefined"));
1808 return;
1809 }
1810
1811 // FIXME: Imaginary root
1812 if (x->sign < 0 && np % 2 == 0) {
1813 mperr(_("nth root of negative number is undefined for even n")gettext ("nth root of negative number is undefined for even n"
));
1814 mp_set_from_integer(0, z);
1815 return;
1816 }
1817
1818 /* DIVIDE EXPONENT BY NP */
1819 ex = x->exponent / np;
1820
1821 /* Get initial approximation */
1822 mp_set_from_mp(x, &t1);
1823 t1.exponent = 0;
1824 approximation = exp(((float)(np * ex - x->exponent) * log((float)MP_BASE10000) -
1825 log((fabs(mp_cast_to_float(&t1))))) / (float)np);
1826 mp_set_from_float(approximation, &t1);
1827 t1.sign = x->sign;
1828 t1.exponent -= ex;
1829
1830 /* MAIN ITERATION LOOP */
1831 it0 = t = 3;
1832 while(1) {
1833 int ts2, ts3;
1834
1835 /* t1 = t1 - ((t1 * ((x * t1^np) - 1)) / np) */
1836 mp_xpowy_integer(&t1, np, &t2);
1837 mp_multiply(x, &t2, &t2);
1838 mp_add_integer(&t2, -1, &t2);
1839 mp_multiply(&t1, &t2, &t2);
1840 mp_divide_integer(&t2, np, &t2);
1841 mp_subtract(&t1, &t2, &t1);
1842
1843 /* FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE BECAUSE
1844 * NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).
1845 */
1846 if (t >= MP_T100)
1847 break;
1848
1849 ts3 = t;
1850 t = MP_T100;
1851 do {
1852 ts2 = t;
1853 t = (t + it0) / 2;
1854 } while (t > ts3);
1855 t = min(ts2, MP_T)((ts2) <= (100) ? (ts2) : (100));
1856 }
1857
1858 /* NOW R(I2) IS X**(-1/NP)
1859 * CHECK THAT NEWTON ITERATION WAS CONVERGING
1860 */
1861 if (t2.sign != 0 && (t1.exponent - t2.exponent) << 1 < MP_T100 - it0) {
1862 /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL,
1863 * OR THAT THE INITIAL APPROXIMATION OBTAINED USING ALOG AND EXP
1864 * IS NOT ACCURATE ENOUGH.
1865 */
1866 mperr("*** ERROR OCCURRED IN MP_ROOT, NEWTON ITERATION NOT CONVERGING PROPERLY ***");
1867 }
1868
1869 if (n >= 0) {
1870 mp_xpowy_integer(&t1, n - 1, &t1);
1871 mp_multiply(x, &t1, z);
1872 return;
1873 }
1874
1875 mp_set_from_mp(&t1, z);
1876}
1877
1878
1879void
1880mp_root(const MPNumber *x, int64_t n, MPNumber *z)
1881{
1882 if (!mp_is_complex(x) && mp_is_negative(x) && n % 2 == 1) {
1883 mp_abs(x, z);
1884 mp_root_real(z, n, z);
1885 mp_invert_sign(z, z);
1886 }
1887 else if (mp_is_complex(x) || mp_is_negative(x)) {
1888 MPNumber r, theta;
1889
1890 mp_abs(x, &r);
1891 mp_arg(x, MP_RADIANS, &theta);
1892
1893 mp_root_real(&r, n, &r);
1894 mp_divide_integer(&theta, n, &theta);
1895 mp_set_from_polar(&r, MP_RADIANS, &theta, z);
1896 }
1897 else
1898 mp_root_real(x, n, z);
1899}
1900
1901
1902void
1903mp_sqrt(const MPNumber *x, MPNumber *z)
1904{
1905 if (mp_is_zero(x))
1906 mp_set_from_integer(0, z);
1907 /* FIXME: Make complex numbers optional */
1908 /*else if (x->sign < 0) {
1909 mperr(_("Square root is undefined for negative values"));
1910 mp_set_from_integer(0, z);
1911 }*/
1912 else {
1913 MPNumber t;
1914
1915 mp_root(x, -2, &t);
1916 mp_multiply(x, &t, z);
1917 mp_ext(t.fraction[0], z->fraction[0], z);
1918 }
1919}
1920
1921
1922void
1923mp_factorial(const MPNumber *x, MPNumber *z)
1924{
1925 int i, value;
1926
1927 /* 0! == 1 */
1928 if (mp_is_zero(x)) {
1929 mp_set_from_integer(1, z);
1930 return;
1931 }
1932 if (!mp_is_natural(x)) {
1933 /* Translators: Error displayed when attempted take the factorial of a fractional number */
1934 mperr(_("Factorial is only defined for natural numbers")gettext ("Factorial is only defined for natural numbers"));
1935 mp_set_from_integer(0, z);
1936 return;
1937 }
1938
1939 /* Convert to integer - if couldn't be converted then the factorial would be too big anyway */
1940 value = mp_cast_to_int(x);
1941 mp_set_from_mp(x, z);
1942 for (i = 2; i < value; i++)
1943 mp_multiply_integer(z, i, z);
1944}
1945
1946
1947void
1948mp_modulus_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
1949{
1950 MPNumber t1, t2;
1951
1952 if (!mp_is_integer(x) || !mp_is_integer(y)) {
1953 /* Translators: Error displayed when attemping to do a modulus division on non-integer numbers */
1954 mperr(_("Modulus division is only defined for integers")gettext ("Modulus division is only defined for integers"));
1955 mp_set_from_integer(0, z);
1956 }
1957
1958 mp_divide(x, y, &t1);
1959 mp_floor(&t1, &t1);
1960 mp_multiply(&t1, y, &t2);
1961 mp_subtract(x, &t2, z);
1962
1963 mp_set_from_integer(0, &t1);
1964 if ((mp_is_less_than(y, &t1) && mp_is_greater_than(z, &t1)) || mp_is_less_than(z, &t1))
1965 mp_add(z, y, z);
1966}
1967
1968
1969void
1970mp_xpowy(const MPNumber *x, const MPNumber *y, MPNumber *z)
1971{
1972 if (mp_is_integer(y)) {
1
Taking false branch
→
1973 mp_xpowy_integer(x, mp_cast_to_int(y), z);
1974 } else {
1975 MPNumber reciprocal;
1976 mp_reciprocal(y, &reciprocal);
1977 if (mp_is_integer(&reciprocal))
2
←
Taking false branch
→
1978 mp_root(x, mp_cast_to_int(&reciprocal), z);
1979 else
1980 mp_pwr(x, y, z);
3
←
Calling 'mp_pwr'
→
1981 }
1982}
1983
1984
1985void
1986mp_xpowy_integer(const MPNumber *x, int64_t n, MPNumber *z)
1987{
1988 int i;
1989 MPNumber t;
1990
1991 /* 0^-n invalid */
1992 if (mp_is_zero(x) && n < 0) {
1993 /* Translators: Error displayed when attempted to raise 0 to a negative exponent */
1994 mperr(_("The power of zero is undefined for a negative exponent")gettext ("The power of zero is undefined for a negative exponent"
));
1995 mp_set_from_integer(0, z);
1996 return;
1997 }
1998
1999 /* x^0 = 1 */
2000 if (n == 0) {
2001 mp_set_from_integer(1, z);
2002 return;
2003 }
2004
2005 /* 0^n = 0 */
2006 if (mp_is_zero(x)) {
2007 mp_set_from_integer(0, z);
2008 return;
2009 }
2010
2011 /* x^1 = x */
2012 if (n == 1) {
2013 mp_set_from_mp(x, z);
2014 return;
2015 }
2016
2017 if (n < 0) {
2018 mp_reciprocal(x, &t);
2019 n = -n;
2020 }
2021 else
2022 mp_set_from_mp(x, &t);
2023
2024 /* Multply x n times */
2025 // FIXME: Can do mp_multiply(z, z, z) until close to answer (each call doubles number of multiples) */
2026 mp_set_from_integer(1, z);
2027 for (i = 0; i < n; i++)
2028 mp_multiply(z, &t, z);
2029}
2030
2031
2032GList*
2033mp_factorize(const MPNumber *x)
2034{
2035 GList *list = NULL((void*)0);
2036 MPNumber *factor = g_slice_alloc0(sizeof(MPNumber));
2037 MPNumber value, tmp, divisor, root;
2038
2039 mp_abs(x, &value);
2040
2041 if (mp_is_zero(&value)) {
2042 mp_set_from_mp(&value, factor);
2043 list = g_list_append(list, factor);
2044 return list;
2045 }
2046
2047 mp_set_from_integer(1, &tmp);
2048 if (mp_is_equal(&value, &tmp)) {
2049 mp_set_from_mp(x, factor);
2050 list = g_list_append(list, factor);
2051 return list;
2052 }
2053
2054 mp_set_from_integer(2, &divisor);
2055 while (TRUE(!(0))) {
2056 mp_divide(&value, &divisor, &tmp);
2057 if (mp_is_integer(&tmp)) {
2058 value = tmp;
2059 mp_set_from_mp(&divisor, factor);
2060 list = g_list_append(list, factor);
2061 factor = g_slice_alloc0(sizeof(MPNumber));
2062 } else {
2063 break;
2064 }
2065 }
2066
2067 mp_set_from_integer(3, &divisor);
2068 mp_sqrt(&value, &root);
2069 while (mp_is_less_equal(&divisor, &root)) {
2070 mp_divide(&value, &divisor, &tmp);
2071 if (mp_is_integer(&tmp)) {
2072 value = tmp;
2073 mp_sqrt(&value, &root);
2074 mp_set_from_mp(&divisor, factor);
2075 list = g_list_append(list, factor);
2076 factor = g_slice_alloc0(sizeof(MPNumber));
2077 } else {
2078 mp_add_integer(&divisor, 2, &tmp);
2079 divisor = tmp;
2080 }
2081 }
2082
2083 mp_set_from_integer(1, &tmp);
2084 if (mp_is_greater_than(&value, &tmp)) {
2085 mp_set_from_mp(&value, factor);
2086 list = g_list_append(list, factor);
2087 } else {
2088 g_slice_free(MPNumber, factor)do { if (1) g_slice_free1 (sizeof (MPNumber), (factor)); else
(void) ((MPNumber*) 0 == (factor)); } while (0);
2089 }
2090
2091 if (mp_is_negative(x)) {
2092 mp_invert_sign(list->data, list->data);
2093 }
2094
2095 return list;
2096}