2 * Copyright (C) 1987-2008 Sun Microsystems, Inc. All Rights Reserved.
3 * Copyright (C) 2008-2011 Robert Ancell
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
20 #include "mp-private.h"
23 mp_set_from_mp(const MPNumber
*x
, MPNumber
*z
)
26 memcpy(z
, x
, sizeof(MPNumber
));
31 mp_set_from_float(float rx
, MPNumber
*z
)
36 mp_set_from_integer(0, z
);
42 } else if (rx
> 0.0f
) {
46 /* IF RX = 0E0 RETURN 0 */
47 mp_set_from_integer(0, z
);
51 /* INCREASE IE AND DIVIDE RJ BY 16. */
57 while (rj
< 0.0625f
) {
62 /* NOW RJ IS DY DIVIDED BY SUITABLE POWER OF 16.
67 /* CONVERSION LOOP (ASSUME SINGLE-PRECISION OPS. EXACT) */
68 for (i
= 0; i
< MP_T
+ 4; i
++) {
69 rj
*= (float) MP_BASE
;
70 z
->fraction
[i
] = (int) rj
;
71 rj
-= (float) z
->fraction
[i
];
74 /* NORMALIZE RESULT */
78 ib
= max(MP_BASE
* 7 * MP_BASE
, 32767) / 16;
81 /* NOW MULTIPLY BY 16**IE */
84 for (i
= 1; i
<= k
; ++i
) {
86 if (tp
<= ib
&& tp
!= MP_BASE
&& i
< k
)
88 mp_divide_integer(z
, tp
, z
);
92 for (i
= 1; i
<= ie
; ++i
) {
94 if (tp
<= ib
&& tp
!= MP_BASE
&& i
< ie
)
96 mp_multiply_integer(z
, tp
, z
);
104 mp_set_from_double(double dx
, MPNumber
*z
)
106 int i
, k
, ib
, ie
, tp
;
109 mp_set_from_integer(0, z
);
115 } else if (dx
> 0.0) {
119 mp_set_from_integer(0, z
);
123 /* INCREASE IE AND DIVIDE DJ BY 16. */
124 for (ie
= 0; dj
>= 1.0; ie
++)
127 for ( ; dj
< 1.0/16.0; ie
--)
130 /* NOW DJ IS DY DIVIDED BY SUITABLE POWER OF 16
135 /* CONVERSION LOOP (ASSUME DOUBLE-PRECISION OPS. EXACT) */
136 for (i
= 0; i
< MP_T
+ 4; i
++) {
137 dj
*= (double) MP_BASE
;
138 z
->fraction
[i
] = (int) dj
;
139 dj
-= (double) z
->fraction
[i
];
142 /* NORMALIZE RESULT */
146 ib
= max(MP_BASE
* 7 * MP_BASE
, 32767) / 16;
149 /* NOW MULTIPLY BY 16**IE */
152 for (i
= 1; i
<= k
; ++i
) {
154 if (tp
<= ib
&& tp
!= MP_BASE
&& i
< k
)
156 mp_divide_integer(z
, tp
, z
);
160 for (i
= 1; i
<= ie
; ++i
) {
162 if (tp
<= ib
&& tp
!= MP_BASE
&& i
< ie
)
164 mp_multiply_integer(z
, tp
, z
);
172 mp_set_from_integer(int64_t x
, MPNumber
*z
)
176 memset(z
, 0, sizeof(MPNumber
));
191 z
->fraction
[z
->exponent
] = x
% MP_BASE
;
195 for (i
= 0; i
< z
->exponent
/ 2; i
++) {
196 int t
= z
->fraction
[i
];
197 z
->fraction
[i
] = z
->fraction
[z
->exponent
- i
- 1];
198 z
->fraction
[z
->exponent
- i
- 1] = t
;
204 mp_set_from_unsigned_integer(uint64_t x
, MPNumber
*z
)
208 mp_set_from_integer(0, z
);
217 z
->fraction
[z
->exponent
] = x
% MP_BASE
;
221 for (i
= 0; i
< z
->exponent
/ 2; i
++) {
222 int t
= z
->fraction
[i
];
223 z
->fraction
[i
] = z
->fraction
[z
->exponent
- i
- 1];
224 z
->fraction
[z
->exponent
- i
- 1] = t
;
230 mp_set_from_fraction(int64_t numerator
, int64_t denominator
, MPNumber
*z
)
232 mp_gcd(&numerator
, &denominator
);
234 if (denominator
== 0) {
235 mperr("*** J == 0 IN CALL TO MP_SET_FROM_FRACTION ***\n");
236 mp_set_from_integer(0, z
);
240 if (denominator
< 0) {
241 numerator
= -numerator
;
242 denominator
= -denominator
;
245 mp_set_from_integer(numerator
, z
);
246 if (denominator
!= 1)
247 mp_divide_integer(z
, denominator
, z
);
252 mp_set_from_polar(const MPNumber
*r
, MPAngleUnit unit
, const MPNumber
*theta
, MPNumber
*z
)
256 mp_cos(theta
, unit
, &x
);
257 mp_multiply(&x
, r
, &x
);
258 mp_sin(theta
, unit
, &y
);
259 mp_multiply(&y
, r
, &y
);
260 mp_set_from_complex(&x
, &y
, z
);
265 mp_set_from_complex(const MPNumber
*x
, const MPNumber
*y
, MPNumber
*z
)
267 /* NOTE: Do imaginary component first as z may be x or y */
268 z
->im_sign
= y
->sign
;
269 z
->im_exponent
= y
->exponent
;
270 memcpy(z
->im_fraction
, y
->fraction
, sizeof(int) * MP_SIZE
);
273 z
->exponent
= x
->exponent
;
275 memcpy(z
->fraction
, x
->fraction
, sizeof(int) * MP_SIZE
);
280 mp_set_from_random(MPNumber
*z
)
282 mp_set_from_double(drand48(), z
);
287 mp_cast_to_int(const MPNumber
*x
)
293 if (x
->sign
== 0 || x
->exponent
<= 0)
296 /* Multiply digits together */
297 for (i
= 0; i
< x
->exponent
; i
++) {
301 z
= z
* MP_BASE
+ x
->fraction
[i
];
303 /* Check for overflow */
308 /* Validate result */
310 for (i
= x
->exponent
- 1; i
>= 0; i
--) {
314 digit
= v
- (v
/ MP_BASE
) * MP_BASE
;
315 if (x
->fraction
[i
] != digit
)
328 mp_cast_to_unsigned_int(const MPNumber
*x
)
334 if (x
->sign
<= 0 || x
->exponent
<= 0)
337 /* Multiply digits together */
338 for (i
= 0; i
< x
->exponent
; i
++) {
342 z
= z
* MP_BASE
+ x
->fraction
[i
];
344 /* Check for overflow */
349 /* Validate result */
351 for (i
= x
->exponent
- 1; i
>= 0; i
--) {
355 digit
= v
- (v
/ MP_BASE
) * MP_BASE
;
356 if (x
->fraction
[i
] != digit
)
369 mppow_ri(float ap
, int bp
)
398 mp_cast_to_float(const MPNumber
*x
)
406 for (i
= 0; i
< MP_T
; i
++) {
407 rz
= (float) MP_BASE
* rz
+ (float)x
->fraction
[i
];
409 /* CHECK IF FULL SINGLE-PRECISION ACCURACY ATTAINED */
414 /* NOW ALLOW FOR EXPONENT */
415 rz
*= mppow_ri((float) MP_BASE
, x
->exponent
- i
- 1);
417 /* CHECK REASONABLENESS OF RESULT */
418 /* LHS SHOULD BE <= 0.5, BUT ALLOW FOR SOME ERROR IN ALOG */
419 if (rz
<= (float)0. ||
420 fabs((float) x
->exponent
- (log(rz
) / log((float) MP_BASE
) + (float).5)) > (float).6) {
421 /* FOLLOWING MESSAGE INDICATES THAT X IS TOO LARGE OR SMALL -
422 * TRY USING MPCMRE INSTEAD.
424 mperr("*** FLOATING-POINT OVER/UNDER-FLOW IN MP_CAST_TO_FLOAT ***\n");
436 mppow_di(double ap
, int bp
)
442 if (ap
== 0) return(pow
);
447 if (bp
& 01) pow
*= ap
;
448 if (bp
>>= 1) ap
*= ap
;
458 mp_cast_to_double(const MPNumber
*x
)
461 double d__1
, dz2
, ret_val
= 0.0;
466 for (i
= 0; i
< MP_T
; i
++) {
467 ret_val
= (double) MP_BASE
* ret_val
+ (double) x
->fraction
[i
];
470 /* CHECK IF FULL DOUBLE-PRECISION ACCURACY ATTAINED */
473 /* TEST BELOW NOT ALWAYS EQUIVALENT TO - IF (DZ2.LE.DZ) GO TO 20,
474 * FOR EXAMPLE ON CYBER 76.
476 if (dz2
- ret_val
<= 0.0)
480 /* NOW ALLOW FOR EXPONENT */
481 ret_val
*= mppow_di((double) MP_BASE
, x
->exponent
- tm
- 1);
483 /* CHECK REASONABLENESS OF RESULT. */
484 /* LHS SHOULD BE .LE. 0.5 BUT ALLOW FOR SOME ERROR IN DLOG */
486 ((d__1
= (double) ((float) x
->exponent
) - (log(ret_val
) / log((double)
487 ((float) MP_BASE
)) + .5), abs(d__1
)) > .6)) {
488 /* FOLLOWING MESSAGE INDICATES THAT X IS TOO LARGE OR SMALL -
489 * TRY USING MPCMDE INSTEAD.
491 mperr("*** FLOATING-POINT OVER/UNDER-FLOW IN MP_CAST_TO_DOUBLE ***\n");
503 char_val(char **c
, int base
)
505 int i
, j
, value
, offset
;
506 const char *digits
[][10] = {{"٠", "١", "٢", "٣", "٤", "٥", "٦", "٧", "٨", "٩"},
507 {"〇", "〡", "〢", "〣", "〤", "〥", "〦", "〧", "〨", "〩"},
508 {"۰", "۱", "۲", "۳", "۴", "۵", "۶", "۷", "۸", "۹"},
509 {"߀", "߁", "߂", "߃", "߄", "߅", "߆", "߇", "߈", "߉"},
510 {"०", "१", "२", "३", "४", "५", "६", "७", "८", "९"},
511 {"০", "১", "২", "৩", "৪", "৫", "৬", "৭", "৮", "৯"},
512 {"੦", "੧", "੨", "੩", "੪", "੫", "੬", "੭", "੮", "੯"},
513 {"૦", "૧", "૨", "૩", "૪", "૫", "૬", "૭", "૮", "૯"},
514 {"୦", "୧", "୨", "୩", "୪", "୫", "୬", "୭", "୮", "୯"},
515 {"௦", "௧", "௨", "௩", "௪", "௫", "௬", "௭", "௮", "௯"},
516 {"౦", "౧", "౨", "౩", "౪", "౫", "౬", "౭", "౮", "౯"},
517 {"೦", "೧", "೨", "೩", "೪", "೫", "೬", "೭", "೮", "೯"},
518 {"൦", "൧", "൨", "൩", "൪", "൫", "൬", "൭", "൮", "൯"},
519 {"๐", "๑", "๒", "๓", "๔", "๕", "๖", "๗", "๘", "๙"},
520 {"໐", "໑", "໒", "໓", "໔", "໕", "໖", "໗", "໘", "໙"},
521 {"༠", "༡", "༢", "༣", "༤", "༥", "༦", "༧", "༨", "༩"},
522 {"၀", "၁", "၂", "၃", "၄", "၅", "၆", "၇", "၈", "၉"},
523 {"႐", "႑", "႒", "႓", "႔", "႕", "႖", "႗", "႘", "႙"},
524 {"០", "១", "២", "៣", "៤", "៥", "៦", "៧", "៨", "៩"},
525 {"᠐", "᠑", "᠒", "᠓", "᠔", "᠕", "᠖", "᠗", "᠘", "᠙"},
526 {"᥆", "᥇", "᥈", "᥉", "᥊", "᥋", "᥌", "᥍", "᥎", "᥏"},
527 {"᧐", "᧑", "᧒", "᧓", "᧔", "᧕", "᧖", "᧗", "᧘", "᧙"},
528 {"᭐", "᭑", "᭒", "᭓", "᭔", "᭕", "᭖", "᭗", "᭘", "᭙"},
529 {"᮰", "᮱", "᮲", "᮳", "᮴", "᮵", "᮶", "᮷", "᮸", "᮹"},
530 {"᱀", "᱁", "᱂", "᱃", "᱄", "᱅", "᱆", "᱇", "᱈", "᱉"},
531 {"᱐", "᱑", "᱒", "᱓", "᱔", "᱕", "᱖", "᱗", "᱘", "᱙"},
532 {"꘠", "꘡", "꘢", "꘣", "꘤", "꘥", "꘦", "꘧", "꘨", "꘩"},
533 {"꣐", "꣑", "꣒", "꣓", "꣔", "꣕", "꣖", "꣗", "꣘", "꣙"},
534 {"꤀", "꤁", "꤂", "꤃", "꤄", "꤅", "꤆", "꤇", "꤈", "꤉"},
535 {"꩐", "꩑", "꩒", "꩓", "꩔", "꩕", "꩖", "꩗", "꩘", "꩙"},
536 {"𐒠", "𐒡", "𐒢", "𐒣", "𐒤", "𐒥", "𐒦", "𐒧", "𐒨", "𐒩"},
537 {NULL
, NULL
, NULL
, NULL
, NULL
, NULL
, NULL
, NULL
, NULL
, NULL
}};
539 if (**c
>= '0' && **c
<= '9') {
542 } else if (**c
>= 'a' && **c
<= 'f') {
543 value
= **c
- 'a' + 10;
545 } else if (**c
>= 'A' && **c
<= 'F') {
546 value
= **c
- 'A' + 10;
549 for (i
= 0; digits
[i
][0]; i
++) {
550 for (j
= 0; j
< 10; j
++) {
551 if (strncmp(*c
, digits
[i
][j
], strlen(digits
[i
][j
])) == 0)
557 if (digits
[i
][0] == NULL
)
560 offset
= strlen(digits
[i
][j
]);
572 ends_with(const char *start
, const char *end
, const char *word
)
574 size_t word_len
= strlen(word
);
576 if (word_len
> end
- start
)
579 return strncmp(end
- word_len
, word
, word_len
) == 0;
583 // FIXME: Doesn't handle errors well (e.g. trailing space)
585 set_from_sexagesimal(const char *str
, int length
, MPNumber
*z
)
587 int degrees
= 0, minutes
= 0;
588 char seconds
[length
+1];
593 n_matched
= sscanf(str
, "%d°%d'%s\"", °rees
, &minutes
, seconds
);
597 mp_set_from_integer(degrees
, z
);
599 mp_set_from_integer(minutes
, &t
);
600 mp_divide_integer(&t
, 60, &t
);
604 mp_set_from_string(seconds
, 10, &t
);
605 mp_divide_integer(&t
, 3600, &t
);
614 mp_set_from_string(const char *str
, int default_base
, MPNumber
*z
)
616 int i
, base
, negate
= 0, multiplier
= 0, base_multiplier
= 1;
618 gboolean has_fraction
= FALSE
;
620 const char *base_digits
[] = {"₀", "₁", "₂", "₃", "₄", "₅", "₆", "₇", "₈", "₉", NULL
};
621 const char *fractions
[] = {"½", "⅓", "⅔", "¼", "¾", "⅕", "⅖", "⅗", "⅘", "⅙", "⅚", "⅛", "⅜", "⅝", "⅞", NULL
};
622 int numerators
[] = { 1, 1, 2, 1, 3, 1, 2, 3, 4, 1, 5, 1, 3, 5, 7};
623 int denominators
[] = { 2, 3, 3, 4, 4, 5, 5, 5, 5, 6, 6, 8, 8, 8, 8};
625 if (strstr(str
, "°"))
626 return set_from_sexagesimal(str
, strlen(str
), z
);
634 for (i
= 0; base_digits
[i
] != NULL
; i
++) {
635 if (ends_with(str
, end
, base_digits
[i
])) {
636 base
+= i
* base_multiplier
;
637 end
-= strlen(base_digits
[i
]);
638 base_multiplier
*= 10;
642 if (base_digits
[i
] == NULL
)
645 if (base_multiplier
== 1)
648 /* Check if this has a sign */
652 } else if (*c
== '-') {
655 } else if (strncmp(c
, "−", strlen("−")) == 0) {
660 /* Convert integer part */
661 mp_set_from_integer(0, z
);
662 while ((i
= char_val((char **)&c
, base
)) >= 0) {
665 mp_multiply_integer(z
, base
, z
);
666 mp_add_integer(z
, i
, z
);
669 /* Look for fraction characters, e.g. ⅚ */
670 for (i
= 0; fractions
[i
] != NULL
; i
++) {
671 if (ends_with(str
, end
, fractions
[i
])) {
672 end
-= strlen(fractions
[i
]);
676 if (fractions
[i
] != NULL
) {
678 mp_set_from_fraction(numerators
[i
], denominators
[i
], &fraction
);
679 mp_add(z
, &fraction
, z
);
687 /* Convert fractional part */
689 MPNumber numerator
, denominator
;
691 mp_set_from_integer(0, &numerator
);
692 mp_set_from_integer(1, &denominator
);
693 while ((i
= char_val((char **)&c
, base
)) >= 0) {
694 mp_multiply_integer(&denominator
, base
, &denominator
);
695 mp_multiply_integer(&numerator
, base
, &numerator
);
696 mp_add_integer(&numerator
, i
, &numerator
);
698 mp_divide(&numerator
, &denominator
, &numerator
);
699 mp_add(z
, &numerator
, z
);
706 if (multiplier
!= 0) {
708 mp_set_from_integer(10, &t
);
709 mp_xpowy_integer(&t
, multiplier
, &t
);
710 mp_multiply(z
, &t
, z
);
714 mp_invert_sign(z
, z
);