Index: basictypes.h =================================================================== --- basictypes.h (revision 4073) +++ basictypes.h (working copy) @@ -24,10 +24,17 @@ // MS Visual C++ does not include stdint.h typedef __int64 int64; typedef unsigned __int64 uint64; +#define INT64_MAX LLONG_MAX +#define UINT64_MAX ULLONG_MAX #else #include +#include typedef int64_t int64; typedef uint64_t uint64; +#define INT64_MAX 9223372036854775807LL +#define UINT64_MAX 0xffffffffffffffffULL +//#define INT64_MAX numeric_limits::max() +//#define UINT64_MAX numeric_limits::max() #endif #endif // _BASICTYPES_H_ Index: bigfix.cpp =================================================================== --- bigfix.cpp (revision 4073) +++ bigfix.cpp (working copy) @@ -1,7 +1,12 @@ // bigfix.cpp // -// Copyright (C) 2001, Chris Laurel +// Copyright (C) 2007-2008, Chris Laurel // +// 128-bit fixed point (64.64) numbers for high-precision celestial +// coordinates. When you need millimeter accurate navigation across a scale +// of thousands of light years, double precision floating point numbers +// are inadequate. +// // This program is free software; 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 2 @@ -12,85 +17,110 @@ #include "bigfix.h" +static const double POW2_31 = 2147483648.0; +static const double POW2_32 = 4294967296.0; +static const double POW2_64 = POW2_32 * POW2_32; + +static const double WORD0_FACTOR = 1.0 / POW2_64; +static const double WORD1_FACTOR = 1.0 / POW2_32; +static const double WORD2_FACTOR = 1.0; +static const double WORD3_FACTOR = POW2_32; + + /*** Constructors ***/ +/* +// Compute the additive inverse of a 128-bit twos complement value +// represented by two 64-bit values. +inline void negate128(uint64& hi, uint64& lo) +{ + // For a twos-complement number, -n = ~n + 1 + hi = ~hi; + lo = ~lo; + lo++; + if (lo == 0) + hi++; +} +*/ + + // Create a BigFix initialized to zero BigFix::BigFix() { - for (int i = 0; i < N_WORDS; i++) - n[i] = 0; + hi = 0; + lo = 0; } -BigFix::BigFix(int i) +BigFix::BigFix(uint64 i) { - n[N_WORDS / 2] = (unsigned short) (i & 0xffff); - n[N_WORDS / 2 + 1] = (unsigned short) ((i >> 16) & 0xffff); - - // Sign extend if negative - if (i < 0) - { - for (int j = N_WORDS / 2 + 2; j < N_WORDS; j++) - n[j] = 0xffff; - } + hi = i; + lo = 0; } BigFix::BigFix(double d) { - int sign = 1; + bool isNegative = false; + // Handle negative values by inverting them before conversion, + // then inverting the converted value. if (d < 0) { - sign = -1; - d = -d; + isNegative = true; + d = -d; } + + // Need to break the number into 32-bit chunks because a 64-bit + // integer has more bits of precision than a double. + double e = floor(d * (1.0 / WORD3_FACTOR)); + if (e < POW2_31) + { + uint32 w3 = (uint32) e; + d -= w3 * WORD3_FACTOR; + uint32 w2 = (uint32) (d * (1.0 / WORD2_FACTOR)); + d -= w2 * WORD2_FACTOR; + uint32 w1 = (uint32) (d * (1.0 / WORD1_FACTOR)); + d -= w1 * WORD1_FACTOR; + uint32 w0 = (uint32) (d * (1.0 / WORD0_FACTOR)); - double e = floor(d / (65536.0 * 65536.0 * 65536.0)); - - // Check for overflow - if (e < 32767) - { - n[7] = (unsigned short) e; - d -= n[7] * 65536.0 * 65536.0 * 65536.0; - n[6] = (unsigned short) (d / (65536.0 * 65536.0)); - d -= n[6] * 65536.0 * 65536.0; - n[5] = (unsigned short) (d / 65536.0); - d -= n[5] * 65536.0; - n[4] = (unsigned short) d; - d -= n[4]; - n[3] = (unsigned short) (d * 65536.0); - d -= n[3] / 65536.0; - n[2] = (unsigned short) (d * 65536.0 * 65536.0); - d -= n[2] / (65536.0 * 65536.0); - n[1] = (unsigned short) (d * 65536.0 * 65536.0 * 65536.0); - d -= n[1] / (65536.0 * 65536.0 * 65536.0); - n[0] = (unsigned short) (d * 65536.0 * 65536.0 * 65536.0 * 65536.0); + hi = ((uint64) w3 << 32) | w2; + lo = ((uint64) w1 << 32) | w0; } - if (sign < 0) - *this = -*this; + if (isNegative) + negate128(hi, lo); } BigFix::operator double() const { - double d = 0; - double e = 1.0 / (65536.0 * 65536.0 * 65536.0 * 65536.0); - BigFix f; + // Handle negative values by inverting them before conversion, + // then inverting the converted value. + int sign = 1; + uint64 l = lo; + uint64 h = hi; - if ((n[N_WORDS - 1] & 0x8000) == 0) - f = *this; - else - f = -*this; - - for (int i = 0; i < N_WORDS; i++) + if (isNegative()) { - d += e * f.n[i]; - e *= 65536; + negate128(h, l); + sign = -1; } - return ((n[N_WORDS - 1] & 0x8000) == 0) ? d : -d; + // Need to break the number into 32-bit chunks because a 64-bit + // integer has more bits of precision than a double. + uint32 w0 = l & 0xffffffff; + uint32 w1 = l >> 32; + uint32 w2 = h & 0xffffffff; + uint32 w3 = h >> 32; + double d; + + d = (w0 * WORD0_FACTOR + + w1 * WORD1_FACTOR + + w2 * WORD2_FACTOR + + w3 * WORD3_FACTOR) * sign; + + return d; } @@ -100,95 +130,188 @@ } -BigFix BigFix::operator-() const +bool operator==(const BigFix& a, const BigFix& b) { - int i; - BigFix f; + return a.hi == b.hi && a.lo == b.lo; +} - for (i = 0; i < N_WORDS; i++) - f.n[i] = ~n[i]; - unsigned int carry = 1; - i = 0; - while (carry != 0 && i < N_WORDS) - { - unsigned int m = f.n[i] + carry; - f.n[i] = (unsigned short) (m & 0xffff); - carry = m >> 16; - i++; - } - - return f; +bool operator!=(const BigFix& a, const BigFix& b) +{ + return a.hi != b.hi || a.lo != b.lo; } -BigFix operator+(BigFix a, BigFix b) +bool operator<(const BigFix& a, const BigFix& b) { - unsigned int carry = 0; - BigFix c; - - for (int i = 0; i < N_WORDS; i++) + if (a.isNegative() == b.isNegative()) { - unsigned int m = a.n[i] + b.n[i] + carry; - c.n[i] = (unsigned short) (m & 0xffff); - carry = m >> 16; + if (a.hi == b.hi) + return a.lo < b.lo; + else + return a.hi < b.hi; } - - return c; + else + { + return a.isNegative(); + } } -BigFix operator-(BigFix a, BigFix b) +bool operator>(const BigFix& a, const BigFix& b) { - return a + -b; + return b < a; } -bool operator==(BigFix a, BigFix b) +// TODO: probably faster to do this by converting the double to fixed +// point and using the fix*fix multiplication. +BigFix operator*(BigFix f, double d) { - for (int i = 0; i < N_WORDS; i++) - if (a.n[i] != b.n[i]) - return false; - return true; + // Need to break the number into 32-bit chunks because a 64-bit + // integer has more bits of precision than a double. + uint32 w0 = f.lo & 0xffffffff; + uint32 w1 = f.lo >> 32; + uint32 w2 = f.hi & 0xffffffff; + uint32 w3 = f.hi >> 32; + + return BigFix(w0 * d * WORD0_FACTOR) + + BigFix(w1 * d * WORD1_FACTOR) + + BigFix(w2 * d * WORD2_FACTOR) + + BigFix(w3 * d * WORD3_FACTOR); } -bool operator!=(BigFix a, BigFix b) +BigFix operator*(const BigFix& a, const BigFix& b) { - return !(a == b); + // Multiply two fixed point values together using partial products. + + uint64 ah = a.hi; + uint64 al = a.lo; + if (a.isNegative()) + BigFix::negate128(ah, al); + + uint64 bh = b.hi; + uint64 bl = b.lo; + if (b.isNegative()) + BigFix::negate128(bh, bl); + + // Break the values down into 32-bit words so that the partial products + // will fit into 64-bit words. + uint64 aw[4]; + aw[0] = al & 0xffffffff; + aw[1] = al >> 32; + aw[2] = ah & 0xffffffff; + aw[3] = ah >> 32; + + uint64 bw[4]; + bw[0] = bl & 0xffffffff; + bw[1] = bl >> 32; + bw[2] = bh & 0xffffffff; + bw[3] = bh >> 32; + + // Set the high and low non-zero words; this will + // speed up multiplicatoin of integers and values + // less than one. + unsigned int hiworda = ah == 0 ? 1 : 3; + unsigned int loworda = al == 0 ? 2 : 0; + unsigned int hiwordb = bh == 0 ? 1 : 3; + unsigned int lowordb = bl == 0 ? 2 : 0; + + uint32 result[8]; + + unsigned int i; + for (i = 0; i < 8; i++) + result[i] = 0; + + for (i = lowordb; i <= hiwordb; i++) + { + uint32 carry = 0; + + unsigned int j; + for (j = loworda; j <= hiworda; j++) + { + uint64 partialProd = aw[j] * bw[i]; + + // This sum will never overflow. Let N = 2^32; + // the max value of the partial product is (N-1)^2. + // The max values of result[i + j] and carry are + // N-1. Thus the max value of the sum is + // (N-1)^2 + 2(N-1) = (N^2 - 2N + 1) + 2(N-1) = N^2-1 + uint64 q = (uint64) result[i + j] + partialProd + (uint64) carry; + carry = (uint32) (q >> 32); + result[i + j] = (uint32) (q & 0xffffffff); + } + + result[i + j] = carry; + } + + // TODO: check for overflow + // (as simple as result[0] != 0 || result[1] != 0 || highbit(result[2])) + BigFix c; + c.lo = (uint64) result[2] + ((uint64) result[3] << 32); + c.hi = (uint64) result[4] + ((uint64) result[5] << 32); + + bool resultNegative = a.isNegative() != b.isNegative(); + if (resultNegative) + return -c; + else + return c; } -int BigFix::sign() +int BigFix::sign() const { - if ((n[N_WORDS - 1] & 0x8000) != 0) + + if (hi == 0 && lo == 0) + return 0; + else if (hi > INT64_MAX) return -1; - - for (int i = 0; i < N_WORDS - 1; i++) - if (n[N_WORDS] != 0) - return 1; - - return 0; + else + return 1; } // For debugging void BigFix::dump() { - for (int i = 7; i >= 0; i--) - printf("%04x ", n[i]); - printf("\n"); + printf("%08x %08x %08x %08x", + (uint32) (hi >> 32), + (uint32) (hi & 0xffffffff), + (uint32) (lo >> 32), + (uint32) (lo & 0xffffffff)); } +#if 0 +int main(int argc, char* argv[]) +{ + if (argc != 3) + return 1; + double a = 0.0; + if (sscanf(argv[1], "%lf", &a) != 1) + return 1; + + double b = 0.0; + if (sscanf(argv[2], "%lf", &b) != 1) + return 1; + + printf(" sum:\n%f\n%f\n", a + b, (double) (BigFix(a) + BigFix(b))); + printf(" diff:\n%f\n%f\n", a - b, (double) (BigFix(a) - BigFix(b))); + printf("product:\n%f\n%f\n", a * b, (double) (BigFix(a) * BigFix(b))); + printf(" lt: %u %u\n", a < b, BigFix(a) < BigFix(b)); + + return 0; +} +#endif + static unsigned char alphabet[65] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/"; BigFix::BigFix(const std::string& val) { static char inalphabet[256], decoder[256]; int i, bits, c, char_count; - /*int errors = 0; Unused*/; for (i = (sizeof alphabet) - 1; i >= 0 ; i--) { @@ -196,6 +319,10 @@ decoder[alphabet[i]] = i; } + uint16 n[8]; + + // Code from original BigFix class to convert base64 string into + // array of 8 16-bit values. for (i = 0; i < 8 ; i++) n[i] = 0; @@ -204,7 +331,7 @@ i = 0; - for (int j = 0; j < (int)val.length(); j++) + for (int j = 0; j < (int) val.length(); j++) { c = val[j]; if (c == '=') @@ -252,12 +379,38 @@ if (i & 1) n[i/2] >>= 8; + + // Now, convert the 8 16-bit values to a 2 64-bit values + lo = ((uint64) n[0] | + ((uint64) n[1] << 16) | + ((uint64) n[2] << 32) | + ((uint64) n[3] << 48)); + hi = ((uint64) n[4] | + ((uint64) n[5] << 16) | + ((uint64) n[6] << 32) | + ((uint64) n[7] << 48)); } std::string BigFix::toString() { + // Old BigFix class used 8 16-bit words. The bulk of this function + // is copied from that class, so first we'll convert from two + // 64-bit words to 8 16-bit words so that the old code can work + // as-is. + unsigned short n[8]; + n[0] = lo & 0xffff; + n[1] = (lo >> 16) & 0xffff; + n[2] = (lo >> 32) & 0xffff; + n[3] = (lo >> 48) & 0xffff; + + n[4] = hi & 0xffff; + n[5] = (hi >> 16) & 0xffff; + n[6] = (hi >> 32) & 0xffff; + n[7] = (hi >> 48) & 0xffff; + + // Conversion using code from the original BigFix class. std::string encoded(""); int bits, c, char_count, started, i, j; @@ -313,5 +466,3 @@ return encoded; } - - Index: bigfix.h =================================================================== --- bigfix.h (revision 4073) +++ bigfix.h (working copy) @@ -1,49 +1,145 @@ // bigfix.h // -// Copyright (C) 2001, Chris Laurel +// Copyright (C) 2007-2008, Chris Laurel // // 128-bit fixed point (64.64) numbers for high-precision celestial // coordinates. When you need millimeter accurate navigation across a scale -// of thousands of light years, doubles just don't do it. +// of thousands of light years, double precision floating point numbers +// are inadequate. // // This program is free software; 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 2 // of the License, or (at your option) any later version. +#ifndef _CELUTIL_BIGFIX64_H_ +#define _CELUTIL_BIGFIX64_H_ -#ifndef _BIGFIX_H_ -#define _BIGFIX_H_ - -#define N_WORDS 8 #include +#include "basictypes.h" +/*! 64.64 signed fixed point numbers. + */ + class BigFix { -public: - BigFix(); - BigFix(int); - BigFix(double); - BigFix(const std::string&); + public: + BigFix(); + BigFix(uint64); + BigFix(double); + BigFix(const std::string&); - operator double() const; - operator float() const; + operator double() const; + operator float() const; - BigFix operator-() const; + BigFix operator-() const; + BigFix operator+=(const BigFix&); + BigFix operator-=(const BigFix&); - friend BigFix operator+(BigFix, BigFix); - friend BigFix operator-(BigFix, BigFix); - friend bool operator==(BigFix, BigFix); - friend bool operator!=(BigFix, BigFix); + friend BigFix operator+(const BigFix&, const BigFix&); + friend BigFix operator-(const BigFix&, const BigFix&); + friend BigFix operator*(const BigFix&, const BigFix&); + friend BigFix operator*(BigFix, double); + friend bool operator==(const BigFix&, const BigFix&); + friend bool operator!=(const BigFix&, const BigFix&); + friend bool operator<(const BigFix&, const BigFix&); + friend bool operator>(const BigFix&, const BigFix&); - int sign(); + int sign() const; - // for debugging - void dump(); - std::string toString(); + // for debugging + void dump(); + std::string toString(); -private: - unsigned short n[N_WORDS]; // high 16-bits is n[N_WORDS - 1] + private: + bool isNegative() const + { + return hi > INT64_MAX; + } + + static void BigFix::negate128(uint64& hi, uint64& lo); + + private: + uint64 hi; + uint64 lo; }; -#endif // _BIGFIX_H_ + +// Compute the additive inverse of a 128-bit twos complement value +// represented by two 64-bit values. +inline void BigFix::negate128(uint64& hi, uint64& lo) +{ + // For a twos-complement number, -n = ~n + 1 + hi = ~hi; + lo = ~lo; + lo++; + if (lo == 0) + hi++; +} + +inline BigFix BigFix::operator-() const +{ + BigFix result = *this; + + negate128(result.hi, result.lo); + + return result; +} + + +inline BigFix BigFix::operator+=(const BigFix& a) +{ + lo += a.lo; + hi += a.hi; + + // carry + if (lo < a.lo) + hi++; + + return *this; +} + + +inline BigFix BigFix::operator-=(const BigFix& a) +{ + lo -= a.lo; + hi -= a.hi; + + // borrow + if (lo > a.lo) + hi--; + + return *this; +} + + +inline BigFix operator+(const BigFix& a, const BigFix& b) +{ + BigFix c; + + c.lo = a.lo + b.lo; + c.hi = a.hi + b.hi; + + // carry + if (c.lo < a.lo) + c.hi++; + + return c; +} + + +inline BigFix operator-(const BigFix& a, const BigFix& b) +{ + BigFix c; + + c.lo = a.lo - b.lo; + c.hi = a.hi - b.hi; + + // borrow + if (c.lo > a.lo) + c.hi--; + + return c; +} + +#endif // _CELUTIL_BIGFIX64_H_