#include <NTL/RR.h>
#include <NTL/new.h>
NTL_START_IMPL
long RR::prec = 150;
void RR::SetPrecision(long p)
{
if (p < 53)
p = 53;
if (NTL_OVERFLOW(p, 1, 0))
Error("RR: precision too high");
prec = p;
}
long RR::oprec = 10;
void RR::SetOutputPrecision(long p)
{
if (p < 1)
p = 1;
if (NTL_OVERFLOW(p, 1, 0))
Error("RR: output precision too high");
oprec = p;
}
static
void normalize1(RR& z, const ZZ& y_x, long y_e, long prec, long residual)
{
long len = NumBits(y_x);
if (len > prec) {
long correction = ZZ_RoundCorrection(y_x, len - prec, residual);
RightShift(z.x, y_x, len - prec);
if (correction)
add(z.x, z.x, correction);
z.e = y_e + len - prec;
}
else if (len == 0) {
clear(z.x);
z.e = 0;
}
else {
z.x = y_x;
z.e = y_e;
}
if (!IsOdd(z.x))
z.e += MakeOdd(z.x);
if (z.e >= NTL_OVFBND)
Error("RR: overflow");
if (z.e <= -NTL_OVFBND)
Error("RR: underflow");
}
void normalize(RR& z, const RR& y, long residual = 0)
{
normalize1(z, y.x, y.e, RR::prec, residual);
}
void MakeRR(RR& z, const ZZ& a, long e)
{
if (e >= NTL_OVFBND)
Error("MakeRR: e too big");
if (e <= -NTL_OVFBND)
Error("MakeRR: e too small");
normalize1(z, a, e, RR::prec, 0);
}
void MakeRRPrec(RR& x, const ZZ& a, long e, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("MakeRRPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
MakeRR(x, a, e);
RR::prec = old_p;
}
void random(RR& z)
{
static RR t;
RandomBits(t.x, RR::prec);
t.e = -RR::prec;
normalize(z, t);
}
static inline
void xcopy(RR& x, const RR& a)
{ normalize(x, a); }
// xcopy emulates old assignment semantics...
// many routines here implicitly assume assignment normalizes,
// but that is no longer the case as of v3.0.
void ConvPrec(RR& x, const RR& a, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("ConvPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
normalize(x, a);
RR::prec = old_p;
}
void RoundToPrecision(RR& x, const RR& a, long p)
{
ConvPrec(x, a, p);
}
void conv(RR& x, const RR& a)
{
normalize(x, a);
}
long IsZero(const RR& a)
{
return IsZero(a.x);
}
long IsOne(const RR& a)
{
return a.e == 0 && IsOne(a.x);
}
long sign(const RR& a)
{
return sign(a.x);
}
void clear(RR& z)
{
z.e = 0;
clear(z.x);
}
void set(RR& z)
{
z.e = 0;
set(z.x);
}
void add(RR& z, const RR& a, const RR& b)
{
static RR t;
if (IsZero(a.x)) {
xcopy(z, b);
return;
}
if (IsZero(b.x)) {
xcopy(z, a);
return;
}
if (a.e > b.e) {
if (a.e-b.e - max(RR::prec-NumBits(a.x),0) >= NumBits(b.x) + 2)
normalize(z, a, sign(b));
else {
LeftShift(t.x, a.x, a.e-b.e);
add(t.x, t.x, b.x);
t.e = b.e;
normalize(z, t);
}
}
else if (a.e < b.e) {
if (b.e-a.e - max(RR::prec-NumBits(b.x),0) >= NumBits(a.x) + 2)
normalize(z, b, sign(a));
else {
LeftShift(t.x, b.x, b.e-a.e);
add(t.x, t.x, a.x);
t.e = a.e;
normalize(z, t);
}
}
else {
add(t.x, a.x, b.x);
t.e = a.e;
normalize(z, t);
}
}
void AddPrec(RR& x, const RR& a, const RR& b, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("AddPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
add(x, a, b);
RR::prec = old_p;
}
void sub(RR& z, const RR& a, const RR& b)
{
static RR t;
if (IsZero(a.x)) {
negate(z, b);
return;
}
if (IsZero(b.x)) {
xcopy(z, a);
return;
}
if (a.e > b.e) {
if (a.e-b.e - max(RR::prec-NumBits(a.x),0) >= NumBits(b.x) + 2)
normalize(z, a, -sign(b));
else {
LeftShift(t.x, a.x, a.e-b.e);
sub(t.x, t.x, b.x);
t.e = b.e;
xcopy(z, t);
}
}
else if (a.e < b.e) {
if (b.e-a.e - max(RR::prec-NumBits(b.x),0) >= NumBits(a.x) + 2) {
normalize(z, b, -sign(a));
negate(z.x, z.x);
}
else {
LeftShift(t.x, b.x, b.e-a.e);
sub(t.x, a.x, t.x);
t.e = a.e;
xcopy(z, t);
}
}
else {
sub(t.x, a.x, b.x);
t.e = a.e;
normalize(z, t);
}
}
void SubPrec(RR& x, const RR& a, const RR& b, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("SubPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
sub(x, a, b);
RR::prec = old_p;
}
void negate(RR& z, const RR& a)
{
xcopy(z, a);
negate(z.x, z.x);
}
void NegatePrec(RR& x, const RR& a, const RR& b, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("NegatePrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
negate(x, a);
RR::prec = old_p;
}
void abs(RR& z, const RR& a)
{
xcopy(z, a);
abs(z.x, z.x);
}
void AbsPrec(RR& x, const RR& a, const RR& b, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("AbsPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
abs(x, a);
RR::prec = old_p;
}
void mul(RR& z, const RR& a, const RR& b)
{
static RR t;
mul(t.x, a.x, b.x);
t.e = a.e + b.e;
xcopy(z, t);
}
void MulPrec(RR& x, const RR& a, const RR& b, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("MulPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
mul(x, a, b);
RR::prec = old_p;
}
void sqr(RR& z, const RR& a)
{
static RR t;
sqr(t.x, a.x);
t.e = a.e + a.e;
xcopy(z, t);
}
void SqrPrec(RR& x, const RR& a, const RR& b, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("SqrPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
sqr(x, a);
RR::prec = old_p;
}
void div(RR& z, const RR& a, const RR& b)
{
if (IsZero(b))
Error("RR: division by zero");
if (IsZero(a)) {
clear(z);
return;
}
long la = NumBits(a.x);
long lb = NumBits(b.x);
long neg = (sign(a) != sign(b));
long k = RR::prec - la + lb + 1;
if (k < 0) k = 0;
static RR t;
static ZZ A, B, R;
abs(A, a.x);
LeftShift(A, A, k);
abs(B, b.x);
DivRem(t.x, R, A, B);
t.e = a.e - b.e - k;
normalize(z, t, !IsZero(R));
if (neg)
negate(z.x, z.x);
}
void DivPrec(RR& x, const RR& a, const RR& b, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("DivPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
div(x, a, b);
RR::prec = old_p;
}
void SqrRoot(RR& z, const RR& a)
{
if (sign(a) < 0)
Error("RR: attempt to take square root of negative number");
if (IsZero(a)) {
clear(z);
return;
}
RR t;
ZZ T1, T2;
long k;
k = 2*RR::prec - NumBits(a.x) + 1;
if (k < 0) k = 0;
if ((a.e - k) & 1) k++;
LeftShift(T1, a.x, k);
// since k >= 2*prec - bits(a) + 1, T1 has at least 2*prec+1 bits,
// thus T1 >= 2^(2*prec)
SqrRoot(t.x, T1); // t.x >= 2^prec thus t.x contains the round bit
t.e = (a.e - k)/2;
sqr(T2, t.x);
// T1-T2 is the (lower part of the) sticky bit
normalize(z, t, T2 < T1);
}
void SqrRootPrec(RR& x, const RR& a, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("SqrRootPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
SqrRoot(x, a);
RR::prec = old_p;
}
void swap(RR& a, RR& b)
{
swap(a.x, b.x);
swap(a.e, b.e);
}
long compare(const RR& a, const RR& b)
{
static RR t;
SubPrec(t, a, b, 1);
return sign(t);
}
long operator==(const RR& a, const RR& b)
{
return a.e == b.e && a.x == b.x;
}
void trunc(RR& z, const RR& a)
{
static RR t;
if (a.e >= 0)
xcopy(z, a);
else {
RightShift(t.x, a.x, -a.e);
t.e = 0;
xcopy(z, t);
}
}
void TruncPrec(RR& x, const RR& a, const RR& b, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("TruncPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
trunc(x, a);
RR::prec = old_p;
}
void floor(RR& z, const RR& a)
{
static RR t;
if (a.e >= 0)
xcopy(z, a);
else {
RightShift(t.x, a.x, -a.e);
if (sign(a.x) < 0)
add(t.x, t.x, -1);
t.e = 0;
xcopy(z, t);
}
}
void FloorPrec(RR& x, const RR& a, const RR& b, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("FloorPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
floor(x, a);
RR::prec = old_p;
}
void ceil(RR& z, const RR& a)
{
static RR t;
if (a.e >= 0)
xcopy(z, a);
else {
RightShift(t.x, a.x, -a.e);
if (sign(a.x) > 0)
add(t.x, t.x, 1);
t.e = 0;
xcopy(z, t);
}
}
void CeilPrec(RR& x, const RR& a, const RR& b, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("CeilPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
ceil(x, a);
RR::prec = old_p;
}
void round(RR& z, const RR& a)
{
if (a.e >= 0) {
xcopy(z, a);
return;
}
long len = NumBits(a.x);
if (-a.e > len) {
z = 0;
return;
}
if (-a.e == len) {
if (len == 1)
z = 0;
else
z = sign(a.x);
return;
}
static RR t;
ConvPrec(t, a, len+a.e);
xcopy(z, t);
}
void RoundPrec(RR& x, const RR& a, const RR& b, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("RoundPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
round(x, a);
RR::prec = old_p;
}
void conv(RR& z, const ZZ& a)
{
normalize1(z, a, 0, RR::prec, 0);
}
void ConvPrec(RR& x, const ZZ& a, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("ConvPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
conv(x, a);
RR::prec = old_p;
}
void conv(RR& z, long a)
{
if (a == 0) {
clear(z);
return;
}
if (a == 1) {
set(z);
return;
}
static ZZ t;
t = a;
conv(z, t);
}
void ConvPrec(RR& x, long a, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("ConvPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
conv(x, a);
RR::prec = old_p;
}
void conv(RR& z, unsigned long a)
{
if (a == 0) {
clear(z);
return;
}
if (a == 1) {
set(z);
return;
}
static ZZ t;
conv(t, a);
conv(z, t);
}
void ConvPrec(RR& x, unsigned long a, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("ConvPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
conv(x, a);
RR::prec = old_p;
}
void conv(RR& z, double a)
{
if (a == 0) {
clear(z);
return;
}
if (a == 1) {
set(z);
return;
}
if (!IsFinite(&a))
Error("RR: conversion of a non-finite double");
int e;
double f;
static RR t;
f = frexp(a, &e);
f = f * NTL_FDOUBLE_PRECISION;
f = f * 4;
conv(t.x, f);
t.e = e - (NTL_DOUBLE_PRECISION + 1);
xcopy(z, t);
}
void ConvPrec(RR& x, double a, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("ConvPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
conv(x, a);
RR::prec = old_p;
}
void conv(ZZ& z, const RR& a)
{
if (a.e >= 0)
LeftShift(z, a.x, a.e);
else {
long sgn = sign(a.x);
RightShift(z, a.x, -a.e);
if (sgn < 0)
sub(z, z, 1);
}
}
void CeilToZZ(ZZ& z, const RR& a)
{
if (a.e >= 0)
LeftShift(z, a.x, a.e);
else {
long sgn = sign(a.x);
RightShift(z, a.x, -a.e);
if (sgn > 0)
add(z, z, 1);
}
}
void TruncToZZ(ZZ& z, const RR& a)
{
if (a.e >= 0)
LeftShift(z, a.x, a.e);
else
RightShift(z, a.x, -a.e);
}
void RoundToZZ(ZZ& z, const RR& a)
{
if (a.e >= 0) {
LeftShift(z, a.x, a.e);
return;
}
long len = NumBits(a.x);
if (-a.e > len) {
z = 0;
return;
}
if (-a.e == len) {
if (len == 1)
z = 0;
else
z = sign(a.x);
return;
}
static RR t;
ConvPrec(t, a, len+a.e);
LeftShift(z, t.x, t.e);
}
void conv(long& z, const RR& a)
{
ZZ t;
conv(t, a);
conv(z, t);
}
void conv(double& z, const RR& aa)
{
double x;
static RR a;
ConvPrec(a, aa, NTL_DOUBLE_PRECISION);
// round to NTL_DOUBLE_PRECISION bits to avoid double overflow
conv(x, a.x);
z = _ntl_ldexp(x, a.e);
}
void add(RR& z, const RR& a, double b)
{
static RR B;
B = b;
add(z, a, B);
}
void sub(RR& z, const RR& a, double b)
{
static RR B;
B = b;
sub(z, a, B);
}
void sub(RR& z, double a, const RR& b)
{
static RR A;
A = a;
sub(z, A, b);
}
void mul(RR& z, const RR& a, double b)
{
static RR B;
B = b;
mul(z, a, B);
}
void div(RR& z, const RR& a, double b)
{
static RR B;
B = b;
div(z, a, B);
}
void div(RR& z, double a, const RR& b)
{
static RR A;
A = a;
div(z, A, b);
}
void inv(RR& z, const RR& a)
{
static RR one = to_RR(1);
div(z, one, a);
}
void InvPrec(RR& x, const RR& a, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("InvPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
inv(x, a);
RR::prec = old_p;
}
long compare(const RR& a, double b)
{
if (b == 0) return sign(a);
static RR B;
B = b;
return compare(a, B);
}
long operator==(const RR& a, double b)
{
if (b == 0) return IsZero(a);
if (b == 1) return IsOne(a);
static RR B;
B = b;
return a == B;
}
void power(RR& z, const RR& a, long e)
{
RR b, res;
long n = NumBits(e);
long p = RR::precision();
RR::SetPrecision(p + n + 10);
xcopy(b, a);
set(res);
long i;
for (i = n-1; i >= 0; i--) {
sqr(res, res);
if (bit(e, i))
mul(res, res, b);
}
RR::SetPrecision(p);
if (e < 0)
inv(z, res);
else
xcopy(z, res);
}
istream& operator>>(istream& s, RR& x)
{
long c;
long cval;
long sign;
ZZ a, b;
if (!s) Error("bad RR input");
c = s.peek();
while (IsWhiteSpace(c)) {
s.get();
c = s.peek();
}
if (c == '-') {
sign = -1;
s.get();
c = s.peek();
}
else
sign = 1;
long got1 = 0;
long got_dot = 0;
long got2 = 0;
a = 0;
b = 1;
cval = CharToIntVal(c);
if (cval >= 0 && cval <= 9) {
got1 = 1;
while (cval >= 0 && cval <= 9) {
mul(a, a, 10);
add(a, a, cval);
s.get();
c = s.peek();
cval = CharToIntVal(c);
}
}
if (c == '.') {
got_dot = 1;
s.get();
c = s.peek();
cval = CharToIntVal(c);
if (cval >= 0 && cval <= 9) {
got2 = 1;
while (cval >= 0 && cval <= 9) {
mul(a, a, 10);
add(a, a, cval);
mul(b, b, 10);
s.get();
c = s.peek();
cval = CharToIntVal(c);
}
}
}
if (got_dot && !got1 && !got2) Error("bad RR input");
ZZ e;
long got_e = 0;
long e_sign;
if (c == 'e' || c == 'E') {
got_e = 1;
s.get();
c = s.peek();
if (c == '-') {
e_sign = -1;
s.get();
c = s.peek();
}
else if (c == '+') {
e_sign = 1;
s.get();
c = s.peek();
}
else
e_sign = 1;
cval = CharToIntVal(c);
if (cval < 0 || cval > 9) Error("bad RR input");
e = 0;
while (cval >= 0 && cval <= 9) {
mul(e, e, 10);
add(e, e, cval);
s.get();
c = s.peek();
cval = CharToIntVal(c);
}
}
if (!got1 && !got2 && !got_e) Error("bad RR input");
RR t1, t2, v;
long old_p = RR::precision();
if (got1 || got2) {
ConvPrec(t1, a, max(NumBits(a), 1));
ConvPrec(t2, b, NumBits(b));
if (got_e)
RR::SetPrecision(old_p + 10);
div(v, t1, t2);
}
else
set(v);
if (sign < 0)
negate(v, v);
if (got_e) {
if (e >= NTL_OVFBND) Error("RR input overflow");
long E;
conv(E, e);
if (e_sign < 0) E = -E;
RR::SetPrecision(old_p + 10);
power(t1, to_RR(10), E);
mul(v, v, t1);
RR::prec = old_p;
}
xcopy(x, v);
return s;
}
void InputPrec(RR& x, istream& s, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("ConvPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
s >> x;
RR::prec = old_p;
}
void conv(RR& z, const xdouble& a)
{
conv(z, a.mantissa());
if (a.exponent() > ((2*NTL_OVFBND)/(2*NTL_XD_HBOUND_LOG)))
Error("RR: overlow");
if (a.exponent() < -((2*NTL_OVFBND)/(2*NTL_XD_HBOUND_LOG)))
Error("RR: underflow");
z.e += a.exponent()*(2*NTL_XD_HBOUND_LOG);
if (z.e >= NTL_OVFBND)
Error("RR: overflow");
if (z.e <= -NTL_OVFBND)
Error("RR: underflow");
}
void ConvPrec(RR& x, const xdouble& a, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("ConvPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
conv(x, a);
RR::prec = old_p;
}
void conv(xdouble& z, const RR& a)
{
xdouble x;
xdouble y;
conv(x, a.x);
power2(y, a.e);
z = x*y;
}
void power2(RR& z, long e)
{
if (e >= NTL_OVFBND)
Error("RR: overflow");
if (e <= -NTL_OVFBND)
Error("RR: underflow");
set(z.x);
z.e = e;
}
void conv(RR& z, const quad_float& a)
{
static RR hi, lo, res;
ConvPrec(hi, a.hi, NTL_DOUBLE_PRECISION);
ConvPrec(lo, a.lo, NTL_DOUBLE_PRECISION);
add(res, hi, lo);
z = res;
}
void ConvPrec(RR& x, const quad_float& a, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("ConvPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
conv(x, a);
RR::prec = old_p;
}
void conv(quad_float& z, const RR& a)
{
long old_p;
static RR a_hi, a_lo;
old_p = RR::prec;
ConvPrec(a_hi, a, NTL_DOUBLE_PRECISION); // high order bits
SubPrec(a_lo, a, a_hi, NTL_DOUBLE_PRECISION); // low order bits
z = to_quad_float(a_hi.x)*power2_quad_float(a_hi.e) +
to_quad_float(a_lo.x)*power2_quad_float(a_lo.e);
}
void conv(RR& x, const char *s)
{
long c;
long cval;
long sign;
ZZ a, b;
long i = 0;
if (!s) Error("bad RR input");
c = s[i];
while (IsWhiteSpace(c)) {
i++;
c = s[i];
}
if (c == '-') {
sign = -1;
i++;
c = s[i];
}
else
sign = 1;
long got1 = 0;
long got_dot = 0;
long got2 = 0;
a = 0;
b = 1;
cval = CharToIntVal(c);
if (cval >= 0 && cval <= 9) {
got1 = 1;
while (cval >= 0 && cval <= 9) {
mul(a, a, 10);
add(a, a, cval);
i++;
c = s[i];
cval = CharToIntVal(c);
}
}
if (c == '.') {
got_dot = 1;
i++;
c = s[i];
cval = CharToIntVal(c);
if (cval >= 0 && cval <= 9) {
got2 = 1;
while (cval >= 0 && cval <= 9) {
mul(a, a, 10);
add(a, a, cval);
mul(b, b, 10);
i++;
c = s[i];
cval = CharToIntVal(c);
}
}
}
if (got_dot && !got1 && !got2) Error("bad RR input");
ZZ e;
long got_e = 0;
long e_sign;
if (c == 'e' || c == 'E') {
got_e = 1;
i++;
c = s[i];
if (c == '-') {
e_sign = -1;
i++;
c = s[i];
}
else if (c == '+') {
e_sign = 1;
i++;
c = s[i];
}
else
e_sign = 1;
cval = CharToIntVal(c);
if (cval < 0 || cval > 9) Error("bad RR input");
e = 0;
while (cval >= 0 && cval <= 9) {
mul(e, e, 10);
add(e, e, cval);
i++;
c = s[i];
cval = CharToIntVal(c);
}
}
if (!got1 && !got2 && !got_e) Error("bad RR input");
RR t1, t2, v;
long old_p = RR::precision();
if (got1 || got2) {
ConvPrec(t1, a, max(NumBits(a), 1));
ConvPrec(t2, b, NumBits(b));
if (got_e)
RR::SetPrecision(old_p + 10);
div(v, t1, t2);
}
else
set(v);
if (sign < 0)
negate(v, v);
if (got_e) {
if (e >= NTL_OVFBND) Error("RR input overflow");
long E;
conv(E, e);
if (e_sign < 0) E = -E;
RR::SetPrecision(old_p + 10);
power(t1, to_RR(10), E);
mul(v, v, t1);
RR::prec = old_p;
}
xcopy(x, v);
}
void ConvPrec(RR& x, const char *s, long p)
{
if (p < 1 || NTL_OVERFLOW(p, 1, 0))
Error("ConvPrec: bad precsion");
long old_p = RR::prec;
RR::prec = p;
conv(x, s);
RR::prec = old_p;
}
void ReallyComputeE(RR& res)
{
long p = RR::precision();
RR::SetPrecision(p + NumBits(p) + 10);
RR s, s1, t;
s = 1;
t = 1;
long i;
for (i = 2; ; i++) {
add(s1, s, t);
if (s == s1) break;
xcopy(s, s1);
div(t, t, i);
}
RR::SetPrecision(p);
xcopy(res, s);
}
void ComputeE(RR& res)
{
static long prec = 0;
static RR e;
long p = RR::precision();
if (prec <= p + 10) {
prec = p + 20;
RR::SetPrecision(prec);
ReallyComputeE(e);
RR::SetPrecision(p);
}
xcopy(res, e);
}
void exp(RR& res, const RR& x)
{
if (x >= NTL_OVFBND || x <= -NTL_OVFBND)
Error("RR: overflow");
long p = RR::precision();
// step 0: write x = n + f, n an integer and |f| <= 1/2
// careful -- we want to compute f to > p bits of precision
RR f, nn;
RR::SetPrecision(NTL_BITS_PER_LONG);
round(nn, x);
RR::SetPrecision(p + 10);
sub(f, x, nn);
long n = to_long(nn);
// step 1: calculate t1 = e^n by repeated squaring
RR::SetPrecision(p + NumBits(n) + 10);
RR e;
ComputeE(e);
RR::SetPrecision(p + 10);
RR t1;
power(t1, e, n);
// step 2: calculate t2 = e^f using Taylor series expansion
RR::SetPrecision(p + NumBits(p) + 10);
RR t2, s, s1, t;
long i;
s = 0;
t = 1;
for (i = 1; ; i++) {
add(s1, s, t);
if (s == s1) break;
xcopy(s, s1);
mul(t, t, f);
div(t, t, i);
}
xcopy(t2, s);
RR::SetPrecision(p);
mul(res, t1, t2);
}
void ReallyComputeLn2(RR& res)
{
long p = RR::precision();
RR::SetPrecision(p + NumBits(p) + 10);
RR s, s1, t, t1;
s = 0;
t = 0.5;
t1 = 0.5;
long i;
for (i = 2; ; i++) {
add(s1, s, t);
if (s == s1) break;
xcopy(s, s1);
mul(t1, t1, 0.5);
div(t, t1, i);
}
RR::SetPrecision(p);
xcopy(res, s);
}
void ComputeLn2(RR& res)
{
static long prec = 0;
static RR ln2;
long p = RR::precision();
if (prec <= p + 10) {
prec = p + 20;
RR::SetPrecision(prec);
ReallyComputeLn2(ln2);
RR::SetPrecision(p);
}
xcopy(res, ln2);
}
long Lg2(const RR& x)
{
return NumBits(x.mantissa()) + x.exponent();
}
void log(RR& res, const RR& x)
{
if (x <= 0) Error("argument to log must be positive");
long p = RR::precision();
RR::SetPrecision(p + NumBits(p) + 10);
RR y;
long n;
// re-write x = 2^n * (1 - y), where -1/2 < y < 1/4 (so 3/4 < 1-y < 3/2)
if (x > 0.75 && x < 1.5) {
n = 0;
sub(y, 1, x);
}
else {
n = Lg2(x) - 1;
RR t;
power2(t, -n);
mul(t, t, x);
while (t > 1.5) {
mul(t, t, 0.5);
n++;
}
sub(y, 1, t);
}
// compute s = - ln(1-y) by power series expansion
RR s, s1, t, t1;
s = 0;
xcopy(t, y);
xcopy(t1, y);
long i;
for (i = 2; ; i++) {
add(s1, s, t);
if (s == s1) break;
xcopy(s, s1);
mul(t1, t1, y);
div(t, t1, i);
}
if (n == 0)
t = 0;
else {
ComputeLn2(t);
mul(t, t, n);
}
RR::SetPrecision(p);
sub(res, t, s);
}
void ComputeLn10(RR& res)
{
static long prec = 0;
static RR ln10;
long p = RR::precision();
if (prec <= p + 10) {
prec = p + 20;
RR::SetPrecision(prec);
log(ln10, to_RR(10));
RR::SetPrecision(p);
}
xcopy(res, ln10);
}
void log10(RR& res, const RR& x)
{
long p = RR::precision();
RR::SetPrecision(p + 10);
RR ln10, t1, t2;
ComputeLn10(ln10);
log(t1, x);
div(t2, t1, ln10);
RR::SetPrecision(p);
xcopy(res, t2);
}
void expm1(RR& res, const RR& x)
{
if (x < -0.5 || x > 0.5) {
RR t;
exp(t, x);
sub(res, t, 1);
return;
}
long p = RR::precision();
RR::SetPrecision(p + NumBits(p) + 10);
RR f;
xcopy(f, x);
RR s, s1, t;
long i;
s = 0;
xcopy(t, f);
for (i = 2; ; i++) {
add(s1, s, t);
if (s == s1) break;
xcopy(s, s1);
mul(t, t, f);
div(t, t, i);
}
RR::SetPrecision(p);
xcopy(res, s);
}
void log1p(RR& res, const RR& x)
{
long p = RR::precision();
RR y;
if (x < -0.5 || x > 0.5) {
RR::SetPrecision(p + 10);
log(y, 1 + x);
RR::SetPrecision(p);
xcopy(res, y);
return;
}
RR::SetPrecision(p + NumBits(p) + 10);
negate(y, x);
// compute s = - ln(1-y) by power series expansion
RR s, s1, t, t1;
s = 0;
xcopy(t, y);
xcopy(t1, y);
long i;
for (i = 2; ; i++) {
add(s1, s, t);
if (s == s1) break;
xcopy(s, s1);
mul(t1, t1, y);
div(t, t1, i);
}
RR::SetPrecision(p);
negate(res, s);
}
void pow(RR& res, const RR& x, const RR& y)
{
if (y == 0) {
res = 1;
return;
}
if (x == 0) {
res = 0;
return;
}
if (x == 1) {
res = 1;
return;
}
if (x < 0) {
Error("pow: sorry...first argument to pow must be nonnegative");
}
long p = RR::precision();
// calculate working precison...one could use p + NTL_BITS_PER_LONG + 10,
// for example, but we want the behaviour to be machine independent.
// so we calculate instead a rough approximation to log |y log(x)|
RR t1, t2;
long k;
if (x > 0.5 && x < 1.5) {
xcopy(t1, x - 1);
k = Lg2(t1);
}
else {
k = NumBits(Lg2(x));
}
k += Lg2(y);
if (k > NTL_BITS_PER_LONG+10) Error("RR: overflow");
if (k < 0) k = 0;
RR::SetPrecision(p + k + 10);
t1 = y*log(x);
RR::SetPrecision(p);
t2 = exp(t1);
res = t2;
}
void ReallyComputePi(RR& res)
{
long p = RR::precision();
RR::SetPrecision(p + NumBits(p) + 10);
RR sum1;
RR s, s1, t, t1;
s = 0;
t = 0.5;
t1 = 0.5;
long i;
for (i = 3; ; i+=2) {
add(s1, s, t);
if (s == s1) break;
xcopy(s, s1);
mul(t1, t1, -0.25);
div(t, t1, i);
}
xcopy(sum1, s);
RR g;
inv(g, to_RR(3)); // g = 1/3
s = 0;
xcopy(t, g);
xcopy(t1, g);
sqr(g, g);
negate(g, g); // g = -1/9
for (i = 3; ; i+=2) {
add(s1, s, t);
if (s == s1) break;
xcopy(s, s1);
mul(t1, t1, g);
div(t, t1, i);
}
add(s, s, sum1);
mul(s, s, 4);
RR::SetPrecision(p);
xcopy(res, s);
}
void ComputePi(RR& res)
{
static long prec = 0;
static RR pi;
long p = RR::precision();
if (prec <= p + 10) {
prec = p + 20;
RR::SetPrecision(prec);
ReallyComputePi(pi);
RR::SetPrecision(p);
}
xcopy(res, pi);
}
void sin(RR& res, const RR& x)
{
if (x == 0) {
res = 0;
return;
}
if (Lg2(x) > 1000)
Error("sin: sorry...argument too large in absolute value");
long p = RR::precision();
RR pi, t1, f;
RR n;
// we want to make x^2 < 3, so that the series for sin(x)
// converges nicely, without any nasty cancellations in the
// first terms of the series.
RR::SetPrecision(p + NumBits(p) + 10);
if (x*x < 3) {
xcopy(f, x);
}
else {
// we want to write x/pi = n + f, |f| < 1/2....
// but we have to do *this* very carefully, so that f is computed
// to precision > p. I know, this is sick!
long p1;
p1 = p + Lg2(x) + 20;
for (;;) {
RR::SetPrecision(p1);
ComputePi(pi);
xcopy(t1, x/pi);
xcopy(n, floor(t1));
xcopy(f, t1 - n);
if (f > 0.5) {
n++;
xcopy(f, t1 - n);
}
if (f == 0 || p1 < p - Lg2(f) + Lg2(n) + 10) {
// we don't have enough bits of f...increase p1 and continue
p1 = p1 + max(20, p1/10);
}
else
break;
}
RR::SetPrecision(p + NumBits(p) + 10);
ComputePi(pi);
xcopy(f, pi * f);
if (n != 0 && n.exponent() == 0) {
// n is odd, so we negate f, which negates sin(f)
xcopy(f, -f);
}
}
// Boy, that was painful, but now its over, and we can simply apply
// the series for sin(f)
RR t2, s, s1, t;
long i;
s = 0;
xcopy(t, f);
for (i = 3; ; i=i+2) {
add(s1, s, t);
if (s == s1) break;
xcopy(s, s1);
mul(t, t, f);
mul(t, t, f);
div(t, t, i-1);
div(t, t, i);
negate(t, t);
}
RR::SetPrecision(p);
xcopy(res, s);
}
void cos(RR& res, const RR& x)
{
if (x == 0) {
res = 1;
return;
}
if (Lg2(x) > 1000)
Error("cos: sorry...argument too large in absolute value");
long p = RR::precision();
RR pi, t1, f;
RR n;
// we want to write x/pi = (n+1/2) + f, |f| < 1/2....
// but we have to do *this* very carefully, so that f is computed
// to precision > p. I know, this is sick!
long p1;
p1 = p + Lg2(x) + 20;
for (;;) {
RR::SetPrecision(p1);
ComputePi(pi);
xcopy(t1, x/pi);
xcopy(n, floor(t1));
xcopy(f, t1 - (n + 0.5));
if (f == 0 || p1 < p - Lg2(f) + Lg2(n) + 10) {
// we don't have enough bits of f...increase p1 and continue
p1 = p1 + max(20, p1/10);
}
else
break;
}
RR::SetPrecision(p + NumBits(p) + 10);
ComputePi(pi);
xcopy(f, pi * f);
if (n == 0 || n.exponent() != 0) {
// n is even, so we negate f, which negates sin(f)
xcopy(f, -f);
}
// Boy, that was painful, but now its over, and we can simply apply
// the series for sin(f)
RR t2, s, s1, t;
long i;
s = 0;
xcopy(t, f);
for (i = 3; ; i=i+2) {
add(s1, s, t);
if (s == s1) break;
xcopy(s, s1);
mul(t, t, f);
mul(t, t, f);
div(t, t, i-1);
div(t, t, i);
negate(t, t);
}
RR::SetPrecision(p);
xcopy(res, s);
}
ostream& operator<<(ostream& s, const RR& a)
{
if (IsZero(a)) {
s << "0";
return s;
}
long old_p = RR::precision();
// we compute new_p and log_10_a precisely using sufficient
// precision---this is necessary to achieve accuracy and
// platform independent behaviour
long temp_p = max(NumBits(RR::OutputPrecision()),
NumBits(Lg2(a))) + 10;
RR::SetPrecision(temp_p);
RR ln2, ln10, log_2_10;
ComputeLn2(ln2);
ComputeLn10(ln10);
log_2_10 = ln10/ln2;
long new_p = to_long(RR::OutputPrecision()*log_2_10) + 20;
long log_10_a = to_long(Lg2(a)/log_2_10);
RR::SetPrecision(new_p);
RR b;
long neg;
if (a < 0) {
negate(b, a);
neg = 1;
}
else {
xcopy(b, a);
neg = 0;
}
long k = RR::OutputPrecision() - log_10_a;
RR c, d;
power(c, to_RR(10), RR::OutputPrecision());
power(d, to_RR(10), log_10_a);
div(b, b, d);
mul(b, b, c);
while (b < c) {
mul(b, b, 10);
k++;
}
while (b >= c) {
div(b, b, 10);
k--;
}
add(b, b, 0.5);
k = -k;
ZZ B;
conv(B, b);
long bp_len = RR::OutputPrecision()+10;
char *bp = NTL_NEW_OP char[bp_len];
if (!bp) Error("RR output: out of memory");
long len, i;
len = 0;
do {
if (len >= bp_len) Error("RR output: buffer overflow");
bp[len] = IntValToChar(DivRem(B, B, 10));
len++;
} while (B > 0);
for (i = 0; i < len/2; i++) {
char tmp;
tmp = bp[i];
bp[i] = bp[len-1-i];
bp[len-1-i] = tmp;
}
i = len-1;
while (bp[i] == '0') i--;
k += (len-1-i);
len = i+1;
bp[len] = '\0';
if (k > 3 || k < -len - 3) {
// use scientific notation
if (neg) s << "-";
s << "0." << bp << "e" << (k + len);
}
else if (k >= 0) {
if (neg) s << "-";
s << bp;
for (i = 0; i < k; i++)
s << "0";
}
else if (k <= -len) {
if (neg) s << "-";
s << "0.";
for (i = 0; i < -len-k; i++)
s << "0";
s << bp;
}
else {
if (neg) s << "-";
for (i = 0; i < len+k; i++)
s << bp[i];
s << ".";
for (i = len+k; i < len; i++)
s << bp[i];
}
RR::SetPrecision(old_p);
delete [] bp;
return s;
}
NTL_END_IMPL
syntax highlighted by Code2HTML, v. 0.9.1