From 93a007616a0ea689db8cb0153122576821e0a75f Mon Sep 17 00:00:00 2001 From: Clement Jamin Date: Mon, 15 Dec 2014 21:34:29 +0100 Subject: [PATCH] Header-only CGAL_Core --- CGAL_Core/include/CGAL/CORE/BigFloat.h | 33 +- CGAL_Core/include/CGAL/CORE/BigFloatRep.h | 13 +- CGAL_Core/include/CGAL/CORE/BigFloat_impl.h | 1325 +++++++++++++++++ CGAL_Core/include/CGAL/CORE/BigInt.h | 12 + CGAL_Core/include/CGAL/CORE/CoreAux.h | 12 + CGAL_Core/include/CGAL/CORE/CoreAux_impl.h | 229 +++ CGAL_Core/include/CGAL/CORE/CoreDefs.h | 125 +- CGAL_Core/include/CGAL/CORE/CoreDefs_impl.h | 170 +++ CGAL_Core/include/CGAL/CORE/CoreIO_impl.h | 488 ++++++ CGAL_Core/include/CGAL/CORE/Expr.h | 50 +- CGAL_Core/include/CGAL/CORE/ExprRep.h | 53 +- CGAL_Core/include/CGAL/CORE/Expr_impl.h | 1213 +++++++++++++++ CGAL_Core/include/CGAL/CORE/Filter.h | 9 +- CGAL_Core/include/CGAL/CORE/Gmp.h | 12 + CGAL_Core/include/CGAL/CORE/Gmp_impl.h | 280 ++++ CGAL_Core/include/CGAL/CORE/Real.h | 29 +- CGAL_Core/include/CGAL/CORE/extLong.h | 12 + CGAL_Core/include/CGAL/CORE/poly/Poly.h | 3 +- CGAL_Core/src/CGAL_Core/BigFloat.cpp | 1246 +--------------- CGAL_Core/src/CGAL_Core/CoreAux.cpp | 180 +-- CGAL_Core/src/CGAL_Core/CoreDefs.cpp | 123 +- CGAL_Core/src/CGAL_Core/CoreIO.cpp | 430 +----- CGAL_Core/src/CGAL_Core/Expr.cpp | 1118 +------------- CGAL_Core/src/CGAL_Core/GmpIO.cpp | 254 +--- CGAL_Core/src/CGAL_Core/Real.cpp | 239 +-- CGAL_Core/src/CGAL_Core/extLong.cpp | 154 +- Number_types/include/CGAL/CORE_BigFloat.h | 22 +- .../include/CGAL/CORE_coercion_traits.h | 6 +- .../include/CGAL/Spatial_lock_grid_3.h | 58 - 29 files changed, 4006 insertions(+), 3892 deletions(-) create mode 100644 CGAL_Core/include/CGAL/CORE/BigFloat_impl.h create mode 100644 CGAL_Core/include/CGAL/CORE/CoreAux_impl.h create mode 100644 CGAL_Core/include/CGAL/CORE/CoreDefs_impl.h create mode 100644 CGAL_Core/include/CGAL/CORE/CoreIO_impl.h create mode 100644 CGAL_Core/include/CGAL/CORE/Expr_impl.h create mode 100644 CGAL_Core/include/CGAL/CORE/Gmp_impl.h diff --git a/CGAL_Core/include/CGAL/CORE/BigFloat.h b/CGAL_Core/include/CGAL/CORE/BigFloat.h index 5ed23446997..4ffaabbb0e9 100644 --- a/CGAL_Core/include/CGAL/CORE/BigFloat.h +++ b/CGAL_Core/include/CGAL/CORE/BigFloat.h @@ -37,6 +37,13 @@ #include +#ifdef CGAL_HEADER_ONLY +#undef CGAL_EXPORT // CJTODO: TEMPORARY +#undef CGAL_CORE_EXPORT +#define CGAL_EXPORT +#define CGAL_CORE_EXPORT +#endif + namespace CORE { class Expr; @@ -83,8 +90,8 @@ public: BigFloat(const BigInt& I) : RCBigFloat(new BigFloatRep(I)) {} /// constructor for BigRat - BigFloat(const BigRat& R, const extLong& r = defRelPrec, - const extLong& a = defAbsPrec) + BigFloat(const BigRat& R, const extLong& r = get_static_defRelPrec(), + const extLong& a = get_static_defAbsPrec()) : RCBigFloat(new BigFloatRep()) { rep->approx(R, r, a); } @@ -93,8 +100,8 @@ public: // know about Expr, but BigFloat has a special role in our system! // =============================== /// constructor for Expr - explicit BigFloat(const Expr& E, const extLong& r = defRelPrec, - const extLong& a = defAbsPrec); + explicit BigFloat(const Expr& E, const extLong& r = get_static_defRelPrec(), + const extLong& a = get_static_defAbsPrec()); //Dummy explicit BigFloat(const BigFloat& E, const extLong& , @@ -161,7 +168,7 @@ public: /// operator/= BigFloat& operator/= (const BigFloat& x) { BigFloat z; - z.rep->div(*rep, *x.rep, defBFdivRelPrec); + z.rep->div(*rep, *x.rep, get_static_defBFdivRelPrec()); *this = z; return *this; } @@ -182,11 +189,11 @@ public: /// \name String Conversion Functions //@{ /// set value from const char* (base = 10) - void fromString(const char* s, const extLong& p=defBigFloatInputDigits) { + void fromString(const char* s, const extLong& p=get_static_defBigFloatInputDigits()) { rep->fromString(s, p); } /// convert to std::string (base = 10) - std::string toString(long prec=defBigFloatOutputDigits, bool sci=false) const { + std::string toString(long prec=get_static_defBigFloatOutputDigits(), bool sci=false) const { return rep->toString(prec, sci); } std::string str() const { @@ -435,7 +442,7 @@ inline BigFloat operator* (const BigFloat& x, const BigFloat& y) { /// operator/ inline BigFloat operator/ (const BigFloat& x, const BigFloat& y) { BigFloat z; - z.getRep().div(x.getRep(),y.getRep(),defBFdivRelPrec); + z.getRep().div(x.getRep(),y.getRep(),get_static_defBFdivRelPrec()); return z; } @@ -490,12 +497,12 @@ inline BigFloat power(const BigFloat& x, unsigned long p) { /// The argument x is an initial approximation. BigFloat root(const BigFloat&, unsigned long k, const extLong&, const BigFloat&); inline BigFloat root(const BigFloat& x, unsigned long k) { - return root(x, k, defBFsqrtAbsPrec, x); + return root(x, k, get_static_defBFsqrtAbsPrec(), x); } /// sqrt to defAbsPrec: inline BigFloat sqrt(const BigFloat& x) { - return x.sqrt(defBFsqrtAbsPrec); + return x.sqrt(get_static_defBFsqrtAbsPrec()); } /// convert an BigFloat Interval to a BigFloat with error bits @@ -625,4 +632,10 @@ inline BigRat::BigRat(const BigFloat& f) : RCBigRat(new BigRatRep()){ *this = f.BigRatValue(); } } //namespace CORE + +#ifdef CGAL_HEADER_ONLY +#include +#include +#endif // CGAL_HEADER_ONLY + #endif // _CORE_BIGFLOAT_H_ diff --git a/CGAL_Core/include/CGAL/CORE/BigFloatRep.h b/CGAL_Core/include/CGAL/CORE/BigFloatRep.h index 7cef72c07dc..61b39131269 100644 --- a/CGAL_Core/include/CGAL/CORE/BigFloatRep.h +++ b/CGAL_Core/include/CGAL/CORE/BigFloatRep.h @@ -40,6 +40,13 @@ #include #include +#ifdef CGAL_HEADER_ONLY +#undef CGAL_EXPORT // CJTODO: TEMPORARY +#undef CGAL_CORE_EXPORT +#define CGAL_EXPORT +#define CGAL_CORE_EXPORT +#endif + namespace CORE { // forward reference @@ -144,11 +151,11 @@ private: // conversion // toString() Joaquin Grech 31/5/2003 - std::string toString(long prec=defBigFloatOutputDigits, bool sci=false) const; + std::string toString(long prec=get_static_defBigFloatOutputDigits(), bool sci=false) const; std::string round(std::string inRep, long& L10, unsigned int width) const; - DecimalOutput toDecimal(unsigned int width=defBigFloatOutputDigits, + DecimalOutput toDecimal(unsigned int width=get_static_defBigFloatOutputDigits(), bool Scientific=false) const; - void fromString(const char *p, const extLong & prec = defBigFloatInputDigits); + void fromString(const char *p, const extLong & prec = get_static_defBigFloatInputDigits()); void dump() const; //inline long adjustE(long E, BigInt M, long e) const; diff --git a/CGAL_Core/include/CGAL/CORE/BigFloat_impl.h b/CGAL_Core/include/CGAL/CORE/BigFloat_impl.h new file mode 100644 index 00000000000..d1750ef6d4d --- /dev/null +++ b/CGAL_Core/include/CGAL/CORE/BigFloat_impl.h @@ -0,0 +1,1325 @@ +/**************************************************************************** + * Core Library Version 1.7, August 2004 + * Copyright (c) 1995-2004 Exact Computation Project + * All rights reserved. + * + * This file is part of CORE (http://cs.nyu.edu/exact/core/). + * You can redistribute it and/or modify it under the terms of the GNU + * General Public License as published by the Free Software Foundation, + * either version 3 of the License, or (at your option) any later version. + * + * Licensees holding a valid commercial license may use this file in + * accordance with the commercial license agreement provided with the + * software. + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * + * + * File: BigFloat.cpp + * Synopsis: + * BigFloat numbers with error bounds + * + * EXACTNESS PROPERTY: + * ================== + * For BigFloats that are exact (i.e., error=0), + * addition/subtraction and multiplication return the + * exact result (i.e., error=0). We also introduce the operation + * div2(), which simply divides a BigFloat by 2, + * but this again preserves exactness. Such exactness + * properties are used in our Newton iteration/Sturm Sequences. + * + * Written by + * Chee Yap + * Chen Li + * Zilin Du + * + * WWW URL: http://cs.nyu.edu/exact/ + * Email: exact@cs.nyu.edu + * + * $URL$ + * $Id$ + ***************************************************************************/ + +#ifdef CGAL_HEADER_ONLY +#define CGAL_INLINE_FUNCTION inline +#else +#define CGAL_INLINE_FUNCTION +#endif + +#include +#include +#include + +namespace CORE { + + +//////////////////////////////////////////////////////////// +// Misc Helper Functions +//////////////////////////////////////////////////////////// + +CGAL_INLINE_FUNCTION +BigInt FiveTo(unsigned long exp) { + if (exp == 0) + return BigInt(1); + else if (exp == 1) + return BigInt(5); + else { + BigInt x = FiveTo(exp / 2); + + x = x * x; + + if (exp & 1) + x *= 5; + + return x; + } +} + +//////////////////////////////////////////////////////////// +// class BigFloat +//////////////////////////////////////////////////////////// + +// STATIC BIGFLOAT CONSTANTS +// ZERO +CGAL_INLINE_FUNCTION +const BigFloat& BigFloat::getZero() { + static BigFloat Zero(0); + return Zero; +} +// ONE +CGAL_INLINE_FUNCTION +const BigFloat& BigFloat::getOne() { + static BigFloat One(1); + return One; +} + +// A special constructor for BigFloat from Expr +// -- this method is somewhat of an anomaly (we normally do not expect +// BigFloats to know about Expr). +CGAL_INLINE_FUNCTION +BigFloat::BigFloat(const Expr& E, const extLong& r, const extLong& a) + : RCBigFloat(new BigFloatRep()) { + *this = E.approx(r, a).BigFloatValue(); // lazy implementaion, any other way? +} + +//////////////////////////////////////////////////////////// +// class BigFloatRep +//////////////////////////////////////////////////////////// + +CGAL_INLINE_FUNCTION +BigFloatRep::BigFloatRep(double d) : m(0), err(0), exp(0) { + if (d != 0.0) { + int isNegative = 0; + + if (d < 0.0) { + isNegative = 1; + d = - d; + } + + int binExp; + double f = frexp(d, &binExp); + + exp = chunkFloor(binExp); + + long s = binExp - bits(exp); + + long stop = 0; + double intPart; + + // convert f into a BigInt + while (f != 0.0 && stop < DBL_MAX_CHUNK) { + f = ldexp(f, (int)CHUNK_BIT); + f = modf(f, &intPart); + m <<= CHUNK_BIT; + m += (long)intPart; + exp--; + stop++; + } +#ifdef CORE_DEBUG + assert (s >= 0); +#endif + + if (s) + m <<= s; + if (isNegative) + negate(m); + } +}//BigFloatRep constructor + +// approximation +CGAL_INLINE_FUNCTION +void BigFloatRep::trunc(const BigInt& I, const extLong& r, const extLong& a) { + if (sign(I)) { + long tr = chunkFloor((- r + bitLength(I)).asLong()); + long ta = chunkFloor(- a.asLong()); + long t; + + if (r.isInfty() || a.isTiny()) + t = ta; + else if (a.isInfty()) + t = tr; + else + t = ta < tr ? tr : ta; + + if (t > 0) { // BigInt remainder; + m = chunkShift(I, - t); + err = 1; + exp = t; + } else { // t <= 0 + m = I; + err = 0; + exp = 0; + } + } else {// I == 0 + m = 0; + err = 0; + exp = 0; + } +} + +CGAL_INLINE_FUNCTION +void BigFloatRep :: truncM(const BigFloatRep& B, const extLong& r, const extLong& a) { + if (sign(B.m)) { + long tr = chunkFloor((- 1 - r + bitLength(B.m)).asLong()); + long ta = chunkFloor(- 1 - a.asLong()) - B.exp; + long t; + + if (r.isInfty() || a.isTiny()) + t = ta; + else if (a.isInfty()) + t = tr; + else + t = ta < tr ? tr : ta; + + if (t >= chunkCeil(clLg(B.err))) { + m = chunkShift(B.m, - t); + err = 2; + exp = B.exp + t; + } else // t < chunkCeil(clLg(B.err)) + core_error(std::string("BigFloat error: truncM called with stricter") + + "precision than current error.", __FILE__, __LINE__, true); + } else {// B.m == 0 + long t = chunkFloor(- a.asLong()) - B.exp; + + if (t >= chunkCeil(clLg(B.err))) { + m = 0; + err = 1; + exp = B.exp + t; + } else // t < chunkCeil(clLg(B.err)) + core_error(std::string("BigFloat error: truncM called with stricter") + + "precision than current error.", __FILE__, __LINE__, true); + } +} + +// This is the main approximation function +// REMARK: would be useful to have a self-modifying version +// of this function (e.g., for Newton). +CGAL_INLINE_FUNCTION +void BigFloatRep::approx(const BigFloatRep& B, + const extLong& r, const extLong& a) { + if (B.err) { + if (1 + clLg(B.err) <= bitLength(B.m)) + truncM(B, r + 1, a); + else // 1 + clLg(B.err) > lg(B.m) + truncM(B, CORE_posInfty, a); + } else {// B.err == 0 + trunc(B.m, r, a - bits(B.exp)); + exp += B.exp; + } + // Call normalization globally -- IP 10/9/98 + normal(); +} + +CGAL_INLINE_FUNCTION +void BigFloatRep::div(const BigInt& N, const BigInt& D, + const extLong& r, const extLong& a) { + if (sign(D)) { + if (sign(N)) { + long tr = chunkFloor((- r + bitLength(N) - bitLength(D) - 1).asLong()); + long ta = chunkFloor(- a.asLong()); + + if (r.isInfty() || a.isTiny()) + exp = ta; + else if (a.isInfty()) + exp = tr; + else + exp = ta < tr ? tr : ta; + + BigInt remainder; + + // divide(chunkShift(N, - exp), D, m, remainder); + div_rem(m, remainder, chunkShift(N, - exp), D); + + if (exp <= 0 && sign(remainder) == 0) + err = 0; + else + err = 1; + } else {// N == 0 + m = 0; + err = 0; + exp = 0; + } + } else // D == 0 + core_error( "BigFloat error: zero divisor.", __FILE__, __LINE__, true); + + // Call normalization globally -- IP 10/9/98 + normal(); +}//div + +// error-normalization +CGAL_INLINE_FUNCTION +void BigFloatRep::normal() { + long le = flrLg(err); + + if (le >= CHUNK_BIT + 2) { // so we do not carry more than 16 = CHUNK_BIT + 2 + // bits of error + long f = chunkFloor(--le); // f is roughly equal to floor(le/CHUNK_BIT) + long bits_f = bits(f); // f chunks will have bits_f many bits +#ifdef CORE_DEBUG + assert (bits_f >= 0); +#endif + + m >>= bits_f; // reduce mantissa by bits_f many bits + err >>= bits_f; // same for err + err += 2; // why 2? + exp += f; + } + if (err == 0) // unlikely, if err += 2 above + eliminateTrailingZeroes(); +} + +// bigNormal(err) +// convert a bigInt error value (=err) into an error that fits into +// a long number. This is done by +// by increasing the exponent, and corresponding decrease +// in the bit lengths of the mantissa and error. +// +CGAL_INLINE_FUNCTION +void BigFloatRep::bigNormal(BigInt& bigErr) { + long le = bitLength(bigErr); + + if (le < CHUNK_BIT + 2) { + err = ulongValue(bigErr); + } else { + long f = chunkFloor(--le); + long bits_f = bits(f); +#ifdef CORE_DEBUG + assert(bits_f >= 0); +#endif + + m >>= bits_f; + bigErr >>= bits_f; + err = ulongValue(bigErr) + 2; // you need to add "2" because "1" comes + // from truncation error in the mantissa, and another + // "1" comes from the truncation error in the bigErr. + // (But there is danger of overflow...) + exp += f; + } + + if (err == 0) + eliminateTrailingZeroes(); +} + +// ARITHMETIC: +// Addition +CGAL_INLINE_FUNCTION +void BigFloatRep::add(const BigFloatRep& x, const BigFloatRep& y) { + long expDiff = x.exp - y.exp; + + if (expDiff > 0) {// x.exp > y.exp + if (!x.err) { + m = chunkShift(x.m, expDiff) + y.m; + err = y.err; + exp = y.exp; + } else {// x.err > 0 + m = x.m + chunkShift(y.m, - expDiff); // negative shift! + err = x.err + 5; // To account for y.err (but why 5?) + exp = x.exp; // + // normal(); + } + } else if (!expDiff) {// x.exp == y.exp + m = x.m + y.m; + err = x.err + y.err; + exp = x.exp; + // normal(); + } else {// x.exp < y.exp + if (!y.err) { + m = x.m + chunkShift(y.m, - expDiff); + err = x.err; + exp = x.exp; + } else {// y.err > 0 + m = chunkShift(x.m, expDiff) + y.m; + err = y.err + 5; + exp = y.exp; + // normal(); + } + } + // Call normalization globally -- IP 10/9/98 + normal(); +} + +// Subtraction +CGAL_INLINE_FUNCTION +void BigFloatRep::sub(const BigFloatRep& x, const BigFloatRep& y) { + long expDiff = x.exp - y.exp; + + if (expDiff > 0) {// x.exp > y.exp + if (!x.err) { + m = chunkShift(x.m, expDiff) - y.m; + err = y.err; + exp = y.exp; + } else {// x.err > 0 + m = x.m - chunkShift(y.m, - expDiff); + err = x.err + 5; + exp = x.exp; + // normal(); + } + } else if (!expDiff) { + m = x.m - y.m; + err = x.err + y.err; + exp = x.exp; + // normal(); + } else { // x.exp < y.exp + if (!y.err) { + m = x.m - chunkShift(y.m, - expDiff); + err = x.err; + exp = x.exp; + } else {// y.err > 0 + m = chunkShift(x.m, expDiff) - y.m; + err = y.err + 5; + exp = y.exp; + // normal(); + } + } + // Call normalization globally -- IP 10/9/98 + normal(); +} + +CGAL_INLINE_FUNCTION +void BigFloatRep::mul(const BigFloatRep& x, const BigFloatRep& y) { + m = x.m * y.m; + exp = x.exp + y.exp; + // compute error (new code, much faster. Zilin Du, Nov 2003) + if (x.err == 0 && y.err == 0) { + err = 0; + eliminateTrailingZeroes(); + } else { + BigInt bigErr(0); + if (y.err != 0) + bigErr += abs(x.m)*y.err; + if (x.err != 0) + bigErr += abs(y.m)*x.err; + if (x.err !=0 && y.err != 0) + bigErr += x.err*y.err; + bigNormal(bigErr); + } +} +// BigFloat div2 will half the value of x, exactly with NO error +// REMARK: should generalize this to dividing by any power of 2 +// We need this in our use of BigFloats to maintain isolation +// intervals (e.g., in Sturm sequences) --Chee/Vikram 4/2003 +// +CGAL_INLINE_FUNCTION +void BigFloatRep :: div2(const BigFloatRep& x) { + if (isEven(x.m)) { + m = (x.m >> 1); + exp = x.exp ; + } else { + m = (x.m << static_cast(CHUNK_BIT-1)); + exp = x.exp -1; + } +} + +// Converts a BigFloat interval into one BigFloat with almost same error bound +// This routine ignores the errors in inputs a and b. +// But you cannot really ignore them since, they are taken into account +// when you compute "r.sub(a,b)"... +CGAL_INLINE_FUNCTION +void BigFloatRep::centerize(const BigFloatRep& a, const BigFloatRep& b) { + if ((a.m == b.m) && (a.err == b.err) && (a.exp == b.exp)) { + m = a.m; + err = a.err; + exp = a.exp; + return; + } + + BigFloatRep r; + r.sub(a, b); + r.div2(r); + + //setup mantissa and exponent, but not error bits + // But this already sets the error bits? Chee + add(a,b); + div2(*this); + // error bits = ceil ( B^{-exp}*|a-b|/2 ) + + // bug fixed: possible overflow on converting + // Zilin & Vikram, 08/24/04 + // err = 1 + longValue(chunkShift(r.m, r.exp - exp)); + BigInt E = chunkShift(r.m, r.exp - exp); + bigNormal(E); +} + +// BigFloat Division, computing x/y: +// Unlike +,-,*, this one takes a relative precision bound R +// Note that R is only used when x and y are error-free! +// (This remark means that we may be less efficient than we could be) +// +// Assert( R>0 && R< CORE_Infty ) +// +CGAL_INLINE_FUNCTION +void BigFloatRep :: div(const BigFloatRep& x, const BigFloatRep& y, + const extLong& R) { + if (!y.isZeroIn()) { // y.m > y.err, so we are not dividing by 0 + if (!x.err && !y.err) { + if (R < 0 || R.isInfty()) //Oct 9, 2002: fixed major bug! [Zilin/Chee] + div(x.m, y.m, get_static_defBFdivRelPrec(), CORE_posInfty); + else + div(x.m, y.m, R, CORE_posInfty); + exp += x.exp - y.exp; // chen: adjust exp. + } else {// x.err > 0 or y.err > 0 + BigInt bigErr, errRemainder; + + if (x.isZeroIn()) { // x.m <= x.err + m = 0; + exp = x.exp - y.exp; + + div_rem(bigErr, errRemainder, abs(x.m) + static_cast(x.err), + abs(y.m) - static_cast(y.err)); + } else { // x.m > x.err + long lx = bitLength(x.m); + long ly = bitLength(y.m); + long r; + + if (!x.err) // x.err == 0 and y.err > 0 + r = ly + 2; + else if(!y.err) // x.err > 0 and y.err == 0 + r = lx + 2; + else // x.err > 0 and y.err > 0 + r = lx < ly ? lx + 2: ly + 2; + + long t = chunkFloor(- r + lx - ly - 1); + BigInt remainder; + + div_rem(m, remainder, chunkShift(x.m, - t), y.m); + exp = t + x.exp - y.exp; + + long delta = ((t > 0) ? 2 : 0); + + // Chen Li: 9/9/99 + // here again, it use ">>" operator with a negative + // right operand. So the result is not well defined. + // Erroneous code: + // divide(abs(remainder) + (static_cast(x.err) >> bits(t)) + // + delta + static_cast(y.err) * abs(m), + // abs(y.m) - static_cast(y.err), + // bigErr, + // errRemainder); + // New code: + BigInt errx_over_Bexp = x.err; + long bits_Bexp = bits(t); + if (bits_Bexp >= 0) { + errx_over_Bexp >>= bits_Bexp; + } else { + errx_over_Bexp <<= (-bits_Bexp); + } + + // divide(abs(remainder) + errx_over_Bexp + // + delta + static_cast(y.err) * abs(m), + // abs(y.m) - static_cast(y.err), + // bigErr, + // errRemainder); + div_rem(bigErr, errRemainder, + abs(remainder) + errx_over_Bexp + delta + static_cast(y.err) * abs(m), + abs(y.m) - static_cast(y.err)); + } + + if (sign(errRemainder)) + ++bigErr; + + bigNormal(bigErr); + } + } else {// y.m <= y.err + core_error("BigFloat error: possible zero divisor.", + __FILE__, __LINE__, true); + } + + // Call normalization globally -- IP 10/9/98 + // normal(); -- Chen: after calling bigNormal, this call is redundant. +}// BigFloatRep::div + +// squareroot for BigInt argument, without initial approximation +// sqrt(x,a) computes sqrt of x to absolute precision a. +// -- this is where Newton is applied +// -- this is called by BigFloatRep::sqrt(BigFloat, extLong) +CGAL_INLINE_FUNCTION +void BigFloatRep::sqrt(const BigInt& x, const extLong& a) { + sqrt(x, a, BigFloat(x, 0, 0)); +} // sqrt(BigInt x, extLong a) , without initial approx + +// sqrt(x,a,A) where +// x = bigInt whose sqrt is to be computed +// a = absolute precision bound +// A = initial approximation in BigFloat +// -- this is where Newton is applied +// -- it is called by BigFloatRep::sqrt(BigFloatRep, extLong, BigFloat) +CGAL_INLINE_FUNCTION +void BigFloatRep::sqrt(const BigInt& x, const extLong& a, const BigFloat& A) { + if (sign(x) == 0) { + m = 0; + err = 0; + exp = 0; + } else if (x == 1) { + m = 1; + err = 0; + exp = 0; + } else {// main case + // here is where we use the initial approximation + m = A.m(); + err = 0; + exp = A.exp(); + + BigFloatRep q, z; + extLong aa; + // need this to make sure that in case the + // initial approximation A is less than sqrt(x) + // then Newton iteration will still proceed at + // least one step. + bool firstTime = true; + for (;;) { + aa = a - bits(exp); + q.div(x, m, CORE_posInfty, aa); + q.err = 0; + q.exp -= exp; + + z.sub(*this, q); // this=current approximation, so z = this - q + /*if (sign(z.m) <= 0 || z.MSB() < - a) // justification: see Koji's + break; // thesis (p. 28) which states + // that we can exit when + // " (*this) <= q + 2**(-a)" + */ + // The preceding code is replaced by what follows: + if (z.MSB() < -a) + break; + if (sign(z.m) <= 0) { + if (firstTime) + firstTime = false; + else + break; + } + + z.add(*this, q); + // Chen Li: a bug fixed here. + // m = z.m >> 1; + // err = 0; + // exp = z.exp; + if ((z.m > 1) && isEven(z.m)) { + m = z.m >> 1; // exact division by 2 + err = 0; + exp = z.exp; + } else { // need to shift left before division by 2 + m = chunkShift(z.m, 1) >> 1; + err = 0; + exp = z.exp - 1; + }//else + }//for + }//else +} // sqrt of BigInt, with initial approx + +// MAIN ENTRY INTO SQRT FUNCTION (BIGFLOAT ARGUMENT, WITHOUT INITIAL APPROX) +CGAL_INLINE_FUNCTION +void BigFloatRep::sqrt(const BigFloatRep& x, const extLong& a) { + sqrt(x, a, BigFloat(x.m, 0, x.exp)); +} //sqrt(BigFloat, extLong a) + +// MAIN ENTRY INTO SQRT FUNCTION (BIGFLOAT ARGUMENT WITH INITIAL APPROXIMATION) +CGAL_INLINE_FUNCTION +void BigFloatRep::sqrt(const BigFloatRep& x, const extLong& a, const BigFloat& A) { + // This computes the sqrt of x to absolute precision a, starting with + // the initial approximation A + if (sign(x.m) >= 0) { // x.m >= 0 + int delta = x.exp & 1; // delta=0 if x.exp is even, otherwise delta=1 + + if (x.isZeroIn()) { // x.m <= x.err + m = 0; + if (!x.err) + err = 0; + else { // x.err > 0 + err = (long)(std::sqrt((double)x.err)); + err++; + err <<= 1; + if (delta) + err <<= HALF_CHUNK_BIT; + } + exp = x.exp >> 1; + normal(); + } else { + long aExp = A.exp() - (x.exp >> 1); + BigFloat AA( chunkShift(A.m(), delta), 0, aExp); + + if (!x.err) { // x.m > x.err = 0 (ERROR FREE CASE) + BigFloatRep z; + extLong ppp; + if (a.isInfty()) //Oct 9, 2002: fixed major bug! [Zilin/Chee] + ppp = get_static_defBFsqrtAbsPrec(); + else + ppp = a + EXTLONG_ONE; + extLong absp = ppp + bits(x.exp >> 1); + + z.sqrt(chunkShift(x.m, delta), absp, AA); // call sqrt(BigInt, a, AA) + + long p = (absp + bits(z.exp)).asLong(); + + // Next, normalize the error: + if (p <= 0) { + m = z.m; + // Chen Li: a bug fixed + // BigInt bigErr = 1 << (-p); + BigInt bigErr(1); + bigErr = bigErr << static_cast(-p); + exp = z.exp + (x.exp >> 1); + bigNormal(bigErr); + } else { // p > 0 + m = chunkShift(z.m, chunkCeil(p)); + long r = CHUNK_BIT - 1 - (p + CHUNK_BIT - 1) % CHUNK_BIT; +#ifdef CORE_DEBUG + assert(r >= 0); +#endif + + err = 1 >> r; + exp = - chunkCeil(ppp.asLong()); + normal(); + } + } else { // x.m > x.err > 0 (mantissa has error) + BigFloatRep z; + extLong absp=-flrLg(x.err)+bitLength(x.m)-(bits(delta) >> 1)+EXTLONG_FOUR; + + z.sqrt(chunkShift(x.m, delta), absp, AA); + + long qqq = - 1 + (bitLength(x.m) >> 1) - delta * HALF_CHUNK_BIT; + long qq = qqq - clLg(x.err); + long q = qq + bits(z.exp); + + if (q <= 0) { + m = z.m; + long qqqq = - qqq - bits(z.exp); + // Chen Li (09/08/99), a bug fixed here: + // BigInt bigErr = x.err << - qqqq; + // when (-qqqq) is negative, the result is not correct. + // how "<<" and ">>" process negative second operand is + // not well defined. Seems it just take it as a unsigned + // integer and extract the last few bits. + // x.err is a long number which easily overflows. + // From page 22 of Koji's paper, I think the exponent is + // wrong here. So I rewrote it as: + BigInt bigErr = x.err; + if (qqqq >= 0) { + bigErr <<= qqqq; + } else { + bigErr >>= (-qqqq); + ++bigErr; // we need to keep its ceiling. + } + + exp = z.exp + (x.exp >> 1); + bigNormal(bigErr); + } else { // q > 0 + m = chunkShift(z.m, chunkCeil(q)); + long r = CHUNK_BIT - 1 - (q + CHUNK_BIT - 1) % CHUNK_BIT; +#ifdef CORE_DEBUG + assert(r >= 0); +#endif + + err = 1 >> r; + exp = (x.exp >> 1) - chunkCeil(qq); + normal(); + } + } // end of case with error in mantissa + }//else + } else + core_error("BigFloat error: squareroot called with negative operand.", + __FILE__, __LINE__, true); +} //sqrt with initial approximation + +// compareMExp(x) +// returns 1 if *this > x +// 0 if *this = x, +// -1 if *this < x, +// +// Main comparison method for BigFloat +// This is called by BigFloat::compare() +// BE CAREFUL: The error bits are ignored! +// Need another version if we want to take care of error bits + +CGAL_INLINE_FUNCTION +int BigFloatRep :: compareMExp(const BigFloatRep& x) const { + int st = sign(m); + int sx = sign(x.m); + + if (st > sx) + return 1; + else if (st == 0 && sx == 0) + return 0; + else if (st < sx) + return - 1; + else { // need to compare m && exp + long expDiff = exp - x.exp; + + if (expDiff > 0) // exp > x.exp + return cmp(chunkShift(m, expDiff), x.m); + else if (!expDiff) + return cmp(m, x.m); + else // exp < x.exp + return cmp(m, chunkShift(x.m, - expDiff)); + } +} + +// 3/6/2000: +// This is a private function used by BigFloatRep::operator<< +// to get the exact value +// of floor(log10(M * 2^ e)) where E is an initial guess. +// We will return the correct E which satisfies +// 10^E <= M * 2^e < 10^{E+1} +// But we convert this into +// mm <= M < 10.mm + +CGAL_INLINE_FUNCTION +long BigFloatRep :: adjustE( long E, BigInt M, long ee) const { + if (M<0) + M=-M; + BigInt mm(1); + if (ee > 0) + M = (M<(ee)); + else + mm = (mm << static_cast(-ee)); + if (E > 0) + mm *= (FiveTo(E)<< static_cast(E)); + else + M *= (FiveTo(-E) << static_cast(-E)); + + if (M < mm) { + do { + E--; + M *= 10; + } while (M < mm); + } else if (M >= 10*mm) { + mm *= 10; + do { + E++; + mm *= 10; + } while (M >= mm); + } + return E; +} + +CGAL_INLINE_FUNCTION +BigFloatRep::DecimalOutput +BigFloatRep::toDecimal(unsigned int width, bool Scientific) const { + BigFloatRep::DecimalOutput decOut; // to be returned + if (err > 0) { + decOut.isExact = false; + } else { // err == 0 + decOut.isExact = true; + } + + if (err > 0 && err >= abs(m)) { + // if err is larger than mantissa, sign and significant values + // can not be determined. + core_error("BigFloat error: Error is too big!", + __FILE__, __LINE__, false); + decOut.rep = "0.0e0"; // error is too big + decOut.isScientific = false; + decOut.noSignificant = 0; + decOut.errorCode = 1; // sign of this number is unknown + return decOut; + } + + decOut.sign = sign(m); + decOut.errorCode = 0; + + BigInt M(m); // temporary mantissa + long lm = bitLength(M); // binary length of mantissa + long e2 = bits(exp); // binary shift length represented by exponent + long le = clLg(err); // binary length of err + if (le == -1) + le = 0; + + long L10 = 0; + if (M != 0) { + L10 = (long)std::floor((lm + e2) / lgTenM); + L10 = adjustE(L10, m, e2); // L10: floor[log10(M 2^(e2))], M != 0 + } else { + L10 = 0; + } + // Convention: in the positional format, when the output is + // the following string of 8 characters: + // (d0, d1, d2, d3, ".", d4, d5, d6, d7) + // then the decimal point is said to be in the 4th position. + // E.g., (d0, ".", d1, d2) has the decimal point in the 1st position. + // The value of L10 says that the decimal point of output should be at + // the (L10 + 1)st position. This is + // true regardingless of whether M = 0 or not. For zero, we output + // {0.0*} so L10=0. In general, the |value| is less than 10 + // if and only if L10 is 0 and the + // decimal point is in the 1st place. Note that L10 is defined even if + // the output is an integer (in which case it does not physically appear + // but conceptually terminates the sequence of digits). + + // First, get the decimal representaion of (m * B^(exp)). + if (e2 < 0) { + M *= FiveTo(-e2); // M = x * 10^(-e2) + } else if (e2 > 0) { + M <<= e2; // M = x * 2^(e2) + } + + std::string decRep = M.get_str(); + // Determine the "significant part" of this string, i.e. the part which + // is guaranteed to be correct in the presence of error, + // except that the last digit which might be subject to +/- 1. + + if (err != 0) { // valid = number of significant digits + unsigned long valid = floorlg10(m) - (long)std::floor(std::log10(float(err))); + if (decRep.length() > valid) { + decRep.erase(valid); + } + } + + // All the digits in decM are correct, except the last one might + // subject to an error +/- 1. + + if ((decRep[0] == '+') || (decRep[0] == '-')) { + decRep.erase(0, 1); + } + + // Second, make choice between positional representation + // and scientific notation. Use scientific notation when: + // 0) if scientific notation flag is on + // 1) err * B^exp >= 1, the error contribute to the integral part. + // 2) (1 + L10) >= width, there is not have enough places to hold the + // positional representation, not including decimal point. + // 3) The distance between the first significant digit and decimal + // point is too large for the width limit. This is equivalent to + // Either ((L10 >= 0 and (L10 + 1) > width)) + // Or ((L10 < 0) and (-L10 + 1) > width). + + if (Scientific || + ((err > 0) && (le + e2) >= 0) || // if err*B^exp >= 1 + ((L10 >= 0) && (L10 + 1 >= (long)width )) || + ((L10 < 0) && (-L10 + 1 > (long)width ))) { + // use scientific notation + decRep = round(decRep, L10, width); + decOut.noSignificant = width; + decRep.insert(1, "."); + if (L10 != 0) { + decRep += 'e'; + if (L10 > 0) { + decRep += '+'; + } else { // L10 < 0 + decRep += '-'; + } + char eBuf[48]; // enought to hold long number L10 + int ne = 0; + if ((ne = sprintf(eBuf, "%ld", labs(L10))) >= 0) { + eBuf[ne] = '\0'; + } else { + //perror("BigFloat.cpp: Problem in outputing the exponent!"); + core_error("BigFloat error: Problem in outputing the exponent", + __FILE__, __LINE__, true); + } + decRep += eBuf; + decOut.isScientific = true; + } + } else { + // use conventional positional notation. + if (L10 >= 0) { // x >= 1 or x == 0 and L10 + 1 <= width + // round when necessary + if (decRep.length() > width ) { + decRep = round(decRep, L10, width ); + if (decRep.length() > width ) { + // overflow happens! use scientific notation + return toDecimal(width, true); + } + } + decOut.noSignificant = decRep.length(); + if (L10 + 1 < (long)width ) { + decRep.insert(L10 + 1, "."); + } else { // L10 + 1 == width + // do nothing + } + } else { // L10 < 0, 0 < x < 1 + // (-L10) leading zeroes, including one to the left of decimal dot + // need to be added in beginning. + decRep = std::string(-L10, '0') + decRep; + // then round when necessary + if (decRep.length() > width ) { + decRep = round(decRep, L10, width ); + // cannot overflow since there are L10 leading zeroes. + } + decOut.noSignificant = decRep.length() - (-L10); + decRep.insert(1, "."); + } + decOut.isScientific = false; + } +#ifdef CORE_DEBUG + assert(decOut.noSignificant >= 0); +#endif + + decOut.rep = decRep; + return decOut; +}//toDecimal + +CGAL_INLINE_FUNCTION +std::string BigFloatRep::round(std::string inRep, long& L10, unsigned int width) const { + // round inRep so that the length would not exceed width. + if (inRep.length() <= width) + return inRep; + + int i = width; // < length + bool carry = false; + + if ((inRep[i] >= '5') && (inRep[i] <= '9')) { + carry = true; + i--; + while ((i >= 0) && carry) { + if (carry) { + inRep[i] ++; + if (inRep[i] > '9') { + inRep[i] = '0'; + carry = true; + } else { + carry = false; + } + } + i-- ; + } + + if ((i < 0) && carry) { // overflow + inRep.insert(inRep.begin(), '1'); + L10 ++; + width ++; + } + } + + return inRep.substr(0, width); +}//round(string,width) + + +// This function fromString(str, prec) is similar to the +// constructor Real(char * str, extLong prec) +// See the file Real.cc for the differences + +CGAL_INLINE_FUNCTION +void BigFloatRep :: fromString(const char *str, const extLong & prec ) { + // NOTE: prec defaults to get_static_defBigFloatInputDigits() (see BigFloat.h) + // check that prec is not INFTY + if (prec.isInfty()) + core_error("BigFloat error: infinite precision not allowed", + __FILE__, __LINE__, true); + + const char *e = strchr(str, 'e'); + int dot = 0; + long e10 = 0; + if (e != NULL) + e10 = atol(e+1); // e10 is decimal precision of the input string + // i.e., input is A/10^{e10}. + else { + e = str + strlen(str); +#ifdef CORE_DEBUG + assert(*e == '\0'); +#endif + + } + + const char *p = str; + if (*p == '-' || *p == '+') + p++; + m = 0; + exp = 0; + + for (; p < e; p++) { + if (*p == '.') { + dot = 1; + continue; + } + m = m * 10 + (*p - '0'); + if (dot) + e10--; + } + + BigInt one = 1; + long t = (e10 < 0) ? -e10 : e10; + BigInt ten = FiveTo(t) * (one << static_cast(t)); + + // HERE IS WHERE WE USE THE SYSTEM CONSTANT + // defBigFloatInputDigits + // Note: this constant is rather similar to defInputDigits which + // is used by Real and Expr for controlling + // input accuracy. The difference is that defInputDigits can + // be CORE_INFTY, but defBigFloatInputDigits must be finite. + + if (e10 < 0) + div(m, ten, CORE_posInfty, 4 * prec); + else + m *= ten; + if (*str == '-') + m = -m; +}//BigFloatRep::fromString + +CGAL_INLINE_FUNCTION +std::istream& BigFloatRep :: operator >>(std::istream& i) { + int size = 20; + char *str = new char[size]; + char *p = str; + char c; + int d = 0, e = 0, s = 0; + // d=1 means dot is found + // e=1 means 'e' or 'E' is found + // int done = 0; + + // Chen Li: fixed a bug, the original statement is + // for (i.get(c); c == ' '; i.get(c)); + // use isspace instead of testing c == ' ', since it must also + // skip tab, catridge/return, etc. + // Change to: + // int status; + do { + c = i.get(); + } while (isspace(c)); /* loop if met end-of-file, or + char read in is white-space. */ + // Chen Li, "if (c == EOF)" is unsafe since c is of char type and + // EOF is of int tyep with a negative value -1 + if (i.eof()) { + i.clear(std::ios::eofbit | std::ios::failbit); + return i; + } + + // the current content in "c" should be the first non-whitespace char + if (c == '-' || c == '+') { + *p++ = c; + i.get(c); + } + + for (; isdigit(c) || (!d && c=='.') || + (!e && ((c=='e') || (c=='E'))) || (!s && (c=='-' || c=='+')); i.get(c)) { + if (!e && (c == '-' || c == '+')) + break; + // Chen Li: put one more rule to prohibite input like + // xxxx.xxxe+xxx.xxx: + if (e && (c == '.')) + break; + if (p - str == size) { + char *t = str; + str = new char[size*2]; + memcpy(str, t, size); + delete [] t; + p = str + size; + size *= 2; + } +#ifdef CORE_DEBUG + assert((p-str) < size); +#endif + + *p++ = c; + if (c == '.') + d = 1; + // Chen Li: fix a bug -- the sign of exponent can not happen before + // the character "e" appears! It must follow the "e' actually. + // if (e || c == '-' || c == '+') s = 1; + if (e) + s = 1; + if ((c == 'e') || (c=='E')) + e = 1; + } + + // chenli: make sure that the p is still in the range + if (p - str >= size) { + int len = p - str; + char *t = str; + str = new char[len + 1]; + memcpy(str, t, len); + delete [] t; + p = str + len; + } + +#ifdef CORE_DEBUG + assert(p - str < size); +#endif + + *p = '\0'; + i.putback(c); + fromString(str); + delete [] str; + return i; +}//operator >> + + +// BigFloatRep::toDouble() +// converts the BigFloat to a machine double +// This is a dangerous function as the method +// is silent when it does not fit into a machine double! +// ToDo: fix this by return a machine NaN, +/- Infinity, +/- 0, +// when appropriate. +// Return NaN when error is larger than mantissa +// Return +/- Infinity if BigFloat is too big +// Return +/- 0 if BigFloat is too small +#ifdef _MSC_VER +#pragma warning(disable: 4723) +#endif +CGAL_INLINE_FUNCTION +double BigFloatRep :: toDouble() const { + if (m == 0) + return (sign(m) * 0.0); + + long e2 = bits(exp); + long le = clLg(err); // if err=0, le will be -1 + if (le == -1) + le = 0; + + BigInt M = m >> static_cast(le);// remove error bits in mantissa + + // Below, we want to return NaN by computing 0.0/0.0. + // To avoid compiler warnings about divide by zero, we do this: + + double foolCompilerZero; + foolCompilerZero = 0.0; + + // COMMENT: we should directly store the + // special IEEE values NaN, +/-Infinity, +/-0 in the code!! + + if (M == 0) + return ( 0.0/foolCompilerZero ) ; // return NaN + + e2 += le; // adjust exponent corresponding to error bits + + int len = bitLength(M) - 53; // this is positive if M is too large + + if (len > 0) { + M >>= len; + e2 += len; + } + + double tt = doubleValue(M); + + int ee = e2 + bitLength(M) - 1; // real exponent. + + if (ee >= 1024) // overflow! + return ( sign(m)/foolCompilerZero ); // return a signed infinity + + if (ee <= -1075) // underflow! + // NOTE: if (-52 < ee <= 0) get denormalized number + return ( sign(m) * 0.0 ); // return signed zero. + + // Execute this loop if e2 < 0; + for (int i = 0; i > e2; i--) + tt /= 2; + + // Execute this loop if e2 > 0; + for (int j = 0; j < e2; j++) + tt *= 2; + + return tt; +}//toDouble +#ifdef _MSC_VER +#pragma warning(default: 4723) +#endif +CGAL_INLINE_FUNCTION +BigInt BigFloatRep::toBigInt() const { + long e2 = bits(exp); + long le = clLg(err); + if (le == -1) + le = 0; +#ifdef CORE_DEBUG + assert (le >= 0); +#endif + + BigInt M = m >> static_cast(le); // discard the contaminated bits. + e2 += le; // adjust the exponent + + if (e2 < 0) + return M >> static_cast(-e2); + else if (e2 > 0) + return M << static_cast(e2); + else + return M; +} + +CGAL_INLINE_FUNCTION +long BigFloatRep :: toLong() const { + // convert a BigFloat to a long integer, rounded toward -\infty. + long e2 = bits(exp); + long le = clLg(err); +#ifdef CORE_DEBUG + assert (le >= 0); +#endif + + BigInt M = m >> static_cast(le); // discard the contaminated bits. + e2 += le; // adjust the exponent + long t; + if (e2 < 0) + t = ulongValue(M >> static_cast(-e2)); + else if (e2 > 0) + t = ulongValue(M << static_cast(e2)); + else + t = ulongValue(M); + // t = M.as_long(); + // Note: as_long() will return LONG_MAX in case of overflow. + + return t; +} + +// pow(r,n) function for BigFloat +// Note: power(r,n) calls pow(r,n) +CGAL_INLINE_FUNCTION +BigFloat pow(const BigFloat& r, unsigned long n) { + if (n == 0) + return BigFloat(1); + else if (n == 1) + return r; + else { + BigFloat x = r; + while ((n % 2) == 0) { // n is even + x *= x; + n >>= 1; + } + BigFloat u = x; + while (true) { + n >>= 1; + if (n == 0) + return u; + x *= x; + if ((n % 2) == 1) // n is odd + u *= x; + } + //return u; // unreachable + } +}//pow + +// experimental +CGAL_INLINE_FUNCTION +BigFloat root(const BigFloat& x, unsigned long k, + const extLong& a, const BigFloat& A) { + if (x.sign() == 0) { + return BigFloat(0); + } else if (x == 1) { + return BigFloat(1); + } else { + BigFloat q, del, zz; + BigFloat z = A; + BigFloat bk = long(k); + for (; ;) { + zz = pow(z, k-1); + q = x.div(zz, a); + q.makeExact(); + del = z - q; + del.makeExact(); + if (del.MSB() < -a) + break; + z = ((bk-1)*z + q).div(bk, a); + // newton's iteration: z_{n+1}=((k-1)z_n+x/z_n^{k-1})/k + z.makeExact(); + } + return z; + } +}//root + +} //namespace CORE diff --git a/CGAL_Core/include/CGAL/CORE/BigInt.h b/CGAL_Core/include/CGAL/CORE/BigInt.h index aa388c2db72..23e89b6aaf1 100644 --- a/CGAL_Core/include/CGAL/CORE/BigInt.h +++ b/CGAL_Core/include/CGAL/CORE/BigInt.h @@ -39,6 +39,13 @@ #include #include +#ifdef CGAL_HEADER_ONLY +#undef CGAL_EXPORT // CJTODO: TEMPORARY +#undef CGAL_CORE_EXPORT +#define CGAL_EXPORT +#define CGAL_CORE_EXPORT +#endif + namespace CORE { @@ -559,4 +566,9 @@ inline BigInt randomize(const BigInt& a) { //@} } //namespace CORE + +#ifdef CGAL_HEADER_ONLY +//#include +#endif // CGAL_HEADER_ONLY + #endif // _CORE_BIGINT_H_ diff --git a/CGAL_Core/include/CGAL/CORE/CoreAux.h b/CGAL_Core/include/CGAL/CORE/CoreAux.h index 82aea0734c5..fcc8cd00799 100644 --- a/CGAL_Core/include/CGAL/CORE/CoreAux.h +++ b/CGAL_Core/include/CGAL/CORE/CoreAux.h @@ -39,6 +39,13 @@ #include #include "CGAL/CORE/Impl.h" +#ifdef CGAL_HEADER_ONLY +#undef CGAL_EXPORT // CJTODO: TEMPORARY +#undef CGAL_CORE_EXPORT +#define CGAL_EXPORT +#define CGAL_CORE_EXPORT +#endif + namespace CORE { #ifndef LONG_BIT // such as in Linux @@ -181,4 +188,9 @@ inline void core_debug(std::string msg){ } //namespace CORE + +#ifdef CGAL_HEADER_ONLY +#include +#endif // CGAL_HEADER_ONLY + #endif // _CORE_COREAUX_H_ diff --git a/CGAL_Core/include/CGAL/CORE/CoreAux_impl.h b/CGAL_Core/include/CGAL/CORE/CoreAux_impl.h new file mode 100644 index 00000000000..98f31a78e43 --- /dev/null +++ b/CGAL_Core/include/CGAL/CORE/CoreAux_impl.h @@ -0,0 +1,229 @@ +/**************************************************************************** + * Core Library Version 1.7, August 2004 + * Copyright (c) 1995-2004 Exact Computation Project + * All rights reserved. + * + * This file is part of CORE (http://cs.nyu.edu/exact/core/). + * You can redistribute it and/or modify it under the terms of the GNU + * General Public License as published by the Free Software Foundation, + * either version 3 of the License, or (at your option) any later version. + * + * Licensees holding a valid commercial license may use this file in + * accordance with the commercial license agreement provided with the + * software. + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * + * + * File: CoreAux.cpp + * Synopsis: + * Auxiliary routines such as ceiling of log_2, etc. + * they are not specific to any Core classes. + * + * Written by + * Chee Yap + * Chen Li + * Zilin Du + * + * WWW URL: http://cs.nyu.edu/exact/ + * Email: exact@cs.nyu.edu + * + * $URL$ + * $Id$ + ***************************************************************************/ + +#ifdef CGAL_HEADER_ONLY +#define CGAL_INLINE_FUNCTION inline +#else +#define CGAL_INLINE_FUNCTION +#endif + +#include +#include + +namespace CORE { + +//////////////////////////////////////////////////////////// +// More useful functions to implement: +// +// To convert digits into bits: +// given X, compute X * log_2(10) +// To convert bits into digits: +// given X, compute X * log_10(2) +// +//////////////////////////////////////////////////////////// + +//////////////////////////////////////////////////////////// +// flrLg(x) +// returns floor log base 2 of abs(x) +// CONVENTION lg(0) = -1 (Slight improvement, Zilin/Chee 8/5/02) +//////////////////////////////////////////////////////////// +CGAL_INLINE_FUNCTION +int flrLg(long x) { + if (x == LONG_MIN) { + // special treatment as -LONG_MIN would be not representable as "long" + return LONG_BIT - 1; + } else { + // 1 <= |x| <= LONG_MAX + if (x < 0) + x = - x; + + int lg = -1; + while (x > 0) { + lg++; + x >>= 1; + } + return lg; + } +} + +//////////////////////////////////////////////////////////// +// floor log base 2 of unsigned long x +// CONVENTION lg(0) = -1 (Slight improvement, Zilin/Chee 8/5/02) +//////////////////////////////////////////////////////////// +CGAL_INLINE_FUNCTION +int flrLg(unsigned long x) { + int lg = -1; + while (x > 0) { + lg++; + x >>= 1; + } + return lg; +} + +//////////////////////////////////////////////////////////// +// ceiling log base 2 of abs(x) +// CONVENTION lg(0) = -1 (Slight improvement, Zilin/Chee 8/5/02) +//////////////////////////////////////////////////////////// +CGAL_INLINE_FUNCTION +int clLg(long x) { + if (x == LONG_MIN) + return LONG_BIT - 1; + if (x < 0) + x = -x; // use absolute value + if (x > (LONG_MAX >> 1)) // the leading data bit is 1 + return (LONG_BIT - 1); // exclude the sign bit + if (x >= 2) + return flrLg((unsigned long)((x << 1) - 1)); + // SINCE ceilLog_2(x) = floorLog_2(2x-1) for x>=2 + if (x == 1) + return 0; + return -1; // x must be 0 here +} + +//////////////////////////////////////////////////////////// +// ceiling log base 2 of unsigned long x +// CONVENTION lg(0) = -1 +//////////////////////////////////////////////////////////// +CGAL_INLINE_FUNCTION +int clLg(unsigned long x) { + if (x > (ULONG_MAX >> 1)) // the leading bit is 1 + return LONG_BIT; + if (x >= 2) + return flrLg((x << 1) - 1); + // SINCE ceilLog_2(x) = floorLog_2(2x-1) for x>=2. + if (x == 1) + return 0; + return -1; // x must be equal to 0 +} + +/// gcd for machine type long +/** This is needed when we construct Polynomials with int or long coefficients */ +CGAL_INLINE_FUNCTION +long gcd(long m, long n) { + if (m == 0) + return core_abs(n); + if (n == 0) + return core_abs(m); + long p = core_abs(m); + long q = core_abs(n); + if (p0) { + long r = p % q; + p = q; + q = r; + } + return p; +} + +// return a gmp_randstate_t structure +CGAL_INLINE_FUNCTION +gmp_randstate_t* getRandstate() { + static gmp_randstate_t rstate; + static bool initialized = false; + if (!initialized) { + gmp_randinit(rstate, GMP_RAND_ALG_DEFAULT, 32L); + initialized = true; + } + return &rstate; +} + +// char* core_itoa(int n, char* buffer) +// returns a pointer to the null-terminated string in buffer +// NOTES: +// 0. Buffer size should be 17 bytes (resp., 33 bytes, 65 bytes) on 16-bit +// (resp., 32-bit, 64-bit) machines. Formula: 1+sizeof(int)*8 bytes. +// 1. itoa(...) is available on some stdlib.h, but it is +// not defined by ANSI-C and so not all compilers support it. +// 2. Our use of sprintf(...) to do the job is known to +// be inefficient, but this is hardly critical for our usage. +// 3. A more general program should take a 3rd argument (the radix of +// output number). We assume radix 10. +CGAL_INLINE_FUNCTION +char * core_itoa(int n, char* buffer) { + std::sprintf(buffer, "%d", n); + return buffer; +} + +/// implements the "integer mantissa" function +// (See CORE_PATH/progs/ieee/frexp.cpp for details) +CGAL_INLINE_FUNCTION +double IntMantissa(double d) { + int e; + return std::ldexp(std::frexp(d, &e), 53); +} + +/// implements the "integer exponent" function +// (See CORE_PATH/progs/ieee/frexp.cpp for details) +CGAL_INLINE_FUNCTION +int IntExponent(double d) { + int e; + std::frexp(d, &e); + return e-53; +} + +/// CORE_DIAGFILE is file name for core_error(..) output. +const char* CORE_DIAGFILE = "Core_Diagnostics"; // global file name + +/// core_error is the method to write Core Library warning or error messages +/** Both warnings and errors are written to a file called CORE_DIAGFILE. + * But errors are also written on std:cerr (similar to std::perror()). + * */ +// Usage: core_error(message, file_with_error, line_number, err_type) +// where err_type=0 means WARNING, error_type=0 means ERROR +CGAL_INLINE_FUNCTION +void core_error(std::string msg, std::string file, int lineno, bool err) { + std::ofstream outFile(CORE_DIAGFILE, std::ios::app); // open to append + if (!outFile) { + // perror("CORE ERROR: cannot open Core Diagnostics file"); + std::cerr << "CORE ERROR: can't open Core Diagnostics file"< +#ifdef CGAL_HEADER_ONLY +#undef CGAL_EXPORT // CJTODO: TEMPORARY +#undef CGAL_CORE_EXPORT +#define CGAL_EXPORT +#define CGAL_CORE_EXPORT +#endif + +#ifdef CGAL_HEADER_ONLY + + #define CGAL_GLOBAL_STATE_VAR(TYPE, NAME, VALUE) \ + inline TYPE & get_static_##NAME() \ + { \ + static TYPE NAME = VALUE; \ + return NAME; \ + } + +#else // CGAL_HEADER_ONLY + + #define CGAL_GLOBAL_STATE_VAR(TYPE, NAME, VALUE) \ + CGAL_EXPORT extern TYPE NAME; \ + inline TYPE& get_static_##NAME() \ + { \ + return NAME; \ + } + +#endif // CGAL_HEADER_ONLY + namespace CORE { ////////////////////////////////////////////////////////////// @@ -60,34 +87,34 @@ namespace CORE { /** The normal behavior is to abort when an invalid expression * is constructed. This flag can be used to turn off this abort. * In any case, an error message will be printed */ -CGAL_CORE_EXPORT extern bool AbortFlag; +CGAL_GLOBAL_STATE_VAR(bool, AbortFlag, true) /// Invalid Flag -- initiallly value is non-negative /** If the Abort Flag is false, then the Invalid flag will be set to * a negative value whenever an invalid expression is constructed. * It is the user's responsibility to check this flag and to make * it non-negative again. */ -CGAL_CORE_EXPORT extern int InvalidFlag; +CGAL_GLOBAL_STATE_VAR(int, InvalidFlag, 0) /// Escape Precision in bits -CGAL_CORE_EXPORT extern extLong EscapePrec; +CGAL_GLOBAL_STATE_VAR(extLong, EscapePrec, CORE_posInfty) /// current ur when EscapePrec triggered /** this flag becomes negative when default EscapePrec is applied */ -CGAL_CORE_EXPORT extern long EscapePrecFlag; +CGAL_GLOBAL_STATE_VAR(long, EscapePrecFlag, 0) /// Escape Precision Warning Flag /** this flag is true by default, and will cause a warning to be printed when EscapePrec is reached */ -CGAL_CORE_EXPORT extern bool EscapePrecWarning; +CGAL_GLOBAL_STATE_VAR(bool, EscapePrecWarning, true) // These following two values determine the precision of computing // approximations in Expr. /// default Relative Precision in bits -CGAL_CORE_EXPORT extern extLong defRelPrec; +CGAL_GLOBAL_STATE_VAR(extLong, defRelPrec, 60) /// default Absolute Precision in bits -CGAL_CORE_EXPORT extern extLong defAbsPrec; +CGAL_GLOBAL_STATE_VAR(extLong, defAbsPrec, CORE_posInfty) /// default # of decimal digits for conversion from a BF to string. /** This value cannot be CORE_INFTY. @@ -97,42 +124,41 @@ CGAL_CORE_EXPORT extern extLong defAbsPrec; "controls the printout precision of std::cout for BigFloat" Perhaps, we should merge defOutputDigits and defBigFloatOutputDigits? */ -CGAL_CORE_EXPORT extern long defBigFloatOutputDigits; +CGAL_GLOBAL_STATE_VAR(long, defBigFloatOutputDigits, 10) /// default input precision in digits for converting a string to a Real or Expr /** This value can be CORE_INFTY */ -CGAL_CORE_EXPORT extern extLong defInputDigits; +CGAL_GLOBAL_STATE_VAR(extLong, defInputDigits, CORE_posInfty) /// controls the printout precision of std::cout for Real and Expr /** This value cannot be CORE_INFTY See also defBigFloatOutputDigits. (it really should be an int, as in std::cout.setprecision(int)). */ -CGAL_CORE_EXPORT extern long defOutputDigits; +CGAL_GLOBAL_STATE_VAR(long, defOutputDigits, get_static_defBigFloatOutputDigits()) /// default input precision in digits for converting a string to a BigFloat /** This value cannot be CORE_INFTY. */ -CGAL_CORE_EXPORT extern long defBigFloatInputDigits; +CGAL_GLOBAL_STATE_VAR(long, defBigFloatInputDigits, 16) /// default BigFloat Division Relative Precision -CGAL_CORE_EXPORT extern extLong defBFdivRelPrec; - +CGAL_GLOBAL_STATE_VAR(extLong, defBFdivRelPrec, 54) /// default BigFloat Sqrt Absolute Precision -CGAL_CORE_EXPORT extern extLong defBFsqrtAbsPrec; +CGAL_GLOBAL_STATE_VAR(extLong, defBFsqrtAbsPrec, 54) ////////////////////////////////////////////////////////////// // Mode parameters: incremental, progressive, filters ////////////////////////////////////////////////////////////// /// floating point filter flag -CGAL_CORE_EXPORT extern bool fpFilterFlag; +CGAL_GLOBAL_STATE_VAR(bool, fpFilterFlag, true) /// if true, evaluation of expressions would be incremental -CGAL_CORE_EXPORT extern bool incrementalEvalFlag; +CGAL_GLOBAL_STATE_VAR(bool, incrementalEvalFlag, true) /// progressive evaluation flag -CGAL_CORE_EXPORT extern bool progressiveEvalFlag; +CGAL_GLOBAL_STATE_VAR(bool, progressiveEvalFlag, true) /// rational reduction flag -CGAL_CORE_EXPORT extern bool rationalReduceFlag; +CGAL_GLOBAL_STATE_VAR(bool, rationalReduceFlag, false) /// default initial (bit) precision for AddSub Progressive Evaluation -CGAL_CORE_EXPORT extern long defInitialProgressivePrec; +CGAL_GLOBAL_STATE_VAR(long, defInitialProgressivePrec, 64) ////////////////////////////////////////////////////////////// // methods for setting global precision parameters @@ -144,87 +170,87 @@ CGAL_CORE_EXPORT extern long defInitialProgressivePrec; /** It determines the precision to which an Expr evaluates its (exact, implicit) constant value. */ inline void setDefaultPrecision(const extLong &r, const extLong &a) { - defRelPrec = r; - defAbsPrec = a; + get_static_defRelPrec() = r; + get_static_defAbsPrec() = a; } /// set default relative precision inline extLong setDefaultRelPrecision(const extLong &r) { - extLong old = defRelPrec; - defRelPrec = r; + extLong old = get_static_defRelPrec(); + get_static_defRelPrec() = r; return old; } /// set default absolute precision inline extLong setDefaultAbsPrecision(const extLong &a) { - extLong old = defAbsPrec; - defAbsPrec = a; + extLong old = get_static_defAbsPrec(); + get_static_defAbsPrec() = a; return old; } /// set default input digits (for Expr, Real) /** it controls the absolute error */ inline extLong setDefaultInputDigits(const extLong &d) { - extLong old = defInputDigits; - defInputDigits = d; + extLong old = get_static_defInputDigits(); + get_static_defInputDigits() = d; return old; } /// set default output digits (for Expr, Real) -inline long setDefaultOutputDigits(long d = defOutputDigits, +inline long setDefaultOutputDigits(long d = get_static_defOutputDigits(), std::ostream& o = std::cout) { - long old = defOutputDigits; - defOutputDigits = d; + long old = get_static_defOutputDigits(); + get_static_defOutputDigits() = d; o.precision(d); return old; } /// set default input digits for BigFloat inline long setDefaultBFInputDigits(long d) { - long old = defBigFloatInputDigits; - defBigFloatInputDigits = d; + long old = get_static_defBigFloatInputDigits(); + get_static_defBigFloatInputDigits() = d; return old; } /// set default output digits for BigFloat inline long setDefaultBFOutputDigits(long d) { - long old = defBigFloatOutputDigits; - defBigFloatOutputDigits = d; + long old = get_static_defBigFloatOutputDigits(); + get_static_defBigFloatOutputDigits() = d; return old; } /// turn floating-point filter on/off inline bool setFpFilterFlag(bool f) { - bool oldf = fpFilterFlag; - fpFilterFlag = f; + bool oldf = get_static_fpFilterFlag(); + get_static_fpFilterFlag() = f; return oldf; } /// turn incremental evaluation flag on/off inline bool setIncrementalEvalFlag(bool f) { - bool oldf = incrementalEvalFlag; - incrementalEvalFlag = f; + bool oldf = get_static_incrementalEvalFlag(); + get_static_incrementalEvalFlag() = f; return oldf; } /// turn progressive evaluation flag on/off inline bool setProgressiveEvalFlag(bool f) { - bool oldf = progressiveEvalFlag; - progressiveEvalFlag = f; + bool oldf = get_static_progressiveEvalFlag(); + get_static_progressiveEvalFlag() = f; return oldf; } /// set initial bit precision for progressive evaluation: inline long setDefInitialProgressivePrec(long n) { - long oldn = defInitialProgressivePrec; - defInitialProgressivePrec = n; + long oldn = get_static_defInitialProgressivePrec(); + get_static_defInitialProgressivePrec() = n; return oldn; } /// turn rational reduction flag on/off inline bool setRationalReduceFlag(bool f) { - bool oldf = rationalReduceFlag; - rationalReduceFlag = f; + bool oldf = get_static_rationalReduceFlag(); + get_static_rationalReduceFlag() = f; return oldf; } @@ -235,9 +261,9 @@ inline bool setRationalReduceFlag(bool f) { e.g., overriding the default std::cout precision (most systems initializes this value to 6) to our own */ inline void CORE_init(long d) { - defAbsPrec = CORE_posInfty; - defOutputDigits = d; - std::setprecision(defOutputDigits); + get_static_defAbsPrec() = CORE_posInfty; + get_static_defOutputDigits() = d; + std::setprecision(get_static_defOutputDigits()); } /// change to scientific output format @@ -251,4 +277,9 @@ inline void setPositionalFormat(std::ostream& o = std::cout) { } } //namespace CORE + +#ifdef CGAL_HEADER_ONLY +#include +#endif // CGAL_HEADER_ONLY + #endif // _CORE_COREDEFS_H_ diff --git a/CGAL_Core/include/CGAL/CORE/CoreDefs_impl.h b/CGAL_Core/include/CGAL/CORE/CoreDefs_impl.h new file mode 100644 index 00000000000..8fdd7b1e10e --- /dev/null +++ b/CGAL_Core/include/CGAL/CORE/CoreDefs_impl.h @@ -0,0 +1,170 @@ +/**************************************************************************** + * Core Library Version 1.7, August 2004 + * Copyright (c) 1995-2004 Exact Computation Project + * All rights reserved. + * + * This file is part of CORE (http://cs.nyu.edu/exact/core/). + * You can redistribute it and/or modify it under the terms of the GNU + * General Public License as published by the Free Software Foundation, + * either version 3 of the License, or (at your option) any later version. + * + * Licensees holding a valid commercial license may use this file in + * accordance with the commercial license agreement provided with the + * software. + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * + * + * File: CoreDefs.cpp + * Synopsis: + * Useful parameters for Core Library which users may change + * + * Written by + * Chee Yap + * Chen Li + * Zilin Du + * + * WWW URL: http://cs.nyu.edu/exact/ + * Email: exact@cs.nyu.edu + * + * $URL$ + * $Id$ + ***************************************************************************/ + +#include "CGAL/CORE/CoreDefs.h" + +#ifdef CGAL_HEADER_ONLY +#undef CGAL_EXPORT // CJTODO: TEMPORARY +#undef CGAL_CORE_EXPORT +#define CGAL_EXPORT +#define CGAL_CORE_EXPORT +#endif + +namespace CORE { + +// Default Values + +/* ************************************************************ + * ERROR FLAGS + * ************************************************************ */ + +#ifndef CGAL_HEADER_ONLY + +/** I/O error flag (Default value 0, indicating no error) + * User's responsibility to check and reset value to 0. */ +// This is currently used in geom2d/points2d.cpp for I/O of points + +// Note from 2014: does not seem to be used anywhere, and it is not declared +// in CoreDefs.h so it is not accessible +// Left here for compatibilty when CGAL_HEADER_ONLY is not defined + +int IOErrorFlag = 0; + +/** + * If AbortFlag is true when invalid expression is constructed, system will abort + */ + +bool AbortFlag = true; + +/** + * InvalidFlag is set to negative whenever an invalid expression is constructed. + * The user has the responsibility to reset to non-negative value. + */ + +int InvalidFlag = 0; + +/* ************************************************************ + * PRECISION PARAMETERS + * ************************************************************ */ + +/** + * Default BigFloat Division Relative Precision + * -- this is used by BigFloat division when the arguments are error-free. + */ + +extLong defBFdivRelPrec = 54; + +/** + * Default BigFloat Sqrt Absolute Precision + * -- this is used by BigFloat sqrt when the argument is error-free. + */ + +extLong defBFsqrtAbsPrec = 54; + +/** + * Escape Precision + * -- we will not compare a number to precision higher than this + * -- if this is infinity, there there is no escape precision */ +extLong EscapePrec = CORE_posInfty; + +/** this flag becomes negative if the EscapePrec is used. */ +long EscapePrecFlag = 0; + +/// Escape Precision Warning Flag +/** this flag is true by default, and will cause a warning to be printed + when EscapePrec is reached */ +bool EscapePrecWarning = true; + +/** The Composite Precision [defAbsPrec, defRelPrec] + * determines the precision to which an Expr evaluates its + * (exact, implicit) constant value. */ + +/** absolute precision = 2^31 - 1 */ +extLong defAbsPrec = CORE_posInfty; +/** default relative precision is 60 relative bits. + * Why 60? We would really like this to be 54, so that the default + * conversion duplicates the IEEE double precision. But it turns out + * (see README file under BUGS), we need 59 to ensure this. + * Chee Yap (7/1/01) */ +extLong defRelPrec = 60; + +/** number of BigFloat digits to print out */ +long defBigFloatOutputDigits = 10; + +/** NORMALLY, we like to make this equal to defBigFloatOutputDigits + * 8/3/01, Chee: re-introduced this parameter */ +long defOutputDigits = defBigFloatOutputDigits; + +/** String Input Precision */ + +/** Set this to 16 if you want machine precision. This controls the + * absolute error in string decimal inputs to Real or Expr. + * If defInputDigits is finite, then the absolute error will be + * at most 10^{-defInputDigits}. Otherwise, the input is exactly + * represented by some BigFloat or BigRat value. */ +extLong defInputDigits = CORE_posInfty; + +/** This controls the absolute error in converting from string to BigFloat + * The absolute error will be at most 10^{-defInputDigits} */ +long defBigFloatInputDigits = 16; + +/* ************************************************************ + * EVALUATION FLAGS + * ************************************************************ */ + +/** Floating Point filter + * true = turn on floating point filter */ +bool fpFilterFlag = true; + +/** IncrementaL evaluation flag + * incremental evaluation is requested, This means, we try to use previous + * approximate values to improve an approximation */ +bool incrementalEvalFlag = true; + +/** Progressive evaluation flag + * true = turn on progressive evaluation flag */ +bool progressiveEvalFlag = true; + +/** Initial progressive evaluation precision + * Used by AddSubRep */ +long defInitialProgressivePrec = 64; + +/** RATIONAL REDUCTION FLAG + * true = turn on rational reduction */ +bool rationalReduceFlag = false; + +#endif // CGAL_HEADER_ONLY + +} //namespace CORE + diff --git a/CGAL_Core/include/CGAL/CORE/CoreIO_impl.h b/CGAL_Core/include/CGAL/CORE/CoreIO_impl.h new file mode 100644 index 00000000000..2e3e10bcd8f --- /dev/null +++ b/CGAL_Core/include/CGAL/CORE/CoreIO_impl.h @@ -0,0 +1,488 @@ +/**************************************************************************** + * Core Library Version 1.7, August 2004 + * Copyright (c) 1995-2004 Exact Computation Project + * All rights reserved. + * + * This file is part of CORE (http://cs.nyu.edu/exact/core/). + * You can redistribute it and/or modify it under the terms of the GNU + * General Public License as published by the Free Software Foundation, + * either version 3 of the License, or (at your option) any later version. + * + * Licensees holding a valid commercial license may use this file in + * accordance with the commercial license agreement provided with the + * software. + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * + * + * File: CoreIO.cpp + * + * Written by + * Zilin Du + * Chee Yap + * + * WWW URL: http://cs.nyu.edu/exact/ + * Email: exact@cs.nyu.edu + * + * $URL$ + * $Id$ + ***************************************************************************/ + +#ifndef _COREIO_IMPL_H_ +#define _COREIO_IMPL_H_ + +#ifdef CGAL_HEADER_ONLY +#define CGAL_INLINE_FUNCTION inline +#else +#define CGAL_INLINE_FUNCTION +#endif + +#ifdef CGAL_HEADER_ONLY +#define CGAL_HEADERONLY_STATIC_FUNCTION static +#else +#define CGAL_HEADERONLY_STATIC_FUNCTION +#endif + +#include +#include +#include + +namespace CORE { + +CGAL_HEADERONLY_STATIC_FUNCTION +void core_io_error_handler(const char *f, const char *m) { + std::cout << "\n error_handler"; + std::cout << "::" << f << "::" << m << "\n"; + std::cout.flush(); + std::abort(); +} + +CGAL_HEADERONLY_STATIC_FUNCTION +void core_io_memory_handler(char *t, const char *f, const char *m) { + if (t == NULL) { + std::cout << "\n memory_handler"; + std::cout << "::" << f << "::" << m; + std::cout << "memory exhausted\n"; + std::cout.flush(); + std::abort(); + } +} + +// s has size old_size and will be resized to new_size. +CGAL_HEADERONLY_STATIC_FUNCTION +void allocate (char * &s, int old_size, int new_size) { + if (old_size > new_size) + old_size = new_size; + + if (s == NULL) + old_size = 0; + + char *t = new char[new_size]; + core_io_memory_handler(t, "CoreIO", "allocate::out of memory error"); + + int i; + for (i = 0; i < old_size; i++) + t[i] = s[i]; + + delete[] s; + s = t; +} + +// appends c to s at position pos. +// sz is the size of s +CGAL_HEADERONLY_STATIC_FUNCTION +void append_char (char * &s, int & sz, int pos, char c) { + if (pos > sz) + core_io_error_handler("CoreIO", "append_char::invalid argument"); + + if (pos == sz) { + allocate(s, sz, 2*sz); + sz *= 2; + } + + s[pos] = c; +} + +// skip blanks, tabs, line breaks and comment lines +CGAL_HEADERONLY_STATIC_FUNCTION +int skip_comment_line (std::istream & in) { + int c; + + do { + c = in.get(); + while ( c == '#' ) { + do { + c = in.get(); + } while ( c != '\n' ); + c = in.get(); + } + } while (c == ' ' || c == '\t' || c == '\n'); + + if (c == EOF) + core_io_error_handler("CoreIO::read_from_file()","unexpected end of file."); + + in.putback(c); + return c; +} + +// skips '\\' followed by '\n' +CGAL_HEADERONLY_STATIC_FUNCTION +int skip_backslash_new_line (std::istream & in) { + int c = in.get(); + + while (c == '\\') { + c = in.get(); + + if (c == '\n') + c = in.get(); + else + core_io_error_handler("CoreIO::operator>>", "\\ must be immediately followed by new line."); + } + + return c; +} + +CGAL_HEADERONLY_STATIC_FUNCTION +void read_string(std::istream& in, char* &buffer, int sz) { + int c, pos=0; + skip_comment_line(in); + + while ( (c = in.get()) != EOF ) { + if ( c == ' ' || c == '\t' || c == '\n' || c == '#') + break; + else + append_char(buffer, sz, pos++, c); + } + append_char(buffer, sz, pos, '\0'); +} + +CGAL_HEADERONLY_STATIC_FUNCTION +void read_base_number(std::istream& in, BigInt& m, long length, long maxBits) { + char *buffer; + int size, offset; + int base; + bool is_negate; + + int c, pos = 0; + skip_comment_line(in); + + // read sign + c = in.get(); + if (c == '-') { + is_negate = true; + c = in.get(); + } else + is_negate = false; + + // read base and compute digits + if (c == '0') { + c = in.get(); + if (c == 'b') { + base = 2; + size = (maxBits == 0 || maxBits > length) ? length : maxBits; + offset = length - size; + } else if (c == 'x') { + base = 16; + size = (maxBits == 0) ? length : (maxBits+3) >> 2; + size = (size > length) ? length : size; + offset = (length - size) << 2; + } else { + base = 8; + size = (maxBits == 0) ? length : (maxBits+2) / 3; + size = (size > length) ? length : size; + offset = (length - size) * 3; + in.putback(c); + } + } else { + base = 10; + size = (maxBits == 0) ? length : (int)std::ceil(maxBits*std::log(2.0)/std::log(10.0)); + size = (size > length) ? length : size; + offset = length - size; + in.putback(c); + } + + buffer = new char[size+2]; + // read digits + for (int i=0; (i 0 && base != 10) { + m <<= offset; + } + + if (is_negate) + negate(m); +} + + +CGAL_HEADERONLY_STATIC_FUNCTION +void write_base_number(std::ostream& out, char* buffer, int length, int base, int charsPerLine) { + // write big number in a format that gmp's mpz_set_str() can + // automatically recognize with argument base = 0. + if (base == 2) + out << "0b"; + else if (base == 16) + out << "0x"; + else if (base == 8) + out << '0'; + + // write big number in charsPerLine. + char* start, *end, c; + for (int i=0; i= length) + out << start; + else { + end = start + charsPerLine; + c = *end; + *end = '\0'; + + out << start << "\\\n"; + *end = c; + } + } +} + +CGAL_INLINE_FUNCTION +void readFromFile(BigInt& z, std::istream& in, long maxLength) { + char *buffer; + long length; + + // check type name whether it is Integer or not. + buffer = new char[8]; + read_string(in, buffer, sizeof(buffer)); + if ( std::strcmp(buffer, "Integer") != 0) + core_io_error_handler("BigInt::read_from_file()","type name expected."); + delete[] buffer; + + // read the bit length field. + buffer = new char[100]; + read_string(in, buffer, sizeof(buffer)); + length = std::atol(buffer); + delete[] buffer; + + // read bigint + read_base_number(in, z, length, maxLength); +} + +CGAL_INLINE_FUNCTION +void writeToFile(const BigInt& z, std::ostream& out, int base, int charsPerLine) { + BigInt c = abs(z); + + // get the absoulte value string + char* buffer = new char[mpz_sizeinbase(c.get_mp(), base) + 2]; + mpz_get_str(buffer, base, c.get_mp()); + int length = std::strlen(buffer); + + // write type name of big number and length + //out << "# This is an experimental big number format.\n"; + out << "Integer " << length << "\n"; + + // if bigint is negative, then write an sign '-'. + if ( sign(z) < 0 ) + out << '-'; + + write_base_number(out, buffer, length, base, charsPerLine); + out << "\n"; + delete[] buffer; +} + +CGAL_INLINE_FUNCTION +void readFromFile(BigFloat& bf, std::istream& in, long maxLength) { + char *buffer; + long length; + long exponent; + BigInt mantissa; + + // check type name whether it is Float + buffer = new char[6]; + read_string(in, buffer, sizeof(buffer)); + if (std::strcmp(buffer, "Float") != 0) + core_io_error_handler("BigFloat::read_from_file()", "type name expected"); + delete[] buffer; + + // read base (default is 16384) + buffer = new char[8]; + read_string(in, buffer, sizeof(buffer)); + if (std::strcmp(buffer, "(16384)") != 0) + core_io_error_handler("BigFloat::read_from_file()", "base expected"); + delete[] buffer; + + // read the bit length field. + buffer = new char[100]; + read_string(in, buffer, sizeof(buffer)); + length = std::atol(buffer); + delete[] buffer; + + // read exponent + buffer = new char[100]; + read_string(in, buffer, sizeof(buffer)); + exponent = std::atol(buffer); + delete[] buffer; + + // read mantissa + read_base_number(in, mantissa, length, maxLength); + + // construct BigFloat + bf = BigFloat(mantissa, 0, exponent); +} + +CGAL_INLINE_FUNCTION +void writeToFile(const BigFloat& bf, std::ostream& out, int base, int charsPerLine) { + BigInt c(CORE::abs(bf.m())); + + // get the absoulte value string + char* buffer = new char[mpz_sizeinbase(c.get_mp(), base) + 2]; + mpz_get_str(buffer, base, c.get_mp()); + int length = std::strlen(buffer); + + + // write type name, base, length + //out << "# This is an experimental Big Float format." << std::endl; + out << "Float (16384) " << length << std::endl; + // write exponent + out << bf.exp() << std::endl; + + // write mantissa + if ( CORE::sign(bf.m()) < 0 ) + out << '-'; + + write_base_number(out, buffer, length, base, charsPerLine); + out << '\n'; + delete[] buffer; +} + +/* Underconstruction ---- +void BigFloat::read_from_file2(std::istream& in, long maxLength) { + long length = 1024; + char *buffer; + + // check type name whether it is Float + buffer = new char[7]; + BigInt::read_string(in, buffer, sizeof(buffer)); + if (strcmp(buffer, "NFloat") != 0) + core_io_error_handler("BigFloat::read_from_file2()", "type name expected"); + delete[] buffer; + + // read base (default is 16) + buffer = new char[5]; + BigInt::read_string(in, buffer, sizeof(buffer)); + if (strcmp(buffer, "(16)") != 0) + core_io_error_handler("BigFloat::read_from_file2()", "base expected"); + delete[] buffer; + + // read length field + buffer = new char[100]; + BigInt::read_string(in, buffer, sizeof(buffer)); + + // get the length field if it is not null. + if (buffer[0] != '\0') { + length = atol(buffer); + if (maxLength > 0 && length >= maxLength) + length = maxLength; + } + delete[] buffer; + + // read exponent + buffer = new char[100]; + BigInt::read_string(in, buffer, sizeof(buffer)); + long exp16 = atol(buffer); + delete[] buffer; + + // read mantissa + buffer = new char[length+2]; + //BigInt::read_base_number(in, buffer, length); + + BigInt m16(buffer); + delete[] buffer; + + // convert to base CHUNK_BIT + exp16 = exp16 - length + 1; + if ( m16.is_negative() ) + exp16 ++; + + long tmp_exp = exp16 * 4; + long q = tmp_exp / CHUNK_BIT; + long r = tmp_exp % CHUNK_BIT; + if ( r < 0 ) { + r += CHUNK_BIT; + q --; + } + + BigInt mantissa = m16 << r; + long exponent = q; + + // construct BigFloat + if (--rep->refCount == 0) + delete rep; + + rep = new BigFloatRep(mantissa, 0, exponent); + rep->refCount++; + +} + +// write normal float +// now it assumed to write in hex base, i.e. B=2^4=16 +// (note: our default base B=2^(CHUNK_BIT)=2^14=16384 +void BigFloat::write_to_file2(std::ostream& out, int base, int charsPerLine) { + // convert to base 16. + long new_base = 4; // 2^4 = 16 + long tmp_exp = (rep->exp) * CHUNK_BIT; + long q = tmp_exp / new_base; + long r = tmp_exp % new_base; + std::cout << "CORE_DEBUG: q=" << q << ", r=" << r << std::endl; + if ( r < 0 ) { + r += new_base; + q--; + } + std::cout << "CORE_DEBUG: q=" << q << ", r=" << r << std::endl; + + BigInt m16 = (rep->m) << r; + + int size = mpz_sizeinbase(m16.I, base) + 2; + std::cout << "size=" << size << std::endl; + char* buffer = new char[size]; + + int length = bigint_to_string(m16, buffer, base); + std::cout << "length=" << length << std::endl; + + long exp16 = q + length - 1; + if ( m16.is_negative() ) + exp16 --; + + // write type name, base, length + out << "# This is an experimental Big Float format." << std::endl; + out << "NFloat (16) " << length << std::endl; + + // write exponent + out << exp16 << std::endl; + + // write mantissa + if ( m16.is_negative() ) { + out << '-'; + buffer ++; + } + + BigInt::write_base_number(out, buffer, length, base, charsPerLine); + out << '\n'; + delete[] buffer; +} +*/ + +} //namespace CORE + +#endif // _COREIO_IMPL_H_ diff --git a/CGAL_Core/include/CGAL/CORE/Expr.h b/CGAL_Core/include/CGAL/CORE/Expr.h index f38795eb720..90f9ca3e784 100644 --- a/CGAL_Core/include/CGAL/CORE/Expr.h +++ b/CGAL_Core/include/CGAL/CORE/Expr.h @@ -36,11 +36,22 @@ * $Id$ ***************************************************************************/ +// We need to include BigFloat.h here because there is a circular dependency +// between Expr and BigFloat. +#include + #ifndef _CORE_EXPR_H_ #define _CORE_EXPR_H_ #include +#ifdef CGAL_HEADER_ONLY +#undef CGAL_EXPORT // CJTODO: TEMPORARY +#undef CGAL_CORE_EXPORT +#define CGAL_EXPORT +#define CGAL_CORE_EXPORT +#endif + namespace CORE { /// \class Expr Expr.h @@ -78,9 +89,9 @@ public: // (i.e., not infinite and not NaN) if (! CGAL_CORE_finite(f)) { std::cerr << " ERROR : constructed an invalid float! " << std::endl; - if (AbortFlag) + if (get_static_AbortFlag()) abort(); - InvalidFlag = -1; + get_static_InvalidFlag() = -1; } rep = new ConstDoubleRep(f); } @@ -89,9 +100,9 @@ public: // (i.e., not infinite and not NaN) if (! CGAL_CORE_finite(d)) { std::cerr << " ERROR : constructed an invalid double! " << std::endl; - if (AbortFlag) + if (get_static_AbortFlag()) abort(); - InvalidFlag = -2; + get_static_InvalidFlag() = -2; } rep = new ConstDoubleRep(d); } @@ -111,11 +122,11 @@ public: * it is generally recommended that the (String) constructor be used in * preference to the (double) constructor. */ - Expr(const char *s, const extLong& p = defInputDigits) + Expr(const char *s, const extLong& p = get_static_defInputDigits()) : RCExpr(new ConstRealRep(Real(s, p))) {} /// constructor for std::string - Expr(const std::string& s, const extLong& p = defInputDigits) + Expr(const std::string& s, const extLong& p = get_static_defInputDigits()) : RCExpr(new ConstRealRep(Real(s, p))) {} /// constructor for Real @@ -179,9 +190,9 @@ public: Expr& operator/=(const Expr& e) { if ((e.rep)->getSign() == 0) { std::cerr << " ERROR : division by zero ! " << std::endl; - if (AbortFlag) + if (get_static_AbortFlag()) abort(); - InvalidFlag = -3; + get_static_InvalidFlag() = -3; } *this = new DivRep(rep, e.rep); return *this; @@ -225,12 +236,12 @@ public: /// \name String Conversion Functions //@{ /// set value from const char* - void fromString(const char* s, const extLong& prec = defInputDigits) { + void fromString(const char* s, const extLong& prec = get_static_defInputDigits()) { *this = Expr(s, prec); } /// convert to std::string /** give decimal string representation */ - std::string toString(long prec=defOutputDigits, bool sci=false) const { + std::string toString(long prec=get_static_defOutputDigits(), bool sci=false) const { return rep->toString(prec, sci); } //@} @@ -282,8 +293,8 @@ public: /** Here is the definition of what this means: If e is the exact value and ee is the approximate value, then |e - ee| <= 2^{-a} or |e - ee| <= 2^{-r} |e|. */ - const Real & approx(const extLong& relPrec = defRelPrec, - const extLong& absPrec = defAbsPrec) const { + const Real & approx(const extLong& relPrec = get_static_defRelPrec(), + const extLong& absPrec = get_static_defAbsPrec()) const { return rep->getAppValue(relPrec, absPrec); } //@} @@ -355,7 +366,7 @@ inline std::ostream& operator<<(std::ostream& o, const Expr& e) { /// I/O Stream operator>> inline std::istream& operator>>(std::istream& i, Expr& e) { Real rVal; - i >> rVal; // precision is = defInputDigits + i >> rVal; // precision is = get_static_defInputDigits() if (i) e = rVal; // only assign when reading is successful. return i; @@ -382,9 +393,9 @@ inline Expr operator*(const Expr& e1, const Expr& e2) { inline Expr operator/(const Expr& e1, const Expr& e2) { if (e2.sign() == 0) { std::cerr << " ERROR : division by zero ! " << std::endl; - if (AbortFlag) + if (get_static_AbortFlag()) abort(); - InvalidFlag = -4; + get_static_InvalidFlag() = -4; } return Expr(new DivRep(e1.Rep(), e2.Rep())); } @@ -485,9 +496,9 @@ inline bool isDivisible(const Expr& e1, const Expr& e2) { inline Expr sqrt(const Expr& e) { if (e.sign() < 0) { std::cerr << " ERROR : sqrt of negative value ! " << std::endl; - if (AbortFlag) + if (get_static_AbortFlag()) abort(); - InvalidFlag = -5; + get_static_InvalidFlag() = -5; } return Expr(new SqrtRep(e.Rep())); } @@ -548,4 +559,9 @@ inline Expr radical(const NT& n, int m) { #include } //namespace CORE + +#ifdef CGAL_HEADER_ONLY +#include +#endif // CGAL_HEADER_ONLY + #endif // _CORE_EXPR_H_ diff --git a/CGAL_Core/include/CGAL/CORE/ExprRep.h b/CGAL_Core/include/CGAL/CORE/ExprRep.h index 19c5bbc84cb..40b228aa1ab 100644 --- a/CGAL_Core/include/CGAL/CORE/ExprRep.h +++ b/CGAL_Core/include/CGAL/CORE/ExprRep.h @@ -43,9 +43,16 @@ #include #include +#ifdef CGAL_HEADER_ONLY +#undef CGAL_EXPORT // CJTODO: TEMPORARY +#undef CGAL_CORE_EXPORT +#define CGAL_EXPORT +#define CGAL_CORE_EXPORT +#endif + namespace CORE { -#ifdef CORE_DEBUG_BOUND +#if defined(CORE_DEBUG_BOUND) && !defined(CGAL_HEADER_ONLY) // These counters are incremented each time each bound is recognized as equal // to the best one in computeBound(). extern unsigned int BFMSS_counter; @@ -189,8 +196,8 @@ public: /// \name Helper Functions //@{ /// Get the approximate value - CGAL_CORE_EXPORT const Real & getAppValue(const extLong& relPrec = defRelPrec, - const extLong& absPrec = defAbsPrec); + CGAL_CORE_EXPORT const Real & getAppValue(const extLong& relPrec = get_static_defRelPrec(), + const extLong& absPrec = get_static_defAbsPrec()); /// Get the sign. int getSign(); int getExactSign(); @@ -389,8 +396,8 @@ public: BigFloat BigFloatValue(); /// represent as a string in decimal value // toString() Joaquin Grech 31/5/2003 - std::string toString(long prec=defOutputDigits, bool sci=false) { - return (getAppValue(defRelPrec, defAbsPrec)).toString(prec,sci); + std::string toString(long prec=get_static_defOutputDigits(), bool sci=false) { + return (getAppValue(get_static_defRelPrec(), get_static_defAbsPrec())).toString(prec,sci); } //@} @@ -677,7 +684,7 @@ protected: tc() = ceilLg(ss.seq[0].getTailCoeff()); // no rational reduction - if (rationalReduceFlag) + if (get_static_rationalReduceFlag()) ratFlag() = -1; flagsComputed() = true; @@ -838,7 +845,12 @@ protected: /// \brief "functor" class used as parameter to AddSubRep<> struct Add { /// name +#ifndef CGAL_HEADER_ONLY CGAL_CORE_EXPORT static const char* name; +#endif + static const char* get_name() { + return "+"; + } /// unary operator template @@ -857,7 +869,12 @@ struct Add { /// \brief "functor" class used as parameter to AddSubRep<> struct Sub { /// name +#ifndef CGAL_HEADER_ONLY CGAL_CORE_EXPORT static const char* name; +#endif + static const char* get_name() { + return "-"; + } /// unary operator template @@ -895,7 +912,7 @@ protected: void computeApproxValue(const extLong&, const extLong&); /// return operator in string const std::string op() const { - return Operator::name; + return Operator::get_name(); } private: static Operator Op; @@ -925,7 +942,7 @@ void AddSubRep::computeExactFlags() { reduceTo(second); sign() = Op(ss); appValue() = Op(appValue()); - if (rationalReduceFlag && ratFlag() > 0) + if (get_static_rationalReduceFlag() && ratFlag() > 0) *(ratValue()) = Op(*(ratValue())); return; } else if (ss == 0) { // second operand is zero @@ -933,7 +950,7 @@ void AddSubRep::computeExactFlags() { return; } // rational node - if (rationalReduceFlag) { + if (get_static_rationalReduceFlag()) { if (first->ratFlag() > 0 && second->ratFlag() > 0) { BigRat val=Op(*(first->ratValue()), *(second->ratValue())); reduceToBigRat(val); @@ -1046,7 +1063,7 @@ void AddSubRep::computeExactFlags() { if (lowBound <= EXTLONG_ZERO) lowBound = EXTLONG_ONE; - if (!progressiveEvalFlag) { + if (!get_static_progressiveEvalFlag()) { // convert the absolute error requirement "lowBound" to // a relative error requirement "ur", s.t. // |x|*2^(-ur) <= 2^(-lowBound). @@ -1085,7 +1102,7 @@ void AddSubRep::computeExactFlags() { // larger than lowBound AND the defaultInitialProgressivePrec, // so that we do at least one iteration of the for-loop. So: // i is the variable for iteration. - extLong i = core_min(defInitialProgressivePrec, lowBound.asLong()); + extLong i = core_min(get_static_defInitialProgressivePrec(), lowBound.asLong()); extLong ua = lowBound.asLong() + EXTLONG_ONE; // NOTE: ua is allowed to be CORE_INFTY @@ -1096,7 +1113,7 @@ void AddSubRep::computeExactFlags() { lMSB() = CORE_negInfty; sign() = 0; - EscapePrecFlag = 0; // reset the Escape Flag + get_static_EscapePrecFlag() = 0; // reset the Escape Flag // Now we try to determine the real lMSB and sign, // in case it is not really zero: @@ -1151,22 +1168,22 @@ void AddSubRep::computeExactFlags() { break; // assert -- this must happen in the loop if nonzero! } //8/9/01, Chee: implement escape precision here: - if (i> EscapePrec) { - EscapePrecFlag = -i.asLong();//negative means EscapePrec is used + if (i> get_static_EscapePrec()) { + get_static_EscapePrecFlag() = -i.asLong();//negative means EscapePrec is used core_error("Escape precision triggered at", __FILE__, __LINE__, false); - if (EscapePrecWarning) + if (get_static_EscapePrecWarning()) std::cout<< "Escape Precision triggered at " - << EscapePrec << " bits" << std::endl; + << get_static_EscapePrec() << " bits" << std::endl; #ifdef CORE_DEBUG - std::cout << "EscapePrecFlags=" << EscapePrecFlag << std::endl; + std::cout << "EscapePrecFlags=" << get_static_EscapePrecFlag() << std::endl; std::cout << "ua =" << ua << ",lowBound=" << lowBound << std::endl; #endif break; }// if }// for (long i=1...) -#ifdef CORE_DEBUG_BOUND +#if defined(CORE_DEBUG_BOUND) && !defined(CGAL_HEADER_ONLY) rootBoundHitCounter++; #endif diff --git a/CGAL_Core/include/CGAL/CORE/Expr_impl.h b/CGAL_Core/include/CGAL/CORE/Expr_impl.h new file mode 100644 index 00000000000..42ef3136c61 --- /dev/null +++ b/CGAL_Core/include/CGAL/CORE/Expr_impl.h @@ -0,0 +1,1213 @@ +/**************************************************************************** + * Core Library Version 1.7, August 2004 + * Copyright (c) 1995-2004 Exact Computation Project + * All rights reserved. + * + * This file is part of CORE (http://cs.nyu.edu/exact/core/). + * You can redistribute it and/or modify it under the terms of the GNU + * General Public License as published by the Free Software Foundation, + * either version 3 of the License, or (at your option) any later version. + * + * Licensees holding a valid commercial license may use this file in + * accordance with the commercial license agreement provided with the + * software. + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * + * + * File: Expr.cpp + * + * Written by + * Koji Ouchi + * Chee Yap + * Igor Pechtchanski + * Vijay Karamcheti + * Chen Li + * Zilin Du + * Sylvain Pion + * + * WWW URL: http://cs.nyu.edu/exact/ + * Email: exact@cs.nyu.edu + * + * $URL$ + * $Id$ + ***************************************************************************/ + +#ifdef CGAL_HEADER_ONLY +#define CGAL_INLINE_FUNCTION inline +#else +#define CGAL_INLINE_FUNCTION +#endif + +#include +#include + +namespace CORE { + +#if defined(CORE_DEBUG_BOUND) && !defined(CGAL_HEADER_ONLY) +unsigned int BFMSS_counter = 0; +unsigned int BFMSS_only_counter = 0; +unsigned int Measure_counter = 0; +unsigned int Measure_only_counter = 0; +unsigned int Cauchy_counter = 0; +unsigned int Cauchy_only_counter = 0; +unsigned int LiYap_counter = 0; +unsigned int LiYap_only_counter = 0; +unsigned int rootBoundHitCounter = 0; +unsigned int computeBoundCallsCounter = 0; +#endif + +#ifndef CGAL_HEADER_ONLY +const char* Add::name = "+"; +const char* Sub::name = "-"; +#endif + +/******************************************************** + * class Expr + ********************************************************/ +CGAL_INLINE_FUNCTION +const Expr& Expr::getZero() { + static Expr Zero(0); + return Zero; +} +CGAL_INLINE_FUNCTION +const Expr& Expr::getOne() { + static Expr One(1); + return One; +} + +// computes an interval comprising a pair of doubles +// Note: +// +// This function returns are two consecutive representable binary +// IEEE double values whichs contain the real value, but when you print out +// them, you might be confused by the decimal represention due to round. +// +CGAL_INLINE_FUNCTION +void Expr::doubleInterval(double & lb, double & ub) const { + double d = doubleValue(); + if (! CGAL_CORE_finite(d)) { // if overflow, underflow or NaN + lb = ub = d; + return; + } + int sign = ((* this) -Expr(d)).sign(); + // Seems like doubleValue() always give a lower bound, + // so sign = 0 or 1 (never -1). + //std::cout << "Sign = " << sign << std::endl; + if (sign == 0) { + lb = ub = d; + return; + } + int exp; + frexp(d, & exp); // get the exponent of d + exp--; // the exp from frexp satisfies + // 2^{exp-1} <= d < 2^{exp} + // But, we want exp to satisfy + // 2^{exp} <= d < 2^{exp+1} + if (sign > 0) { + lb = d; + ub = d + ldexp(1.0, -52+exp); + return; + } else { + ub = d; + lb = d - ldexp(1.0, -52+exp); + return; + } +} + +// floor(e, sub) returns the floor(e), and puts the +// remainder into sub. +CGAL_INLINE_FUNCTION +BigInt floor(const Expr& e, Expr &sub) { + if (e==0) { + return 0; + } + BigInt f = e.approx(CORE_INFTY, 2).BigIntValue(); + sub = e-f; + // Adjustment + if (sub<0) + ++sub, --f; + if (sub>=1) + --sub, ++f; + assert(sub >=0 && sub<1); // got an assertion error? (Chee 3/24/04) + return f; +} + +// Chenli: implemented algorithm from Goldberg's article. +// 7/01: Thanks to Serge Pashkov for fixing infinite loop when n=0. +CGAL_INLINE_FUNCTION +Expr pow(const Expr& e, unsigned long n) { + if (n == 0) + return Expr(1); + else if (n == 1) + return e; + else { + Expr x = e; + while ((n % 2) == 0) { // n is even + x *= x; + n >>= 1; + } + Expr u = x; + while (true) { + n >>= 1; + if (n == 0) + return u; + x *= x; + if ((n % 2) == 1) // n is odd + u *= x; + } + //return u; // unreachable + } +}//pow + +CGAL_INLINE_FUNCTION +NodeInfo::NodeInfo() : appValue(CORE_REAL_ZERO), appComputed(false), + flagsComputed(false), knownPrecision(CORE_negInfty), +#ifdef CORE_DEBUG + relPrecision(EXTLONG_ZERO), absPrecision(CORE_negInfty), numNodes(0), +#endif + // Most of the following data members don't need to be + // initialized here. + d_e(EXTLONG_ZERO), visited(false), sign(0), + uMSB(CORE_negInfty), lMSB(CORE_negInfty), + // length(0), + measure(EXTLONG_ZERO), high(EXTLONG_ZERO), low(EXTLONG_ONE), + lc(EXTLONG_ZERO), tc(EXTLONG_ZERO), + v2p(EXTLONG_ZERO), v2m(EXTLONG_ZERO), + v5p(EXTLONG_ZERO), v5m(EXTLONG_ZERO), + u25(EXTLONG_ZERO), l25(EXTLONG_ZERO), + ratFlag(0), ratValue(NULL) { } + +/******************************************************** + * class ExprRep + ********************************************************/ +// constructor +CGAL_INLINE_FUNCTION +ExprRep::ExprRep() : refCount(1), nodeInfo(NULL), ffVal(0.0) { } + +// Computes the root bit bound of the expression. +// In effect, computeBound() returns the current value of low. +CGAL_INLINE_FUNCTION +extLong ExprRep::computeBound() { + extLong measureBd = measure(); + // extLong cauchyBd = length(); + extLong ourBd = (d_e() - EXTLONG_ONE) * high() + lc(); + // BFMSS[2,5] bound. + extLong bfmsskBd; + if (v2p().isInfty() || v2m().isInfty()) + bfmsskBd = CORE_INFTY; + else + bfmsskBd = l25() + u25() * (d_e() - EXTLONG_ONE) - v2() - ceilLg5(v5()); + + // since we might compute \infty - \infty for this bound + if (bfmsskBd.isNaN()) + bfmsskBd = CORE_INFTY; + + extLong bd = core_min(measureBd, + // core_min(cauchyBd, + core_min(bfmsskBd, ourBd)); +#ifdef CORE_SHOW_BOUNDS + std::cout << "Bounds (" << measureBd << + "," << bfmsskBd << ", " << ourBd << "), "; + std::cout << "MIN = " << bd << std::endl; + std::cout << "d_e= " << d_e() << std::endl; +#endif + +#if defined(CORE_DEBUG_BOUND) && !defined(CGAL_HEADER_ONLY) + // Some statistics about which one is/are the winner[s]. + computeBoundCallsCounter++; + int number_of_winners = 0; + std::cerr << " New contest " << std::endl; + if (bd == bfmsskBd) { + BFMSS_counter++; + number_of_winners++; + std::cerr << " BFMSS is the winner " << std::endl; + } + if (bd == measureBd) { + Measure_counter++; + number_of_winners++; + std::cerr << " measureBd is the winner " << std::endl; + } + /* if (bd == cauchyBd) { + Cauchy_counter++; + number_of_winners++; + std::cerr << " cauchyBd is the winner " << std::endl; + } + */ + if (bd == ourBd) { + LiYap_counter++; + number_of_winners++; + std::cerr << " ourBd is the winner " << std::endl; + } + + assert(number_of_winners >= 1); + + if (number_of_winners == 1) { + if (bd == bfmsskBd) { + BFMSS_only_counter++; + std::cerr << " BFMSSBd is the only winner " << std::endl; + } else if (bd == measureBd) { + Measure_only_counter++; + std::cerr << " measureBd is the only winner " << std::endl; + } + /* else if (bd == cauchyBd) { + Cauchy_only_counter++; + std::cerr << " cauchyBd is the only winner " << std::endl; + } */ + else if (bd == ourBd) { + LiYap_only_counter++; + std::cerr << " ourBd is the only winner " << std::endl; + } + } +#endif + + return bd; +}//computeBound() + +CGAL_INLINE_FUNCTION +void ExprRep::reduceToBigRat(const BigRat& rat) { + Real value(rat); + + //appValue() = value; + appComputed() = false; // since appValue is not assigned until approx() is called + flagsComputed() = true; + knownPrecision() = CORE_negInfty; + +#ifdef CORE_DEBUG + relPrecision() = EXTLONG_ZERO; + absPrecision() = CORE_negInfty; + //numNodes() = numNodes(); +#endif + + d_e() = EXTLONG_ONE; + //visited() = e->visited(); + sign() = value.sign(); + uMSB() = value.MSB(); + lMSB() = value.MSB(); + // length() = value.length(); // fixed? original = 1 + measure() = value.height(); // measure <= height for rational value + + // BFMSS[2,5] bound. + value.ULV_E(u25(), l25(), v2p(), v2m(), v5p(), v5m()); + + extLong u_e = u25() + v2p(); + extLong l_e = l25() + v2m(); + + u_e = u_e + ceilLg5(v5p()); + l_e = l_e + ceilLg5(v5m()); + + if (l_e == EXTLONG_ZERO) { // no divisions introduced + high() = u_e; + low() = EXTLONG_ONE - u_e; // - (u_e - 1) + } else { + high() = u_e - l_e + EXTLONG_ONE; + low() = 2 - high(); + } + + lc() = l_e; + tc() = u_e; + + if (ratValue() == NULL) + ratValue() = new BigRat(rat); + else + *(ratValue()) = rat; +} + +// This only copies the current information of the argument e to +// *this ExprRep. +CGAL_INLINE_FUNCTION +void ExprRep::reduceTo(const ExprRep *e) { + if (e->appComputed()) { + appValue() = e->appValue(); + appComputed() = true; + flagsComputed() = true; + knownPrecision() = e->knownPrecision(); +#ifdef CORE_DEBUG + relPrecision() = e->relPrecision(); + absPrecision() = e->absPrecision(); + numNodes() = e->numNodes(); +#endif + + } + d_e() = e->d_e(); + //visited() = e->visited(); + sign() = e->sign(); + uMSB() = e->uMSB(); + lMSB() = e->lMSB(); + // length() = e->length(); // fixed? original = 1 + measure() = e->measure(); + + // BFMSS[2,5] bound. + u25() = e->u25(); + l25() = e->l25(); + v2p() = e->v2p(); + v2m() = e->v2m(); + v5p() = e->v5p(); + v5m() = e->v5m(); + + high() = e->high(); + low() = e->low(); // fixed? original = 0 + lc() = e->lc(); + tc() = e->tc(); + + // Chee (Mar 23, 2004), Notes on ratFlag(): + // =============================================================== + // For more information on the use of this flag, see progs/pentagon. + // This is an integer valued member of the NodeInfo class. + // Its value is used to determine whether + // we can ``reduce'' an Expression to a single node containing + // a BigRat value. This reduction is done if the global variable + // get_static_rationalReduceFlag()=true. The default value is false. + // This is the intepretation of ratFlag: + // ratFlag < 0 means irrational + // ratFlag = 0 means not initialized + // ratFlag > 0 means rational + // Currently, ratFlag>0 is an upper bound on the size of the expression, + // since we recursively compute + // ratFlag(v) = ratFlag(v.lchild)+ratFlag(v.rchild) + 1. + // PROPOSAL: if ratFlag() > RAT_REDUCE_THRESHHOLD + // then we automatically do a reduction. We must determine + // an empirical value for RAT_REDUCE_THRESHOLD + + if (get_static_rationalReduceFlag()) { + ratFlag() = e->ratFlag(); + + if (e->ratFlag() > 0 && e->ratValue() != NULL) { + ratFlag() ++; + if (ratValue() == NULL) + ratValue() = new BigRat(*(e->ratValue())); + else + *(ratValue()) = *(e->ratValue()); + } else + ratFlag() = -1; + } +} + +CGAL_INLINE_FUNCTION +void ExprRep::reduceToZero() { + appValue() = CORE_REAL_ZERO; + appComputed() = true; + flagsComputed() = true; + knownPrecision() = CORE_negInfty; +#ifdef CORE_DEBUG + relPrecision() = EXTLONG_ZERO; + absPrecision() = CORE_negInfty; + // numNodes() = 0; +#endif + + d_e() = EXTLONG_ONE; + visited() = false; + sign() = 0; + uMSB() = CORE_negInfty; + lMSB() = CORE_negInfty; + // length() = 0; // fixed? original = 1 + measure() = EXTLONG_ZERO; + + // BFMSS[2,5] bound. + u25() = l25() = v2p() = v2m() = v5p() = v5m() = EXTLONG_ZERO; + + low() = EXTLONG_ONE; // fixed? original = 0 + high() = lc() = tc() = EXTLONG_ZERO; + + if (get_static_rationalReduceFlag()) { + if (ratFlag() > 0) { + ratFlag() ++; + if (ratValue() == NULL) + ratValue() = new BigRat(0); + else + *(ratValue()) = 0; + } else + ratFlag() = 1; + } +} + +//////////////////////////////////////////////////////////// +// Test whether the current approximate value satisfies +// the composite precision requirements [relPrec, absPrec]. +//////////////////////////////////////////////////////////// + +CGAL_INLINE_FUNCTION +bool ExprRep::withinKnownPrecision(const extLong& relPrec, + const extLong& absPrec) { + if (appComputed()) { // an approximate value has been evaluated. + if (appValue().isExact()) { + return true; + } else { // approximation has certain error. + // decide to which position it is required to compute correctly. + extLong required = core_max(-absPrec, appValue().lMSB()-relPrec); + // see whether the existing error is smaller than the requirement. + return (knownPrecision() <= required); + } + } else + return false; +}//withinKnownPrecision(a, r) + +// approximate the expression to certain precisions when +// necessary (either no approximate value available or +// the existing one is not accurate enough). +CGAL_INLINE_FUNCTION +void ExprRep::approx(const extLong& relPrec = get_static_defRelPrec(), + const extLong& absPrec = get_static_defAbsPrec()) { + if (!getSign()) + return; // if it is exactly zero... + + // NOTE: The Filter might give a precise enough approximation already. + if (!getExactSign()) + return; + + if (!appComputed() || (!withinKnownPrecision(relPrec, absPrec))) { + // it's necessary to re-evaluate. + // to avoid huge lMSB which would cause long time and problems. + + // if it is a rational node + if (get_static_rationalReduceFlag() && ratFlag() > 0 && ratValue() != NULL) + appValue() = Real(*(ratValue())).approx(relPrec, absPrec); //< shouldn't + // this case be done by computeApproxValue()? + else + computeApproxValue(relPrec, absPrec); + + // update flags + appComputed() = true; + knownPrecision() = appValue().clLgErr(); +#ifdef CORE_DEBUG + if (relPrecision() < relPrec) + relPrecision() = relPrec; + if (absPrecision() < absPrec) + absPrecision() = absPrec; +#endif + + } +} + +// return an approximate value to certain precision. +CGAL_INLINE_FUNCTION +const Real& ExprRep::getAppValue(const extLong& relPrec, + const extLong& absPrec) { + if (getSign()) { + approx(relPrec, absPrec); + return appValue(); + } else + return CORE_REAL_ZERO; +} + +CGAL_INLINE_FUNCTION +std::ostream& operator<<(std::ostream& o, ExprRep& rep) { + if (rep.getSign()) { + rep.approx(get_static_defRelPrec(), get_static_defAbsPrec()); + o << rep.appValue(); + } else { + o << "0"; + } + return o; +} + +// Chee, Zilin: July 17, 2002 +// Original algorithm is wrongly implemented, and can take time +// exponential in the size of the dag. +// +// METHOD: +// Inductively assume that all "visited" flags are false. +// This calls for a reimplementation of "count()" and "clearFlag()". +// Actually, we did not have to fix the count() function. +// +// (1) First recursively compute d_e for each node +// by calling the count() function. +// Important thing is count() will turn the "visited" flags +// to be true, so that there is no double counting. +// Furthermore, if d_e had already been computed, the +// arithmetic for d_e can be avoided (in this case, +// it is only the setting of "visited" flag that we +// are interested in! +// (2) At the end of count(), we have set all reachable nodes +// to "visited", and their d_e have been computed. +// (3) Now, call clearFlag() to recursively clear all reachable +// nodes. NOTE THAT PREVIOUSLY, clearFlag() was called +// first! This obvious is wrong + +CGAL_INLINE_FUNCTION +extLong ExprRep::degreeBound() { + if (d_e() == EXTLONG_ONE) // no radical nodes below + return EXTLONG_ONE; + count(); + clearFlag(); + return d_e(); +} +// class ConstRealRep +// constructor +CGAL_INLINE_FUNCTION +ConstRealRep::ConstRealRep(const Real & r) : value(r) { + if (!value.isExact()) { + // clone the BigFloat and set its error to zero. + value = value.BigFloatValue().makeExact(); + } + ffVal = filteredFp(value); +} + +// initialize nodeInfo +CGAL_INLINE_FUNCTION +void ConstRep::initNodeInfo() { + nodeInfo = new NodeInfo(); + d_e() = EXTLONG_ONE; +} +CGAL_INLINE_FUNCTION +void UnaryOpRep::initNodeInfo() { + if (child->nodeInfo == NULL) + child->initNodeInfo(); + nodeInfo = new NodeInfo(); +} +CGAL_INLINE_FUNCTION +void BinOpRep::initNodeInfo() { + if (first->nodeInfo == NULL) + first->initNodeInfo(); + if (second->nodeInfo == NULL) + second->initNodeInfo(); + nodeInfo = new NodeInfo(); +} + +#ifdef CORE_DEBUG +CGAL_INLINE_FUNCTION +unsigned long ConstRep::dagSize() { + if (!visited()) { + visited() = true; + numNodes() = 1; + } else + numNodes() = 0; + return numNodes(); +} + +CGAL_INLINE_FUNCTION +unsigned long UnaryOpRep::dagSize() { + if (!visited()) { + visited() = true; + numNodes() = child->dagSize() + 1; + } else + numNodes() = 0; + return numNodes(); +} + +CGAL_INLINE_FUNCTION +unsigned long BinOpRep::dagSize() { + if (!visited()) { + visited() = true; + numNodes() = first->dagSize() + second->dagSize() + 1; + } else + numNodes() = 0; + return numNodes(); +} + +CGAL_INLINE_FUNCTION +void ConstRep::fullClearFlag() { + if (visited()) + visited() = false; +} + +CGAL_INLINE_FUNCTION +void UnaryOpRep::fullClearFlag() { + if (visited()) { + child->fullClearFlag(); + visited() = false; + } +} + +CGAL_INLINE_FUNCTION +void BinOpRep::fullClearFlag() { + if (visited()) { + first->fullClearFlag(); + second->fullClearFlag(); + visited() = false; + } +} +#endif + +// +// clear visited flag +// +/* see Expr.h + void ConstRep::clearFlag() + { visited = false; } +*/ +CGAL_INLINE_FUNCTION +void UnaryOpRep::clearFlag() { + if (d_e() == EXTLONG_ONE) + return; // no radicals below. + if (visited()) { + visited() = false; + child->clearFlag(); + } +} +// class BinOpRep +CGAL_INLINE_FUNCTION +void BinOpRep::clearFlag() { + if (d_e() == EXTLONG_ONE) + return; // rational below + if (visited()) { + visited() = false; + first->clearFlag(); + second->clearFlag(); + } +} + +// +// count # of squareroot +// +CGAL_INLINE_FUNCTION +extLong ConstRep::count() { + if (visited()) + return EXTLONG_ONE; + visited() = true; + return d_e(); +} + +CGAL_INLINE_FUNCTION +extLong NegRep::count() { + if (d_e() == EXTLONG_ONE) + return EXTLONG_ONE; + if (visited()) + return EXTLONG_ONE; + visited() = true; + d_e() = child->count(); + return d_e(); +} + +CGAL_INLINE_FUNCTION +extLong SqrtRep::count() { + if (d_e() == EXTLONG_ONE) + return EXTLONG_ONE; + if (visited()) + return EXTLONG_ONE; + visited() = true; + d_e() = child->count() * EXTLONG_TWO; + return d_e(); +} + +CGAL_INLINE_FUNCTION +extLong BinOpRep::count() { + if (d_e() == EXTLONG_ONE) + return EXTLONG_ONE; + if (visited()) + return EXTLONG_ONE; + visited() = true; + d_e() = first->count() * second->count(); + return d_e(); +} + +// +// compute exact flags functions +// +// exact value + +CGAL_INLINE_FUNCTION +void computeExactFlags_temp(ConstRep* t, const Real &value) { + // Chen Li: the following is incorrect: + // uMSB = lMSB = value.MSB(); + // because the value could be a BigFloat which is an interval. + if (value.isExact()) { + t->uMSB() = t->lMSB() = value.MSB(); + } else { + t->uMSB() = value.uMSB(); + t->lMSB() = value.lMSB(); + core_error("Leaves in DAG is not exact!", __FILE__, __LINE__, true); + } + + t->sign() = value.sign(); + // t->length() = value.length(); + t->measure() = value.height(); // for rationals and integers, + // measure = height. + + // BFMSS[2,5] bound. + value.ULV_E(t->u25(), t->l25(), t->v2p(), t->v2m(), t->v5p(), t->v5m()); + + // The original BFMSS parameters can be set from the BFMSS[2,5] parameters. + // Here we just need them locally. + extLong u_e = t->u25() + t->v2p() + ceilLg5(t->v5p()); + extLong l_e = t->l25() + t->v2m() + ceilLg5(t->v5m()); + +#ifdef ORIGINAL_BFMSS + // To go back to the original BFMSS : + t->u25() = u_e; + t->l25() = l_e; + t->v2p() = t->v2m() = t->v5p() = t->v5m() = EXTLONG_ZERO; +#elif defined BFMSS_2_ONLY + // To go back to BFMSS[2] only : + t->u25() = t->u25() + ceilLg5(t->v5p()); + t->l25() = t->l25() + ceilLg5(t->v5m()); + t->v5p() = t->v5m() = EXTLONG_ZERO; +#endif + + if (l_e == EXTLONG_ZERO) { // no divisions introduced + t->high() = u_e; + t->low() = EXTLONG_ONE - u_e; // - (u_e - 1) + } else { + t->high() = u_e - l_e + EXTLONG_ONE; + t->low() = EXTLONG_TWO - t->high(); + } + + t->lc() = l_e; + t->tc() = u_e; + + // set BigRat value + if (get_static_rationalReduceFlag()) { + t->ratFlag() = 1; + t->ratValue() = new BigRat(value.BigRatValue()); + } + + t->flagsComputed() = true; +} + +CGAL_INLINE_FUNCTION +void ConstDoubleRep::computeExactFlags() {// can be made more efficient + computeExactFlags_temp(this, Real(ffVal.getValue())); +} + +CGAL_INLINE_FUNCTION +void ConstRealRep::computeExactFlags() { + computeExactFlags_temp(this, value); +} + +CGAL_INLINE_FUNCTION +void NegRep::computeExactFlags() { + if (!child->flagsComputed()) + child->computeExactFlags(); + + if (child->sign() == 0) { + reduceToZero(); + return; + } + + if (get_static_rationalReduceFlag()) { + if (child->ratFlag()>0 && child->ratValue() != NULL) { + BigRat val = -(*(child->ratValue())); + reduceToBigRat(val); + ratFlag() = child->ratFlag()+1; + return; + } else + ratFlag() = -1; + } + + sign() = -child->sign(); + uMSB() = child->uMSB(); + lMSB() = child->lMSB(); + + // length() = child->length(); + measure() = child->measure(); + u25() = child->u25(); + l25() = child->l25(); + v2p() = child->v2p(); + v2m() = child->v2m(); + v5p() = child->v5p(); + v5m() = child->v5m(); + high() = child->high(); + low() = child->low(); + lc() = child->lc(); + tc() = child->tc(); + flagsComputed() = true; +}//NegRep::computeExactFlags + +CGAL_INLINE_FUNCTION +void SqrtRep::computeExactFlags() { + if (!child->flagsComputed()) + child->computeExactFlags(); + + if (get_static_rationalReduceFlag()) + ratFlag() = -1; + + sign() = child->sign(); + if (sign() < 0) + core_error("squareroot is called with negative operand.", + __FILE__, __LINE__, true); + + uMSB() = child->uMSB() / EXTLONG_TWO; + lMSB() = child->lMSB() / EXTLONG_TWO; + + // length() = child->length(); + measure() = child->measure(); + + // BFMSS[2,5] bound. + if (child->v2p() + ceilLg5(child->v5p()) + child->u25() >= + child->v2m() + ceilLg5(child->v5m()) + child->l25()) { + extLong vtilda2 = child->v2p() + child->v2m(); + v2p() = vtilda2 / EXTLONG_TWO; + v2m() = child->v2m(); + extLong vmod2; + if (v2p().isInfty()) + vmod2 = CORE_INFTY; + else + vmod2 = vtilda2 - EXTLONG_TWO*v2p(); // == vtilda2 % 2 + extLong vtilda5 = child->v5p() + child->v5m(); + v5p() = vtilda5 / EXTLONG_TWO; + v5m() = child->v5m(); + extLong vmod5; + if (v5p().isInfty()) + vmod5 = CORE_INFTY; + else + vmod5 = vtilda5 - EXTLONG_TWO*v5p(); // == vtilda5 % 2 + u25() = (child->u25() + child->l25() + vmod2 + ceilLg5(vmod5) + EXTLONG_ONE) / EXTLONG_TWO; + l25() = child->l25(); + } else { + extLong vtilda2 = child->v2p() + child->v2m(); + v2p() = child->v2p(); + v2m() = vtilda2 / EXTLONG_TWO; + extLong vmod2; + if (v2m().isInfty()) + vmod2 = CORE_INFTY; + else + vmod2 = vtilda2 - EXTLONG_TWO*v2m(); // == vtilda2 % 2 + extLong vtilda5 = child->v5p() + child->v5m(); + v5p() = child->v5p(); + v5m() = vtilda5 / EXTLONG_TWO; + u25() = child->u25(); + extLong vmod5; + if (v5m().isInfty()) + vmod5 = CORE_INFTY; + else + vmod5 = vtilda5 - EXTLONG_TWO*v5m(); // == vtilda5 % 2 + l25() = (child->u25() + child->l25() + vmod2 + ceilLg5(vmod5) + EXTLONG_ONE) / EXTLONG_TWO; + } + + high() = (child->high() +EXTLONG_ONE)/EXTLONG_TWO; + low() = child->low() / EXTLONG_TWO; + lc() = child->lc(); + tc() = child->tc(); + flagsComputed() = true; +}// SqrtRep::computeExactFlags + +CGAL_INLINE_FUNCTION +void MultRep::computeExactFlags() { + if (!first->flagsComputed()) + first->computeExactFlags(); + if (!second->flagsComputed()) + second->computeExactFlags(); + + if ((!first->sign()) || (!second->sign())) { + // value must be exactly zero. + reduceToZero(); + return; + } + // rational node + if (get_static_rationalReduceFlag()) { + if (first->ratFlag() > 0 && second->ratFlag() > 0) { + BigRat val = (*(first->ratValue()))*(*(second->ratValue())); + reduceToBigRat(val); + ratFlag() = first->ratFlag() + second->ratFlag(); + return; + } else + ratFlag() = -1; + } + + // value is irrational. + uMSB() = first->uMSB() + second->uMSB() + EXTLONG_ONE; + lMSB() = first->lMSB() + second->lMSB(); + sign() = first->sign() * second->sign(); + + extLong df = first->d_e(); + extLong ds = second->d_e(); + // extLong lf = first->length(); + // extLong ls = second->length(); + + // length() = df * ls + ds * lf; + measure() = (first->measure()) * ds+(second->measure()) * df; + + // BFMSS[2,5] bound. + v2p() = first->v2p() + second->v2p(); + v2m() = first->v2m() + second->v2m(); + v5p() = first->v5p() + second->v5p(); + v5m() = first->v5m() + second->v5m(); + u25() = first->u25() + second->u25(); + l25() = first->l25() + second->l25(); + + high() = first->high() + second->high(); + low() = first->low() + second->low(); + + lc() = ds * first->lc() + df * second->lc(); + tc() = core_min(ds * first->tc() + df * second->tc(), measure()); + + flagsComputed() = true; +}// MultRep::computeExactFlags + +CGAL_INLINE_FUNCTION +void DivRep::computeExactFlags() { + if (!first->flagsComputed()) + first->computeExactFlags(); + if (!second->flagsComputed()) + second->computeExactFlags(); + + if (!second->sign()) + core_error("zero divisor.", __FILE__, __LINE__, true); + + if (!first->sign()) {// value must be exactly zero. + reduceToZero(); + return; + } + + // rational node + if (get_static_rationalReduceFlag()) { + if (first->ratFlag() > 0 && second->ratFlag() > 0) { + BigRat val = (*(first->ratValue()))/(*(second->ratValue())); + reduceToBigRat(val); + ratFlag() = first->ratFlag() + second->ratFlag(); + return; + } else + ratFlag() = -1; + } + + // value is irrational. + uMSB() = first->uMSB() - second->lMSB(); + lMSB() = first->lMSB() - second->uMSB() - EXTLONG_ONE; + sign() = first->sign() * second->sign(); + + extLong df = first->d_e(); + extLong ds = second->d_e(); + // extLong lf = first->length(); + // extLong ls = second->length(); + + // length() = df * ls + ds * lf; + measure() = (first->measure())*ds + (second->measure())*df; + + // BFMSS[2,5] bound. + v2p() = first->v2p() + second->v2m(); + v2m() = first->v2m() + second->v2p(); + v5p() = first->v5p() + second->v5m(); + v5m() = first->v5m() + second->v5p(); + u25() = first->u25() + second->l25(); + l25() = first->l25() + second->u25(); + + high() = first->high() + second->low(); + low() = first->low() + second->high(); + + lc() = ds * first->lc() + df * second->tc(); + tc() = core_min(ds * first->tc() + df * second->lc(), measure()); + + flagsComputed() = true; +} + +// +// approximation functions +// +CGAL_INLINE_FUNCTION +void ConstDoubleRep::computeApproxValue(const extLong& /*relPrec*/, + const extLong& /*absPrec*/) +// can ignore precision bounds since ffVal.getValue() returns exact value +{ + appValue() = Real(ffVal.getValue()); +} + +CGAL_INLINE_FUNCTION +void ConstRealRep::computeApproxValue(const extLong& relPrec, + const extLong& absPrec) { + appValue() = value.approx(relPrec, absPrec); +} + +CGAL_INLINE_FUNCTION +void NegRep::computeApproxValue(const extLong& relPrec, + const extLong& absPrec) { + appValue() = -child->getAppValue(relPrec, absPrec); +} + +CGAL_INLINE_FUNCTION +void SqrtRep::computeApproxValue(const extLong& relPrec, + const extLong& absPrec) { + extLong r = relPrec + relPrec + EXTLONG_EIGHT; // chenli: ??? + extLong a = absPrec + absPrec + EXTLONG_EIGHT; + extLong pr = - lMSB() + r; + extLong p = pr < a ? pr : a; + + Real val = child->getAppValue(r, a); + if (get_static_incrementalEvalFlag()) { + if (appValue() == CORE_REAL_ZERO) + appValue() = val; + appValue() = val.sqrt(p, appValue().BigFloatValue()); + } else + appValue() = val.sqrt(p); +} + +CGAL_INLINE_FUNCTION +void MultRep::computeApproxValue(const extLong& relPrec, + const extLong& absPrec) { + if (lMSB() < EXTLONG_BIG && lMSB() > EXTLONG_SMALL) { + extLong r = relPrec + EXTLONG_FOUR; + extLong afr = - first->lMSB() + EXTLONG_ONE; + extLong afa = second->uMSB() + absPrec + EXTLONG_FIVE; + extLong af = afr > afa ? afr : afa; + extLong asr = - second->lMSB() + EXTLONG_ONE; + extLong asa = first->uMSB() + absPrec + EXTLONG_FIVE; + extLong as = asr > asa ? asr : asa; + appValue() = first->getAppValue(r, af)*second->getAppValue(r, as); + } else { + std::cerr << "lMSB = " << lMSB() << std::endl; + core_error("a huge lMSB in MulRep", __FILE__, __LINE__, false); + } +} + +CGAL_INLINE_FUNCTION +void DivRep::computeApproxValue(const extLong& relPrec, + const extLong& absPrec) { + if (lMSB() < EXTLONG_BIG && lMSB() > EXTLONG_SMALL) { + extLong rr = relPrec + EXTLONG_SEVEN; // These rules come from + extLong ra = uMSB() + absPrec + EXTLONG_EIGHT; // Koji's Master Thesis, page 65 + extLong ra2 = core_max(ra, EXTLONG_TWO); + extLong r = core_min(rr, ra2); + extLong af = - first->lMSB() + r; + extLong as = - second->lMSB() + r; + + extLong pr = relPrec + EXTLONG_SIX; + extLong pa = uMSB() + absPrec + EXTLONG_SEVEN; + extLong p = core_min(pr, pa); // Seems to be an error: + // p can be negative here! + // Also, this does not conform to + // Koji's thesis which has a default + // relative precision (p.65). + + appValue() = first->getAppValue(r, af).div(second->getAppValue(r, as), p); + } else { + std::cerr << "lMSB = " << lMSB() << std::endl; + core_error("a huge lMSB in DivRep", __FILE__, __LINE__, false); + } +} + +// +// Debug Help Functions +// +CGAL_INLINE_FUNCTION +void Expr::debug(int mode, int level, int depthLimit) const { + std::cout << "-------- Expr debug() -----------" << std::endl; + std::cout << "rep = " << rep << std::endl; + if (mode == Expr::LIST_MODE) + rep->debugList(level, depthLimit); + else if (mode == Expr::TREE_MODE) + rep->debugTree(level, 0, depthLimit); + else + core_error("unknown debugging mode", __FILE__, __LINE__, false); + std::cout << "---- End Expr debug(): " << std::endl; +} + + +CGAL_INLINE_FUNCTION +const std::string ExprRep::dump(int level) const { + std::ostringstream ost; + if (level == OPERATOR_ONLY) { + ost << op(); + } else if (level == VALUE_ONLY) { + ost << appValue(); + } else if (level == OPERATOR_VALUE) { + ost << op() << "[val: " << appValue() << "]"; + } else if (level == FULL_DUMP) { + ost << op() + << "[val: " << appValue() << "; " + << "kp: " << knownPrecision() << "; " +#ifdef CORE_DEBUG + << "r: " << relPrecision() << "; " + << "a: " << absPrecision() << "; " +#endif + << "lMSB: " << lMSB() << "; " + << "uMSB: " << uMSB() << "; " + << "sign: " << sign() << "; " + // << "length: " << length() << "; " + << "measure: " << measure() << "; " + << "d_e: " << d_e() << "; " + << "u25: " << u25() << "; " + << "l25: " << l25() << "; " + << "v2p: " << v2p() << "; " + << "v2m: " << v2m() << "; " + << "v5p: " << v5p() << "; " + << "v5m: " << v5m() << "; " + << "high: " << high() << "; " + << "low: " << low() << "; " + << "lc: " << lc() << "; " + << "tc: " << tc() + << "]"; + } + return std::string(ost.str()); + // note that str() return an array not properly terminated! +} + + +CGAL_INLINE_FUNCTION +void UnaryOpRep::debugList(int level, int depthLimit) const { + if (depthLimit <= 0) + return; + if (level == Expr::SIMPLE_LEVEL) { + std::cout << "(" << dump(OPERATOR_VALUE); + child->debugList(level, depthLimit - 1); + std::cout << ")"; + } else if (level == Expr::DETAIL_LEVEL) { + std::cout << "(" << dump(FULL_DUMP); + child->debugList(level, depthLimit - 1); + std::cout << ")"; + } +} + +CGAL_INLINE_FUNCTION +void UnaryOpRep::debugTree(int level, int indent, int depthLimit) const { + if (depthLimit <= 0) + return; + for (int i = 0; idebugTree(level, indent + 2, depthLimit - 1); +} + +CGAL_INLINE_FUNCTION +void ConstRep::debugList(int level, int depthLimit) const { + if (depthLimit <= 0) + return; + if (level == Expr::SIMPLE_LEVEL) { + std::cout << "(" << dump(OPERATOR_VALUE) << ")"; + } else if (level == Expr::DETAIL_LEVEL) { + std::cout << "(" << dump(FULL_DUMP) << ")"; + } +} + +CGAL_INLINE_FUNCTION +void ConstRep::debugTree(int level, int indent, int depthLimit) const { + if (depthLimit <= 0) + return; + for (int i=0; idebugList(level, depthLimit - 1); + std::cout << ", "; + second->debugList(level, depthLimit - 1); + std::cout << ")" ; +} + +CGAL_INLINE_FUNCTION +void BinOpRep::debugTree(int level, int indent, int depthLimit) const { + if (depthLimit <= 0) + return; + for (int i=0; idebugTree(level, indent + 2, depthLimit - 1); + second->debugTree(level, indent + 2, depthLimit - 1); +} + +} //namespace CORE diff --git a/CGAL_Core/include/CGAL/CORE/Filter.h b/CGAL_Core/include/CGAL/CORE/Filter.h index 0f1f5901228..42a3af162e0 100644 --- a/CGAL_Core/include/CGAL/CORE/Filter.h +++ b/CGAL_Core/include/CGAL/CORE/Filter.h @@ -42,6 +42,13 @@ #include #include +#ifdef CGAL_HEADER_ONLY +#undef CGAL_EXPORT // CJTODO: TEMPORARY +#undef CGAL_CORE_EXPORT +#define CGAL_EXPORT +#define CGAL_CORE_EXPORT +#endif + #if !defined CGAL_CFG_NO_CPP0X_ISFINITE #define CGAL_CORE_finite(x) std::isfinite(x) #elif defined (_MSC_VER) || defined (__MINGW32__) // add support for MinGW @@ -102,7 +109,7 @@ public: } /// check whether the sign (!) of the filtered value is OK bool isOK() const { - return (fpFilterFlag && // To disable filter + return (get_static_fpFilterFlag() && // To disable filter CGAL_CORE_finite(fpVal) && // Test for infinite and NaNs (core_abs(fpVal) >= maxAbs*ind*CORE_EPS)); } diff --git a/CGAL_Core/include/CGAL/CORE/Gmp.h b/CGAL_Core/include/CGAL/CORE/Gmp.h index 164095fd5d9..4b13adfa029 100644 --- a/CGAL_Core/include/CGAL/CORE/Gmp.h +++ b/CGAL_Core/include/CGAL/CORE/Gmp.h @@ -27,6 +27,13 @@ #include #include +#ifdef CGAL_HEADER_ONLY +#undef CGAL_EXPORT // CJTODO: TEMPORARY +#undef CGAL_CORE_EXPORT +#define CGAL_EXPORT +#define CGAL_CORE_EXPORT +#endif + namespace CORE { CGAL_CORE_EXPORT std::ostream& io_write (std::ostream &, mpz_srcptr); @@ -39,4 +46,9 @@ CGAL_CORE_EXPORT std::istream& io_read (std::istream &, mpq_ptr); //std::istream& operator>> (std::istream &, mpq_ptr); } //namespace CORE + +#ifdef CGAL_HEADER_ONLY +#include +#endif // CGAL_HEADER_ONLY + #endif // _CORE_GMP_H_ diff --git a/CGAL_Core/include/CGAL/CORE/Gmp_impl.h b/CGAL_Core/include/CGAL/CORE/Gmp_impl.h new file mode 100644 index 00000000000..32cec37e7ed --- /dev/null +++ b/CGAL_Core/include/CGAL/CORE/Gmp_impl.h @@ -0,0 +1,280 @@ +/**************************************************************************** + * Core Library Version 1.7, August 2004 + * Copyright (c) 1995-2004 Exact Computation Project + * All rights reserved. + * + * file: GmpIO.cpp + * Adapted from multi-files under /cxx in GMP's source distribution + * + * Zilin Du, 2003 + * + * $URL$ + * $Id$ + ***************************************************************************/ + +/* Auxiliary functions for C++-style input of GMP types. + +Copyright 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +The GNU MP Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, +MA 02110-1301, USA. */ + +#ifdef CGAL_HEADER_ONLY +#define CGAL_INLINE_FUNCTION inline +#else +#define CGAL_INLINE_FUNCTION +#endif + +#include +#include +#include +#include +#include + +using namespace std; + +namespace CORE { + +CGAL_INLINE_FUNCTION +int +__gmp_istream_set_base (istream &i, char &c, bool &zero, bool &showbase) +{ + int base; + + zero = showbase = false; + switch (i.flags() & ios::basefield) + { + case ios::dec: + base = 10; + break; + case ios::hex: + base = 16; + break; + case ios::oct: + base = 8; + break; + default: + showbase = true; // look for initial "0" or "0x" or "0X" + if (c == '0') + { + if (! i.get(c)) + c = 0; // reset or we might loop indefinitely + + if (c == 'x' || c == 'X') + { + base = 16; + i.get(c); + } + else + { + base = 8; + zero = true; // if no other digit is read, the "0" counts + } + } + else + base = 10; + break; + } + + return base; +} + +CGAL_INLINE_FUNCTION +void +__gmp_istream_set_digits (string &s, istream &i, char &c, bool &ok, int base) +{ + switch (base) + { + case 10: + while (isdigit(c)) + { + ok = true; // at least a valid digit was read + s += c; + if (! i.get(c)) + break; + } + break; + case 8: + while (isdigit(c) && c != '8' && c != '9') + { + ok = true; // at least a valid digit was read + s += c; + if (! i.get(c)) + break; + } + break; + case 16: + while (isxdigit(c)) + { + ok = true; // at least a valid digit was read + s += c; + if (! i.get(c)) + break; + } + break; + } +} + +CGAL_INLINE_FUNCTION +istream & +//operator>> (istream &i, mpz_ptr z) +io_read (istream &i, mpz_ptr z) +{ + int base; + char c = 0; + string s; + bool ok = false, zero, showbase; + + i.get(c); // start reading + + if (i.flags() & ios::skipws) // skip initial whitespace + while (isspace(c) && i.get(c)) + ; + + if (c == '-' || c == '+') // sign + { + if (c == '-') // mpz_set_str doesn't accept '+' + s = "-"; + i.get(c); + } + + while (isspace(c) && i.get(c)) // skip whitespace + ; + + base = __gmp_istream_set_base(i, c, zero, showbase); // select the base + __gmp_istream_set_digits(s, i, c, ok, base); // read the number + + if (i.good()) // last character read was non-numeric + i.putback(c); + else if (i.eof() && (ok || zero)) // stopped just before eof + i.clear(); + + if (ok) + mpz_set_str(z, s.c_str(), base); // extract the number + else if (zero) + mpz_set_ui(z, 0); + else + i.setstate(ios::failbit); // read failed + + return i; +} + +CGAL_INLINE_FUNCTION +istream & +//operator>> (istream &i, mpq_ptr q) +io_read (istream &i, mpq_ptr q) +{ + int base; + char c = 0; + string s; + bool ok = false, zero, showbase; + + i.get(c); // start reading + + if (i.flags() & ios::skipws) // skip initial whitespace + while (isspace(c) && i.get(c)) + ; + + if (c == '-' || c == '+') // sign + { + if (c == '-') + s = "-"; + i.get(c); + } + + while (isspace(c) && i.get(c)) // skip whitespace + ; + + base = __gmp_istream_set_base(i, c, zero, showbase); // select the base + __gmp_istream_set_digits(s, i, c, ok, base); // read the numerator + + if (! ok && zero) // the only digit read was "0" + { + base = 10; + s += '0'; + ok = true; + } + + if (i.flags() & ios::skipws) + while (isspace(c) && i.get(c)) // skip whitespace + ; + + if (c == '/') // there's a denominator + { + bool zero2 = false; + int base2 = base; + + s += '/'; + ok = false; // denominator is mandatory + i.get(c); + + while (isspace(c) && i.get(c)) // skip whitespace + ; + + if (showbase) // check base of denominator + base2 = __gmp_istream_set_base(i, c, zero2, showbase); + + if (base2 == base || base2 == 10) // read the denominator + __gmp_istream_set_digits(s, i, c, ok, base); + + if (! ok && zero2) // the only digit read was "0" + { // denominator is 0, but that's your business + s += '0'; + ok = true; + } + } + + if (i.good()) // last character read was non-numeric + i.putback(c); + else if (i.eof() && ok) // stopped just before eof + i.clear(); + + if (ok) + mpq_set_str(q, s.c_str(), base); // extract the number + else + i.setstate(ios::failbit); // read failed + + return i; +} + +CGAL_INLINE_FUNCTION +ostream& +//operator<< (ostream &o, mpz_srcptr z) +io_write (ostream &o, mpz_srcptr z) +{ + char *str = new char [mpz_sizeinbase(z,10) + 2]; + str = mpz_get_str(str, 10, z); + o << str ; + delete[] str; + return o; +} + +CGAL_INLINE_FUNCTION +ostream& +//operator<< (ostream &o, mpq_srcptr q) +io_write (ostream &o, mpq_srcptr q) +{ + // size according to GMP documentation + char *str = new char [mpz_sizeinbase(mpq_numref(q), 10) + + mpz_sizeinbase (mpq_denref(q), 10) + 3]; + str = mpq_get_str(str, 10, q); + o << str ; + delete[] str; + return o; +} + +} //namespace CORE diff --git a/CGAL_Core/include/CGAL/CORE/Real.h b/CGAL_Core/include/CGAL/CORE/Real.h index 52c802fed5b..66cc4718626 100644 --- a/CGAL_Core/include/CGAL/CORE/Real.h +++ b/CGAL_Core/include/CGAL/CORE/Real.h @@ -39,6 +39,13 @@ #define _CORE_REAL_H_ #include "RealRep.h" +#ifdef CGAL_HEADER_ONLY +#undef CGAL_EXPORT // CJTODO: TEMPORARY +#undef CGAL_CORE_EXPORT +#define CGAL_EXPORT +#define CGAL_CORE_EXPORT +#endif + namespace CORE { // class Real typedef RCImpl RCReal; @@ -57,10 +64,10 @@ public: Real(const BigInt& I) : RCReal(new RealBigInt(I)) {} Real(const BigRat& R) : RCReal(new RealBigRat(R)) {} Real(const BigFloat& F) : RCReal(new RealBigFloat(F)) {} - Real(const char* s, const extLong& prec=defInputDigits) : RCReal(NULL) { + Real(const char* s, const extLong& prec=get_static_defInputDigits()) : RCReal(NULL) { constructFromString(s, prec); } - Real(const std::string& s, const extLong& prec=defInputDigits) : RCReal(NULL){ + Real(const std::string& s, const extLong& prec=get_static_defInputDigits()) : RCReal(NULL){ constructFromString(s.c_str(), prec); } @@ -134,12 +141,12 @@ public: /// \name String Conversion Functions //@{ /// set value from const char* - void fromString(const char* s, const extLong& prec = defInputDigits) { + void fromString(const char* s, const extLong& prec = get_static_defInputDigits()) { *this = Real(s, prec); } /// convert to std::string /** give decimal string representation */ - std::string toString(long prec=defOutputDigits, bool sci=false) const { + std::string toString(long prec=get_static_defOutputDigits(), bool sci=false) const { return rep->toString(prec, sci); } //@} @@ -179,7 +186,8 @@ public: /// \name Aprroximation Function //@{ /// approximation - Real approx(const extLong& r=defRelPrec, const extLong& a=defAbsPrec) const { + Real approx(const extLong& r=get_static_defRelPrec(), + const extLong& a=get_static_defAbsPrec()) const { return rep->approx(r, a); } //@} @@ -395,7 +403,7 @@ inline Real& Real::operator*=(const Real& rhs) { return *this; } inline Real& Real::operator/=(const Real& rhs) { - *this = real_div::eval(getRep(), rhs.getRep(), defRelPrec); + *this = real_div::eval(getRep(), rhs.getRep(), get_static_defRelPrec()); return *this; } @@ -413,7 +421,7 @@ inline Real operator*(const Real& x, const Real& y) { } // operator/ inline Real operator/(const Real& x, const Real& y) { - return real_div::eval(x.getRep(), y.getRep(), defRelPrec); + return real_div::eval(x.getRep(), y.getRep(), get_static_defRelPrec()); } // div w/ precision inline Real Real::div(const Real& x, const extLong& r) const { @@ -478,7 +486,7 @@ inline Real power(const Real& r, unsigned long p) { } /// square root inline Real sqrt(const Real& x) { - return x.sqrt(defAbsPrec); + return x.sqrt(get_static_defAbsPrec()); } // class Realbase_for (need defined after Real) @@ -493,4 +501,9 @@ inline Real RealLong::operator-() const { } } //namespace CORE + +#ifdef CGAL_HEADER_ONLY +#include +#endif // CGAL_HEADER_ONLY + #endif // _CORE_REAL_H_ diff --git a/CGAL_Core/include/CGAL/CORE/extLong.h b/CGAL_Core/include/CGAL/CORE/extLong.h index 3091bb740b9..5060cb51d92 100644 --- a/CGAL_Core/include/CGAL/CORE/extLong.h +++ b/CGAL_Core/include/CGAL/CORE/extLong.h @@ -41,6 +41,13 @@ #include #include +#ifdef CGAL_HEADER_ONLY +#undef CGAL_EXPORT // CJTODO: TEMPORARY +#undef CGAL_CORE_EXPORT +#define CGAL_EXPORT +#define CGAL_CORE_EXPORT +#endif + namespace CORE { #ifndef LONG_MAX @@ -294,4 +301,9 @@ inline bool extLong::isNaN() const { } } //namespace CORE + +#ifdef CGAL_HEADER_ONLY +#include +#endif // CGAL_HEADER_ONLY + #endif // _CORE_EXTLONG_H_ diff --git a/CGAL_Core/include/CGAL/CORE/poly/Poly.h b/CGAL_Core/include/CGAL/CORE/poly/Poly.h index 3bb3f5b9065..ff9afa0cf7b 100644 --- a/CGAL_Core/include/CGAL/CORE/poly/Poly.h +++ b/CGAL_Core/include/CGAL/CORE/poly/Poly.h @@ -205,7 +205,8 @@ public: /// Polynomial evaluation where the coefficients are approximated first /// Returns a BigFloat with error that contains the value BigFloat evalApprox(const BigFloat& f, - const extLong& r=defRelPrec, const extLong& a=defAbsPrec) const; + const extLong& r=get_static_defRelPrec(), + const extLong& a=get_static_defAbsPrec()) const; /// Polynomial evaluation at a BigFloat value. /// The returned BigFloat (with error) has the exact sign. /// In particular, if the value is 0, we return 0. diff --git a/CGAL_Core/src/CGAL_Core/BigFloat.cpp b/CGAL_Core/src/CGAL_Core/BigFloat.cpp index 6cd7da52618..301f7b49b13 100644 --- a/CGAL_Core/src/CGAL_Core/BigFloat.cpp +++ b/CGAL_Core/src/CGAL_Core/BigFloat.cpp @@ -41,1247 +41,9 @@ * $Id$ ***************************************************************************/ -#include +#ifndef CGAL_HEADER_ONLY + #include -#include +#include -namespace CORE { - - -//////////////////////////////////////////////////////////// -// Misc Helper Functions -//////////////////////////////////////////////////////////// - -BigInt FiveTo(unsigned long exp) { - if (exp == 0) - return BigInt(1); - else if (exp == 1) - return BigInt(5); - else { - BigInt x = FiveTo(exp / 2); - - x = x * x; - - if (exp & 1) - x *= 5; - - return x; - } -} - -//////////////////////////////////////////////////////////// -// class BigFloat -//////////////////////////////////////////////////////////// - -// STATIC BIGFLOAT CONSTANTS -// ZERO -const BigFloat& BigFloat::getZero() { - static BigFloat Zero(0); - return Zero; -} -// ONE -const BigFloat& BigFloat::getOne() { - static BigFloat One(1); - return One; -} - -// A special constructor for BigFloat from Expr -// -- this method is somewhat of an anomaly (we normally do not expect -// BigFloats to know about Expr). -BigFloat::BigFloat(const Expr& E, const extLong& r, const extLong& a) - : RCBigFloat(new BigFloatRep()) { - *this = E.approx(r, a).BigFloatValue(); // lazy implementaion, any other way? -} - -//////////////////////////////////////////////////////////// -// class BigFloatRep -//////////////////////////////////////////////////////////// - -BigFloatRep::BigFloatRep(double d) : m(0), err(0), exp(0) { - if (d != 0.0) { - int isNegative = 0; - - if (d < 0.0) { - isNegative = 1; - d = - d; - } - - int binExp; - double f = frexp(d, &binExp); - - exp = chunkFloor(binExp); - - long s = binExp - bits(exp); - - long stop = 0; - double intPart; - - // convert f into a BigInt - while (f != 0.0 && stop < DBL_MAX_CHUNK) { - f = ldexp(f, (int)CHUNK_BIT); - f = modf(f, &intPart); - m <<= CHUNK_BIT; - m += (long)intPart; - exp--; - stop++; - } -#ifdef CORE_DEBUG - assert (s >= 0); -#endif - - if (s) - m <<= s; - if (isNegative) - negate(m); - } -}//BigFloatRep constructor - -// approximation -void BigFloatRep::trunc(const BigInt& I, const extLong& r, const extLong& a) { - if (sign(I)) { - long tr = chunkFloor((- r + bitLength(I)).asLong()); - long ta = chunkFloor(- a.asLong()); - long t; - - if (r.isInfty() || a.isTiny()) - t = ta; - else if (a.isInfty()) - t = tr; - else - t = ta < tr ? tr : ta; - - if (t > 0) { // BigInt remainder; - m = chunkShift(I, - t); - err = 1; - exp = t; - } else { // t <= 0 - m = I; - err = 0; - exp = 0; - } - } else {// I == 0 - m = 0; - err = 0; - exp = 0; - } -} - -void BigFloatRep :: truncM(const BigFloatRep& B, const extLong& r, const extLong& a) { - if (sign(B.m)) { - long tr = chunkFloor((- 1 - r + bitLength(B.m)).asLong()); - long ta = chunkFloor(- 1 - a.asLong()) - B.exp; - long t; - - if (r.isInfty() || a.isTiny()) - t = ta; - else if (a.isInfty()) - t = tr; - else - t = ta < tr ? tr : ta; - - if (t >= chunkCeil(clLg(B.err))) { - m = chunkShift(B.m, - t); - err = 2; - exp = B.exp + t; - } else // t < chunkCeil(clLg(B.err)) - core_error(std::string("BigFloat error: truncM called with stricter") - + "precision than current error.", __FILE__, __LINE__, true); - } else {// B.m == 0 - long t = chunkFloor(- a.asLong()) - B.exp; - - if (t >= chunkCeil(clLg(B.err))) { - m = 0; - err = 1; - exp = B.exp + t; - } else // t < chunkCeil(clLg(B.err)) - core_error(std::string("BigFloat error: truncM called with stricter") - + "precision than current error.", __FILE__, __LINE__, true); - } -} - -// This is the main approximation function -// REMARK: would be useful to have a self-modifying version -// of this function (e.g., for Newton). -void BigFloatRep::approx(const BigFloatRep& B, - const extLong& r, const extLong& a) { - if (B.err) { - if (1 + clLg(B.err) <= bitLength(B.m)) - truncM(B, r + 1, a); - else // 1 + clLg(B.err) > lg(B.m) - truncM(B, CORE_posInfty, a); - } else {// B.err == 0 - trunc(B.m, r, a - bits(B.exp)); - exp += B.exp; - } - // Call normalization globally -- IP 10/9/98 - normal(); -} - -void BigFloatRep::div(const BigInt& N, const BigInt& D, - const extLong& r, const extLong& a) { - if (sign(D)) { - if (sign(N)) { - long tr = chunkFloor((- r + bitLength(N) - bitLength(D) - 1).asLong()); - long ta = chunkFloor(- a.asLong()); - - if (r.isInfty() || a.isTiny()) - exp = ta; - else if (a.isInfty()) - exp = tr; - else - exp = ta < tr ? tr : ta; - - BigInt remainder; - - // divide(chunkShift(N, - exp), D, m, remainder); - div_rem(m, remainder, chunkShift(N, - exp), D); - - if (exp <= 0 && sign(remainder) == 0) - err = 0; - else - err = 1; - } else {// N == 0 - m = 0; - err = 0; - exp = 0; - } - } else // D == 0 - core_error( "BigFloat error: zero divisor.", __FILE__, __LINE__, true); - - // Call normalization globally -- IP 10/9/98 - normal(); -}//div - -// error-normalization -void BigFloatRep::normal() { - long le = flrLg(err); - - if (le >= CHUNK_BIT + 2) { // so we do not carry more than 16 = CHUNK_BIT + 2 - // bits of error - long f = chunkFloor(--le); // f is roughly equal to floor(le/CHUNK_BIT) - long bits_f = bits(f); // f chunks will have bits_f many bits -#ifdef CORE_DEBUG - assert (bits_f >= 0); -#endif - - m >>= bits_f; // reduce mantissa by bits_f many bits - err >>= bits_f; // same for err - err += 2; // why 2? - exp += f; - } - if (err == 0) // unlikely, if err += 2 above - eliminateTrailingZeroes(); -} - -// bigNormal(err) -// convert a bigInt error value (=err) into an error that fits into -// a long number. This is done by -// by increasing the exponent, and corresponding decrease -// in the bit lengths of the mantissa and error. -// -void BigFloatRep::bigNormal(BigInt& bigErr) { - long le = bitLength(bigErr); - - if (le < CHUNK_BIT + 2) { - err = ulongValue(bigErr); - } else { - long f = chunkFloor(--le); - long bits_f = bits(f); -#ifdef CORE_DEBUG - assert(bits_f >= 0); -#endif - - m >>= bits_f; - bigErr >>= bits_f; - err = ulongValue(bigErr) + 2; // you need to add "2" because "1" comes - // from truncation error in the mantissa, and another - // "1" comes from the truncation error in the bigErr. - // (But there is danger of overflow...) - exp += f; - } - - if (err == 0) - eliminateTrailingZeroes(); -} - -// ARITHMETIC: -// Addition -void BigFloatRep::add(const BigFloatRep& x, const BigFloatRep& y) { - long expDiff = x.exp - y.exp; - - if (expDiff > 0) {// x.exp > y.exp - if (!x.err) { - m = chunkShift(x.m, expDiff) + y.m; - err = y.err; - exp = y.exp; - } else {// x.err > 0 - m = x.m + chunkShift(y.m, - expDiff); // negative shift! - err = x.err + 5; // To account for y.err (but why 5?) - exp = x.exp; // - // normal(); - } - } else if (!expDiff) {// x.exp == y.exp - m = x.m + y.m; - err = x.err + y.err; - exp = x.exp; - // normal(); - } else {// x.exp < y.exp - if (!y.err) { - m = x.m + chunkShift(y.m, - expDiff); - err = x.err; - exp = x.exp; - } else {// y.err > 0 - m = chunkShift(x.m, expDiff) + y.m; - err = y.err + 5; - exp = y.exp; - // normal(); - } - } - // Call normalization globally -- IP 10/9/98 - normal(); -} - -// Subtraction -void BigFloatRep::sub(const BigFloatRep& x, const BigFloatRep& y) { - long expDiff = x.exp - y.exp; - - if (expDiff > 0) {// x.exp > y.exp - if (!x.err) { - m = chunkShift(x.m, expDiff) - y.m; - err = y.err; - exp = y.exp; - } else {// x.err > 0 - m = x.m - chunkShift(y.m, - expDiff); - err = x.err + 5; - exp = x.exp; - // normal(); - } - } else if (!expDiff) { - m = x.m - y.m; - err = x.err + y.err; - exp = x.exp; - // normal(); - } else { // x.exp < y.exp - if (!y.err) { - m = x.m - chunkShift(y.m, - expDiff); - err = x.err; - exp = x.exp; - } else {// y.err > 0 - m = chunkShift(x.m, expDiff) - y.m; - err = y.err + 5; - exp = y.exp; - // normal(); - } - } - // Call normalization globally -- IP 10/9/98 - normal(); -} - -void BigFloatRep::mul(const BigFloatRep& x, const BigFloatRep& y) { - m = x.m * y.m; - exp = x.exp + y.exp; - // compute error (new code, much faster. Zilin Du, Nov 2003) - if (x.err == 0 && y.err == 0) { - err = 0; - eliminateTrailingZeroes(); - } else { - BigInt bigErr(0); - if (y.err != 0) - bigErr += abs(x.m)*y.err; - if (x.err != 0) - bigErr += abs(y.m)*x.err; - if (x.err !=0 && y.err != 0) - bigErr += x.err*y.err; - bigNormal(bigErr); - } -} -// BigFloat div2 will half the value of x, exactly with NO error -// REMARK: should generalize this to dividing by any power of 2 -// We need this in our use of BigFloats to maintain isolation -// intervals (e.g., in Sturm sequences) --Chee/Vikram 4/2003 -// -void BigFloatRep :: div2(const BigFloatRep& x) { - if (isEven(x.m)) { - m = (x.m >> 1); - exp = x.exp ; - } else { - m = (x.m << static_cast(CHUNK_BIT-1)); - exp = x.exp -1; - } -} - -// Converts a BigFloat interval into one BigFloat with almost same error bound -// This routine ignores the errors in inputs a and b. -// But you cannot really ignore them since, they are taken into account -// when you compute "r.sub(a,b)"... -void BigFloatRep::centerize(const BigFloatRep& a, const BigFloatRep& b) { - if ((a.m == b.m) && (a.err == b.err) && (a.exp == b.exp)) { - m = a.m; - err = a.err; - exp = a.exp; - return; - } - - BigFloatRep r; - r.sub(a, b); - r.div2(r); - - //setup mantissa and exponent, but not error bits - // But this already sets the error bits? Chee - add(a,b); - div2(*this); - // error bits = ceil ( B^{-exp}*|a-b|/2 ) - - // bug fixed: possible overflow on converting - // Zilin & Vikram, 08/24/04 - // err = 1 + longValue(chunkShift(r.m, r.exp - exp)); - BigInt E = chunkShift(r.m, r.exp - exp); - bigNormal(E); -} - -// BigFloat Division, computing x/y: -// Unlike +,-,*, this one takes a relative precision bound R -// Note that R is only used when x and y are error-free! -// (This remark means that we may be less efficient than we could be) -// -// Assert( R>0 && R< CORE_Infty ) -// -void BigFloatRep :: div(const BigFloatRep& x, const BigFloatRep& y, - const extLong& R) { - if (!y.isZeroIn()) { // y.m > y.err, so we are not dividing by 0 - if (!x.err && !y.err) { - if (R < 0 || R.isInfty()) //Oct 9, 2002: fixed major bug! [Zilin/Chee] - div(x.m, y.m, defBFdivRelPrec, CORE_posInfty); - else - div(x.m, y.m, R, CORE_posInfty); - exp += x.exp - y.exp; // chen: adjust exp. - } else {// x.err > 0 or y.err > 0 - BigInt bigErr, errRemainder; - - if (x.isZeroIn()) { // x.m <= x.err - m = 0; - exp = x.exp - y.exp; - - div_rem(bigErr, errRemainder, abs(x.m) + static_cast(x.err), - abs(y.m) - static_cast(y.err)); - } else { // x.m > x.err - long lx = bitLength(x.m); - long ly = bitLength(y.m); - long r; - - if (!x.err) // x.err == 0 and y.err > 0 - r = ly + 2; - else if(!y.err) // x.err > 0 and y.err == 0 - r = lx + 2; - else // x.err > 0 and y.err > 0 - r = lx < ly ? lx + 2: ly + 2; - - long t = chunkFloor(- r + lx - ly - 1); - BigInt remainder; - - div_rem(m, remainder, chunkShift(x.m, - t), y.m); - exp = t + x.exp - y.exp; - - long delta = ((t > 0) ? 2 : 0); - - // Chen Li: 9/9/99 - // here again, it use ">>" operator with a negative - // right operand. So the result is not well defined. - // Erroneous code: - // divide(abs(remainder) + (static_cast(x.err) >> bits(t)) - // + delta + static_cast(y.err) * abs(m), - // abs(y.m) - static_cast(y.err), - // bigErr, - // errRemainder); - // New code: - BigInt errx_over_Bexp = x.err; - long bits_Bexp = bits(t); - if (bits_Bexp >= 0) { - errx_over_Bexp >>= bits_Bexp; - } else { - errx_over_Bexp <<= (-bits_Bexp); - } - - // divide(abs(remainder) + errx_over_Bexp - // + delta + static_cast(y.err) * abs(m), - // abs(y.m) - static_cast(y.err), - // bigErr, - // errRemainder); - div_rem(bigErr, errRemainder, - abs(remainder) + errx_over_Bexp + delta + static_cast(y.err) * abs(m), - abs(y.m) - static_cast(y.err)); - } - - if (sign(errRemainder)) - ++bigErr; - - bigNormal(bigErr); - } - } else {// y.m <= y.err - core_error("BigFloat error: possible zero divisor.", - __FILE__, __LINE__, true); - } - - // Call normalization globally -- IP 10/9/98 - // normal(); -- Chen: after calling bigNormal, this call is redundant. -}// BigFloatRep::div - -// squareroot for BigInt argument, without initial approximation -// sqrt(x,a) computes sqrt of x to absolute precision a. -// -- this is where Newton is applied -// -- this is called by BigFloatRep::sqrt(BigFloat, extLong) -void BigFloatRep::sqrt(const BigInt& x, const extLong& a) { - sqrt(x, a, BigFloat(x, 0, 0)); -} // sqrt(BigInt x, extLong a) , without initial approx - -// sqrt(x,a,A) where -// x = bigInt whose sqrt is to be computed -// a = absolute precision bound -// A = initial approximation in BigFloat -// -- this is where Newton is applied -// -- it is called by BigFloatRep::sqrt(BigFloatRep, extLong, BigFloat) -void BigFloatRep::sqrt(const BigInt& x, const extLong& a, const BigFloat& A) { - if (sign(x) == 0) { - m = 0; - err = 0; - exp = 0; - } else if (x == 1) { - m = 1; - err = 0; - exp = 0; - } else {// main case - // here is where we use the initial approximation - m = A.m(); - err = 0; - exp = A.exp(); - - BigFloatRep q, z; - extLong aa; - // need this to make sure that in case the - // initial approximation A is less than sqrt(x) - // then Newton iteration will still proceed at - // least one step. - bool firstTime = true; - for (;;) { - aa = a - bits(exp); - q.div(x, m, CORE_posInfty, aa); - q.err = 0; - q.exp -= exp; - - z.sub(*this, q); // this=current approximation, so z = this - q - /*if (sign(z.m) <= 0 || z.MSB() < - a) // justification: see Koji's - break; // thesis (p. 28) which states - // that we can exit when - // " (*this) <= q + 2**(-a)" - */ - // The preceding code is replaced by what follows: - if (z.MSB() < -a) - break; - if (sign(z.m) <= 0) { - if (firstTime) - firstTime = false; - else - break; - } - - z.add(*this, q); - // Chen Li: a bug fixed here. - // m = z.m >> 1; - // err = 0; - // exp = z.exp; - if ((z.m > 1) && isEven(z.m)) { - m = z.m >> 1; // exact division by 2 - err = 0; - exp = z.exp; - } else { // need to shift left before division by 2 - m = chunkShift(z.m, 1) >> 1; - err = 0; - exp = z.exp - 1; - }//else - }//for - }//else -} // sqrt of BigInt, with initial approx - -// MAIN ENTRY INTO SQRT FUNCTION (BIGFLOAT ARGUMENT, WITHOUT INITIAL APPROX) -void BigFloatRep::sqrt(const BigFloatRep& x, const extLong& a) { - sqrt(x, a, BigFloat(x.m, 0, x.exp)); -} //sqrt(BigFloat, extLong a) - -// MAIN ENTRY INTO SQRT FUNCTION (BIGFLOAT ARGUMENT WITH INITIAL APPROXIMATION) -void BigFloatRep::sqrt(const BigFloatRep& x, const extLong& a, const BigFloat& A) { - // This computes the sqrt of x to absolute precision a, starting with - // the initial approximation A - if (sign(x.m) >= 0) { // x.m >= 0 - int delta = x.exp & 1; // delta=0 if x.exp is even, otherwise delta=1 - - if (x.isZeroIn()) { // x.m <= x.err - m = 0; - if (!x.err) - err = 0; - else { // x.err > 0 - err = (long)(std::sqrt((double)x.err)); - err++; - err <<= 1; - if (delta) - err <<= HALF_CHUNK_BIT; - } - exp = x.exp >> 1; - normal(); - } else { - long aExp = A.exp() - (x.exp >> 1); - BigFloat AA( chunkShift(A.m(), delta), 0, aExp); - - if (!x.err) { // x.m > x.err = 0 (ERROR FREE CASE) - BigFloatRep z; - extLong ppp; - if (a.isInfty()) //Oct 9, 2002: fixed major bug! [Zilin/Chee] - ppp = defBFsqrtAbsPrec; - else - ppp = a + EXTLONG_ONE; - extLong absp = ppp + bits(x.exp >> 1); - - z.sqrt(chunkShift(x.m, delta), absp, AA); // call sqrt(BigInt, a, AA) - - long p = (absp + bits(z.exp)).asLong(); - - // Next, normalize the error: - if (p <= 0) { - m = z.m; - // Chen Li: a bug fixed - // BigInt bigErr = 1 << (-p); - BigInt bigErr(1); - bigErr = bigErr << static_cast(-p); - exp = z.exp + (x.exp >> 1); - bigNormal(bigErr); - } else { // p > 0 - m = chunkShift(z.m, chunkCeil(p)); - long r = CHUNK_BIT - 1 - (p + CHUNK_BIT - 1) % CHUNK_BIT; -#ifdef CORE_DEBUG - assert(r >= 0); -#endif - - err = 1 >> r; - exp = - chunkCeil(ppp.asLong()); - normal(); - } - } else { // x.m > x.err > 0 (mantissa has error) - BigFloatRep z; - extLong absp=-flrLg(x.err)+bitLength(x.m)-(bits(delta) >> 1)+EXTLONG_FOUR; - - z.sqrt(chunkShift(x.m, delta), absp, AA); - - long qqq = - 1 + (bitLength(x.m) >> 1) - delta * HALF_CHUNK_BIT; - long qq = qqq - clLg(x.err); - long q = qq + bits(z.exp); - - if (q <= 0) { - m = z.m; - long qqqq = - qqq - bits(z.exp); - // Chen Li (09/08/99), a bug fixed here: - // BigInt bigErr = x.err << - qqqq; - // when (-qqqq) is negative, the result is not correct. - // how "<<" and ">>" process negative second operand is - // not well defined. Seems it just take it as a unsigned - // integer and extract the last few bits. - // x.err is a long number which easily overflows. - // From page 22 of Koji's paper, I think the exponent is - // wrong here. So I rewrote it as: - BigInt bigErr = x.err; - if (qqqq >= 0) { - bigErr <<= qqqq; - } else { - bigErr >>= (-qqqq); - ++bigErr; // we need to keep its ceiling. - } - - exp = z.exp + (x.exp >> 1); - bigNormal(bigErr); - } else { // q > 0 - m = chunkShift(z.m, chunkCeil(q)); - long r = CHUNK_BIT - 1 - (q + CHUNK_BIT - 1) % CHUNK_BIT; -#ifdef CORE_DEBUG - assert(r >= 0); -#endif - - err = 1 >> r; - exp = (x.exp >> 1) - chunkCeil(qq); - normal(); - } - } // end of case with error in mantissa - }//else - } else - core_error("BigFloat error: squareroot called with negative operand.", - __FILE__, __LINE__, true); -} //sqrt with initial approximation - -// compareMExp(x) -// returns 1 if *this > x -// 0 if *this = x, -// -1 if *this < x, -// -// Main comparison method for BigFloat -// This is called by BigFloat::compare() -// BE CAREFUL: The error bits are ignored! -// Need another version if we want to take care of error bits - -int BigFloatRep :: compareMExp(const BigFloatRep& x) const { - int st = sign(m); - int sx = sign(x.m); - - if (st > sx) - return 1; - else if (st == 0 && sx == 0) - return 0; - else if (st < sx) - return - 1; - else { // need to compare m && exp - long expDiff = exp - x.exp; - - if (expDiff > 0) // exp > x.exp - return cmp(chunkShift(m, expDiff), x.m); - else if (!expDiff) - return cmp(m, x.m); - else // exp < x.exp - return cmp(m, chunkShift(x.m, - expDiff)); - } -} - -// 3/6/2000: -// This is a private function used by BigFloatRep::operator<< -// to get the exact value -// of floor(log10(M * 2^ e)) where E is an initial guess. -// We will return the correct E which satisfies -// 10^E <= M * 2^e < 10^{E+1} -// But we convert this into -// mm <= M < 10.mm - -long BigFloatRep :: adjustE( long E, BigInt M, long ee) const { - if (M<0) - M=-M; - BigInt mm(1); - if (ee > 0) - M = (M<(ee)); - else - mm = (mm << static_cast(-ee)); - if (E > 0) - mm *= (FiveTo(E)<< static_cast(E)); - else - M *= (FiveTo(-E) << static_cast(-E)); - - if (M < mm) { - do { - E--; - M *= 10; - } while (M < mm); - } else if (M >= 10*mm) { - mm *= 10; - do { - E++; - mm *= 10; - } while (M >= mm); - } - return E; -} - -BigFloatRep::DecimalOutput -BigFloatRep::toDecimal(unsigned int width, bool Scientific) const { - BigFloatRep::DecimalOutput decOut; // to be returned - if (err > 0) { - decOut.isExact = false; - } else { // err == 0 - decOut.isExact = true; - } - - if (err > 0 && err >= abs(m)) { - // if err is larger than mantissa, sign and significant values - // can not be determined. - core_error("BigFloat error: Error is too big!", - __FILE__, __LINE__, false); - decOut.rep = "0.0e0"; // error is too big - decOut.isScientific = false; - decOut.noSignificant = 0; - decOut.errorCode = 1; // sign of this number is unknown - return decOut; - } - - decOut.sign = sign(m); - decOut.errorCode = 0; - - BigInt M(m); // temporary mantissa - long lm = bitLength(M); // binary length of mantissa - long e2 = bits(exp); // binary shift length represented by exponent - long le = clLg(err); // binary length of err - if (le == -1) - le = 0; - - long L10 = 0; - if (M != 0) { - L10 = (long)std::floor((lm + e2) / lgTenM); - L10 = adjustE(L10, m, e2); // L10: floor[log10(M 2^(e2))], M != 0 - } else { - L10 = 0; - } - // Convention: in the positional format, when the output is - // the following string of 8 characters: - // (d0, d1, d2, d3, ".", d4, d5, d6, d7) - // then the decimal point is said to be in the 4th position. - // E.g., (d0, ".", d1, d2) has the decimal point in the 1st position. - // The value of L10 says that the decimal point of output should be at - // the (L10 + 1)st position. This is - // true regardingless of whether M = 0 or not. For zero, we output - // {0.0*} so L10=0. In general, the |value| is less than 10 - // if and only if L10 is 0 and the - // decimal point is in the 1st place. Note that L10 is defined even if - // the output is an integer (in which case it does not physically appear - // but conceptually terminates the sequence of digits). - - // First, get the decimal representaion of (m * B^(exp)). - if (e2 < 0) { - M *= FiveTo(-e2); // M = x * 10^(-e2) - } else if (e2 > 0) { - M <<= e2; // M = x * 2^(e2) - } - - std::string decRep = M.get_str(); - // Determine the "significant part" of this string, i.e. the part which - // is guaranteed to be correct in the presence of error, - // except that the last digit which might be subject to +/- 1. - - if (err != 0) { // valid = number of significant digits - unsigned long valid = floorlg10(m) - (long)std::floor(std::log10(float(err))); - if (decRep.length() > valid) { - decRep.erase(valid); - } - } - - // All the digits in decM are correct, except the last one might - // subject to an error +/- 1. - - if ((decRep[0] == '+') || (decRep[0] == '-')) { - decRep.erase(0, 1); - } - - // Second, make choice between positional representation - // and scientific notation. Use scientific notation when: - // 0) if scientific notation flag is on - // 1) err * B^exp >= 1, the error contribute to the integral part. - // 2) (1 + L10) >= width, there is not have enough places to hold the - // positional representation, not including decimal point. - // 3) The distance between the first significant digit and decimal - // point is too large for the width limit. This is equivalent to - // Either ((L10 >= 0 and (L10 + 1) > width)) - // Or ((L10 < 0) and (-L10 + 1) > width). - - if (Scientific || - ((err > 0) && (le + e2) >= 0) || // if err*B^exp >= 1 - ((L10 >= 0) && (L10 + 1 >= (long)width )) || - ((L10 < 0) && (-L10 + 1 > (long)width ))) { - // use scientific notation - decRep = round(decRep, L10, width); - decOut.noSignificant = width; - decRep.insert(1, "."); - if (L10 != 0) { - decRep += 'e'; - if (L10 > 0) { - decRep += '+'; - } else { // L10 < 0 - decRep += '-'; - } - char eBuf[48]; // enought to hold long number L10 - int ne = 0; - if ((ne = sprintf(eBuf, "%ld", labs(L10))) >= 0) { - eBuf[ne] = '\0'; - } else { - //perror("BigFloat.cpp: Problem in outputing the exponent!"); - core_error("BigFloat error: Problem in outputing the exponent", - __FILE__, __LINE__, true); - } - decRep += eBuf; - decOut.isScientific = true; - } - } else { - // use conventional positional notation. - if (L10 >= 0) { // x >= 1 or x == 0 and L10 + 1 <= width - // round when necessary - if (decRep.length() > width ) { - decRep = round(decRep, L10, width ); - if (decRep.length() > width ) { - // overflow happens! use scientific notation - return toDecimal(width, true); - } - } - decOut.noSignificant = decRep.length(); - if (L10 + 1 < (long)width ) { - decRep.insert(L10 + 1, "."); - } else { // L10 + 1 == width - // do nothing - } - } else { // L10 < 0, 0 < x < 1 - // (-L10) leading zeroes, including one to the left of decimal dot - // need to be added in beginning. - decRep = std::string(-L10, '0') + decRep; - // then round when necessary - if (decRep.length() > width ) { - decRep = round(decRep, L10, width ); - // cannot overflow since there are L10 leading zeroes. - } - decOut.noSignificant = decRep.length() - (-L10); - decRep.insert(1, "."); - } - decOut.isScientific = false; - } -#ifdef CORE_DEBUG - assert(decOut.noSignificant >= 0); -#endif - - decOut.rep = decRep; - return decOut; -}//toDecimal - -std::string BigFloatRep::round(std::string inRep, long& L10, unsigned int width) const { - // round inRep so that the length would not exceed width. - if (inRep.length() <= width) - return inRep; - - int i = width; // < length - bool carry = false; - - if ((inRep[i] >= '5') && (inRep[i] <= '9')) { - carry = true; - i--; - while ((i >= 0) && carry) { - if (carry) { - inRep[i] ++; - if (inRep[i] > '9') { - inRep[i] = '0'; - carry = true; - } else { - carry = false; - } - } - i-- ; - } - - if ((i < 0) && carry) { // overflow - inRep.insert(inRep.begin(), '1'); - L10 ++; - width ++; - } - } - - return inRep.substr(0, width); -}//round(string,width) - - -// This function fromString(str, prec) is similar to the -// constructor Real(char * str, extLong prec) -// See the file Real.cc for the differences - -void BigFloatRep :: fromString(const char *str, const extLong & prec ) { - // NOTE: prec defaults to defBigFloatInputDigits (see BigFloat.h) - // check that prec is not INFTY - if (prec.isInfty()) - core_error("BigFloat error: infinite precision not allowed", - __FILE__, __LINE__, true); - - const char *e = strchr(str, 'e'); - int dot = 0; - long e10 = 0; - if (e != NULL) - e10 = atol(e+1); // e10 is decimal precision of the input string - // i.e., input is A/10^{e10}. - else { - e = str + strlen(str); -#ifdef CORE_DEBUG - assert(*e == '\0'); -#endif - - } - - const char *p = str; - if (*p == '-' || *p == '+') - p++; - m = 0; - exp = 0; - - for (; p < e; p++) { - if (*p == '.') { - dot = 1; - continue; - } - m = m * 10 + (*p - '0'); - if (dot) - e10--; - } - - BigInt one = 1; - long t = (e10 < 0) ? -e10 : e10; - BigInt ten = FiveTo(t) * (one << static_cast(t)); - - // HERE IS WHERE WE USE THE SYSTEM CONSTANT - // defBigFloatInputDigits - // Note: this constant is rather similar to defInputDigits which - // is used by Real and Expr for controlling - // input accuracy. The difference is that defInputDigits can - // be CORE_INFTY, but defBigFloatInputDigits must be finite. - - if (e10 < 0) - div(m, ten, CORE_posInfty, 4 * prec); - else - m *= ten; - if (*str == '-') - m = -m; -}//BigFloatRep::fromString - -std::istream& BigFloatRep :: operator >>(std::istream& i) { - int size = 20; - char *str = new char[size]; - char *p = str; - char c; - int d = 0, e = 0, s = 0; - // d=1 means dot is found - // e=1 means 'e' or 'E' is found - // int done = 0; - - // Chen Li: fixed a bug, the original statement is - // for (i.get(c); c == ' '; i.get(c)); - // use isspace instead of testing c == ' ', since it must also - // skip tab, catridge/return, etc. - // Change to: - // int status; - do { - c = i.get(); - } while (isspace(c)); /* loop if met end-of-file, or - char read in is white-space. */ - // Chen Li, "if (c == EOF)" is unsafe since c is of char type and - // EOF is of int tyep with a negative value -1 - if (i.eof()) { - i.clear(std::ios::eofbit | std::ios::failbit); - return i; - } - - // the current content in "c" should be the first non-whitespace char - if (c == '-' || c == '+') { - *p++ = c; - i.get(c); - } - - for (; isdigit(c) || (!d && c=='.') || - (!e && ((c=='e') || (c=='E'))) || (!s && (c=='-' || c=='+')); i.get(c)) { - if (!e && (c == '-' || c == '+')) - break; - // Chen Li: put one more rule to prohibite input like - // xxxx.xxxe+xxx.xxx: - if (e && (c == '.')) - break; - if (p - str == size) { - char *t = str; - str = new char[size*2]; - memcpy(str, t, size); - delete [] t; - p = str + size; - size *= 2; - } -#ifdef CORE_DEBUG - assert((p-str) < size); -#endif - - *p++ = c; - if (c == '.') - d = 1; - // Chen Li: fix a bug -- the sign of exponent can not happen before - // the character "e" appears! It must follow the "e' actually. - // if (e || c == '-' || c == '+') s = 1; - if (e) - s = 1; - if ((c == 'e') || (c=='E')) - e = 1; - } - - // chenli: make sure that the p is still in the range - if (p - str >= size) { - int len = p - str; - char *t = str; - str = new char[len + 1]; - memcpy(str, t, len); - delete [] t; - p = str + len; - } - -#ifdef CORE_DEBUG - assert(p - str < size); -#endif - - *p = '\0'; - i.putback(c); - fromString(str); - delete [] str; - return i; -}//operator >> - - -// BigFloatRep::toDouble() -// converts the BigFloat to a machine double -// This is a dangerous function as the method -// is silent when it does not fit into a machine double! -// ToDo: fix this by return a machine NaN, +/- Infinity, +/- 0, -// when appropriate. -// Return NaN when error is larger than mantissa -// Return +/- Infinity if BigFloat is too big -// Return +/- 0 if BigFloat is too small -#ifdef _MSC_VER -#pragma warning(disable: 4723) -#endif -double BigFloatRep :: toDouble() const { - if (m == 0) - return (sign(m) * 0.0); - - long e2 = bits(exp); - long le = clLg(err); // if err=0, le will be -1 - if (le == -1) - le = 0; - - BigInt M = m >> static_cast(le);// remove error bits in mantissa - - // Below, we want to return NaN by computing 0.0/0.0. - // To avoid compiler warnings about divide by zero, we do this: - - double foolCompilerZero; - foolCompilerZero = 0.0; - - // COMMENT: we should directly store the - // special IEEE values NaN, +/-Infinity, +/-0 in the code!! - - if (M == 0) - return ( 0.0/foolCompilerZero ) ; // return NaN - - e2 += le; // adjust exponent corresponding to error bits - - int len = bitLength(M) - 53; // this is positive if M is too large - - if (len > 0) { - M >>= len; - e2 += len; - } - - double tt = doubleValue(M); - - int ee = e2 + bitLength(M) - 1; // real exponent. - - if (ee >= 1024) // overflow! - return ( sign(m)/foolCompilerZero ); // return a signed infinity - - if (ee <= -1075) // underflow! - // NOTE: if (-52 < ee <= 0) get denormalized number - return ( sign(m) * 0.0 ); // return signed zero. - - // Execute this loop if e2 < 0; - for (int i = 0; i > e2; i--) - tt /= 2; - - // Execute this loop if e2 > 0; - for (int j = 0; j < e2; j++) - tt *= 2; - - return tt; -}//toDouble -#ifdef _MSC_VER -#pragma warning(default: 4723) -#endif -BigInt BigFloatRep::toBigInt() const { - long e2 = bits(exp); - long le = clLg(err); - if (le == -1) - le = 0; -#ifdef CORE_DEBUG - assert (le >= 0); -#endif - - BigInt M = m >> static_cast(le); // discard the contaminated bits. - e2 += le; // adjust the exponent - - if (e2 < 0) - return M >> static_cast(-e2); - else if (e2 > 0) - return M << static_cast(e2); - else - return M; -} - -long BigFloatRep :: toLong() const { - // convert a BigFloat to a long integer, rounded toward -\infty. - long e2 = bits(exp); - long le = clLg(err); -#ifdef CORE_DEBUG - assert (le >= 0); -#endif - - BigInt M = m >> static_cast(le); // discard the contaminated bits. - e2 += le; // adjust the exponent - long t; - if (e2 < 0) - t = ulongValue(M >> static_cast(-e2)); - else if (e2 > 0) - t = ulongValue(M << static_cast(e2)); - else - t = ulongValue(M); - // t = M.as_long(); - // Note: as_long() will return LONG_MAX in case of overflow. - - return t; -} - -// pow(r,n) function for BigFloat -// Note: power(r,n) calls pow(r,n) -BigFloat pow(const BigFloat& r, unsigned long n) { - if (n == 0) - return BigFloat(1); - else if (n == 1) - return r; - else { - BigFloat x = r; - while ((n % 2) == 0) { // n is even - x *= x; - n >>= 1; - } - BigFloat u = x; - while (true) { - n >>= 1; - if (n == 0) - return u; - x *= x; - if ((n % 2) == 1) // n is odd - u *= x; - } - //return u; // unreachable - } -}//pow - -// experimental -BigFloat root(const BigFloat& x, unsigned long k, - const extLong& a, const BigFloat& A) { - if (x.sign() == 0) { - return BigFloat(0); - } else if (x == 1) { - return BigFloat(1); - } else { - BigFloat q, del, zz; - BigFloat z = A; - BigFloat bk = long(k); - for (; ;) { - zz = pow(z, k-1); - q = x.div(zz, a); - q.makeExact(); - del = z - q; - del.makeExact(); - if (del.MSB() < -a) - break; - z = ((bk-1)*z + q).div(bk, a); - // newton's iteration: z_{n+1}=((k-1)z_n+x/z_n^{k-1})/k - z.makeExact(); - } - return z; - } -}//root - -} //namespace CORE +#endif \ No newline at end of file diff --git a/CGAL_Core/src/CGAL_Core/CoreAux.cpp b/CGAL_Core/src/CGAL_Core/CoreAux.cpp index 8d0c801cb46..adc3998cb57 100644 --- a/CGAL_Core/src/CGAL_Core/CoreAux.cpp +++ b/CGAL_Core/src/CGAL_Core/CoreAux.cpp @@ -33,181 +33,9 @@ * $Id$ ***************************************************************************/ -#include "CGAL/CORE/CoreAux.h" -#include +#ifndef CGAL_HEADER_ONLY -namespace CORE { +#include +#include -//////////////////////////////////////////////////////////// -// More useful functions to implement: -// -// To convert digits into bits: -// given X, compute X * log_2(10) -// To convert bits into digits: -// given X, compute X * log_10(2) -// -//////////////////////////////////////////////////////////// - -//////////////////////////////////////////////////////////// -// flrLg(x) -// returns floor log base 2 of abs(x) -// CONVENTION lg(0) = -1 (Slight improvement, Zilin/Chee 8/5/02) -//////////////////////////////////////////////////////////// -int flrLg(long x) { - if (x == LONG_MIN) { - // special treatment as -LONG_MIN would be not representable as "long" - return LONG_BIT - 1; - } else { - // 1 <= |x| <= LONG_MAX - if (x < 0) - x = - x; - - int lg = -1; - while (x > 0) { - lg++; - x >>= 1; - } - return lg; - } -} - -//////////////////////////////////////////////////////////// -// floor log base 2 of unsigned long x -// CONVENTION lg(0) = -1 (Slight improvement, Zilin/Chee 8/5/02) -//////////////////////////////////////////////////////////// -int flrLg(unsigned long x) { - int lg = -1; - while (x > 0) { - lg++; - x >>= 1; - } - return lg; -} - -//////////////////////////////////////////////////////////// -// ceiling log base 2 of abs(x) -// CONVENTION lg(0) = -1 (Slight improvement, Zilin/Chee 8/5/02) -//////////////////////////////////////////////////////////// -int clLg(long x) { - if (x == LONG_MIN) - return LONG_BIT - 1; - if (x < 0) - x = -x; // use absolute value - if (x > (LONG_MAX >> 1)) // the leading data bit is 1 - return (LONG_BIT - 1); // exclude the sign bit - if (x >= 2) - return flrLg((unsigned long)((x << 1) - 1)); - // SINCE ceilLog_2(x) = floorLog_2(2x-1) for x>=2 - if (x == 1) - return 0; - return -1; // x must be 0 here -} - -//////////////////////////////////////////////////////////// -// ceiling log base 2 of unsigned long x -// CONVENTION lg(0) = -1 -//////////////////////////////////////////////////////////// -int clLg(unsigned long x) { - if (x > (ULONG_MAX >> 1)) // the leading bit is 1 - return LONG_BIT; - if (x >= 2) - return flrLg((x << 1) - 1); - // SINCE ceilLog_2(x) = floorLog_2(2x-1) for x>=2. - if (x == 1) - return 0; - return -1; // x must be equal to 0 -} - -/// gcd for machine type long -/** This is needed when we construct Polynomials with int or long coefficients */ -long gcd(long m, long n) { - if (m == 0) - return core_abs(n); - if (n == 0) - return core_abs(m); - long p = core_abs(m); - long q = core_abs(n); - if (p0) { - long r = p % q; - p = q; - q = r; - } - return p; -} - -// return a gmp_randstate_t structure -gmp_randstate_t* getRandstate() { - static gmp_randstate_t rstate; - static bool initialized = false; - if (!initialized) { - gmp_randinit(rstate, GMP_RAND_ALG_DEFAULT, 32L); - initialized = true; - } - return &rstate; -} - -// char* core_itoa(int n, char* buffer) -// returns a pointer to the null-terminated string in buffer -// NOTES: -// 0. Buffer size should be 17 bytes (resp., 33 bytes, 65 bytes) on 16-bit -// (resp., 32-bit, 64-bit) machines. Formula: 1+sizeof(int)*8 bytes. -// 1. itoa(...) is available on some stdlib.h, but it is -// not defined by ANSI-C and so not all compilers support it. -// 2. Our use of sprintf(...) to do the job is known to -// be inefficient, but this is hardly critical for our usage. -// 3. A more general program should take a 3rd argument (the radix of -// output number). We assume radix 10. -char * core_itoa(int n, char* buffer) { - std::sprintf(buffer, "%d", n); - return buffer; -} - -/// implements the "integer mantissa" function -// (See CORE_PATH/progs/ieee/frexp.cpp for details) -double IntMantissa(double d) { - int e; - return std::ldexp(std::frexp(d, &e), 53); -} - -/// implements the "integer exponent" function -// (See CORE_PATH/progs/ieee/frexp.cpp for details) -int IntExponent(double d) { - int e; - std::frexp(d, &e); - return e-53; -} - -/// CORE_DIAGFILE is file name for core_error(..) output. -const char* CORE_DIAGFILE = "Core_Diagnostics"; // global file name - -/// core_error is the method to write Core Library warning or error messages -/** Both warnings and errors are written to a file called CORE_DIAGFILE. - * But errors are also written on std:cerr (similar to std::perror()). - * */ -// Usage: core_error(message, file_with_error, line_number, err_type) -// where err_type=0 means WARNING, error_type=0 means ERROR -void core_error(std::string msg, std::string file, int lineno, bool err) { - std::ofstream outFile(CORE_DIAGFILE, std::ios::app); // open to append - if (!outFile) { - // perror("CORE ERROR: cannot open Core Diagnostics file"); - std::cerr << "CORE ERROR: can't open Core Diagnostics file"< +#include +#endif diff --git a/CGAL_Core/src/CGAL_Core/CoreIO.cpp b/CGAL_Core/src/CGAL_Core/CoreIO.cpp index a0ae4003802..a03cb47cd03 100644 --- a/CGAL_Core/src/CGAL_Core/CoreIO.cpp +++ b/CGAL_Core/src/CGAL_Core/CoreIO.cpp @@ -29,430 +29,6 @@ * $Id$ ***************************************************************************/ -#include -#include - -namespace CORE { - -void core_io_error_handler(const char *f, const char *m) { - std::cout << "\n error_handler"; - std::cout << "::" << f << "::" << m << "\n"; - std::cout.flush(); - std::abort(); -} - -void core_io_memory_handler(char *t, const char *f, const char *m) { - if (t == NULL) { - std::cout << "\n memory_handler"; - std::cout << "::" << f << "::" << m; - std::cout << "memory exhausted\n"; - std::cout.flush(); - std::abort(); - } -} - -// s has size old_size and will be resized to new_size. -void allocate (char * &s, int old_size, int new_size) { - if (old_size > new_size) - old_size = new_size; - - if (s == NULL) - old_size = 0; - - char *t = new char[new_size]; - core_io_memory_handler(t, "CoreIO", "allocate::out of memory error"); - - int i; - for (i = 0; i < old_size; i++) - t[i] = s[i]; - - delete[] s; - s = t; -} - -// appends c to s at position pos. -// sz is the size of s -void append_char (char * &s, int & sz, int pos, char c) { - if (pos > sz) - core_io_error_handler("CoreIO", "append_char::invalid argument"); - - if (pos == sz) { - allocate(s, sz, 2*sz); - sz *= 2; - } - - s[pos] = c; -} - -// skip blanks, tabs, line breaks and comment lines -int skip_comment_line (std::istream & in) { - int c; - - do { - c = in.get(); - while ( c == '#' ) { - do { - c = in.get(); - } while ( c != '\n' ); - c = in.get(); - } - } while (c == ' ' || c == '\t' || c == '\n'); - - if (c == EOF) - core_io_error_handler("CoreIO::read_from_file()","unexpected end of file."); - - in.putback(c); - return c; -} - -// skips '\\' followed by '\n' -int skip_backslash_new_line (std::istream & in) { - int c = in.get(); - - while (c == '\\') { - c = in.get(); - - if (c == '\n') - c = in.get(); - else - core_io_error_handler("CoreIO::operator>>", "\\ must be immediately followed by new line."); - } - - return c; -} - -void read_string(std::istream& in, char* &buffer, int sz) { - int c, pos=0; - skip_comment_line(in); - - while ( (c = in.get()) != EOF ) { - if ( c == ' ' || c == '\t' || c == '\n' || c == '#') - break; - else - append_char(buffer, sz, pos++, c); - } - append_char(buffer, sz, pos, '\0'); -} - -void read_base_number(std::istream& in, BigInt& m, long length, long maxBits) { - char *buffer; - int size, offset; - int base; - bool is_negate; - - int c, pos = 0; - skip_comment_line(in); - - // read sign - c = in.get(); - if (c == '-') { - is_negate = true; - c = in.get(); - } else - is_negate = false; - - // read base and compute digits - if (c == '0') { - c = in.get(); - if (c == 'b') { - base = 2; - size = (maxBits == 0 || maxBits > length) ? length : maxBits; - offset = length - size; - } else if (c == 'x') { - base = 16; - size = (maxBits == 0) ? length : (maxBits+3) >> 2; - size = (size > length) ? length : size; - offset = (length - size) << 2; - } else { - base = 8; - size = (maxBits == 0) ? length : (maxBits+2) / 3; - size = (size > length) ? length : size; - offset = (length - size) * 3; - in.putback(c); - } - } else { - base = 10; - size = (maxBits == 0) ? length : (int)std::ceil(maxBits*std::log(2.0)/std::log(10.0)); - size = (size > length) ? length : size; - offset = length - size; - in.putback(c); - } - - buffer = new char[size+2]; - // read digits - for (int i=0; (i 0 && base != 10) { - m <<= offset; - } - - if (is_negate) - negate(m); -} - - -void write_base_number(std::ostream& out, char* buffer, int length, int base, int charsPerLine) { - // write big number in a format that gmp's mpz_set_str() can - // automatically recognize with argument base = 0. - if (base == 2) - out << "0b"; - else if (base == 16) - out << "0x"; - else if (base == 8) - out << '0'; - - // write big number in charsPerLine. - char* start, *end, c; - for (int i=0; i= length) - out << start; - else { - end = start + charsPerLine; - c = *end; - *end = '\0'; - - out << start << "\\\n"; - *end = c; - } - } -} - -void readFromFile(BigInt& z, std::istream& in, long maxLength) { - char *buffer; - long length; - - // check type name whether it is Integer or not. - buffer = new char[8]; - read_string(in, buffer, sizeof(buffer)); - if ( std::strcmp(buffer, "Integer") != 0) - core_io_error_handler("BigInt::read_from_file()","type name expected."); - delete[] buffer; - - // read the bit length field. - buffer = new char[100]; - read_string(in, buffer, sizeof(buffer)); - length = std::atol(buffer); - delete[] buffer; - - // read bigint - read_base_number(in, z, length, maxLength); -} - -void writeToFile(const BigInt& z, std::ostream& out, int base, int charsPerLine) { - BigInt c = abs(z); - - // get the absoulte value string - char* buffer = new char[mpz_sizeinbase(c.get_mp(), base) + 2]; - mpz_get_str(buffer, base, c.get_mp()); - int length = std::strlen(buffer); - - // write type name of big number and length - //out << "# This is an experimental big number format.\n"; - out << "Integer " << length << "\n"; - - // if bigint is negative, then write an sign '-'. - if ( sign(z) < 0 ) - out << '-'; - - write_base_number(out, buffer, length, base, charsPerLine); - out << "\n"; - delete[] buffer; -} - -void readFromFile(BigFloat& bf, std::istream& in, long maxLength) { - char *buffer; - long length; - long exponent; - BigInt mantissa; - - // check type name whether it is Float - buffer = new char[6]; - read_string(in, buffer, sizeof(buffer)); - if (std::strcmp(buffer, "Float") != 0) - core_io_error_handler("BigFloat::read_from_file()", "type name expected"); - delete[] buffer; - - // read base (default is 16384) - buffer = new char[8]; - read_string(in, buffer, sizeof(buffer)); - if (std::strcmp(buffer, "(16384)") != 0) - core_io_error_handler("BigFloat::read_from_file()", "base expected"); - delete[] buffer; - - // read the bit length field. - buffer = new char[100]; - read_string(in, buffer, sizeof(buffer)); - length = std::atol(buffer); - delete[] buffer; - - // read exponent - buffer = new char[100]; - read_string(in, buffer, sizeof(buffer)); - exponent = std::atol(buffer); - delete[] buffer; - - // read mantissa - read_base_number(in, mantissa, length, maxLength); - - // construct BigFloat - bf = BigFloat(mantissa, 0, exponent); -} - -void writeToFile(const BigFloat& bf, std::ostream& out, int base, int charsPerLine) { - BigInt c(CORE::abs(bf.m())); - - // get the absoulte value string - char* buffer = new char[mpz_sizeinbase(c.get_mp(), base) + 2]; - mpz_get_str(buffer, base, c.get_mp()); - int length = std::strlen(buffer); - - - // write type name, base, length - //out << "# This is an experimental Big Float format." << std::endl; - out << "Float (16384) " << length << std::endl; - // write exponent - out << bf.exp() << std::endl; - - // write mantissa - if ( CORE::sign(bf.m()) < 0 ) - out << '-'; - - write_base_number(out, buffer, length, base, charsPerLine); - out << '\n'; - delete[] buffer; -} - -/* Underconstruction ---- -void BigFloat::read_from_file2(std::istream& in, long maxLength) { - long length = 1024; - char *buffer; - - // check type name whether it is Float - buffer = new char[7]; - BigInt::read_string(in, buffer, sizeof(buffer)); - if (strcmp(buffer, "NFloat") != 0) - core_io_error_handler("BigFloat::read_from_file2()", "type name expected"); - delete[] buffer; - - // read base (default is 16) - buffer = new char[5]; - BigInt::read_string(in, buffer, sizeof(buffer)); - if (strcmp(buffer, "(16)") != 0) - core_io_error_handler("BigFloat::read_from_file2()", "base expected"); - delete[] buffer; - - // read length field - buffer = new char[100]; - BigInt::read_string(in, buffer, sizeof(buffer)); - - // get the length field if it is not null. - if (buffer[0] != '\0') { - length = atol(buffer); - if (maxLength > 0 && length >= maxLength) - length = maxLength; - } - delete[] buffer; - - // read exponent - buffer = new char[100]; - BigInt::read_string(in, buffer, sizeof(buffer)); - long exp16 = atol(buffer); - delete[] buffer; - - // read mantissa - buffer = new char[length+2]; - //BigInt::read_base_number(in, buffer, length); - - BigInt m16(buffer); - delete[] buffer; - - // convert to base CHUNK_BIT - exp16 = exp16 - length + 1; - if ( m16.is_negative() ) - exp16 ++; - - long tmp_exp = exp16 * 4; - long q = tmp_exp / CHUNK_BIT; - long r = tmp_exp % CHUNK_BIT; - if ( r < 0 ) { - r += CHUNK_BIT; - q --; - } - - BigInt mantissa = m16 << r; - long exponent = q; - - // construct BigFloat - if (--rep->refCount == 0) - delete rep; - - rep = new BigFloatRep(mantissa, 0, exponent); - rep->refCount++; - -} - -// write normal float -// now it assumed to write in hex base, i.e. B=2^4=16 -// (note: our default base B=2^(CHUNK_BIT)=2^14=16384 -void BigFloat::write_to_file2(std::ostream& out, int base, int charsPerLine) { - // convert to base 16. - long new_base = 4; // 2^4 = 16 - long tmp_exp = (rep->exp) * CHUNK_BIT; - long q = tmp_exp / new_base; - long r = tmp_exp % new_base; - std::cout << "CORE_DEBUG: q=" << q << ", r=" << r << std::endl; - if ( r < 0 ) { - r += new_base; - q--; - } - std::cout << "CORE_DEBUG: q=" << q << ", r=" << r << std::endl; - - BigInt m16 = (rep->m) << r; - - int size = mpz_sizeinbase(m16.I, base) + 2; - std::cout << "size=" << size << std::endl; - char* buffer = new char[size]; - - int length = bigint_to_string(m16, buffer, base); - std::cout << "length=" << length << std::endl; - - long exp16 = q + length - 1; - if ( m16.is_negative() ) - exp16 --; - - // write type name, base, length - out << "# This is an experimental Big Float format." << std::endl; - out << "NFloat (16) " << length << std::endl; - - // write exponent - out << exp16 << std::endl; - - // write mantissa - if ( m16.is_negative() ) { - out << '-'; - buffer ++; - } - - BigInt::write_base_number(out, buffer, length, base, charsPerLine); - out << '\n'; - delete[] buffer; -} -*/ - -} //namespace CORE - +#ifndef CGAL_HEADER_ONLY +#include +#endif \ No newline at end of file diff --git a/CGAL_Core/src/CGAL_Core/Expr.cpp b/CGAL_Core/src/CGAL_Core/Expr.cpp index 89751cd62f7..09f0108d0d1 100644 --- a/CGAL_Core/src/CGAL_Core/Expr.cpp +++ b/CGAL_Core/src/CGAL_Core/Expr.cpp @@ -34,1119 +34,9 @@ * $Id$ ***************************************************************************/ +#ifndef CGAL_HEADER_ONLY + #include -#include +#include -namespace CORE { - -#ifdef CORE_DEBUG_BOUND -unsigned int BFMSS_counter = 0; -unsigned int BFMSS_only_counter = 0; -unsigned int Measure_counter = 0; -unsigned int Measure_only_counter = 0; -unsigned int Cauchy_counter = 0; -unsigned int Cauchy_only_counter = 0; -unsigned int LiYap_counter = 0; -unsigned int LiYap_only_counter = 0; -unsigned int rootBoundHitCounter = 0; -unsigned int computeBoundCallsCounter = 0; -#endif - -const char* Add::name = "+"; -const char* Sub::name = "-"; - -/******************************************************** - * class Expr - ********************************************************/ -const Expr& Expr::getZero() { - static Expr Zero(0); - return Zero; -} -const Expr& Expr::getOne() { - static Expr One(1); - return One; -} - -// computes an interval comprising a pair of doubles -// Note: -// -// This function returns are two consecutive representable binary -// IEEE double values whichs contain the real value, but when you print out -// them, you might be confused by the decimal represention due to round. -// -void Expr::doubleInterval(double & lb, double & ub) const { - double d = doubleValue(); - if (! CGAL_CORE_finite(d)) { // if overflow, underflow or NaN - lb = ub = d; - return; - } - int sign = ((* this) -Expr(d)).sign(); - // Seems like doubleValue() always give a lower bound, - // so sign = 0 or 1 (never -1). - //std::cout << "Sign = " << sign << std::endl; - if (sign == 0) { - lb = ub = d; - return; - } - int exp; - frexp(d, & exp); // get the exponent of d - exp--; // the exp from frexp satisfies - // 2^{exp-1} <= d < 2^{exp} - // But, we want exp to satisfy - // 2^{exp} <= d < 2^{exp+1} - if (sign > 0) { - lb = d; - ub = d + ldexp(1.0, -52+exp); - return; - } else { - ub = d; - lb = d - ldexp(1.0, -52+exp); - return; - } -} - -// floor(e, sub) returns the floor(e), and puts the -// remainder into sub. -BigInt floor(const Expr& e, Expr &sub) { - if (e==0) { - return 0; - } - BigInt f = e.approx(CORE_INFTY, 2).BigIntValue(); - sub = e-f; - // Adjustment - if (sub<0) - ++sub, --f; - if (sub>=1) - --sub, ++f; - assert(sub >=0 && sub<1); // got an assertion error? (Chee 3/24/04) - return f; -} - -// Chenli: implemented algorithm from Goldberg's article. -// 7/01: Thanks to Serge Pashkov for fixing infinite loop when n=0. -Expr pow(const Expr& e, unsigned long n) { - if (n == 0) - return Expr(1); - else if (n == 1) - return e; - else { - Expr x = e; - while ((n % 2) == 0) { // n is even - x *= x; - n >>= 1; - } - Expr u = x; - while (true) { - n >>= 1; - if (n == 0) - return u; - x *= x; - if ((n % 2) == 1) // n is odd - u *= x; - } - //return u; // unreachable - } -}//pow - -NodeInfo::NodeInfo() : appValue(CORE_REAL_ZERO), appComputed(false), - flagsComputed(false), knownPrecision(CORE_negInfty), -#ifdef CORE_DEBUG - relPrecision(EXTLONG_ZERO), absPrecision(CORE_negInfty), numNodes(0), -#endif - // Most of the following data members don't need to be - // initialized here. - d_e(EXTLONG_ZERO), visited(false), sign(0), - uMSB(CORE_negInfty), lMSB(CORE_negInfty), - // length(0), - measure(EXTLONG_ZERO), high(EXTLONG_ZERO), low(EXTLONG_ONE), - lc(EXTLONG_ZERO), tc(EXTLONG_ZERO), - v2p(EXTLONG_ZERO), v2m(EXTLONG_ZERO), - v5p(EXTLONG_ZERO), v5m(EXTLONG_ZERO), - u25(EXTLONG_ZERO), l25(EXTLONG_ZERO), - ratFlag(0), ratValue(NULL) { } - -/******************************************************** - * class ExprRep - ********************************************************/ -// constructor -ExprRep::ExprRep() : refCount(1), nodeInfo(NULL), ffVal(0.0) { } - -// Computes the root bit bound of the expression. -// In effect, computeBound() returns the current value of low. -extLong ExprRep::computeBound() { - extLong measureBd = measure(); - // extLong cauchyBd = length(); - extLong ourBd = (d_e() - EXTLONG_ONE) * high() + lc(); - // BFMSS[2,5] bound. - extLong bfmsskBd; - if (v2p().isInfty() || v2m().isInfty()) - bfmsskBd = CORE_INFTY; - else - bfmsskBd = l25() + u25() * (d_e() - EXTLONG_ONE) - v2() - ceilLg5(v5()); - - // since we might compute \infty - \infty for this bound - if (bfmsskBd.isNaN()) - bfmsskBd = CORE_INFTY; - - extLong bd = core_min(measureBd, - // core_min(cauchyBd, - core_min(bfmsskBd, ourBd)); -#ifdef CORE_SHOW_BOUNDS - std::cout << "Bounds (" << measureBd << - "," << bfmsskBd << ", " << ourBd << "), "; - std::cout << "MIN = " << bd << std::endl; - std::cout << "d_e= " << d_e() << std::endl; -#endif - -#ifdef CORE_DEBUG_BOUND - // Some statistics about which one is/are the winner[s]. - computeBoundCallsCounter++; - int number_of_winners = 0; - std::cerr << " New contest " << std::endl; - if (bd == bfmsskBd) { - BFMSS_counter++; - number_of_winners++; - std::cerr << " BFMSS is the winner " << std::endl; - } - if (bd == measureBd) { - Measure_counter++; - number_of_winners++; - std::cerr << " measureBd is the winner " << std::endl; - } - /* if (bd == cauchyBd) { - Cauchy_counter++; - number_of_winners++; - std::cerr << " cauchyBd is the winner " << std::endl; - } - */ - if (bd == ourBd) { - LiYap_counter++; - number_of_winners++; - std::cerr << " ourBd is the winner " << std::endl; - } - - assert(number_of_winners >= 1); - - if (number_of_winners == 1) { - if (bd == bfmsskBd) { - BFMSS_only_counter++; - std::cerr << " BFMSSBd is the only winner " << std::endl; - } else if (bd == measureBd) { - Measure_only_counter++; - std::cerr << " measureBd is the only winner " << std::endl; - } - /* else if (bd == cauchyBd) { - Cauchy_only_counter++; - std::cerr << " cauchyBd is the only winner " << std::endl; - } */ - else if (bd == ourBd) { - LiYap_only_counter++; - std::cerr << " ourBd is the only winner " << std::endl; - } - } -#endif - - return bd; -}//computeBound() - -void ExprRep::reduceToBigRat(const BigRat& rat) { - Real value(rat); - - //appValue() = value; - appComputed() = false; // since appValue is not assigned until approx() is called - flagsComputed() = true; - knownPrecision() = CORE_negInfty; - -#ifdef CORE_DEBUG - relPrecision() = EXTLONG_ZERO; - absPrecision() = CORE_negInfty; - //numNodes() = numNodes(); -#endif - - d_e() = EXTLONG_ONE; - //visited() = e->visited(); - sign() = value.sign(); - uMSB() = value.MSB(); - lMSB() = value.MSB(); - // length() = value.length(); // fixed? original = 1 - measure() = value.height(); // measure <= height for rational value - - // BFMSS[2,5] bound. - value.ULV_E(u25(), l25(), v2p(), v2m(), v5p(), v5m()); - - extLong u_e = u25() + v2p(); - extLong l_e = l25() + v2m(); - - u_e = u_e + ceilLg5(v5p()); - l_e = l_e + ceilLg5(v5m()); - - if (l_e == EXTLONG_ZERO) { // no divisions introduced - high() = u_e; - low() = EXTLONG_ONE - u_e; // - (u_e - 1) - } else { - high() = u_e - l_e + EXTLONG_ONE; - low() = 2 - high(); - } - - lc() = l_e; - tc() = u_e; - - if (ratValue() == NULL) - ratValue() = new BigRat(rat); - else - *(ratValue()) = rat; -} - -// This only copies the current information of the argument e to -// *this ExprRep. -void ExprRep::reduceTo(const ExprRep *e) { - if (e->appComputed()) { - appValue() = e->appValue(); - appComputed() = true; - flagsComputed() = true; - knownPrecision() = e->knownPrecision(); -#ifdef CORE_DEBUG - relPrecision() = e->relPrecision(); - absPrecision() = e->absPrecision(); - numNodes() = e->numNodes(); -#endif - - } - d_e() = e->d_e(); - //visited() = e->visited(); - sign() = e->sign(); - uMSB() = e->uMSB(); - lMSB() = e->lMSB(); - // length() = e->length(); // fixed? original = 1 - measure() = e->measure(); - - // BFMSS[2,5] bound. - u25() = e->u25(); - l25() = e->l25(); - v2p() = e->v2p(); - v2m() = e->v2m(); - v5p() = e->v5p(); - v5m() = e->v5m(); - - high() = e->high(); - low() = e->low(); // fixed? original = 0 - lc() = e->lc(); - tc() = e->tc(); - - // Chee (Mar 23, 2004), Notes on ratFlag(): - // =============================================================== - // For more information on the use of this flag, see progs/pentagon. - // This is an integer valued member of the NodeInfo class. - // Its value is used to determine whether - // we can ``reduce'' an Expression to a single node containing - // a BigRat value. This reduction is done if the global variable - // rationalReduceFlag=true. The default value is false. - // This is the intepretation of ratFlag: - // ratFlag < 0 means irrational - // ratFlag = 0 means not initialized - // ratFlag > 0 means rational - // Currently, ratFlag>0 is an upper bound on the size of the expression, - // since we recursively compute - // ratFlag(v) = ratFlag(v.lchild)+ratFlag(v.rchild) + 1. - // PROPOSAL: if ratFlag() > RAT_REDUCE_THRESHHOLD - // then we automatically do a reduction. We must determine - // an empirical value for RAT_REDUCE_THRESHOLD - - if (rationalReduceFlag) { - ratFlag() = e->ratFlag(); - - if (e->ratFlag() > 0 && e->ratValue() != NULL) { - ratFlag() ++; - if (ratValue() == NULL) - ratValue() = new BigRat(*(e->ratValue())); - else - *(ratValue()) = *(e->ratValue()); - } else - ratFlag() = -1; - } -} - -void ExprRep::reduceToZero() { - appValue() = CORE_REAL_ZERO; - appComputed() = true; - flagsComputed() = true; - knownPrecision() = CORE_negInfty; -#ifdef CORE_DEBUG - relPrecision() = EXTLONG_ZERO; - absPrecision() = CORE_negInfty; - // numNodes() = 0; -#endif - - d_e() = EXTLONG_ONE; - visited() = false; - sign() = 0; - uMSB() = CORE_negInfty; - lMSB() = CORE_negInfty; - // length() = 0; // fixed? original = 1 - measure() = EXTLONG_ZERO; - - // BFMSS[2,5] bound. - u25() = l25() = v2p() = v2m() = v5p() = v5m() = EXTLONG_ZERO; - - low() = EXTLONG_ONE; // fixed? original = 0 - high() = lc() = tc() = EXTLONG_ZERO; - - if (rationalReduceFlag) { - if (ratFlag() > 0) { - ratFlag() ++; - if (ratValue() == NULL) - ratValue() = new BigRat(0); - else - *(ratValue()) = 0; - } else - ratFlag() = 1; - } -} - -//////////////////////////////////////////////////////////// -// Test whether the current approximate value satisfies -// the composite precision requirements [relPrec, absPrec]. -//////////////////////////////////////////////////////////// - -bool ExprRep::withinKnownPrecision(const extLong& relPrec, - const extLong& absPrec) { - if (appComputed()) { // an approximate value has been evaluated. - if (appValue().isExact()) { - return true; - } else { // approximation has certain error. - // decide to which position it is required to compute correctly. - extLong required = core_max(-absPrec, appValue().lMSB()-relPrec); - // see whether the existing error is smaller than the requirement. - return (knownPrecision() <= required); - } - } else - return false; -}//withinKnownPrecision(a, r) - -// approximate the expression to certain precisions when -// necessary (either no approximate value available or -// the existing one is not accurate enough). -void ExprRep::approx(const extLong& relPrec = defRelPrec, - const extLong& absPrec = defAbsPrec) { - if (!getSign()) - return; // if it is exactly zero... - - // NOTE: The Filter might give a precise enough approximation already. - if (!getExactSign()) - return; - - if (!appComputed() || (!withinKnownPrecision(relPrec, absPrec))) { - // it's necessary to re-evaluate. - // to avoid huge lMSB which would cause long time and problems. - - // if it is a rational node - if (rationalReduceFlag && ratFlag() > 0 && ratValue() != NULL) - appValue() = Real(*(ratValue())).approx(relPrec, absPrec); //< shouldn't - // this case be done by computeApproxValue()? - else - computeApproxValue(relPrec, absPrec); - - // update flags - appComputed() = true; - knownPrecision() = appValue().clLgErr(); -#ifdef CORE_DEBUG - if (relPrecision() < relPrec) - relPrecision() = relPrec; - if (absPrecision() < absPrec) - absPrecision() = absPrec; -#endif - - } -} - -// return an approximate value to certain precision. -const Real& ExprRep::getAppValue(const extLong& relPrec, - const extLong& absPrec) { - if (getSign()) { - approx(relPrec, absPrec); - return appValue(); - } else - return CORE_REAL_ZERO; -} - -std::ostream& operator<<(std::ostream& o, ExprRep& rep) { - if (rep.getSign()) { - rep.approx(defRelPrec, defAbsPrec); - o << rep.appValue(); - } else { - o << "0"; - } - return o; -} - -// Chee, Zilin: July 17, 2002 -// Original algorithm is wrongly implemented, and can take time -// exponential in the size of the dag. -// -// METHOD: -// Inductively assume that all "visited" flags are false. -// This calls for a reimplementation of "count()" and "clearFlag()". -// Actually, we did not have to fix the count() function. -// -// (1) First recursively compute d_e for each node -// by calling the count() function. -// Important thing is count() will turn the "visited" flags -// to be true, so that there is no double counting. -// Furthermore, if d_e had already been computed, the -// arithmetic for d_e can be avoided (in this case, -// it is only the setting of "visited" flag that we -// are interested in! -// (2) At the end of count(), we have set all reachable nodes -// to "visited", and their d_e have been computed. -// (3) Now, call clearFlag() to recursively clear all reachable -// nodes. NOTE THAT PREVIOUSLY, clearFlag() was called -// first! This obvious is wrong - -extLong ExprRep::degreeBound() { - if (d_e() == EXTLONG_ONE) // no radical nodes below - return EXTLONG_ONE; - count(); - clearFlag(); - return d_e(); -} -// class ConstRealRep -// constructor -ConstRealRep::ConstRealRep(const Real & r) : value(r) { - if (!value.isExact()) { - // clone the BigFloat and set its error to zero. - value = value.BigFloatValue().makeExact(); - } - ffVal = filteredFp(value); -} - -// initialize nodeInfo -void ConstRep::initNodeInfo() { - nodeInfo = new NodeInfo(); - d_e() = EXTLONG_ONE; -} -void UnaryOpRep::initNodeInfo() { - if (child->nodeInfo == NULL) - child->initNodeInfo(); - nodeInfo = new NodeInfo(); -} -void BinOpRep::initNodeInfo() { - if (first->nodeInfo == NULL) - first->initNodeInfo(); - if (second->nodeInfo == NULL) - second->initNodeInfo(); - nodeInfo = new NodeInfo(); -} - -#ifdef CORE_DEBUG -unsigned long ConstRep::dagSize() { - if (!visited()) { - visited() = true; - numNodes() = 1; - } else - numNodes() = 0; - return numNodes(); -} - -unsigned long UnaryOpRep::dagSize() { - if (!visited()) { - visited() = true; - numNodes() = child->dagSize() + 1; - } else - numNodes() = 0; - return numNodes(); -} - -unsigned long BinOpRep::dagSize() { - if (!visited()) { - visited() = true; - numNodes() = first->dagSize() + second->dagSize() + 1; - } else - numNodes() = 0; - return numNodes(); -} - -void ConstRep::fullClearFlag() { - if (visited()) - visited() = false; -} - -void UnaryOpRep::fullClearFlag() { - if (visited()) { - child->fullClearFlag(); - visited() = false; - } -} - -void BinOpRep::fullClearFlag() { - if (visited()) { - first->fullClearFlag(); - second->fullClearFlag(); - visited() = false; - } -} -#endif - -// -// clear visited flag -// -/* see Expr.h - void ConstRep::clearFlag() - { visited = false; } -*/ -void UnaryOpRep::clearFlag() { - if (d_e() == EXTLONG_ONE) - return; // no radicals below. - if (visited()) { - visited() = false; - child->clearFlag(); - } -} -// class BinOpRep -void BinOpRep::clearFlag() { - if (d_e() == EXTLONG_ONE) - return; // rational below - if (visited()) { - visited() = false; - first->clearFlag(); - second->clearFlag(); - } -} - -// -// count # of squareroot -// -extLong ConstRep::count() { - if (visited()) - return EXTLONG_ONE; - visited() = true; - return d_e(); -} - -extLong NegRep::count() { - if (d_e() == EXTLONG_ONE) - return EXTLONG_ONE; - if (visited()) - return EXTLONG_ONE; - visited() = true; - d_e() = child->count(); - return d_e(); -} - -extLong SqrtRep::count() { - if (d_e() == EXTLONG_ONE) - return EXTLONG_ONE; - if (visited()) - return EXTLONG_ONE; - visited() = true; - d_e() = child->count() * EXTLONG_TWO; - return d_e(); -} - -extLong BinOpRep::count() { - if (d_e() == EXTLONG_ONE) - return EXTLONG_ONE; - if (visited()) - return EXTLONG_ONE; - visited() = true; - d_e() = first->count() * second->count(); - return d_e(); -} - -// -// compute exact flags functions -// -// exact value - -void computeExactFlags_temp(ConstRep* t, const Real &value) { - // Chen Li: the following is incorrect: - // uMSB = lMSB = value.MSB(); - // because the value could be a BigFloat which is an interval. - if (value.isExact()) { - t->uMSB() = t->lMSB() = value.MSB(); - } else { - t->uMSB() = value.uMSB(); - t->lMSB() = value.lMSB(); - core_error("Leaves in DAG is not exact!", __FILE__, __LINE__, true); - } - - t->sign() = value.sign(); - // t->length() = value.length(); - t->measure() = value.height(); // for rationals and integers, - // measure = height. - - // BFMSS[2,5] bound. - value.ULV_E(t->u25(), t->l25(), t->v2p(), t->v2m(), t->v5p(), t->v5m()); - - // The original BFMSS parameters can be set from the BFMSS[2,5] parameters. - // Here we just need them locally. - extLong u_e = t->u25() + t->v2p() + ceilLg5(t->v5p()); - extLong l_e = t->l25() + t->v2m() + ceilLg5(t->v5m()); - -#ifdef ORIGINAL_BFMSS - // To go back to the original BFMSS : - t->u25() = u_e; - t->l25() = l_e; - t->v2p() = t->v2m() = t->v5p() = t->v5m() = EXTLONG_ZERO; -#elif defined BFMSS_2_ONLY - // To go back to BFMSS[2] only : - t->u25() = t->u25() + ceilLg5(t->v5p()); - t->l25() = t->l25() + ceilLg5(t->v5m()); - t->v5p() = t->v5m() = EXTLONG_ZERO; -#endif - - if (l_e == EXTLONG_ZERO) { // no divisions introduced - t->high() = u_e; - t->low() = EXTLONG_ONE - u_e; // - (u_e - 1) - } else { - t->high() = u_e - l_e + EXTLONG_ONE; - t->low() = EXTLONG_TWO - t->high(); - } - - t->lc() = l_e; - t->tc() = u_e; - - // set BigRat value - if (rationalReduceFlag) { - t->ratFlag() = 1; - t->ratValue() = new BigRat(value.BigRatValue()); - } - - t->flagsComputed() = true; -} - -void ConstDoubleRep::computeExactFlags() {// can be made more efficient - computeExactFlags_temp(this, Real(ffVal.getValue())); -} - -void ConstRealRep::computeExactFlags() { - computeExactFlags_temp(this, value); -} - -void NegRep::computeExactFlags() { - if (!child->flagsComputed()) - child->computeExactFlags(); - - if (child->sign() == 0) { - reduceToZero(); - return; - } - - if (rationalReduceFlag) { - if (child->ratFlag()>0 && child->ratValue() != NULL) { - BigRat val = -(*(child->ratValue())); - reduceToBigRat(val); - ratFlag() = child->ratFlag()+1; - return; - } else - ratFlag() = -1; - } - - sign() = -child->sign(); - uMSB() = child->uMSB(); - lMSB() = child->lMSB(); - - // length() = child->length(); - measure() = child->measure(); - u25() = child->u25(); - l25() = child->l25(); - v2p() = child->v2p(); - v2m() = child->v2m(); - v5p() = child->v5p(); - v5m() = child->v5m(); - high() = child->high(); - low() = child->low(); - lc() = child->lc(); - tc() = child->tc(); - flagsComputed() = true; -}//NegRep::computeExactFlags - -void SqrtRep::computeExactFlags() { - if (!child->flagsComputed()) - child->computeExactFlags(); - - if (rationalReduceFlag) - ratFlag() = -1; - - sign() = child->sign(); - if (sign() < 0) - core_error("squareroot is called with negative operand.", - __FILE__, __LINE__, true); - - uMSB() = child->uMSB() / EXTLONG_TWO; - lMSB() = child->lMSB() / EXTLONG_TWO; - - // length() = child->length(); - measure() = child->measure(); - - // BFMSS[2,5] bound. - if (child->v2p() + ceilLg5(child->v5p()) + child->u25() >= - child->v2m() + ceilLg5(child->v5m()) + child->l25()) { - extLong vtilda2 = child->v2p() + child->v2m(); - v2p() = vtilda2 / EXTLONG_TWO; - v2m() = child->v2m(); - extLong vmod2; - if (v2p().isInfty()) - vmod2 = CORE_INFTY; - else - vmod2 = vtilda2 - EXTLONG_TWO*v2p(); // == vtilda2 % 2 - extLong vtilda5 = child->v5p() + child->v5m(); - v5p() = vtilda5 / EXTLONG_TWO; - v5m() = child->v5m(); - extLong vmod5; - if (v5p().isInfty()) - vmod5 = CORE_INFTY; - else - vmod5 = vtilda5 - EXTLONG_TWO*v5p(); // == vtilda5 % 2 - u25() = (child->u25() + child->l25() + vmod2 + ceilLg5(vmod5) + EXTLONG_ONE) / EXTLONG_TWO; - l25() = child->l25(); - } else { - extLong vtilda2 = child->v2p() + child->v2m(); - v2p() = child->v2p(); - v2m() = vtilda2 / EXTLONG_TWO; - extLong vmod2; - if (v2m().isInfty()) - vmod2 = CORE_INFTY; - else - vmod2 = vtilda2 - EXTLONG_TWO*v2m(); // == vtilda2 % 2 - extLong vtilda5 = child->v5p() + child->v5m(); - v5p() = child->v5p(); - v5m() = vtilda5 / EXTLONG_TWO; - u25() = child->u25(); - extLong vmod5; - if (v5m().isInfty()) - vmod5 = CORE_INFTY; - else - vmod5 = vtilda5 - EXTLONG_TWO*v5m(); // == vtilda5 % 2 - l25() = (child->u25() + child->l25() + vmod2 + ceilLg5(vmod5) + EXTLONG_ONE) / EXTLONG_TWO; - } - - high() = (child->high() +EXTLONG_ONE)/EXTLONG_TWO; - low() = child->low() / EXTLONG_TWO; - lc() = child->lc(); - tc() = child->tc(); - flagsComputed() = true; -}// SqrtRep::computeExactFlags - -void MultRep::computeExactFlags() { - if (!first->flagsComputed()) - first->computeExactFlags(); - if (!second->flagsComputed()) - second->computeExactFlags(); - - if ((!first->sign()) || (!second->sign())) { - // value must be exactly zero. - reduceToZero(); - return; - } - // rational node - if (rationalReduceFlag) { - if (first->ratFlag() > 0 && second->ratFlag() > 0) { - BigRat val = (*(first->ratValue()))*(*(second->ratValue())); - reduceToBigRat(val); - ratFlag() = first->ratFlag() + second->ratFlag(); - return; - } else - ratFlag() = -1; - } - - // value is irrational. - uMSB() = first->uMSB() + second->uMSB() + EXTLONG_ONE; - lMSB() = first->lMSB() + second->lMSB(); - sign() = first->sign() * second->sign(); - - extLong df = first->d_e(); - extLong ds = second->d_e(); - // extLong lf = first->length(); - // extLong ls = second->length(); - - // length() = df * ls + ds * lf; - measure() = (first->measure()) * ds+(second->measure()) * df; - - // BFMSS[2,5] bound. - v2p() = first->v2p() + second->v2p(); - v2m() = first->v2m() + second->v2m(); - v5p() = first->v5p() + second->v5p(); - v5m() = first->v5m() + second->v5m(); - u25() = first->u25() + second->u25(); - l25() = first->l25() + second->l25(); - - high() = first->high() + second->high(); - low() = first->low() + second->low(); - - lc() = ds * first->lc() + df * second->lc(); - tc() = core_min(ds * first->tc() + df * second->tc(), measure()); - - flagsComputed() = true; -}// MultRep::computeExactFlags - -void DivRep::computeExactFlags() { - if (!first->flagsComputed()) - first->computeExactFlags(); - if (!second->flagsComputed()) - second->computeExactFlags(); - - if (!second->sign()) - core_error("zero divisor.", __FILE__, __LINE__, true); - - if (!first->sign()) {// value must be exactly zero. - reduceToZero(); - return; - } - - // rational node - if (rationalReduceFlag) { - if (first->ratFlag() > 0 && second->ratFlag() > 0) { - BigRat val = (*(first->ratValue()))/(*(second->ratValue())); - reduceToBigRat(val); - ratFlag() = first->ratFlag() + second->ratFlag(); - return; - } else - ratFlag() = -1; - } - - // value is irrational. - uMSB() = first->uMSB() - second->lMSB(); - lMSB() = first->lMSB() - second->uMSB() - EXTLONG_ONE; - sign() = first->sign() * second->sign(); - - extLong df = first->d_e(); - extLong ds = second->d_e(); - // extLong lf = first->length(); - // extLong ls = second->length(); - - // length() = df * ls + ds * lf; - measure() = (first->measure())*ds + (second->measure())*df; - - // BFMSS[2,5] bound. - v2p() = first->v2p() + second->v2m(); - v2m() = first->v2m() + second->v2p(); - v5p() = first->v5p() + second->v5m(); - v5m() = first->v5m() + second->v5p(); - u25() = first->u25() + second->l25(); - l25() = first->l25() + second->u25(); - - high() = first->high() + second->low(); - low() = first->low() + second->high(); - - lc() = ds * first->lc() + df * second->tc(); - tc() = core_min(ds * first->tc() + df * second->lc(), measure()); - - flagsComputed() = true; -} - -// -// approximation functions -// -void ConstDoubleRep::computeApproxValue(const extLong& /*relPrec*/, - const extLong& /*absPrec*/) -// can ignore precision bounds since ffVal.getValue() returns exact value -{ - appValue() = Real(ffVal.getValue()); -} - -void ConstRealRep::computeApproxValue(const extLong& relPrec, - const extLong& absPrec) { - appValue() = value.approx(relPrec, absPrec); -} - -void NegRep::computeApproxValue(const extLong& relPrec, - const extLong& absPrec) { - appValue() = -child->getAppValue(relPrec, absPrec); -} - -void SqrtRep::computeApproxValue(const extLong& relPrec, - const extLong& absPrec) { - extLong r = relPrec + relPrec + EXTLONG_EIGHT; // chenli: ??? - extLong a = absPrec + absPrec + EXTLONG_EIGHT; - extLong pr = - lMSB() + r; - extLong p = pr < a ? pr : a; - - Real val = child->getAppValue(r, a); - if (incrementalEvalFlag) { - if (appValue() == CORE_REAL_ZERO) - appValue() = val; - appValue() = val.sqrt(p, appValue().BigFloatValue()); - } else - appValue() = val.sqrt(p); -} - -void MultRep::computeApproxValue(const extLong& relPrec, - const extLong& absPrec) { - if (lMSB() < EXTLONG_BIG && lMSB() > EXTLONG_SMALL) { - extLong r = relPrec + EXTLONG_FOUR; - extLong afr = - first->lMSB() + EXTLONG_ONE; - extLong afa = second->uMSB() + absPrec + EXTLONG_FIVE; - extLong af = afr > afa ? afr : afa; - extLong asr = - second->lMSB() + EXTLONG_ONE; - extLong asa = first->uMSB() + absPrec + EXTLONG_FIVE; - extLong as = asr > asa ? asr : asa; - appValue() = first->getAppValue(r, af)*second->getAppValue(r, as); - } else { - std::cerr << "lMSB = " << lMSB() << std::endl; - core_error("a huge lMSB in MulRep", __FILE__, __LINE__, false); - } -} - -void DivRep::computeApproxValue(const extLong& relPrec, - const extLong& absPrec) { - if (lMSB() < EXTLONG_BIG && lMSB() > EXTLONG_SMALL) { - extLong rr = relPrec + EXTLONG_SEVEN; // These rules come from - extLong ra = uMSB() + absPrec + EXTLONG_EIGHT; // Koji's Master Thesis, page 65 - extLong ra2 = core_max(ra, EXTLONG_TWO); - extLong r = core_min(rr, ra2); - extLong af = - first->lMSB() + r; - extLong as = - second->lMSB() + r; - - extLong pr = relPrec + EXTLONG_SIX; - extLong pa = uMSB() + absPrec + EXTLONG_SEVEN; - extLong p = core_min(pr, pa); // Seems to be an error: - // p can be negative here! - // Also, this does not conform to - // Koji's thesis which has a default - // relative precision (p.65). - - appValue() = first->getAppValue(r, af).div(second->getAppValue(r, as), p); - } else { - std::cerr << "lMSB = " << lMSB() << std::endl; - core_error("a huge lMSB in DivRep", __FILE__, __LINE__, false); - } -} - -// -// Debug Help Functions -// -void Expr::debug(int mode, int level, int depthLimit) const { - std::cout << "-------- Expr debug() -----------" << std::endl; - std::cout << "rep = " << rep << std::endl; - if (mode == Expr::LIST_MODE) - rep->debugList(level, depthLimit); - else if (mode == Expr::TREE_MODE) - rep->debugTree(level, 0, depthLimit); - else - core_error("unknown debugging mode", __FILE__, __LINE__, false); - std::cout << "---- End Expr debug(): " << std::endl; -} - - -const std::string ExprRep::dump(int level) const { - std::ostringstream ost; - if (level == OPERATOR_ONLY) { - ost << op(); - } else if (level == VALUE_ONLY) { - ost << appValue(); - } else if (level == OPERATOR_VALUE) { - ost << op() << "[val: " << appValue() << "]"; - } else if (level == FULL_DUMP) { - ost << op() - << "[val: " << appValue() << "; " - << "kp: " << knownPrecision() << "; " -#ifdef CORE_DEBUG - << "r: " << relPrecision() << "; " - << "a: " << absPrecision() << "; " -#endif - << "lMSB: " << lMSB() << "; " - << "uMSB: " << uMSB() << "; " - << "sign: " << sign() << "; " - // << "length: " << length() << "; " - << "measure: " << measure() << "; " - << "d_e: " << d_e() << "; " - << "u25: " << u25() << "; " - << "l25: " << l25() << "; " - << "v2p: " << v2p() << "; " - << "v2m: " << v2m() << "; " - << "v5p: " << v5p() << "; " - << "v5m: " << v5m() << "; " - << "high: " << high() << "; " - << "low: " << low() << "; " - << "lc: " << lc() << "; " - << "tc: " << tc() - << "]"; - } - return std::string(ost.str()); - // note that str() return an array not properly terminated! -} - - -void UnaryOpRep::debugList(int level, int depthLimit) const { - if (depthLimit <= 0) - return; - if (level == Expr::SIMPLE_LEVEL) { - std::cout << "(" << dump(OPERATOR_VALUE); - child->debugList(level, depthLimit - 1); - std::cout << ")"; - } else if (level == Expr::DETAIL_LEVEL) { - std::cout << "(" << dump(FULL_DUMP); - child->debugList(level, depthLimit - 1); - std::cout << ")"; - } -} - -void UnaryOpRep::debugTree(int level, int indent, int depthLimit) const { - if (depthLimit <= 0) - return; - for (int i = 0; idebugTree(level, indent + 2, depthLimit - 1); -} - -void ConstRep::debugList(int level, int depthLimit) const { - if (depthLimit <= 0) - return; - if (level == Expr::SIMPLE_LEVEL) { - std::cout << "(" << dump(OPERATOR_VALUE) << ")"; - } else if (level == Expr::DETAIL_LEVEL) { - std::cout << "(" << dump(FULL_DUMP) << ")"; - } -} - -void ConstRep::debugTree(int level, int indent, int depthLimit) const { - if (depthLimit <= 0) - return; - for (int i=0; idebugList(level, depthLimit - 1); - std::cout << ", "; - second->debugList(level, depthLimit - 1); - std::cout << ")" ; -} - -void BinOpRep::debugTree(int level, int indent, int depthLimit) const { - if (depthLimit <= 0) - return; - for (int i=0; idebugTree(level, indent + 2, depthLimit - 1); - second->debugTree(level, indent + 2, depthLimit - 1); -} - -} //namespace CORE +#endif \ No newline at end of file diff --git a/CGAL_Core/src/CGAL_Core/GmpIO.cpp b/CGAL_Core/src/CGAL_Core/GmpIO.cpp index 3d665366a08..9b74ba289e7 100644 --- a/CGAL_Core/src/CGAL_Core/GmpIO.cpp +++ b/CGAL_Core/src/CGAL_Core/GmpIO.cpp @@ -12,257 +12,9 @@ * $Id$ ***************************************************************************/ -/* Auxiliary functions for C++-style input of GMP types. - -Copyright 2001 Free Software Foundation, Inc. - -This file is part of the GNU MP Library. - -The GNU MP Library is free software; you can redistribute it and/or modify -it under the terms of the GNU Lesser General Public License as published by -the Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -The GNU MP Library is distributed in the hope that it will be useful, but -WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public -License for more details. - -You should have received a copy of the GNU Lesser General Public License -along with the GNU MP Library; see the file COPYING.LIB. If not, write to -the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, -MA 02110-1301, USA. */ +#ifndef CGAL_HEADER_ONLY #include -#include -#include -#include -#include +#include -using namespace std; - -namespace CORE { - -int -__gmp_istream_set_base (istream &i, char &c, bool &zero, bool &showbase) -{ - int base; - - zero = showbase = false; - switch (i.flags() & ios::basefield) - { - case ios::dec: - base = 10; - break; - case ios::hex: - base = 16; - break; - case ios::oct: - base = 8; - break; - default: - showbase = true; // look for initial "0" or "0x" or "0X" - if (c == '0') - { - if (! i.get(c)) - c = 0; // reset or we might loop indefinitely - - if (c == 'x' || c == 'X') - { - base = 16; - i.get(c); - } - else - { - base = 8; - zero = true; // if no other digit is read, the "0" counts - } - } - else - base = 10; - break; - } - - return base; -} - -void -__gmp_istream_set_digits (string &s, istream &i, char &c, bool &ok, int base) -{ - switch (base) - { - case 10: - while (isdigit(c)) - { - ok = true; // at least a valid digit was read - s += c; - if (! i.get(c)) - break; - } - break; - case 8: - while (isdigit(c) && c != '8' && c != '9') - { - ok = true; // at least a valid digit was read - s += c; - if (! i.get(c)) - break; - } - break; - case 16: - while (isxdigit(c)) - { - ok = true; // at least a valid digit was read - s += c; - if (! i.get(c)) - break; - } - break; - } -} - -istream & -//operator>> (istream &i, mpz_ptr z) -io_read (istream &i, mpz_ptr z) -{ - int base; - char c = 0; - string s; - bool ok = false, zero, showbase; - - i.get(c); // start reading - - if (i.flags() & ios::skipws) // skip initial whitespace - while (isspace(c) && i.get(c)) - ; - - if (c == '-' || c == '+') // sign - { - if (c == '-') // mpz_set_str doesn't accept '+' - s = "-"; - i.get(c); - } - - while (isspace(c) && i.get(c)) // skip whitespace - ; - - base = __gmp_istream_set_base(i, c, zero, showbase); // select the base - __gmp_istream_set_digits(s, i, c, ok, base); // read the number - - if (i.good()) // last character read was non-numeric - i.putback(c); - else if (i.eof() && (ok || zero)) // stopped just before eof - i.clear(); - - if (ok) - mpz_set_str(z, s.c_str(), base); // extract the number - else if (zero) - mpz_set_ui(z, 0); - else - i.setstate(ios::failbit); // read failed - - return i; -} - -istream & -//operator>> (istream &i, mpq_ptr q) -io_read (istream &i, mpq_ptr q) -{ - int base; - char c = 0; - string s; - bool ok = false, zero, showbase; - - i.get(c); // start reading - - if (i.flags() & ios::skipws) // skip initial whitespace - while (isspace(c) && i.get(c)) - ; - - if (c == '-' || c == '+') // sign - { - if (c == '-') - s = "-"; - i.get(c); - } - - while (isspace(c) && i.get(c)) // skip whitespace - ; - - base = __gmp_istream_set_base(i, c, zero, showbase); // select the base - __gmp_istream_set_digits(s, i, c, ok, base); // read the numerator - - if (! ok && zero) // the only digit read was "0" - { - base = 10; - s += '0'; - ok = true; - } - - if (i.flags() & ios::skipws) - while (isspace(c) && i.get(c)) // skip whitespace - ; - - if (c == '/') // there's a denominator - { - bool zero2 = false; - int base2 = base; - - s += '/'; - ok = false; // denominator is mandatory - i.get(c); - - while (isspace(c) && i.get(c)) // skip whitespace - ; - - if (showbase) // check base of denominator - base2 = __gmp_istream_set_base(i, c, zero2, showbase); - - if (base2 == base || base2 == 10) // read the denominator - __gmp_istream_set_digits(s, i, c, ok, base); - - if (! ok && zero2) // the only digit read was "0" - { // denominator is 0, but that's your business - s += '0'; - ok = true; - } - } - - if (i.good()) // last character read was non-numeric - i.putback(c); - else if (i.eof() && ok) // stopped just before eof - i.clear(); - - if (ok) - mpq_set_str(q, s.c_str(), base); // extract the number - else - i.setstate(ios::failbit); // read failed - - return i; -} - -ostream& -//operator<< (ostream &o, mpz_srcptr z) -io_write (ostream &o, mpz_srcptr z) -{ - char *str = new char [mpz_sizeinbase(z,10) + 2]; - str = mpz_get_str(str, 10, z); - o << str ; - delete[] str; - return o; -} - -ostream& -//operator<< (ostream &o, mpq_srcptr q) -io_write (ostream &o, mpq_srcptr q) -{ - // size according to GMP documentation - char *str = new char [mpz_sizeinbase(mpq_numref(q), 10) + - mpz_sizeinbase (mpq_denref(q), 10) + 3]; - str = mpq_get_str(str, 10, q); - o << str ; - delete[] str; - return o; -} - -} //namespace CORE +#endif diff --git a/CGAL_Core/src/CGAL_Core/Real.cpp b/CGAL_Core/src/CGAL_Core/Real.cpp index 477c02a0cab..cd4604d14c0 100644 --- a/CGAL_Core/src/CGAL_Core/Real.cpp +++ b/CGAL_Core/src/CGAL_Core/Real.cpp @@ -35,242 +35,9 @@ * $Id$ ***************************************************************************/ -#include +#ifndef CGAL_HEADER_ONLY + #include +#include -namespace CORE { - -const Real& Real::getZero() { - static Real Zero(0); - return Zero; -} - -BigInt floor(const Real& r, Real &sub) { - BigInt f = r.approx(CORE_INFTY, 2).BigIntValue(); - sub = r-f; - // Adjustment - if (sub<0) - ++sub, --f; - if (sub>=1) - --sub, ++f; - assert(sub >=0 && sub<1); - return f; -} - -// pow(r,n) and power(r, n) are the same function -// -Real pow(const Real& r, unsigned long n) { - if (n == 0) - return Real(1); - else if (n == 1) - return r; - else { - Real x = r; - while ((n % 2) == 0) { // n is even - x *= x; - n >>= 1; - } - Real u = x; - while (true) { - n >>= 1; - if (n == 0) - return u; - x *= x; - if ((n % 2) == 1) // n is odd - u *= x; - } - //return u; // unreachable - } -}//pow - -extern BigInt FiveTo(unsigned long exp); - -// Construct Real from String -// Note: -// -- Zilin Du: 06/03/2003 -// -- Original it is the code for Real's constructor for "const char*". -// I change it to a function so that two constrcutors can share the code. -// now it is private and no default value. -// -// --Default value of the argument "prec" is defInputDigits -// --If prec is CORE_posInfty, then the input is -// read in exactly. Otherwise, we convert to a RealBigFloat -// with absolute error at most 10^{-prec} - -// Constructor Real( char *str, extLong& prec) -// is very similar to -// BigFloatRep::fromString( char *str, extLong& prec); -// Differences: -// In BigFloat(str, prec), the value of prec cannot be infinity, and -// it defaults to defBigFloatInputDigits; -// In Real(str, prec), the value of prec is allowed to be infinity, and -// it defaults to defInputDigits. -// -// Why do we have the two versions? It allows us to use the BigFloat class -// directly, without relying on Real class. - -void Real::constructFromString(const char *str, const extLong& prec ) -// NOTE: prec defaults to defInputDigits (see Real.h) -{ - // 8/8/01, Chee and Zilin: add a new rational string format: - // this format is indicated by the presence of a slash "/" - // Moreover, the value of prec is ignored (basically - // assumed to be infinity). - - if (std::strchr(str, '/') != NULL) { // this is a rational number - rep = new RealBigRat(BigRat(str)); - return; - } - - const char *e = std::strchr(str, 'e'); - int dot = 0; - long e10 = 0; - if (e != NULL) - e10 = std::atol(e+1); // e10 is decimal precision of the input string - // i.e., input is A/10^{e10}. - else { - e = str + std::strlen(str); -#ifdef CORE_DEBUG - assert(*e == '\0'); #endif - } - - const char *p = str; - if (*p == '-' || *p == '+') - p++; - BigInt m(0); - - for (; p < e; p++) { - if (*p == '.') { - dot = 1; - continue; - } - m = m * 10 + (*p - '0'); - if (dot) - e10--; - } - - long t = (e10 < 0) ? -e10 : e10; - BigInt one(1); - BigInt ten = FiveTo(t) * (one << static_cast(t)); - if (*str == '-') - m = -m; - if (e10 >= 0) { - // convert exactly from integer numbers - m *= ten; - rep = new RealBigInt(m); - } else { // e10 < 0, fractional numbers - // HERE IS WHERE WE USE THE SYSTEM CONSTANT - // defInputDigits - // Note: defInputDigits should be at least log_2(10). - // We default defInputDigits to 4. - //std::cout << "(m,ten)=" << m << "," << ten << std::endl; - BigRat r(m, ten); - if (prec.isInfty()) { // convert exactly! to a big rational - rep = new RealBigRat(r); - } else { - // convert approximately, to a BigFloat within the - // specified precision: - // BigFloat bf(r, CORE_posInfty, prec * lgTenM) ; - BigFloat bf(r, CORE_posInfty, prec * 4) ; - rep = new RealBigFloat(bf); - } - } -}// Real(str, prec) - -// The operator >>(i,x) calls the constructor Real(char*) -std::istream& operator >>(std::istream& i, Real& x) { - int size = 20; - char *str = new char[size]; - char *p = str; - char c; - int d = 0, e = 0, s = 0; - // int done = 0; - - // Chen Li: fixed a bug, the original statement is - // for (i.get(c); c == ' '; i.get(c)); - // use isspace instead of testing c == ' ', since it must also - // skip tab, catridge/return, etc. - // Change to: - // int status; - do { - i.get(c); - } while (!i.eof() && isspace(c)); /* loop if met end-of-file, or - char read in is white-space. */ - // Chen Li, - // original "if (c == EOF) ..." is unsafe since c is of char type and - // EOF is of int tyep with a negative value -1 - - if (i.eof()) { - i.clear(std::ios::eofbit | std::ios::failbit); - delete [] str; - return i; - } - - // the current content in "c" should be the first non-whitespace char - if (c == '-' || c == '+') { - *p++ = c; - i.get(c); - } - - for (; isdigit(c) || (!d && c=='.') || - (!e && c=='e') || (!s && (c=='-' || c=='+')); i.get(c)) { - if (!i) break; - if (!e && (c == '-' || c == '+')) - break; - // Chen Li: put one more rule to prohibite input like - // xxxx.xxxe+xxx.xxx: - if (e && (c == '.')) - break; - if (p - str == size) { - char *t = str; - str = new char[size*2]; - std::memcpy(str, t, size); - delete [] t; - p = str + size; - size *= 2; - } -#ifdef CORE_DEBUG - assert((p-str) < size); -#endif - - *p++ = c; - if (c == '.') - d = 1; - // Chen Li: fix a bug -- the sign of exponent can not happen before - // the character "e" appears! It must follow the "e' actually. - // if (e || c == '-' || c == '+') s = 1; - if (e) - s = 1; - if (c == 'e') - e = 1; - } - - if (!i && !i.eof()) { - delete [] str; - return i; - } - // chenli: make sure that the p is still in the range - if (p - str >= size) { - int len = p - str; - char *t = str; - str = new char[len + 1]; - std::memcpy(str, t, len); - delete [] t; - p = str + len; - } - -#ifdef CORE_DEBUG - assert(p - str < size); -#endif - - *p = '\0'; - i.putback(c); - i.clear(); - // old: x = Real(str, i.precision()); // use precision of input stream. - x = Real(str); // default precision = defInputDigits - delete [] str; - return i; -}//operator >> (std::istream&, Real&) - -} //namespace CORE diff --git a/CGAL_Core/src/CGAL_Core/extLong.cpp b/CGAL_Core/src/CGAL_Core/extLong.cpp index a4fdccbfca5..901b758de58 100644 --- a/CGAL_Core/src/CGAL_Core/extLong.cpp +++ b/CGAL_Core/src/CGAL_Core/extLong.cpp @@ -39,157 +39,9 @@ * $Id$ ***************************************************************************/ +#ifndef CGAL_HEADER_ONLY + #include +#include -namespace CORE { - -const extLong& extLong::getNaNLong() { - static extLong NaNLong(true); - return NaNLong; -} - -const extLong& extLong::getPosInfty() { - static extLong posInfty(EXTLONG_MAX); - return posInfty; -} - -const extLong& extLong::getNegInfty() { - static extLong negInfty(EXTLONG_MIN); - return negInfty; -} - -void extLong::add(extLong& z, long x, long y) { - if (x > 0 && y > 0 && x >= EXTLONG_MAX - y) { - z.val = EXTLONG_MAX; - z.flag = 1; - } else if (x < 0 && y < 0 && x <= EXTLONG_MIN - y) { - z.val = EXTLONG_MIN; - z.flag = -1; - } else { - z.val = x + y; - z.flag = 0; - } -} - -// arithmetic and assignment operators -extLong& extLong::operator+= (const extLong& y) { - if (flag == 2 || y.flag == 2 || (flag * y.flag < 0)) { -#ifdef CORE_DEBUG - if (flag * y.flag < 0) //want a message at the first creation of NaN - core_error("extLong NaN Error in addition.", __FILE__, __LINE__, false); #endif - - *this = CORE_NaNLong; - } else if (flag == 1 || y.flag == 1) { // one of them is +Inf - *this = CORE_posInfty; - } else if (flag == -1 || y.flag == -1) { // one of them is -Inf - *this = CORE_negInfty; - } else { // x and y are normal now - add(*this, val, y.val); - } - return *this; -} - -extLong& extLong::operator-= (const extLong& y) { - if (flag == 2 || y.flag == 2 || (flag * y.flag > 0)) { -#ifdef CORE_DEBUG - if (flag * y.flag > 0) //want a message at the first creation of NaN - core_error("extLong NaN Error in subtraction.", __FILE__, __LINE__, false); -#endif - - *this = CORE_NaNLong; - } else if (flag == 1 || y.flag == -1) { - *this = CORE_posInfty; - } else if (flag == -1 || y.flag == 1) { - *this = CORE_negInfty; - } else { - add(*this, val, -y.val); - } - return *this; -} - -extLong& extLong::operator*= (const extLong& y) { - if (flag == 2 || y.flag == 2) { - *this = CORE_NaNLong; - } else if ((flag != 0) || (y.flag != 0)) { - if (sign() * y.sign() > 0) - *this = CORE_posInfty; - else - *this = CORE_negInfty; - } else { // flag == 0 and y.flag == 0 - double d = double(val) * double(y.val); - long p = val * y.val; - if (std::fabs(d - p) <= std::fabs(d) * relEps) { - val = p; - flag = 0; - } else if (d > EXTLONG_MAX) { - *this = CORE_posInfty; - } else if (d < EXTLONG_MIN) { - *this = CORE_negInfty; - } else { -#ifdef CORE_DEBUG - core_error("extLong NaN Error in multiplication.",__FILE__,__LINE__,false); -#endif - *this = CORE_NaNLong; - } - } - return *this; -} - -extLong& extLong::operator/= (const extLong& y) { - if (flag==2 || y.flag==2 || ((flag != 0) && (y.flag != 0)) || (y.val == 0)) { -#ifdef CORE_DEBUG - if (y.val == 0) - core_error("extLong NaN Error, Divide by Zero.", __FILE__, __LINE__, false); - else if ((flag !=0) && (y.flag !=0)) - core_error("extLong NaN Error, +/-Inf/Inf.", __FILE__, __LINE__, false); -#endif - - *this = CORE_NaNLong; - } else if ((flag != 0) || (y.flag != 0)) { // y.flag == 0 now and y != 0 - if (sign() * y.sign() > 0) - *this = CORE_posInfty; - else - *this = CORE_negInfty; - } else { // flag == 0 and y.flag == 0 - val /= y.val; // no overflow in divisions - flag = 0; - } - return *this; -} - -// unary minus -extLong extLong::operator- () const { - if (flag == 0) - return extLong(-val); - else if (flag == 1) - return CORE_negInfty; - else if (flag == -1) - return CORE_posInfty; - else // NaN - return CORE_NaNLong; -} - -// sign -// You should check "flag" before calling this, otherwise -// you cannot interprete the returned value! -int extLong::sign() const { - if (flag == 2) - core_error("NaN Sign can not be determined!", __FILE__, __LINE__, false); - return ((val == 0) ? 0 : ((val > 0) ? 1 : -1)); -} - -// stream operators -std::ostream& operator<< (std::ostream& o, const extLong& x) { - if (x.flag == 1) - o << " infty "; - else if (x.flag == - 1) - o << " tiny "; - else if (x.flag == 2) - o << " NaN "; - else - o << x.val; - return o; -} - -} //namespace CORE diff --git a/Number_types/include/CGAL/CORE_BigFloat.h b/Number_types/include/CGAL/CORE_BigFloat.h index 48e2dffabd7..2488af7d00d 100644 --- a/Number_types/include/CGAL/CORE_BigFloat.h +++ b/Number_types/include/CGAL/CORE_BigFloat.h @@ -251,7 +251,7 @@ namespace internal{ CORE::BigFloat inline -round(const CORE::BigFloat& x, long rel_prec = CORE::defRelPrec.toLong() ){ +round(const CORE::BigFloat& x, long rel_prec = CORE::get_static_defRelPrec().toLong() ){ CGAL_postcondition(rel_prec >= 0); // since there is not rel prec defined if Zero_in(x) @@ -344,9 +344,9 @@ public: typedef long result_type; long operator() ( long prec ) const { - long result = ::CORE::defRelPrec.toLong(); - ::CORE::defRelPrec = prec; - ::CORE::defBFdivRelPrec = prec; + long result = ::CORE::get_static_defRelPrec().toLong(); + ::CORE::get_static_defRelPrec() = prec; + ::CORE::get_static_defBFdivRelPrec() = prec; return result; } }; @@ -356,7 +356,7 @@ public: typedef long result_type; long operator() () const { - return ::CORE::defRelPrec.toLong(); + return ::CORE::get_static_defRelPrec().toLong(); } }; }; @@ -378,17 +378,19 @@ template <> class Algebraic_structure_traits< CORE::BigFloat > : public std::unary_function< Type, Type > { public: Type operator()( const Type& x ) const { - // What I want is a sqrt computed with ::CORE::defRelPrec bits. + // What I want is a sqrt computed with + // ::CORE::get_static_defRelPrec() bits. // And not ::CORE::defBFsqrtAbsPrec as CORE does. - CGAL_precondition(::CORE::defRelPrec.toLong() > 0); + CGAL_precondition(::CORE::get_static_defRelPrec().toLong() > 0); CGAL_precondition(x > 0); - Type a = CGAL::internal::round(x, ::CORE::defRelPrec.toLong()*2); + Type a = CGAL::internal::round( + x, ::CORE::get_static_defRelPrec().toLong()*2); CGAL_postcondition(a > 0); - Type tmp1 = - CORE::BigFloat(a.m(),0,0).sqrt(::CORE::defRelPrec.toLong()); + Type tmp1 = CORE::BigFloat( + a.m(),0,0).sqrt(::CORE::get_static_defRelPrec().toLong()); Type err = Type(0,long(std::sqrt(double(a.err()))),0) * CORE::BigFloat::exp2(a.exp()*7); diff --git a/Number_types/include/CGAL/CORE_coercion_traits.h b/Number_types/include/CGAL/CORE_coercion_traits.h index d097ffcc6f1..d1ebbc69b72 100644 --- a/Number_types/include/CGAL/CORE_coercion_traits.h +++ b/Number_types/include/CGAL/CORE_coercion_traits.h @@ -84,7 +84,7 @@ struct Coercion_traits{ Type operator()(const CORE::BigFloat& x) const { return x;} Type operator()(const ::CORE::BigInt x) const { CORE::BigFloat result; - result.approx(x,CORE::defRelPrec.toLong(),LONG_MAX); + result.approx(x,CORE::get_static_defRelPrec().toLong(),LONG_MAX); // Do not use MakeFloorExact as it changes the Bigfloat CGAL_postcondition( ::CORE::BigRat(::CORE::BigFloat(result.m()-result.err(),0,result.exp())) <= x ); CGAL_postcondition( ::CORE::BigRat(::CORE::BigFloat(result.m()+result.err(),0,result.exp())) >= x ); @@ -104,7 +104,7 @@ struct Coercion_traits{ Type operator()(const CORE::BigFloat& x) const { return x;} Type operator()(const ::CORE::BigRat x) const { - CORE::BigFloat result(x,CORE::defRelPrec.toLong(),LONG_MAX); + CORE::BigFloat result(x,CORE::get_static_defRelPrec().toLong(),LONG_MAX); // Do not use MakeFloorExact as it changes the Bigfloat CGAL_postcondition( ::CORE::BigRat(::CORE::BigFloat(result.m()-result.err(),0,result.exp())) <= x ); CGAL_postcondition( ::CORE::BigRat(::CORE::BigFloat(result.m()+result.err(),0,result.exp())) >= x ); @@ -123,7 +123,7 @@ struct Coercion_traits{ typedef Type result_type; Type operator()(const CORE::BigFloat& x) const { return x;} Type operator()(const ::CORE::Expr x) const { - CORE::BigFloat result(x, CORE::defRelPrec.toLong(),LONG_MAX); + CORE::BigFloat result(x, CORE::get_static_defRelPrec().toLong(),LONG_MAX); // Do not use MakeFloorExact as it changes the Bigfloat CGAL_postcondition( ::CORE::BigRat(::CORE::BigFloat(result.m()-result.err(),0,result.exp())) <= x ); CGAL_postcondition( ::CORE::BigRat(::CORE::BigFloat(result.m()+result.err(),0,result.exp())) >= x ); diff --git a/STL_Extension/include/CGAL/Spatial_lock_grid_3.h b/STL_Extension/include/CGAL/Spatial_lock_grid_3.h index c2a389361b0..14a5e869790 100644 --- a/STL_Extension/include/CGAL/Spatial_lock_grid_3.h +++ b/STL_Extension/include/CGAL/Spatial_lock_grid_3.h @@ -606,64 +606,6 @@ protected: TLS_thread_uint_ids m_tls_thread_ids; }; -//***************************************************************************** -// class Spatial_lock_grid_3 -// Note: undocumented, for testing only... -//***************************************************************************** - -template <> -class Spatial_lock_grid_3 - : public Spatial_lock_grid_base_3< - Spatial_lock_grid_3 > -{ - typedef Spatial_lock_grid_base_3< - Spatial_lock_grid_3 > Base; - -public: - // Constructors - Spatial_lock_grid_3(const Bbox_3 &bbox, int num_grid_cells_per_axis) - : Base(bbox, num_grid_cells_per_axis) - { - int num_cells = - num_grid_cells_per_axis*num_grid_cells_per_axis*num_grid_cells_per_axis; - m_grid.resize(num_cells); - } - - /// Destructor - ~Spatial_lock_grid_3() - { - } - - bool is_cell_locked_impl(int cell_index) - { - bool locked = m_grid[cell_index].try_lock(); - if (locked) - m_grid[cell_index].unlock(); - return !locked; - } - - template - bool try_lock_cell_impl(int cell_index) - { - bool success = m_grid[cell_index].try_lock(); - if (success) - { - get_thread_local_grid()[cell_index] = true; - m_tls_locked_cells.local().push_back(cell_index); - } - return success; - } - - void unlock_cell_impl(int cell_index) - { - m_grid[cell_index].unlock(); - } - -protected: - - std::vector m_grid; -}; - } //namespace CGAL #else // !CGAL_LINKED_WITH_TBB