5 #define MP_MAX_FREE 500
7 double convert_rational_to_double(U
*);
8 double convert_bignum_to_double(unsigned int *);
9 int ge(unsigned int *, unsigned int *, int);
11 int mtotal
, mfreecount
;
12 unsigned int **free_stack
;
20 if (n
== MP_MIN_SIZE
&& mfreecount
)
21 p
= free_stack
[--mfreecount
];
23 p
= (unsigned int *) malloc((n
+ 3) * sizeof (int));
25 stop("malloc failure");
33 mfree(unsigned int *p
)
37 if (p
[0] == MP_MIN_SIZE
&& mfreecount
< MP_MAX_FREE
)
38 free_stack
[mfreecount
++] = p
;
43 // convert int to bignum
48 unsigned int *p
= mnew(1);
61 mcopy(unsigned int *a
)
69 MLENGTH(b
) = MLENGTH(a
);
71 for (i
= 0; i
< MLENGTH(a
); i
++)
80 ge(unsigned int *a
, unsigned int *b
, int len
)
83 for (i
= len
- 1; i
> 0; i
--)
99 if (isrational(stack
[tos
- 1]) && isrational(stack
[tos
- 2])) {
112 a
= convert_rational_to_double(p1
);
117 b
= convert_rational_to_double(p2
);
125 subtract_numbers(void)
129 if (isrational(stack
[tos
- 1]) && isrational(stack
[tos
- 2])) {
142 a
= convert_rational_to_double(p1
);
147 b
= convert_rational_to_double(p2
);
155 multiply_numbers(void)
159 if (isrational(stack
[tos
- 1]) && isrational(stack
[tos
- 2])) {
172 a
= convert_rational_to_double(p1
);
177 b
= convert_rational_to_double(p2
);
189 if (isrational(stack
[tos
- 1]) && isrational(stack
[tos
- 2])) {
200 stop("divide by zero");
205 a
= convert_rational_to_double(p1
);
210 b
= convert_rational_to_double(p2
);
227 stop("divide by zero");
230 push_double(1 / p1
->u
.d
);
235 a
= mcopy(p1
->u
.q
.a
);
236 b
= mcopy(p1
->u
.q
.b
);
254 compare_rationals(U
*a
, U
*b
)
257 unsigned int *ab
, *ba
;
258 ab
= mmul(a
->u
.q
.a
, b
->u
.q
.b
);
259 ba
= mmul(a
->u
.q
.b
, b
->u
.q
.a
);
267 compare_numbers(U
*a
, U
*b
)
270 if (isrational(a
) && isrational(b
))
271 return compare_rationals(a
, b
);
275 x
= convert_rational_to_double(a
);
279 y
= convert_rational_to_double(b
);
301 p2
->u
.q
.a
= mcopy(p1
->u
.q
.a
);
302 p2
->u
.q
.b
= mcopy(p1
->u
.q
.b
);
303 MSIGN(p2
->u
.q
.a
) *= -1;
307 push_double(-p1
->u
.d
);
310 stop("bug caught in mp_negate_number");
317 bignum_truncate(void)
325 a
= mdiv(p1
->u
.q
.a
, p1
->u
.q
.b
);
356 p2
->u
.q
.a
= mcopy(p1
->u
.q
.a
);
381 p2
->u
.q
.a
= mcopy(p1
->u
.q
.b
);
390 bignum_power_number(int expo
)
392 unsigned int *a
, *b
, *t
;
398 a
= mpow(p1
->u
.q
.a
, abs(expo
));
399 b
= mpow(p1
->u
.q
.b
, abs(expo
));
422 convert_bignum_to_double(unsigned int *p
)
426 for (i
= MLENGTH(p
) - 1; i
>= 0; i
--)
427 d
= 4294967296.0 * d
+ p
[i
];
434 convert_rational_to_double(U
*p
)
437 double a
= 0.0, b
= 0.0;
439 na
= MLENGTH(p
->u
.q
.a
);
440 nb
= MLENGTH(p
->u
.q
.b
);
447 for (i
= 0; i
< n
; i
++) {
448 a
= a
/ 4294967296.0 + p
->u
.q
.a
[i
];
449 b
= b
/ 4294967296.0 + p
->u
.q
.b
[i
];
453 for (i
= nb
; i
< na
; i
++) {
454 a
= a
/ 4294967296.0 + p
->u
.q
.a
[i
];
455 b
= b
/ 4294967296.0;
459 for (i
= na
; i
< nb
; i
++) {
460 a
= a
/ 4294967296.0;
461 b
= b
/ 4294967296.0 + p
->u
.q
.b
[i
];
464 if (MSIGN(p
->u
.q
.a
) == -1)
483 push_double(double d
)
494 push_rational(int a
, int b
)
501 /* FIXME -- normalize */
518 if (isinteger(p1
) && MLENGTH(p1
->u
.q
.a
) == 1) {
523 n
*= MSIGN(p1
->u
.q
.a
);
530 if ((double) n
!= p1
->u
.d
)
544 print_double(U
*p
, int flag
)
547 sprintf(buf
, "%g", p
->u
.d
);
548 if (flag
== 1 && *buf
== '-')
555 bignum_scan_integer(char *s
)
564 if (sign
== '+' || sign
== '-')
584 #define FLT_MIN (1e-999)
585 #define FLT_MAX (9.999999999999999e999)
586 // the following strtod isn't used. but if I remove it, GCC's linker starts giving errors somewhere inside libm
587 // *magic* DO NOT TOUCH *magic*
588 // this strtod was deprecated because it didn't have enough precision
589 double amystrtod(char *s
, char **str_end
) {
590 // TODO handle exponential format, hex format, inf, nan
593 union { double rr
; unsigned long long dd
; } raw
;
595 while (isspace(*s
)) s
++; //skip leading whitespace
596 if ((s
[0] == 'i' || s
[0] == 'I') && (s
[1] == 'n' || s
[1] == 'N') && (s
[2] == 'f' || s
[2] == 'F')) {
599 raw
.dd
= 0x7FF0000000000000; //positive infinity
602 if ((s
[0] == 'n' || s
[0] == 'N') && (s
[1] == 'a' || s
[1] == 'A') && (s
[2] == 'n' || s
[2] == 'N')) {
605 raw
.dd
= 0x7FFFFFFFFFFFFFFF; //QNaN
608 if (!isdigit(*s
) && *s
!= '-' && *s
!= '+' && *s
!= '.') {
610 *str_end
= (char *)s
;
618 // Deliberate fallthrough
641 if (*s
== 'e' || *s
== 'E') {
644 if (*(s
+1) == '-' || *(s
+1) == '+') {
645 if (*(s
+1) == '-') sign
= true;
646 if (isdigit(*(s
+2))) {
649 exponent
= (exponent
*10) + (*s
- '0');
652 r
= pow(r
,((sign
)?-exponent
:exponent
));
654 } else if (isdigit(*(s
+1))) {
657 exponent
= (exponent
*10) + (*s
- '0');
665 *str_end
= (char *)s
;
667 // Portable? Nope. Fast? Yup.
669 raw
.dd
|= (unsigned long long)negative
<< 63;
671 if(raw
.rr
>= FLT_MIN
|| raw
.rr
<= -FLT_MIN
)
675 #define DBL_MAX_EXP 999
676 #define DBL_MIN_EXP (-999)
680 // Convert string to double
682 // Copyright (C) 2002 Michael Ringgaard. All rights reserved.
684 // Redistribution and use in source and binary forms, with or without
685 // modification, are permitted provided that the following conditions
688 // 1. Redistributions of source code must retain the above copyright
689 // notice, this list of conditions and the following disclaimer.
690 // 2. Redistributions in binary form must reproduce the above copyright
691 // notice, this list of conditions and the following disclaimer in the
692 // documentation and/or other materials provided with the distribution.
693 // 3. Neither the name of the project nor the names of its contributors
694 // may be used to endorse or promote products derived from this software
695 // without specific prior written permission.
697 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
698 // ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
699 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
700 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
701 // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
702 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
703 // OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
704 // HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
705 // LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
706 // OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
709 double mystrtod(const char *str
, char **endptr
) {
713 char *p
= (char *) str
;
719 // Skip leading whitespace
720 while (isspace(*p
)) p
++;
722 // Handle optional sign
725 case '-': negative
= 1; // Fall through to increment position
734 // Process string of digits
735 while (isdigit(*p
)) {
736 number
= number
* 10. + (*p
- '0');
741 // Process decimal part
745 while (isdigit(*p
)) {
746 number
= number
* 10. + (*p
- '0');
752 exponent
-= num_decimals
;
755 if (num_digits
== 0) {
761 if (negative
) number
= -number
;
763 // Process an exponent string
764 if (*p
== 'e' || *p
== 'E') {
765 // Handle optional sign
768 case '-': negative
= 1; // Fall through to increment pos
772 // Process string of digits
774 while (isdigit(*p
)) {
775 n
= n
* 10 + (*p
- '0');
786 if (exponent
< DBL_MIN_EXP
|| exponent
> DBL_MAX_EXP
) {
807 if (number
== HUGE_VAL
) errno
= ERANGE
;
808 if (endptr
) *endptr
= p
;
813 bignum_scan_float(char *s
)
815 //push_double(atof(s));
816 push_double(mystrtod((char*)s
, NULL
));
825 static char buf
[100];
829 if (*s
== '+' || *s
== '-')
839 sprintf(buf
, "%.10g", p
->u
.d
);
840 if (*buf
== '+' || *buf
== '-')
858 // if (!isinteger(p1) || !isinteger(p2))
859 // stop("integer args expected for gcd");
865 p3
->u
.q
.a
= mgcd(p1
->u
.q
.a
, p2
->u
.q
.a
);
866 p3
->u
.q
.b
= mgcd(p1
->u
.q
.b
, p2
->u
.q
.b
);
868 MSIGN(p3
->u
.q
.a
) = 1;
883 d
= convert_rational_to_double(p1
);
900 d
= convert_rational_to_double(pop());
904 static unsigned int *__factorial(int);
907 bignum_factorial(int n
)
912 p1
->u
.q
.a
= __factorial(n
);
918 static unsigned int *
922 unsigned int *a
, *b
, *t
;
924 if (n
== 0 || n
== 1) {
933 for (i
= 3; i
<= n
; i
++) {
934 b
[0] = (unsigned int) i
;
945 static unsigned int mask
[32] = {
981 mp_set_bit(unsigned int *x
, unsigned int k
)
983 x
[k
/ 32] |= mask
[k
% 32];
987 mp_clr_bit(unsigned int *x
, unsigned int k
)
989 x
[k
/ 32] &= ~mask
[k
% 32];
993 mshiftright(unsigned int *a
)
998 for (i
= n
- 1; i
>= 0; i
--)
1000 a
[i
] = (a
[i
] >> 1) | c
;
1003 a
[i
] = (a
[i
] >> 1) | c
;
1006 if (n
> 1 && a
[n
- 1] == 0)