Fix 0^n generating error for fractional n
[gcalctool.git] / src / mp-convert.c
blob1052351a9ce7cec9f5a817b9e0b876d20e2ca0d0
1 /*
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
9 * license.
12 #include <stdlib.h>
13 #include <stdio.h>
14 #include <string.h>
15 #include <ctype.h>
16 #include <math.h>
17 #include <langinfo.h>
19 #include "mp.h"
20 #include "mp-private.h"
22 void
23 mp_set_from_mp(const MPNumber *x, MPNumber *z)
25 if (z != x)
26 memcpy(z, x, sizeof(MPNumber));
30 void
31 mp_set_from_float(float rx, MPNumber *z)
33 int i, k, ib, ie, tp;
34 float rj;
36 mp_set_from_integer(0, z);
38 /* CHECK SIGN */
39 if (rx < 0.0f) {
40 z->sign = -1;
41 rj = -(double)(rx);
42 } else if (rx > 0.0f) {
43 z->sign = 1;
44 rj = rx;
45 } else {
46 /* IF RX = 0E0 RETURN 0 */
47 mp_set_from_integer(0, z);
48 return;
51 /* INCREASE IE AND DIVIDE RJ BY 16. */
52 ie = 0;
53 while (rj >= 1.0f) {
54 ++ie;
55 rj *= 0.0625f;
57 while (rj < 0.0625f) {
58 --ie;
59 rj *= 16.0f;
62 /* NOW RJ IS DY DIVIDED BY SUITABLE POWER OF 16.
63 * SET EXPONENT TO 0
65 z->exponent = 0;
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 */
75 mp_normalize(z);
77 /* Computing MAX */
78 ib = max(MP_BASE * 7 * MP_BASE, 32767) / 16;
79 tp = 1;
81 /* NOW MULTIPLY BY 16**IE */
82 if (ie < 0) {
83 k = -ie;
84 for (i = 1; i <= k; ++i) {
85 tp <<= 4;
86 if (tp <= ib && tp != MP_BASE && i < k)
87 continue;
88 mp_divide_integer(z, tp, z);
89 tp = 1;
91 } else if (ie > 0) {
92 for (i = 1; i <= ie; ++i) {
93 tp <<= 4;
94 if (tp <= ib && tp != MP_BASE && i < ie)
95 continue;
96 mp_multiply_integer(z, tp, z);
97 tp = 1;
103 void
104 mp_set_from_double(double dx, MPNumber *z)
106 int i, k, ib, ie, tp;
107 double dj;
109 mp_set_from_integer(0, z);
111 /* CHECK SIGN */
112 if (dx < 0.0) {
113 z->sign = -1;
114 dj = -dx;
115 } else if (dx > 0.0) {
116 z->sign = 1;
117 dj = dx;
118 } else {
119 mp_set_from_integer(0, z);
120 return;
123 /* INCREASE IE AND DIVIDE DJ BY 16. */
124 for (ie = 0; dj >= 1.0; ie++)
125 dj *= 1.0/16.0;
127 for ( ; dj < 1.0/16.0; ie--)
128 dj *= 16.0;
130 /* NOW DJ IS DY DIVIDED BY SUITABLE POWER OF 16
131 * SET EXPONENT TO 0
133 z->exponent = 0;
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 */
143 mp_normalize(z);
145 /* Computing MAX */
146 ib = max(MP_BASE * 7 * MP_BASE, 32767) / 16;
147 tp = 1;
149 /* NOW MULTIPLY BY 16**IE */
150 if (ie < 0) {
151 k = -ie;
152 for (i = 1; i <= k; ++i) {
153 tp <<= 4;
154 if (tp <= ib && tp != MP_BASE && i < k)
155 continue;
156 mp_divide_integer(z, tp, z);
157 tp = 1;
159 } else if (ie > 0) {
160 for (i = 1; i <= ie; ++i) {
161 tp <<= 4;
162 if (tp <= ib && tp != MP_BASE && i < ie)
163 continue;
164 mp_multiply_integer(z, tp, z);
165 tp = 1;
171 void
172 mp_set_from_integer(int64_t x, MPNumber *z)
174 int i;
176 memset(z, 0, sizeof(MPNumber));
178 if (x == 0) {
179 z->sign = 0;
180 return;
183 if (x < 0) {
184 x = -x;
185 z->sign = -1;
187 else if (x > 0)
188 z->sign = 1;
190 while (x != 0) {
191 z->fraction[z->exponent] = x % MP_BASE;
192 z->exponent++;
193 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;
203 void
204 mp_set_from_unsigned_integer(uint64_t x, MPNumber *z)
206 int i;
208 mp_set_from_integer(0, z);
210 if (x == 0) {
211 z->sign = 0;
212 return;
214 z->sign = 1;
216 while (x != 0) {
217 z->fraction[z->exponent] = x % MP_BASE;
218 x = x / MP_BASE;
219 z->exponent++;
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;
229 void
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);
237 return;
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);
251 void
252 mp_set_from_polar(const MPNumber *r, MPAngleUnit unit, const MPNumber *theta, MPNumber *z)
254 MPNumber x, y;
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);
264 void
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);
272 z->sign = x->sign;
273 z->exponent = x->exponent;
274 if (z != x)
275 memcpy(z->fraction, x->fraction, sizeof(int) * MP_SIZE);
279 void
280 mp_set_from_random(MPNumber *z)
282 mp_set_from_double(drand48(), z);
286 int64_t
287 mp_cast_to_int(const MPNumber *x)
289 int i;
290 int64_t z = 0, v;
292 /* |x| <= 1 */
293 if (x->sign == 0 || x->exponent <= 0)
294 return 0;
296 /* Multiply digits together */
297 for (i = 0; i < x->exponent; i++) {
298 int64_t t;
300 t = z;
301 z = z * MP_BASE + x->fraction[i];
303 /* Check for overflow */
304 if (z <= t)
305 return 0;
308 /* Validate result */
309 v = z;
310 for (i = x->exponent - 1; i >= 0; i--) {
311 int64_t digit;
313 /* Get last digit */
314 digit = v - (v / MP_BASE) * MP_BASE;
315 if (x->fraction[i] != digit)
316 return 0;
318 v /= MP_BASE;
320 if (v != 0)
321 return 0;
323 return x->sign * z;
327 uint64_t
328 mp_cast_to_unsigned_int(const MPNumber *x)
330 int i;
331 uint64_t z = 0, v;
333 /* x <= 1 */
334 if (x->sign <= 0 || x->exponent <= 0)
335 return 0;
337 /* Multiply digits together */
338 for (i = 0; i < x->exponent; i++) {
339 uint64_t t;
341 t = z;
342 z = z * MP_BASE + x->fraction[i];
344 /* Check for overflow */
345 if (z <= t)
346 return 0;
349 /* Validate result */
350 v = z;
351 for (i = x->exponent - 1; i >= 0; i--) {
352 uint64_t digit;
354 /* Get last digit */
355 digit = v - (v / MP_BASE) * MP_BASE;
356 if (x->fraction[i] != digit)
357 return 0;
359 v /= MP_BASE;
361 if (v != 0)
362 return 0;
364 return z;
368 static double
369 mppow_ri(float ap, int bp)
371 double pow;
373 if (bp == 0)
374 return 1.0;
376 if (bp < 0) {
377 if (ap == 0)
378 return 1.0;
379 bp = -bp;
380 ap = 1 / ap;
383 pow = 1.0;
384 for (;;) {
385 if (bp & 01)
386 pow *= ap;
387 if (bp >>= 1)
388 ap *= ap;
389 else
390 break;
393 return pow;
397 float
398 mp_cast_to_float(const MPNumber *x)
400 int i;
401 float rz = 0.0;
403 if (mp_is_zero(x))
404 return 0.0;
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 */
410 if (rz + 1.0f <= rz)
411 break;
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");
425 return 0.0;
428 if (x->sign < 0)
429 rz = -(double)(rz);
431 return rz;
435 static double
436 mppow_di(double ap, int bp)
438 double pow = 1.0;
440 if (bp != 0) {
441 if (bp < 0) {
442 if (ap == 0) return(pow);
443 bp = -bp;
444 ap = 1/ap;
446 for (;;) {
447 if (bp & 01) pow *= ap;
448 if (bp >>= 1) ap *= ap;
449 else break;
453 return(pow);
457 double
458 mp_cast_to_double(const MPNumber *x)
460 int i, tm = 0;
461 double d__1, dz2, ret_val = 0.0;
463 if (mp_is_zero(x))
464 return 0.0;
466 for (i = 0; i < MP_T; i++) {
467 ret_val = (double) MP_BASE * ret_val + (double) x->fraction[i];
468 tm = i;
470 /* CHECK IF FULL DOUBLE-PRECISION ACCURACY ATTAINED */
471 dz2 = ret_val + 1.0;
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)
477 break;
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 */
485 if (ret_val <= 0. ||
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");
492 return 0.0;
494 else
496 if (x->sign < 0)
497 ret_val = -ret_val;
498 return ret_val;
502 static int
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') {
540 value = **c - '0';
541 offset = 1;
542 } else if (**c >= 'a' && **c <= 'f') {
543 value = **c - 'a' + 10;
544 offset = 1;
545 } else if (**c >= 'A' && **c <= 'F') {
546 value = **c - 'A' + 10;
547 offset = 1;
548 } else {
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)
552 break;
554 if (j != 10)
555 break;
557 if (digits[i][0] == NULL)
558 return -1;
559 value = j;
560 offset = strlen(digits[i][j]);
562 if (value >= base)
563 return -1;
565 *c += offset;
567 return value;
571 static int
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)
577 return 0;
579 return strncmp(end - word_len, word, word_len) == 0;
583 // FIXME: Doesn't handle errors well (e.g. trailing space)
584 static bool
585 set_from_sexagesimal(const char *str, int length, MPNumber *z)
587 int degrees = 0, minutes = 0;
588 char seconds[length+1];
589 MPNumber t;
590 int n_matched;
592 seconds[0] = '\0';
593 n_matched = sscanf(str, "%d°%d'%s\"", &degrees, &minutes, seconds);
595 if (n_matched < 1)
596 return true;
597 mp_set_from_integer(degrees, z);
598 if (n_matched > 1) {
599 mp_set_from_integer(minutes, &t);
600 mp_divide_integer(&t, 60, &t);
601 mp_add(z, &t, z);
603 if (n_matched > 2) {
604 mp_set_from_string(seconds, 10, &t);
605 mp_divide_integer(&t, 3600, &t);
606 mp_add(z, &t, z);
609 return false;
613 bool
614 mp_set_from_string(const char *str, int default_base, MPNumber *z)
616 int i, base, negate = 0, multiplier = 0, base_multiplier = 1;
617 const char *c, *end;
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);
628 /* Find the base */
629 end = str;
630 while (*end != '\0')
631 end++;
632 base = 0;
633 while (1) {
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;
639 break;
642 if (base_digits[i] == NULL)
643 break;
645 if (base_multiplier == 1)
646 base = default_base;
648 /* Check if this has a sign */
649 c = str;
650 if (*c == '+') {
651 c++;
652 } else if (*c == '-') {
653 negate = 1;
654 c++;
655 } else if (strncmp(c, "−", strlen("−")) == 0) {
656 negate = 1;
657 c += strlen("−");
660 /* Convert integer part */
661 mp_set_from_integer(0, z);
662 while ((i = char_val((char **)&c, base)) >= 0) {
663 if (i > base)
664 return true;
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]);
673 break;
676 if (fractions[i] != NULL) {
677 MPNumber fraction;
678 mp_set_from_fraction(numerators[i], denominators[i], &fraction);
679 mp_add(z, &fraction, z);
682 if (*c == '.') {
683 has_fraction = TRUE;
684 c++;
687 /* Convert fractional part */
688 if (has_fraction) {
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);
702 if (c != end) {
703 return true;
706 if (multiplier != 0) {
707 MPNumber t;
708 mp_set_from_integer(10, &t);
709 mp_xpowy_integer(&t, multiplier, &t);
710 mp_multiply(z, &t, z);
713 if (negate == 1)
714 mp_invert_sign(z, z);
716 return false;