1009 lines
16 KiB
C++
1009 lines
16 KiB
C++
#include "stdafx.h"
|
|
#include "defs.h"
|
|
|
|
#define MP_MIN_SIZE 2
|
|
#define MP_MAX_FREE 500
|
|
|
|
double convert_rational_to_double(U *);
|
|
double convert_bignum_to_double(unsigned int *);
|
|
int ge(unsigned int *, unsigned int *, int);
|
|
|
|
int mtotal, mfreecount;
|
|
unsigned int **free_stack;
|
|
|
|
unsigned int *
|
|
mnew(int n)
|
|
{
|
|
unsigned int *p;
|
|
if (n < MP_MIN_SIZE)
|
|
n = MP_MIN_SIZE;
|
|
if (n == MP_MIN_SIZE && mfreecount)
|
|
p = free_stack[--mfreecount];
|
|
else {
|
|
p = (unsigned int *) malloc((n + 3) * sizeof (int));
|
|
if (p == 0)
|
|
stop("malloc failure");
|
|
}
|
|
p[0] = n;
|
|
mtotal += n;
|
|
return p + 3;
|
|
}
|
|
|
|
void
|
|
mfree(unsigned int *p)
|
|
{
|
|
p -= 3;
|
|
mtotal -= p[0];
|
|
if (p[0] == MP_MIN_SIZE && mfreecount < MP_MAX_FREE)
|
|
free_stack[mfreecount++] = p;
|
|
else
|
|
free(p);
|
|
}
|
|
|
|
// convert int to bignum
|
|
|
|
unsigned int *
|
|
mint(int n)
|
|
{
|
|
unsigned int *p = mnew(1);
|
|
if (n < 0)
|
|
MSIGN(p) = -1;
|
|
else
|
|
MSIGN(p) = 1;
|
|
MLENGTH(p) = 1;
|
|
p[0] = abs(n);
|
|
return p;
|
|
}
|
|
|
|
// copy bignum
|
|
|
|
unsigned int *
|
|
mcopy(unsigned int *a)
|
|
{
|
|
int i;
|
|
unsigned int *b;
|
|
|
|
b = mnew(MLENGTH(a));
|
|
|
|
MSIGN(b) = MSIGN(a);
|
|
MLENGTH(b) = MLENGTH(a);
|
|
|
|
for (i = 0; i < MLENGTH(a); i++)
|
|
b[i] = a[i];
|
|
|
|
return b;
|
|
}
|
|
|
|
// a >= b ?
|
|
|
|
int
|
|
ge(unsigned int *a, unsigned int *b, int len)
|
|
{
|
|
int i;
|
|
for (i = len - 1; i > 0; i--)
|
|
if (a[i] == b[i])
|
|
continue;
|
|
else
|
|
break;
|
|
if (a[i] >= b[i])
|
|
return 1;
|
|
else
|
|
return 0;
|
|
}
|
|
|
|
void
|
|
add_numbers(void)
|
|
{
|
|
double a, b;
|
|
|
|
if (isrational(stack[tos - 1]) && isrational(stack[tos - 2])) {
|
|
qadd();
|
|
return;
|
|
}
|
|
|
|
save();
|
|
|
|
p2 = pop();
|
|
p1 = pop();
|
|
|
|
if (isdouble(p1))
|
|
a = p1->u.d;
|
|
else
|
|
a = convert_rational_to_double(p1);
|
|
|
|
if (isdouble(p2))
|
|
b = p2->u.d;
|
|
else
|
|
b = convert_rational_to_double(p2);
|
|
|
|
push_double(a + b);
|
|
|
|
restore();
|
|
}
|
|
|
|
void
|
|
subtract_numbers(void)
|
|
{
|
|
double a, b;
|
|
|
|
if (isrational(stack[tos - 1]) && isrational(stack[tos - 2])) {
|
|
qsub();
|
|
return;
|
|
}
|
|
|
|
save();
|
|
|
|
p2 = pop();
|
|
p1 = pop();
|
|
|
|
if (isdouble(p1))
|
|
a = p1->u.d;
|
|
else
|
|
a = convert_rational_to_double(p1);
|
|
|
|
if (isdouble(p2))
|
|
b = p2->u.d;
|
|
else
|
|
b = convert_rational_to_double(p2);
|
|
|
|
push_double(a - b);
|
|
|
|
restore();
|
|
}
|
|
|
|
void
|
|
multiply_numbers(void)
|
|
{
|
|
double a, b;
|
|
|
|
if (isrational(stack[tos - 1]) && isrational(stack[tos - 2])) {
|
|
qmul();
|
|
return;
|
|
}
|
|
|
|
save();
|
|
|
|
p2 = pop();
|
|
p1 = pop();
|
|
|
|
if (isdouble(p1))
|
|
a = p1->u.d;
|
|
else
|
|
a = convert_rational_to_double(p1);
|
|
|
|
if (isdouble(p2))
|
|
b = p2->u.d;
|
|
else
|
|
b = convert_rational_to_double(p2);
|
|
|
|
push_double(a * b);
|
|
|
|
restore();
|
|
}
|
|
|
|
void
|
|
divide_numbers(void)
|
|
{
|
|
double a, b;
|
|
|
|
if (isrational(stack[tos - 1]) && isrational(stack[tos - 2])) {
|
|
qdiv();
|
|
return;
|
|
}
|
|
|
|
save();
|
|
|
|
p2 = pop();
|
|
p1 = pop();
|
|
|
|
if (iszero(p2))
|
|
stop("divide by zero");
|
|
|
|
if (isdouble(p1))
|
|
a = p1->u.d;
|
|
else
|
|
a = convert_rational_to_double(p1);
|
|
|
|
if (isdouble(p2))
|
|
b = p2->u.d;
|
|
else
|
|
b = convert_rational_to_double(p2);
|
|
|
|
push_double(a / b);
|
|
|
|
restore();
|
|
}
|
|
|
|
void
|
|
invert_number(void)
|
|
{
|
|
unsigned int *a, *b;
|
|
|
|
save();
|
|
|
|
p1 = pop();
|
|
|
|
if (iszero(p1))
|
|
stop("divide by zero");
|
|
|
|
if (isdouble(p1)) {
|
|
push_double(1 / p1->u.d);
|
|
restore();
|
|
return;
|
|
}
|
|
|
|
a = mcopy(p1->u.q.a);
|
|
b = mcopy(p1->u.q.b);
|
|
|
|
MSIGN(b) = MSIGN(a);
|
|
MSIGN(a) = 1;
|
|
|
|
p1 = alloc();
|
|
|
|
p1->k = NUM;
|
|
|
|
p1->u.q.a = b;
|
|
p1->u.q.b = a;
|
|
|
|
push(p1);
|
|
|
|
restore();
|
|
}
|
|
|
|
int
|
|
compare_rationals(U *a, U *b)
|
|
{
|
|
int t;
|
|
unsigned int *ab, *ba;
|
|
ab = mmul(a->u.q.a, b->u.q.b);
|
|
ba = mmul(a->u.q.b, b->u.q.a);
|
|
t = mcmp(ab, ba);
|
|
mfree(ab);
|
|
mfree(ba);
|
|
return t;
|
|
}
|
|
|
|
int
|
|
compare_numbers(U *a, U *b)
|
|
{
|
|
double x, y;
|
|
if (isrational(a) && isrational(b))
|
|
return compare_rationals(a, b);
|
|
if (isdouble(a))
|
|
x = a->u.d;
|
|
else
|
|
x = convert_rational_to_double(a);
|
|
if (isdouble(b))
|
|
y = b->u.d;
|
|
else
|
|
y = convert_rational_to_double(b);
|
|
if (x < y)
|
|
return -1;
|
|
if (x > y)
|
|
return 1;
|
|
return 0;
|
|
}
|
|
|
|
void
|
|
negate_number(void)
|
|
{
|
|
save();
|
|
p1 = pop();
|
|
if (iszero(p1)) {
|
|
push(p1);
|
|
restore();
|
|
return;
|
|
}
|
|
switch (p1->k) {
|
|
case NUM:
|
|
p2 = alloc();
|
|
p2->k = NUM;
|
|
p2->u.q.a = mcopy(p1->u.q.a);
|
|
p2->u.q.b = mcopy(p1->u.q.b);
|
|
MSIGN(p2->u.q.a) *= -1;
|
|
push(p2);
|
|
break;
|
|
case DOUBLE:
|
|
push_double(-p1->u.d);
|
|
break;
|
|
default:
|
|
stop("bug caught in mp_negate_number");
|
|
break;
|
|
}
|
|
restore();
|
|
}
|
|
|
|
void
|
|
bignum_truncate(void)
|
|
{
|
|
unsigned int *a;
|
|
|
|
save();
|
|
|
|
p1 = pop();
|
|
|
|
a = mdiv(p1->u.q.a, p1->u.q.b);
|
|
|
|
p1 = alloc();
|
|
|
|
p1->k = NUM;
|
|
|
|
p1->u.q.a = a;
|
|
p1->u.q.b = mint(1);
|
|
|
|
push(p1);
|
|
|
|
restore();
|
|
}
|
|
|
|
void
|
|
mp_numerator(void)
|
|
{
|
|
save();
|
|
|
|
p1 = pop();
|
|
|
|
if (p1->k != NUM) {
|
|
push(one);
|
|
restore();
|
|
return;
|
|
}
|
|
|
|
p2 = alloc();
|
|
|
|
p2->k = NUM;
|
|
|
|
p2->u.q.a = mcopy(p1->u.q.a);
|
|
p2->u.q.b = mint(1);
|
|
|
|
push(p2);
|
|
|
|
restore();
|
|
}
|
|
|
|
void
|
|
mp_denominator(void)
|
|
{
|
|
save();
|
|
|
|
p1 = pop();
|
|
|
|
if (p1->k != NUM) {
|
|
push(one);
|
|
restore();
|
|
return;
|
|
}
|
|
|
|
p2 = alloc();
|
|
|
|
p2->k = NUM;
|
|
|
|
p2->u.q.a = mcopy(p1->u.q.b);
|
|
p2->u.q.b = mint(1);
|
|
|
|
push(p2);
|
|
|
|
restore();
|
|
}
|
|
|
|
void
|
|
bignum_power_number(int expo)
|
|
{
|
|
unsigned int *a, *b, *t;
|
|
|
|
save();
|
|
|
|
p1 = pop();
|
|
|
|
a = mpow(p1->u.q.a, abs(expo));
|
|
b = mpow(p1->u.q.b, abs(expo));
|
|
|
|
if (expo < 0) {
|
|
t = a;
|
|
a = b;
|
|
b = t;
|
|
MSIGN(a) = MSIGN(b);
|
|
MSIGN(b) = 1;
|
|
}
|
|
|
|
p1 = alloc();
|
|
|
|
p1->k = NUM;
|
|
|
|
p1->u.q.a = a;
|
|
p1->u.q.b = b;
|
|
|
|
push(p1);
|
|
|
|
restore();
|
|
}
|
|
|
|
double
|
|
convert_bignum_to_double(unsigned int *p)
|
|
{
|
|
int i;
|
|
double d = 0.0;
|
|
for (i = MLENGTH(p) - 1; i >= 0; i--)
|
|
d = 4294967296.0 * d + p[i];
|
|
if (MSIGN(p) == -1)
|
|
d = -d;
|
|
return d;
|
|
}
|
|
|
|
double
|
|
convert_rational_to_double(U *p)
|
|
{
|
|
int i, n, na, nb;
|
|
double a = 0.0, b = 0.0;
|
|
|
|
na = MLENGTH(p->u.q.a);
|
|
nb = MLENGTH(p->u.q.b);
|
|
|
|
if (na < nb)
|
|
n = na;
|
|
else
|
|
n = nb;
|
|
|
|
for (i = 0; i < n; i++) {
|
|
a = a / 4294967296.0 + p->u.q.a[i];
|
|
b = b / 4294967296.0 + p->u.q.b[i];
|
|
}
|
|
|
|
if (na > nb)
|
|
for (i = nb; i < na; i++) {
|
|
a = a / 4294967296.0 + p->u.q.a[i];
|
|
b = b / 4294967296.0;
|
|
}
|
|
|
|
if (na < nb)
|
|
for (i = na; i < nb; i++) {
|
|
a = a / 4294967296.0;
|
|
b = b / 4294967296.0 + p->u.q.b[i];
|
|
}
|
|
|
|
if (MSIGN(p->u.q.a) == -1)
|
|
a = -a;
|
|
|
|
return a / b;
|
|
}
|
|
|
|
void
|
|
push_integer(int n)
|
|
{
|
|
save();
|
|
p1 = alloc();
|
|
p1->k = NUM;
|
|
p1->u.q.a = mint(n);
|
|
p1->u.q.b = mint(1);
|
|
push(p1);
|
|
restore();
|
|
}
|
|
|
|
void
|
|
push_double(double d)
|
|
{
|
|
save();
|
|
p1 = alloc();
|
|
p1->k = DOUBLE;
|
|
p1->u.d = d;
|
|
push(p1);
|
|
restore();
|
|
}
|
|
|
|
void
|
|
push_rational(int a, int b)
|
|
{
|
|
save();
|
|
p1 = alloc();
|
|
p1->k = NUM;
|
|
p1->u.q.a = mint(a);
|
|
p1->u.q.b = mint(b);
|
|
/* FIXME -- normalize */
|
|
push(p1);
|
|
restore();
|
|
}
|
|
|
|
int
|
|
pop_integer(void)
|
|
{
|
|
int n;
|
|
|
|
save();
|
|
|
|
p1 = pop();
|
|
|
|
switch (p1->k) {
|
|
|
|
case NUM:
|
|
if (isinteger(p1) && MLENGTH(p1->u.q.a) == 1) {
|
|
n = p1->u.q.a[0];
|
|
if (n & 0x80000000)
|
|
n = 0x80000000;
|
|
else
|
|
n *= MSIGN(p1->u.q.a);
|
|
} else
|
|
n = 0x80000000;
|
|
break;
|
|
|
|
case DOUBLE:
|
|
n = (int) p1->u.d;
|
|
if ((double) n != p1->u.d)
|
|
n = 0x80000000;
|
|
break;
|
|
|
|
default:
|
|
n = 0x80000000;
|
|
break;
|
|
}
|
|
|
|
restore();
|
|
return n;
|
|
}
|
|
|
|
void
|
|
print_double(U *p, int flag)
|
|
{
|
|
static char buf[80];
|
|
sprintf(buf, "%g", p->u.d);
|
|
if (flag == 1 && *buf == '-')
|
|
print_str(buf + 1);
|
|
else
|
|
print_str(buf);
|
|
}
|
|
|
|
void
|
|
bignum_scan_integer(char *s)
|
|
{
|
|
unsigned int *a;
|
|
char sign;
|
|
|
|
save();
|
|
|
|
sign = *s;
|
|
|
|
if (sign == '+' || sign == '-')
|
|
s++;
|
|
|
|
a = mscan(s);
|
|
|
|
p1 = alloc();
|
|
|
|
p1->k = NUM;
|
|
|
|
p1->u.q.a = a;
|
|
p1->u.q.b = mint(1);
|
|
|
|
push(p1);
|
|
|
|
if (sign == '-')
|
|
negate();
|
|
|
|
restore();
|
|
}
|
|
|
|
#define FLT_MIN (1e-999)
|
|
#define FLT_MAX (9.999999999999999e999)
|
|
// the following strtod isn't used. but if I remove it, GCC's linker starts giving errors somewhere inside libm
|
|
// *magic* DO NOT TOUCH *magic*
|
|
// this strtod was deprecated because it didn't have enough precision
|
|
double amystrtod(char *s, char **str_end) {
|
|
// TODO handle exponential format, hex format, inf, nan
|
|
double r = 0.0;
|
|
int negative = 0;
|
|
union { double rr; unsigned long long dd; } raw;
|
|
|
|
while (isspace(*s)) s++; //skip leading whitespace
|
|
if ((s[0] == 'i' || s[0] == 'I') && (s[1] == 'n' || s[1] == 'N') && (s[2] == 'f' || s[2] == 'F')) {
|
|
if (str_end != NULL)
|
|
*str_end = s+3;
|
|
raw.dd = 0x7FF0000000000000; //positive infinity
|
|
return raw.rr;
|
|
}
|
|
if ((s[0] == 'n' || s[0] == 'N') && (s[1] == 'a' || s[1] == 'A') && (s[2] == 'n' || s[2] == 'N')) {
|
|
if (str_end != NULL)
|
|
*str_end = s+3;
|
|
raw.dd = 0x7FFFFFFFFFFFFFFF; //QNaN
|
|
return raw.rr;
|
|
}
|
|
if (!isdigit(*s) && *s != '-' && *s != '+' && *s != '.') {
|
|
if (str_end != NULL)
|
|
*str_end = (char *)s;
|
|
return 0;
|
|
}
|
|
|
|
switch (*s)
|
|
{
|
|
case '-':
|
|
negative = 1;
|
|
// Deliberate fallthrough
|
|
case '+':
|
|
s++;
|
|
break;
|
|
}
|
|
|
|
while (isdigit(*s))
|
|
{
|
|
r *= 10.0;
|
|
r += *s++ - '0';
|
|
}
|
|
if (*s == '.')
|
|
{
|
|
float f = 10.0f;
|
|
s++;
|
|
while (isdigit(*s))
|
|
{
|
|
r += (*s - '0')/f;
|
|
f *= 10.0f;
|
|
s++;
|
|
}
|
|
}
|
|
|
|
if (*s == 'e' || *s == 'E') {
|
|
int exponent = 0;
|
|
bool sign = false;
|
|
if (*(s+1) == '-' || *(s+1) == '+') {
|
|
if (*(s+1) == '-') sign = true;
|
|
if (isdigit(*(s+2))) {
|
|
s+=2;
|
|
while(isdigit(*s)) {
|
|
exponent = (exponent*10) + (*s - '0');
|
|
s++;
|
|
}
|
|
r = pow(r,((sign)?-exponent:exponent));
|
|
}
|
|
} else if (isdigit(*(s+1))) {
|
|
s++;
|
|
while(isdigit(*s)) {
|
|
exponent = (exponent*10) + (*s - '0');
|
|
s++;
|
|
}
|
|
r = pow(r,exponent);
|
|
}
|
|
}
|
|
|
|
if (str_end != NULL)
|
|
*str_end = (char *)s;
|
|
|
|
// Portable? Nope. Fast? Yup.
|
|
raw.rr = r;
|
|
raw.dd |= (unsigned long long)negative << 63;
|
|
|
|
if(raw.rr >= FLT_MIN || raw.rr <= -FLT_MIN)
|
|
return raw.rr;
|
|
else return 0;
|
|
}
|
|
#define DBL_MAX_EXP 999
|
|
#define DBL_MIN_EXP (-999)
|
|
//
|
|
// strtod.c
|
|
//
|
|
// Convert string to double
|
|
//
|
|
// Copyright (C) 2002 Michael Ringgaard. All rights reserved.
|
|
//
|
|
// Redistribution and use in source and binary forms, with or without
|
|
// modification, are permitted provided that the following conditions
|
|
// are met:
|
|
//
|
|
// 1. Redistributions of source code must retain the above copyright
|
|
// notice, this list of conditions and the following disclaimer.
|
|
// 2. Redistributions in binary form must reproduce the above copyright
|
|
// notice, this list of conditions and the following disclaimer in the
|
|
// documentation and/or other materials provided with the distribution.
|
|
// 3. Neither the name of the project nor the names of its contributors
|
|
// may be used to endorse or promote products derived from this software
|
|
// without specific prior written permission.
|
|
//
|
|
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
|
|
// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
|
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
|
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
|
|
// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
|
|
// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
|
|
// OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
|
|
// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
|
|
// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
|
|
// OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
|
|
// SUCH DAMAGE.
|
|
//
|
|
double mystrtod(const char *str, char **endptr) {
|
|
double number;
|
|
int exponent;
|
|
int negative;
|
|
char *p = (char *) str;
|
|
double p10;
|
|
int n;
|
|
int num_digits;
|
|
int num_decimals;
|
|
|
|
// Skip leading whitespace
|
|
while (isspace(*p)) p++;
|
|
|
|
// Handle optional sign
|
|
negative = 0;
|
|
switch (*p) {
|
|
case '-': negative = 1; // Fall through to increment position
|
|
case '+': p++;
|
|
}
|
|
|
|
number = 0.;
|
|
exponent = 0;
|
|
num_digits = 0;
|
|
num_decimals = 0;
|
|
|
|
// Process string of digits
|
|
while (isdigit(*p)) {
|
|
number = number * 10. + (*p - '0');
|
|
p++;
|
|
num_digits++;
|
|
}
|
|
|
|
// Process decimal part
|
|
if (*p == '.') {
|
|
p++;
|
|
|
|
while (isdigit(*p)) {
|
|
number = number * 10. + (*p - '0');
|
|
p++;
|
|
num_digits++;
|
|
num_decimals++;
|
|
}
|
|
|
|
exponent -= num_decimals;
|
|
}
|
|
|
|
if (num_digits == 0) {
|
|
errno = ERANGE;
|
|
return 0.0;
|
|
}
|
|
|
|
// Correct for sign
|
|
if (negative) number = -number;
|
|
|
|
// Process an exponent string
|
|
if (*p == 'e' || *p == 'E') {
|
|
// Handle optional sign
|
|
negative = 0;
|
|
switch (*++p) {
|
|
case '-': negative = 1; // Fall through to increment pos
|
|
case '+': p++;
|
|
}
|
|
|
|
// Process string of digits
|
|
n = 0;
|
|
while (isdigit(*p)) {
|
|
n = n * 10 + (*p - '0');
|
|
p++;
|
|
}
|
|
|
|
if (negative) {
|
|
exponent -= n;
|
|
} else {
|
|
exponent += n;
|
|
}
|
|
}
|
|
|
|
if (exponent < DBL_MIN_EXP || exponent > DBL_MAX_EXP) {
|
|
errno = ERANGE;
|
|
return HUGE_VAL;
|
|
}
|
|
|
|
// Scale the result
|
|
p10 = 10.;
|
|
n = exponent;
|
|
if (n < 0) n = -n;
|
|
while (n) {
|
|
if (n & 1) {
|
|
if (exponent < 0) {
|
|
number /= p10;
|
|
} else {
|
|
number *= p10;
|
|
}
|
|
}
|
|
n >>= 1;
|
|
p10 *= p10;
|
|
}
|
|
|
|
if (number == HUGE_VAL) errno = ERANGE;
|
|
if (endptr) *endptr = p;
|
|
|
|
return number;
|
|
}
|
|
void
|
|
bignum_scan_float(char *s)
|
|
{
|
|
//push_double(atof(s));
|
|
push_double(mystrtod((char*)s, NULL));
|
|
}
|
|
|
|
// print as unsigned
|
|
|
|
void
|
|
print_number(U *p)
|
|
{
|
|
char *s;
|
|
static char buf[100];
|
|
switch (p->k) {
|
|
case NUM:
|
|
s = mstr(p->u.q.a);
|
|
if (*s == '+' || *s == '-')
|
|
s++;
|
|
print_str(s);
|
|
if (isfraction(p)) {
|
|
print_str("/");
|
|
s = mstr(p->u.q.b);
|
|
print_str(s);
|
|
}
|
|
break;
|
|
case DOUBLE:
|
|
sprintf(buf, "%.10g", p->u.d);
|
|
if (*buf == '+' || *buf == '-')
|
|
print_str(buf + 1);
|
|
else
|
|
print_str(buf);
|
|
break;
|
|
default:
|
|
break;
|
|
}
|
|
}
|
|
|
|
void
|
|
gcd_numbers(void)
|
|
{
|
|
save();
|
|
|
|
p2 = pop();
|
|
p1 = pop();
|
|
|
|
// if (!isinteger(p1) || !isinteger(p2))
|
|
// stop("integer args expected for gcd");
|
|
|
|
p3 = alloc();
|
|
|
|
p3->k = NUM;
|
|
|
|
p3->u.q.a = mgcd(p1->u.q.a, p2->u.q.a);
|
|
p3->u.q.b = mgcd(p1->u.q.b, p2->u.q.b);
|
|
|
|
MSIGN(p3->u.q.a) = 1;
|
|
|
|
push(p3);
|
|
|
|
restore();
|
|
}
|
|
|
|
double
|
|
pop_double(void)
|
|
{
|
|
double d;
|
|
save();
|
|
p1 = pop();
|
|
switch (p1->k) {
|
|
case NUM:
|
|
d = convert_rational_to_double(p1);
|
|
break;
|
|
case DOUBLE:
|
|
d = p1->u.d;
|
|
break;
|
|
default:
|
|
d = 0.0;
|
|
break;
|
|
}
|
|
restore();
|
|
return d;
|
|
}
|
|
|
|
void
|
|
bignum_float(void)
|
|
{
|
|
double d;
|
|
d = convert_rational_to_double(pop());
|
|
push_double(d);
|
|
}
|
|
|
|
static unsigned int *__factorial(int);
|
|
|
|
void
|
|
bignum_factorial(int n)
|
|
{
|
|
save();
|
|
p1 = alloc();
|
|
p1->k = NUM;
|
|
p1->u.q.a = __factorial(n);
|
|
p1->u.q.b = mint(1);
|
|
push(p1);
|
|
restore();
|
|
}
|
|
|
|
static unsigned int *
|
|
__factorial(int n)
|
|
{
|
|
int i;
|
|
unsigned int *a, *b, *t;
|
|
|
|
if (n == 0 || n == 1) {
|
|
a = mint(1);
|
|
return a;
|
|
}
|
|
|
|
a = mint(2);
|
|
|
|
b = mint(0);
|
|
|
|
for (i = 3; i <= n; i++) {
|
|
b[0] = (unsigned int) i;
|
|
t = mmul(a, b);
|
|
mfree(a);
|
|
a = t;
|
|
}
|
|
|
|
mfree(b);
|
|
|
|
return a;
|
|
}
|
|
|
|
static unsigned int mask[32] = {
|
|
0x00000001,
|
|
0x00000002,
|
|
0x00000004,
|
|
0x00000008,
|
|
0x00000010,
|
|
0x00000020,
|
|
0x00000040,
|
|
0x00000080,
|
|
0x00000100,
|
|
0x00000200,
|
|
0x00000400,
|
|
0x00000800,
|
|
0x00001000,
|
|
0x00002000,
|
|
0x00004000,
|
|
0x00008000,
|
|
0x00010000,
|
|
0x00020000,
|
|
0x00040000,
|
|
0x00080000,
|
|
0x00100000,
|
|
0x00200000,
|
|
0x00400000,
|
|
0x00800000,
|
|
0x01000000,
|
|
0x02000000,
|
|
0x04000000,
|
|
0x08000000,
|
|
0x10000000,
|
|
0x20000000,
|
|
0x40000000,
|
|
0x80000000,
|
|
};
|
|
|
|
void
|
|
mp_set_bit(unsigned int *x, unsigned int k)
|
|
{
|
|
x[k / 32] |= mask[k % 32];
|
|
}
|
|
|
|
void
|
|
mp_clr_bit(unsigned int *x, unsigned int k)
|
|
{
|
|
x[k / 32] &= ~mask[k % 32];
|
|
}
|
|
|
|
void
|
|
mshiftright(unsigned int *a)
|
|
{
|
|
int c, i, n;
|
|
n = MLENGTH(a);
|
|
c = 0;
|
|
for (i = n - 1; i >= 0; i--)
|
|
if (a[i] & 1) {
|
|
a[i] = (a[i] >> 1) | c;
|
|
c = 0x80000000;
|
|
} else {
|
|
a[i] = (a[i] >> 1) | c;
|
|
c = 0;
|
|
}
|
|
if (n > 1 && a[n - 1] == 0)
|
|
MLENGTH(a) = n - 1;
|
|
}
|