1 // bigfix.cpp
2 //
3 // Copyright (C) 2007-2008, Chris Laurel <claurel@shatters.net>
4 //
5 // 128-bit fixed point (64.64) numbers for high-precision celestial
6 // coordinates. When you need millimeter accurate navigation across a scale
7 // of thousands of light years, double precision floating point numbers
8 // are inadequate.
9 //
10 // This program is free software; you can redistribute it and/or
11 // modify it under the terms of the GNU General Public License
12 // as published by the Free Software Foundation; either version 2
13 // of the License, or (at your option) any later version.
14
15 #include <math.h>
16 #include <stdio.h>
17 #include "bigfix.h"
18
19
20 static const double POW2_31 = 2147483648.0;
21 static const double POW2_32 = 4294967296.0;
22 static const double POW2_64 = POW2_32 * POW2_32;
23
24 static const double WORD0_FACTOR = 1.0 / POW2_64;
25 static const double WORD1_FACTOR = 1.0 / POW2_32;
26 static const double WORD2_FACTOR = 1.0;
27 static const double WORD3_FACTOR = POW2_32;
28
29
30 /*** Constructors ***/
31
32 /*
33 // Compute the additive inverse of a 128-bit twos complement value
34 // represented by two 64-bit values.
35 inline void negate128(uint64& hi, uint64& lo)
36 {
37 // For a twos-complement number, -n = ~n + 1
38 hi = ~hi;
39 lo = ~lo;
40 lo++;
41 if (lo == 0)
42 hi++;
43 }
44 */
45
46
47 // Create a BigFix initialized to zero
BigFix()48 BigFix::BigFix()
49 {
50 hi = 0;
51 lo = 0;
52 }
53
54
BigFix(uint64 i)55 BigFix::BigFix(uint64 i)
56 {
57 hi = i;
58 lo = 0;
59 }
60
61
BigFix(double d)62 BigFix::BigFix(double d)
63 {
64 bool isNegative = false;
65
66 // Handle negative values by inverting them before conversion,
67 // then inverting the converted value.
68 if (d < 0)
69 {
70 isNegative = true;
71 d = -d;
72 }
73
74 // Need to break the number into 32-bit chunks because a 64-bit
75 // integer has more bits of precision than a double.
76 double e = floor(d * (1.0 / WORD3_FACTOR));
77 if (e < POW2_31)
78 {
79 uint32 w3 = (uint32) e;
80 d -= w3 * WORD3_FACTOR;
81 uint32 w2 = (uint32) (d * (1.0 / WORD2_FACTOR));
82 d -= w2 * WORD2_FACTOR;
83 uint32 w1 = (uint32) (d * (1.0 / WORD1_FACTOR));
84 d -= w1 * WORD1_FACTOR;
85 uint32 w0 = (uint32) (d * (1.0 / WORD0_FACTOR));
86
87 hi = ((uint64) w3 << 32) | w2;
88 lo = ((uint64) w1 << 32) | w0;
89 }
90
91 if (isNegative)
92 negate128(hi, lo);
93 }
94
95
operator double() const96 BigFix::operator double() const
97 {
98 // Handle negative values by inverting them before conversion,
99 // then inverting the converted value.
100 int sign = 1;
101 uint64 l = lo;
102 uint64 h = hi;
103
104 if (isNegative())
105 {
106 negate128(h, l);
107 sign = -1;
108 }
109
110 // Need to break the number into 32-bit chunks because a 64-bit
111 // integer has more bits of precision than a double.
112 uint32 w0 = l & 0xffffffff;
113 uint32 w1 = l >> 32;
114 uint32 w2 = h & 0xffffffff;
115 uint32 w3 = h >> 32;
116 double d;
117
118 d = (w0 * WORD0_FACTOR +
119 w1 * WORD1_FACTOR +
120 w2 * WORD2_FACTOR +
121 w3 * WORD3_FACTOR) * sign;
122
123 return d;
124 }
125
126
operator float() const127 BigFix::operator float() const
128 {
129 return (float) (double) *this;
130 }
131
132
operator ==(const BigFix & a,const BigFix & b)133 bool operator==(const BigFix& a, const BigFix& b)
134 {
135 return a.hi == b.hi && a.lo == b.lo;
136 }
137
138
operator !=(const BigFix & a,const BigFix & b)139 bool operator!=(const BigFix& a, const BigFix& b)
140 {
141 return a.hi != b.hi || a.lo != b.lo;
142 }
143
144
operator <(const BigFix & a,const BigFix & b)145 bool operator<(const BigFix& a, const BigFix& b)
146 {
147 if (a.isNegative() == b.isNegative())
148 {
149 if (a.hi == b.hi)
150 return a.lo < b.lo;
151 else
152 return a.hi < b.hi;
153 }
154 else
155 {
156 return a.isNegative();
157 }
158 }
159
160
operator >(const BigFix & a,const BigFix & b)161 bool operator>(const BigFix& a, const BigFix& b)
162 {
163 return b < a;
164 }
165
166
167 // TODO: probably faster to do this by converting the double to fixed
168 // point and using the fix*fix multiplication.
operator *(BigFix f,double d)169 BigFix operator*(BigFix f, double d)
170 {
171 // Need to break the number into 32-bit chunks because a 64-bit
172 // integer has more bits of precision than a double.
173 uint32 w0 = f.lo & 0xffffffff;
174 uint32 w1 = f.lo >> 32;
175 uint32 w2 = f.hi & 0xffffffff;
176 uint32 w3 = f.hi >> 32;
177
178 return BigFix(w0 * d * WORD0_FACTOR) +
179 BigFix(w1 * d * WORD1_FACTOR) +
180 BigFix(w2 * d * WORD2_FACTOR) +
181 BigFix(w3 * d * WORD3_FACTOR);
182 }
183
184
185 /*! Multiply two BigFix values together. This function does not check for
186 * overflow. This is not a problem in Celestia, where it is used exclusively
187 * in multiplications where one multiplicand has absolute value <= 1.0.
188 */
operator *(const BigFix & a,const BigFix & b)189 BigFix operator*(const BigFix& a, const BigFix& b)
190 {
191 // Multiply two fixed point values together using partial products.
192
193 uint64 ah = a.hi;
194 uint64 al = a.lo;
195 if (a.isNegative())
196 BigFix::negate128(ah, al);
197
198 uint64 bh = b.hi;
199 uint64 bl = b.lo;
200 if (b.isNegative())
201 BigFix::negate128(bh, bl);
202
203 // Break the values down into 32-bit words so that the partial products
204 // will fit into 64-bit words.
205 uint64 aw[4];
206 aw[0] = al & 0xffffffff;
207 aw[1] = al >> 32;
208 aw[2] = ah & 0xffffffff;
209 aw[3] = ah >> 32;
210
211 uint64 bw[4];
212 bw[0] = bl & 0xffffffff;
213 bw[1] = bl >> 32;
214 bw[2] = bh & 0xffffffff;
215 bw[3] = bh >> 32;
216
217 // Set the high and low non-zero words; this will
218 // speed up multiplicatoin of integers and values
219 // less than one.
220 unsigned int hiworda = ah == 0 ? 1 : 3;
221 unsigned int loworda = al == 0 ? 2 : 0;
222 unsigned int hiwordb = bh == 0 ? 1 : 3;
223 unsigned int lowordb = bl == 0 ? 2 : 0;
224
225 uint32 result[8];
226
227 unsigned int i;
228 for (i = 0; i < 8; i++)
229 result[i] = 0;
230
231 for (i = lowordb; i <= hiwordb; i++)
232 {
233 uint32 carry = 0;
234
235 unsigned int j;
236 for (j = loworda; j <= hiworda; j++)
237 {
238 uint64 partialProd = aw[j] * bw[i];
239
240 // This sum will never overflow. Let N = 2^32;
241 // the max value of the partial product is (N-1)^2.
242 // The max values of result[i + j] and carry are
243 // N-1. Thus the max value of the sum is
244 // (N-1)^2 + 2(N-1) = (N^2 - 2N + 1) + 2(N-1) = N^2-1
245 uint64 q = (uint64) result[i + j] + partialProd + (uint64) carry;
246 carry = (uint32) (q >> 32);
247 result[i + j] = (uint32) (q & 0xffffffff);
248 }
249
250 result[i + j] = carry;
251 }
252
253 // TODO: check for overflow
254 // (as simple as result[0] != 0 || result[1] != 0 || highbit(result[2]))
255 BigFix c;
256 c.lo = (uint64) result[2] + ((uint64) result[3] << 32);
257 c.hi = (uint64) result[4] + ((uint64) result[5] << 32);
258
259 bool resultNegative = a.isNegative() != b.isNegative();
260 if (resultNegative)
261 return -c;
262 else
263 return c;
264 }
265
266
sign() const267 int BigFix::sign() const
268 {
269
270 if (hi == 0 && lo == 0)
271 return 0;
272 else if (hi > INT64_MAX)
273 return -1;
274 else
275 return 1;
276 }
277
278
279 // For debugging
dump()280 void BigFix::dump()
281 {
282 printf("%08x %08x %08x %08x",
283 (uint32) (hi >> 32),
284 (uint32) (hi & 0xffffffff),
285 (uint32) (lo >> 32),
286 (uint32) (lo & 0xffffffff));
287 }
288
289
290 #if 0
291 int main(int argc, char* argv[])
292 {
293 if (argc != 3)
294 return 1;
295
296 double a = 0.0;
297 if (sscanf(argv[1], "%lf", &a) != 1)
298 return 1;
299
300 double b = 0.0;
301 if (sscanf(argv[2], "%lf", &b) != 1)
302 return 1;
303
304 printf(" sum:\n%f\n%f\n", a + b, (double) (BigFix(a) + BigFix(b)));
305 printf(" diff:\n%f\n%f\n", a - b, (double) (BigFix(a) - BigFix(b)));
306 printf("product:\n%f\n%f\n", a * b, (double) (BigFix(a) * BigFix(b)));
307 printf(" lt: %u %u\n", a < b, BigFix(a) < BigFix(b));
308
309 return 0;
310 }
311 #endif
312
313 static unsigned char alphabet[65] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/";
314
BigFix(const std::string & val)315 BigFix::BigFix(const std::string& val)
316 {
317 static char inalphabet[256], decoder[256];
318 int i, bits, c, char_count;
319
320 for (i = (sizeof alphabet) - 1; i >= 0 ; i--)
321 {
322 inalphabet[alphabet[i]] = 1;
323 decoder[alphabet[i]] = i;
324 }
325
326 uint16 n[8];
327
328 // Code from original BigFix class to convert base64 string into
329 // array of 8 16-bit values.
330 for (i = 0; i < 8 ; i++)
331 n[i] = 0;
332
333 char_count = 0;
334 bits = 0;
335
336 i = 0;
337
338 for (int j = 0; j < (int) val.length(); j++)
339 {
340 c = val[j];
341 if (c == '=')
342 break;
343 if (c > 255 || !inalphabet[c])
344 continue;
345 bits += decoder[c];
346 char_count++;
347 if (char_count == 4)
348 {
349 n[i/2] >>= 8;
350 n[i/2] += (bits >> 8) & 0xff00;
351 i++;
352 n[i/2] >>= 8;
353 n[i/2] += bits & 0xff00;
354 i++;
355 n[i/2] >>= 8;
356 n[i/2] += (bits << 8) & 0xff00;
357 i++;
358 bits = 0;
359 char_count = 0;
360 }
361 else
362 {
363 bits <<= 6;
364 }
365 }
366
367 switch (char_count)
368 {
369 case 2:
370 n[i/2] >>= 8;
371 n[i/2] += (bits >> 2) & 0xff00;
372 i++;
373 break;
374 case 3:
375 n[i/2] >>= 8;
376 n[i/2] += (bits >> 8) & 0xff00;
377 i++;
378 n[i/2] >>= 8;
379 n[i/2] += bits & 0xff00;
380 i++;
381 break;
382 }
383
384 if (i & 1)
385 n[i/2] >>= 8;
386
387 // Now, convert the 8 16-bit values to a 2 64-bit values
388 lo = ((uint64) n[0] |
389 ((uint64) n[1] << 16) |
390 ((uint64) n[2] << 32) |
391 ((uint64) n[3] << 48));
392 hi = ((uint64) n[4] |
393 ((uint64) n[5] << 16) |
394 ((uint64) n[6] << 32) |
395 ((uint64) n[7] << 48));
396 }
397
398
toString()399 std::string BigFix::toString()
400 {
401 // Old BigFix class used 8 16-bit words. The bulk of this function
402 // is copied from that class, so first we'll convert from two
403 // 64-bit words to 8 16-bit words so that the old code can work
404 // as-is.
405 unsigned short n[8];
406
407 n[0] = lo & 0xffff;
408 n[1] = (lo >> 16) & 0xffff;
409 n[2] = (lo >> 32) & 0xffff;
410 n[3] = (lo >> 48) & 0xffff;
411
412 n[4] = hi & 0xffff;
413 n[5] = (hi >> 16) & 0xffff;
414 n[6] = (hi >> 32) & 0xffff;
415 n[7] = (hi >> 48) & 0xffff;
416
417 // Conversion using code from the original BigFix class.
418 std::string encoded("");
419 int bits, c, char_count, started, i, j;
420
421 char_count = 0;
422 bits = 0;
423 started = 0;
424
425 // Find first significant (non null) byte
426 i = 16;
427 do {
428 i--;
429 c = n[i/2];
430 if (i & 1) c >>= 8;
431 c &= 0xff;
432 } while ((c == 0) && (i != 0));
433
434 if (i == 0)
435 return encoded;
436
437 // Then we encode starting by the LSB (i+1 bytes to encode)
438 j = 0;
439 while (j <= i)
440 {
441 c = n[j/2];
442 if ( j & 1 ) c >>= 8;
443 c &= 0xff;
444 j++;
445 bits += c;
446 char_count++;
447 if (char_count == 3)
448 {
449 encoded += alphabet[bits >> 18];
450 encoded += alphabet[(bits >> 12) & 0x3f];
451 encoded += alphabet[(bits >> 6) & 0x3f];
452 encoded += alphabet[bits & 0x3f];
453 bits = 0;
454 char_count = 0;
455 }
456 else
457 {
458 bits <<= 8;
459 }
460 }
461
462 if (char_count != 0)
463 {
464 bits <<= 16 - (8 * char_count);
465 encoded += alphabet[bits >> 18];
466 encoded += alphabet[(bits >> 12) & 0x3f];
467 if (char_count != 1)
468 encoded += alphabet[(bits >> 6) & 0x3f];
469 }
470
471 return encoded;
472 }
473