/* $Header$ */
/*
* Copyright © 1988-2004 Keith Packard and Bart Massey.
* All Rights Reserved. See the file COPYING in this directory
* for licensing information.
*/
/*
* rational.c
*
* operationalns on rationals
*/
#include "nickle.h"
#include <math.h>
int
RationalInit (void)
{
return 1;
}
#if 0
static Value
natural_to_rational (Natural *n)
{
ENTER ();
RETURN (NewRational (Positive, n, one_natural));
}
#endif
static Value
RationalPlusHelper (Sign sign, Rational *a, Rational *b)
{
ENTER ();
RETURN (NewRational (sign,
NaturalPlus (NaturalTimes (a->num, b->den),
NaturalTimes (b->num, a->den)),
NaturalTimes (a->den, b->den)));
}
static Value
RationalMinusHelper (Rational *a, Rational *b)
{
ENTER ();
Natural *ra, *rb, *t;
Sign sign = Positive;
ra = NaturalTimes (a->num, b->den);
rb = NaturalTimes (b->num, a->den);
if (NaturalLess (ra, rb))
{
sign = Negative;
t = ra;
ra = rb;
rb = t;
}
RETURN (NewRational (sign, NaturalMinus (ra, rb),
NaturalTimes (a->den, b->den)));
}
static Value
RationalPlus (Value av, Value bv, int expandOk)
{
ENTER ();
Rational *a = &av->rational, *b = &bv->rational;
Value ret;
switch (catagorize_signs(a->sign, b->sign)) {
case BothPositive:
case BothNegative:
ret = RationalPlusHelper (a->sign, a, b);
break;
case FirstPositive:
ret = RationalMinusHelper (a, b);
break;
case SecondPositive:
ret = RationalMinusHelper (b, a);
break;
default:
abort();
}
RETURN (ret);
}
static Value
RationalMinus (Value av, Value bv, int expandOk)
{
ENTER ();
Rational *a = &av->rational, *b = &bv->rational;
Value ret;
switch (catagorize_signs(a->sign, b->sign)) {
case BothPositive:
ret = RationalMinusHelper (a, b);
break;
case FirstPositive:
case SecondPositive:
ret = RationalPlusHelper (a->sign, a, b);
break;
case BothNegative:
ret = RationalMinusHelper (b, a);
break;
default:
abort();
}
RETURN (ret);
}
static Value
RationalTimes (Value av, Value bv, int expandOk)
{
ENTER ();
Rational *a = &av->rational, *b = &bv->rational;
Sign sign;
sign = Positive;
if (a->sign != b->sign)
sign = Negative;
RETURN (NewRational (sign,
NaturalTimes (a->num, b->num),
NaturalTimes (a->den, b->den)));
}
static Value
RationalDivide (Value av, Value bv, int expandOk)
{
ENTER ();
Rational *a = &av->rational, *b = &bv->rational;
Sign sign;
if (NaturalZero (b->num))
{
RaiseStandardException (exception_divide_by_zero,
"rational divide by zero",
2, av, bv);
RETURN (Void);
}
sign = Positive;
if (a->sign != b->sign)
sign = Negative;
RETURN (NewRational (sign,
NaturalTimes (a->num, b->den),
NaturalTimes (a->den, b->num)));
}
/*
* Modulus for rational values.
*
* Sorta like for integers:
*
* c/d * (a/b | c/d) + a/b % c/d = a/b
*
* 0 <= a/b % c/d < abs (c/d)
* a/b | c/d is an integer
*
* To calculate modulus (e/f):
*
* c/d * n + e/f = a/b
* e/f = a/b - c/d * n
* (e * b * d) / f = a * d - c * b * n
*
* therefore (e * b * d) / f is integer
*
* c * b * n + (e * b * d) / f = a * d
* (e * b * d) / f = (a * d) % (c * b)
* e / f = ((a * d) % (c * b)) / (b * d)
*/
static Value
RationalMod (Value av, Value bv, int expandOk)
{
ENTER ();
Rational *a = &av->rational, *b = &bv->rational;
Natural *rem, *quo, *div;
if (NaturalZero (b->num))
{
RaiseStandardException (exception_divide_by_zero,
"rational modulus by zero",
2, av, bv);
RETURN (Void);
}
div = NaturalTimes (b->num, a->den);
quo = NaturalDivide (NaturalTimes (a->num, b->den), div, &rem);
if (a->sign == Negative && !NaturalZero (rem))
rem = NaturalMinus (div, rem);
RETURN (NewRational (Positive, rem, NaturalTimes (a->den, b->den)));
}
static Value
RationalLess (Value av, Value bv, int expandOk)
{
ENTER ();
Rational *a = &av->rational, *b = &bv->rational;
Rational *t;
int ret;
switch (catagorize_signs (a->sign, b->sign)) {
case BothNegative:
t = a;
a = b;
b = t;
case BothPositive:
if (!NaturalEqual (a->den, b->den))
ret = NaturalLess (NaturalTimes (a->num, b->den),
NaturalTimes (b->num, a->den));
else
ret = NaturalLess (a->num, b->num);
break;
case FirstPositive:
ret = 0;
break;
case SecondPositive:
ret = 1;
break;
default:
abort();
}
RETURN (ret ? TrueVal : FalseVal);
}
static Value
RationalEqual (Value av, Value bv, int expandOk)
{
Rational *a = &av->rational, *b = &bv->rational;
if (a->sign == b->sign &&
NaturalEqual (a->num, b->num) &&
NaturalEqual (a->den, b->den))
{
return TrueVal;
}
return FalseVal;
}
static Value
RationalNegate (Value av, int expandOk)
{
ENTER ();
Rational *a = &av->rational;
RETURN (NewRational (SignNegate (a->sign), a->num, a->den));
}
static Value
RationalFloor (Value av, int expandOk)
{
ENTER ();
Rational *a = &av->rational;
Natural *quo, *rem;
quo = NaturalDivide (a->num, a->den, &rem);
if (!NaturalZero (rem) && a->sign == Negative)
quo = NaturalPlus (quo, one_natural);
RETURN (NewInteger (a->sign, quo));
}
static Value
RationalCeil (Value av, int expandOk)
{
ENTER ();
Rational *a = &av->rational;
Natural *quo, *rem;
quo = NaturalDivide (a->num, a->den, &rem);
if (!NaturalZero (rem) && a->sign == Positive)
quo = NaturalPlus (quo, one_natural);
RETURN (NewInteger (a->sign, quo));
}
static Value
RationalPromote (Value av, Value bv)
{
ENTER ();
switch (ValueTag(av)) {
case rep_int:
av = NewIntRational (ValueInt(av));
break;
case rep_integer:
av = NewIntegerRational (&av->integer);
break;
default:
break;
}
RETURN (av);
}
static Value
RationalReduce (Value av)
{
ENTER ();
Rational *a = &av->rational;
if (NaturalEqual (a->den, one_natural))
av = Reduce (NewInteger (a->sign, a->num));
RETURN (av);
}
static HashValue
RationalHash (Value av)
{
Rational *a = &av->rational;
return NaturalHash (a->den) ^ NaturalHash(a->num) ^ a->sign;
}
extern ValueRep IntegerRep;
extern Natural *NaturalFactor (Natural *, Natural *);
extern Natural *NaturalSqrt (Natural *);
extern Natural *NaturalIntPow (Natural *, int);
extern Natural *NaturalPow (Natural *, Natural *);
extern Natural *NaturalPowMod (Natural *, Natural *, Natural *);
extern Natural *two_natural;
static Natural *
NaturalPsi(Natural *a, Natural *max)
{
ENTER ();
Natural *p;
int n;
Natural *ret;
Natural *rem;
Natural *next;
Natural *pow;
Natural *fact;
ret = one_natural;
while (!NaturalEqual (a, one_natural))
{
p = NaturalFactor (a, max);
if (!p)
{
ret = 0;
break;
}
n = 0;
for (;;)
{
next = NaturalDivide (a, p, &rem);
if (!NaturalZero (rem))
break;
a = next;
n++;
}
pow = NaturalIntPow (p, n-1);
fact = NaturalMinus (NaturalTimes (pow, p), pow);
ret = NaturalTimes (ret, fact);
if (max && NaturalLess (max, fact))
break;
}
RETURN (ret);
}
#if 0
static int
IntSqrt (int a)
{
int l, h, m;
l = 2;
h = a/2;
while ((h-l) > 1)
{
m = (h+l) >> 1;
if (m * m < a)
l = m;
else
h = m;
}
return h;
}
static int
IntFactor (int a)
{
int v, lim;
if (!a)
return 0;
if ((a & 1) == 0)
return 2;
lim = IntSqrt (a);
for (v = 3; v <= lim; v += 2)
{
if (a % v == 0)
return v;
}
return a;
}
static int
IntPow (int a, int p)
{
int result;
result = 1;
while (p)
{
if (p & 1)
result = result * a;
p >>= 1;
if (p)
a = a * a;
}
return result;
}
static int
IntPowMod (int a, int p, int m)
{
int result;
if (m >= 32767)
{
#if DIGITBITS == 32
signed_digit la = a, lm = m, lr;
lr = 1;
while (p)
{
if (p & 1)
lr = (lr * la) % lm;
p >>= 1;
if (p)
la = (la * la) % lm;
}
result = (int) lr;
#else
ENTER ();
result = NaturalToInt (NaturalPowMod (NewNatural (a),
NewNatural (p),
NewNatural (m)));
EXIT ();
#endif
}
else
{
result = 1;
while (p)
{
if (p & 1)
result = (result * a) % m;
p >>= 1;
if (p)
a = (a * a) % m;
}
}
return result;
}
static int
IntPsi (int a)
{
int p;
int n;
int ret;
ret = 1;
while (a != 1)
{
p = IntFactor (a);
n = 0;
do
{
n++;
a /= p;
} while (a % p == 0);
ret = ret * (IntPow (p, n-1) * (p - 1));
}
return ret;
}
#endif
typedef struct _partial {
DataType *data;
struct _partial *down;
Natural *partial;
int power;
} Partial, *PartialPtr;
static void PartialMark (void *object)
{
PartialPtr p = object;
MemReference (p->partial);
MemReference (p->down);
}
DataType PartialType = { PartialMark, 0, "PartialType" };
static PartialPtr
NewPartial (Natural *partial)
{
ENTER ();
PartialPtr p;
if (!partial)
RETURN (0);
p = ALLOCATE (&PartialType, sizeof (Partial));
p->down = 0;
p->partial = partial;
p->power = 0;
RETURN (p);
}
typedef struct _factor {
DataType *data;
struct _factor *next;
Natural *prime;
int power;
PartialPtr partials;
} Factor, *FactorPtr;
static void FactorMark (void *object)
{
FactorPtr f = object;
MemReference (f->prime);
MemReference (f->next);
MemReference (f->partials);
}
DataType FactorType = { FactorMark, 0, "FactorType" };
static FactorPtr
NewFactor (Natural *prime, int power, FactorPtr next)
{
ENTER ();
FactorPtr f;
f = ALLOCATE (&FactorType, sizeof (Factor));
f->next = next;
f->prime = prime;
f->power = power;
f->partials = 0;
f->partials = NewPartial (prime);
f->partials->power = 0;
RETURN (f);
}
static FactorPtr
GenerateFactors (Natural *n, Natural *max)
{
ENTER ();
FactorPtr f = 0;
Natural *p;
Natural *largest;
Natural *d, *rem;
p = 0;
largest = NaturalSqrt (n);
while (!NaturalEqual (n, one_natural))
{
int power = 1;
for (;;)
{
if (!p)
p = two_natural;
else if (NaturalEqual (p, two_natural))
p = NewNatural (3);
else
p = NaturalPlus (p, two_natural);
d = NaturalDivide (n, p, &rem);
if (NaturalZero (rem))
break;
if (max && NaturalLess (max, p))
RETURN(f);
if (NaturalLess (largest, p))
RETURN (NewFactor (n, 1, f));
}
n = d;
for (;;)
{
d = NaturalDivide (n, p, &rem);
if (!NaturalZero (rem))
break;
n = d;
power++;
}
f = NewFactor (p, power, f);
largest = NaturalSqrt (n);
}
RETURN (f);
}
static Natural *
FactorBump (FactorPtr f)
{
PartialPtr p, minp;
Natural *factor;
ENTER ();
if (!f)
RETURN(0);
p = f->partials;
if (!p)
RETURN(0);
minp = p;
while (p->power)
{
if (!p->down)
p->down = NewPartial (FactorBump (f->next));
p = p->down;
if (!p)
break;
if (NaturalLess (p->partial, minp->partial))
minp = p;
}
if (!minp)
RETURN(0);
factor = minp->partial;
if (minp->power < f->power)
{
minp->partial = NaturalTimes (minp->partial, f->prime);
minp->power++;
}
else
{
f->partials = minp->down;
}
RETURN (factor);
}
static int
RationalRepeatLength (int prec, Natural *nden, int ibase)
{
ENTER ();
Natural *nbase;
Natural *ndigits;
FactorPtr factors;
Natural *factor;
int digits;
Natural *max = 0;
if (NaturalEqual (nden, one_natural))
return 0;
if (prec > 0)
max = NewNatural (prec);
nbase = NewNatural (ibase);
ndigits = NaturalPsi (nden, max);
if (!ndigits)
{
factor = one_natural;
for (factor = one_natural;;
factor = NaturalPlus (factor, one_natural))
{
if (NaturalEqual (NaturalPowMod (nbase, factor, nden),
one_natural))
break;
if (aborting)
break;
if (NaturalLess (max, factor))
{
EXIT ();
return -1;
}
}
}
else
{
factors = GenerateFactors (ndigits, max);
if (aborting)
return 0;
factor = one_natural;
while (factor)
{
if (NaturalEqual (NaturalPowMod (nbase, factor, nden),
one_natural))
break;
if (aborting)
break;
factor = FactorBump (factors);
if (max && factor && NaturalLess (max, factor))
{
EXIT ();
return -1;
}
}
}
if (!factor)
factor = ndigits;
if (NaturalLess (max_int_natural, factor))
factor = max_int_natural;
digits = NaturalToInt (factor);
EXIT ();
return digits;
}
static void
CheckDecimalLength (int prec, Natural *nden, int ibase, int *initial, int *repeat)
{
ENTER ();
Natural *rem;
Natural *nbase;
Natural *g;
int offset;
int rep;
nbase = NewNatural (ibase);
offset = 0;
while (!NaturalEqual ((g = NaturalGcd (nden, nbase)), one_natural))
{
if (aborting)
{
EXIT ();
return;
}
offset++;
if (prec >= 0 && offset > prec)
break;
nden = NaturalDivide (nden, g, &rem);
}
if (prec >= 0 && offset >= prec)
{
if (offset > prec)
offset = -prec;
else
offset = prec;
rep = 0;
}
else if (NaturalEqual (nden, one_natural))
{
rep = 0;
}
else
{
if (prec >= 0)
prec -= offset;
rep = RationalRepeatLength (prec, nden, ibase);
}
*initial = offset;
*repeat = rep;
EXIT ();
}
static Bool
RationalDecimalPrint (Value f, Value rv, char format, int base, int width, int prec, int fill)
{
ENTER ();
Rational *r = &rv->rational;
Natural *quo;
Natural *partial;
Natural *rep, *init;
Natural *dig;
int exponent = 0;
int exponent_width = 0;
char *initial = 0, *in;
char *repeat = 0, *re;
char *whole;
int initial_width, repeat_width = 0;
int frac_width;
int rep_width, brace_width = 0, dot_width = 0;
int whole_width;
int fraction_width;
int print_width;
int min_prec;
Bool use_braces = True;
min_prec = 0;
if (format == 'f' || format == 'e')
{
min_prec = prec;
use_braces = False;
}
if (prec == DEFAULT_OUTPUT_PRECISION)
{
min_prec = 0;
prec = 15;
}
else if (prec == INFINITE_OUTPUT_PRECISION)
prec = -1;
dig = NewNatural (base);
/*
* Check for small numbers for 'e' format
*/
if (NaturalLess (r->num, r->den))
{
Natural *quo, *rem;
Natural *mag;
int bits;
if (format == 'e' || (format == 'g' && prec > 0))
{
quo = NaturalDivide (r->den, r->num, &rem);
bits = NaturalWidth (quo);
exponent = (int) ((double) bits / (log ((double) base) / log (2.0)));
if (exponent < 0)
exponent = 0;
mag = NaturalIntPow (dig, exponent);
while (NaturalLess (mag, quo))
{
mag = NaturalTimes (mag, dig);
exponent++;
}
if (format == 'g' && prec > 0)
if (prec - exponent < 3)
format = 'e';
if (format == 'e')
{
int ev;
rv = RationalTimes (rv, NewRational (Positive, mag, one_natural), True);
r = &rv->rational;
exponent_width = 3;
ev = exponent;
while (ev >= base)
{
exponent_width++;
ev /= base;
}
exponent = -exponent;
}
else
exponent = 0;
}
else
exponent = 0;
}
CheckDecimalLength (prec, r->den, base, &initial_width, &repeat_width);
if (aborting)
{
EXIT ();
return False;
}
if ((rep_width = repeat_width))
{
/*
* When using %f format, just fill the
* result with digits
*/
if (!use_braces && prec != -1)
{
initial_width = -prec;
repeat_width = 0;
rep_width = 0;
}
else
{
if (repeat_width < 0)
rep_width = prec - initial_width;
}
}
if (initial_width)
{
Natural *half_digit;
if (initial_width < 0)
{
initial_width = -initial_width;
half_digit = NaturalTimes (NaturalIntPow (dig, initial_width),
two_natural);
rv = RationalPlusHelper (r->sign,
r,
&NewRational (Positive,
one_natural,
half_digit)->rational);
r = &rv->rational;
}
else
{
if (!repeat_width && initial_width < min_prec)
initial_width = min_prec;
}
initial = malloc (initial_width + 1);
if (!initial)
{
EXIT ();
return False;
}
}
quo = NaturalDivide (r->num, r->den, &partial);
whole = NaturalSprint (0, quo, base, &whole_width);
brace_width = 0;
if (repeat_width)
{
brace_width++;
if (repeat_width > 0)
brace_width++;
}
dot_width = 0;
if (initial_width + rep_width)
dot_width = 1;
/*
* Compute how much space is available for the fractional part
*/
if (width)
{
if (width < 0)
fraction_width = -width;
else
fraction_width = width;
fraction_width = fraction_width - (whole_width + exponent_width);
if (fraction_width < 0)
fraction_width = 0;
if (prec > 0 && fraction_width > prec + dot_width)
fraction_width = prec + dot_width;
}
else if (prec > 0)
fraction_width = prec + dot_width;
else
fraction_width = -1;
/*
* Start paring down parts of the output to fit the desired size
*/
while (fraction_width >= 0 &&
(frac_width = dot_width + initial_width + rep_width + brace_width)
&& frac_width > fraction_width)
{
if (rep_width)
{
if (brace_width > 1)
{
brace_width = 1;
repeat_width = -repeat_width;
}
rep_width = fraction_width - (dot_width + initial_width +
brace_width);
if (rep_width < 0)
{
rep_width = 0;
}
}
else if (brace_width)
brace_width = 0;
else if (initial_width)
{
initial_width = fraction_width - dot_width;
if (initial_width < 0)
initial_width = 0;
}
else
dot_width = 0;
}
if (initial_width)
{
init = NaturalDivide (NaturalTimes (partial,
NaturalIntPow (dig, initial_width)),
r->den,
&partial);
if (aborting)
{
free (initial);
EXIT ();
return False;
}
in = NaturalSprint (initial + initial_width + 1,
init, base, &initial_width);
if (!in)
{
free (initial);
EXIT ();
return False;
}
while (in > initial)
{
*--in = '0';
++initial_width;
}
}
if (rep_width)
{
#define MAX_SENSIBLE 10000000
if (rep_width > MAX_SENSIBLE)
{
repeat_width = -1;
rep_width = MAX_SENSIBLE;
}
/*
* allocate the output buffer; keep trying until this works
*/
while (!(repeat = malloc (rep_width + 1)))
{
repeat_width = -1;
rep_width >>= 1;
}
rep = NaturalDivide (NaturalTimes (partial,
NaturalIntPow (dig, rep_width)),
r->den,
&partial);
if (aborting)
{
free (initial);
free (repeat);
EXIT ();
return False;
}
re = NaturalSprint (repeat + rep_width + 1,
rep, base, &rep_width);
if (!re)
{
free (initial);
free (repeat);
EXIT ();
return False;
}
while (re > repeat)
{
*--re = '0';
++rep_width;
}
if (use_braces)
{
rep_width++; /* open { */
if (repeat_width > 0)
rep_width++; /* close } */
}
}
fraction_width = initial_width + rep_width;
print_width = whole_width + 1 + fraction_width + exponent_width;
if (r->sign == Negative)
print_width = print_width + 1;
while (width > print_width)
{
FileOutchar (f, fill);
width--;
}
if (r->sign == Negative)
FileOutput (f, '-');
FilePuts (f, whole);
FileOutput (f, '.');
if (initial_width)
{
FilePuts (f, initial);
free (initial);
}
if (rep_width)
{
if (use_braces)
FileOutput (f, '{');
FilePuts (f, repeat);
free (repeat);
if (use_braces && repeat_width > 0)
FileOutput (f, '}');
}
if (exponent)
{
FilePrintf (f, "e%d", exponent);
}
while (-width > print_width)
{
FileOutchar (f, fill);
width++;
}
EXIT ();
return True;
}
static Bool
RationalPrint (Value f, Value rv, char format, int base, int width, int prec, int fill)
{
Rational *r = &rv->rational;
char *num, *num_base, *den, *den_base;
int num_width, den_width;
int print_width;
Bool ret = True;
if (base == 0)
base = 10;
switch (format) {
case 'v':
num_width = NaturalEstimateLength (r->num, base);
num_base = malloc (num_width);
num = NaturalSprint (num_base + num_width, r->num, base, &num_width);
if (!num)
{
free (num_base);
ret = False;
break;
}
den_width = NaturalEstimateLength (r->den, base);
den_base = malloc (den_width);
den = NaturalSprint (den_base + den_width, r->den, base, &den_width);
if (!den)
{
free (num_base);
free (den_base);
ret = False;
break;
}
print_width = 1 + num_width + 1 + den_width + 1;
if (r->sign == Negative)
print_width++;
while (width > print_width)
{
FileOutchar (f, fill);
width--;
}
FileOutput (f, '(');
if (r->sign == Negative)
FileOutput (f, '-');
FilePuts (f, num);
FileOutput (f, '/');
FilePuts (f, den);
FileOutput (f, ')');
free (num_base);
free (den_base);
while (-width > print_width)
{
FileOutchar (f, fill);
width++;
}
break;
default:
ret = RationalDecimalPrint (f, rv, format, base, width, prec, fill);
break;
}
return ret;
}
static void
RationalMark (void *object)
{
Rational *rational = object;
MemReference (rational->num);
MemReference (rational->den);
}
ValueRep RationalRep = {
{ RationalMark, 0, "RationalRep" }, /* base */
rep_rational, /* tag */
{ /* binary */
RationalPlus,
RationalMinus,
RationalTimes,
RationalDivide,
NumericDiv,
RationalMod,
RationalLess,
RationalEqual,
0,
0,
},
{ /* unary */
RationalNegate,
RationalFloor,
RationalCeil,
},
RationalPromote,
RationalReduce,
RationalPrint,
0,
RationalHash,
};
Value
NewRational (Sign sign, Natural *num, Natural *den)
{
ENTER ();
Value ret;
Natural *g;
Natural *rem;
if (NaturalZero (num))
den = one_natural;
else
{
if (NaturalLength(den) != 1 || NaturalDigits(den)[0] != 1)
{
g = NaturalGcd (num, den);
if (NaturalLength (g) != 1 || NaturalDigits(g)[0] != 1)
{
num = NaturalDivide (num, g, &rem);
den = NaturalDivide (den, g, &rem);
}
}
}
ret = ALLOCATE (&RationalRep.data, sizeof (Rational));
ret->rational.sign = sign;
ret->rational.num = num;
ret->rational.den = den;
RETURN (ret);
}
Value
NewIntRational (int i)
{
ENTER ();
if (i < 0)
RETURN (NewRational (Negative, NewNatural ((unsigned) -i), one_natural));
else
RETURN (NewRational (Positive, NewNatural ((unsigned) i), one_natural));
}
Value
NewIntegerRational (Integer *i)
{
ENTER ();
RETURN (NewRational (IntegerSign((Value) i), IntegerMag((Value) i),
one_natural));
}
syntax highlighted by Code2HTML, v. 0.9.1