source: trunk/src/real.cpp @ 988

Last change on this file since 988 was 988, checked in by sam, 11 years ago

core: add exp() for real numbers and fix the == operator.

File size: 13.9 KB
Line 
1//
2// Lol Engine
3//
4// Copyright: (c) 2010-2011 Sam Hocevar <sam@hocevar.net>
5//   This program is free software; you can redistribute it and/or
6//   modify it under the terms of the Do What The Fuck You Want To
7//   Public License, Version 2, as published by Sam Hocevar. See
8//   http://sam.zoy.org/projects/COPYING.WTFPL for more details.
9//
10
11#if defined HAVE_CONFIG_H
12#   include "config.h"
13#endif
14
15#include <cstring>
16#include <cstdio>
17
18#include "core.h"
19
20using namespace std;
21
22namespace lol
23{
24
25real::real(float f) { *this = (double)f; }
26real::real(int i) { *this = (double)i; }
27real::real(unsigned int i) { *this = (double)i; }
28
29real::real(double d)
30{
31    union { double d; uint64_t x; } u = { d };
32
33    uint32_t sign = (u.x >> 63) << 31;
34    uint32_t exponent = (u.x << 1) >> 53;
35
36    switch (exponent)
37    {
38    case 0x00:
39        m_signexp = sign;
40        break;
41    case 0x7ff:
42        m_signexp = sign | 0x7fffffffu;
43        break;
44    default:
45        m_signexp = sign | (exponent + (1 << 30) - (1 << 10));
46        break;
47    }
48
49    m_mantissa[0] = u.x >> 36;
50    m_mantissa[1] = u.x >> 20;
51    m_mantissa[2] = u.x >> 4;
52    m_mantissa[3] = u.x << 12;
53    memset(m_mantissa + 4, 0, sizeof(m_mantissa) - 4 * sizeof(m_mantissa[0]));
54}
55
56real::operator float() const { return (float)(double)(*this); }
57real::operator int() const { return (int)(double)(*this); }
58real::operator unsigned int() const { return (unsigned int)(double)(*this); }
59
60real::operator double() const
61{
62    union { double d; uint64_t x; } u;
63
64    /* Get sign */
65    u.x = m_signexp >> 31;
66    u.x <<= 11;
67
68    /* Compute new exponent */
69    uint32_t exponent = (m_signexp << 1) >> 1;
70    int e = (int)exponent - (1 << 30) + (1 << 10);
71
72    if (e < 0)
73        u.x <<= 52;
74    else if (e >= 0x7ff)
75    {
76        u.x |= 0x7ff;
77        u.x <<= 52;
78    }
79    else
80    {
81        u.x |= e;
82
83        /* Store mantissa if necessary */
84        u.x <<= 16;
85        u.x |= m_mantissa[0];
86        u.x <<= 16;
87        u.x |= m_mantissa[1];
88        u.x <<= 16;
89        u.x |= m_mantissa[2];
90        u.x <<= 4;
91        u.x |= m_mantissa[3] >> 12;
92        /* Rounding */
93        u.x += (m_mantissa[3] >> 11) & 1;
94    }
95
96    return u.d;
97}
98
99real real::operator -() const
100{
101    real ret = *this;
102    ret.m_signexp ^= 0x80000000u;
103    return ret;
104}
105
106real real::operator +(real const &x) const
107{
108    if (x.m_signexp << 1 == 0)
109        return *this;
110
111    /* Ensure both arguments are positive. Otherwise, switch signs,
112     * or replace + with -. */
113    if (m_signexp >> 31)
114        return -(-*this + -x);
115
116    if (x.m_signexp >> 31)
117        return *this - (-x);
118
119    /* Ensure *this is the larger exponent (no need to be strictly larger,
120     * as in subtraction). Otherwise, switch. */
121    if ((m_signexp << 1) < (x.m_signexp << 1))
122        return x + *this;
123
124    real ret;
125
126    int e1 = m_signexp - (1 << 30) + 1;
127    int e2 = x.m_signexp - (1 << 30) + 1;
128
129    int bigoff = (e1 - e2) / (sizeof(uint16_t) * 8);
130    int off = e1 - e2 - bigoff * (sizeof(uint16_t) * 8);
131
132    ret.m_signexp = m_signexp;
133
134    uint32_t carry = 0;
135    for (int i = BIGITS; i--; )
136    {
137        carry += m_mantissa[i];
138        if (i - bigoff >= 0)
139            carry += x.m_mantissa[i - bigoff] >> off;
140
141        if (i - bigoff > 0)
142            carry += (x.m_mantissa[i - bigoff - 1] << (16 - off)) & 0xffffu;
143        else if (i - bigoff == 0)
144            carry += 0x0001u << (16 - off);
145
146        ret.m_mantissa[i] = carry;
147        carry >>= 16;
148    }
149
150    /* Renormalise in case we overflowed the mantissa */
151    if (carry)
152    {
153        carry--;
154        for (int i = 0; i < BIGITS; i++)
155        {
156            uint16_t tmp = ret.m_mantissa[i];
157            ret.m_mantissa[i] = (carry << 15) | (tmp >> 1);
158            carry = tmp & 0x0001u;
159        }
160        ret.m_signexp++;
161    }
162
163    return ret;
164}
165
166real real::operator -(real const &x) const
167{
168    if (x.m_signexp << 1 == 0)
169        return *this;
170
171    /* Ensure both arguments are positive. Otherwise, switch signs,
172     * or replace - with +. */
173    if (m_signexp >> 31)
174        return -(-*this + x);
175
176    if (x.m_signexp >> 31)
177        return (*this) + (-x);
178
179    /* Ensure *this is larger than x */
180    if (*this < x)
181        return -(x - *this);
182
183    real ret;
184
185    int e1 = m_signexp - (1 << 30) + 1;
186    int e2 = x.m_signexp - (1 << 30) + 1;
187
188    int bigoff = (e1 - e2) / (sizeof(uint16_t) * 8);
189    int off = e1 - e2 - bigoff * (sizeof(uint16_t) * 8);
190
191    ret.m_signexp = m_signexp;
192
193    int32_t carry = 0;
194    for (int i = 0; i < bigoff; i++)
195    {
196        carry -= x.m_mantissa[BIGITS - i];
197        carry = (carry & 0xffff0000u) | (carry >> 16);
198    }
199    carry -= x.m_mantissa[BIGITS - 1 - bigoff] & ((1 << off) - 1);
200    carry /= (1 << off);
201
202    for (int i = BIGITS; i--; )
203    {
204        carry += m_mantissa[i];
205        if (i - bigoff >= 0)
206            carry -= x.m_mantissa[i - bigoff] >> off;
207
208        if (i - bigoff > 0)
209            carry -= (x.m_mantissa[i - bigoff - 1] << (16 - off)) & 0xffffu;
210        else if (i - bigoff == 0)
211            carry -= 0x0001u << (16 - off);
212
213        ret.m_mantissa[i] = carry;
214        carry = (carry & 0xffff0000u) | (carry >> 16);
215    }
216
217    carry += 1;
218
219    /* Renormalise if we underflowed the mantissa */
220    if (carry == 0)
221    {
222        /* How much do we need to shift the mantissa? FIXME: this could
223         * be computed above */
224        off = 0;
225        for (int i = 0; i < BIGITS; i++)
226        {
227            if (!ret.m_mantissa[i])
228            {
229                off += sizeof(uint16_t) * 8;
230                continue;
231            }
232
233            for (uint16_t tmp = ret.m_mantissa[i]; tmp < 0x8000u; tmp <<= 1)
234                off++;
235            break;
236        }
237        if (off == BIGITS * sizeof(uint16_t) * 8)
238            ret.m_signexp &= 0x80000000u;
239        else
240        {
241            off++; /* Shift one more to get rid of the leading one */
242            ret.m_signexp -= off;
243
244            bigoff = off / (sizeof(uint16_t) * 8);
245            off -= bigoff * sizeof(uint16_t) * 8;
246
247            for (int i = 0; i < BIGITS; i++)
248            {
249                uint16_t tmp = 0;
250                if (i + bigoff < BIGITS)
251                    tmp |= ret.m_mantissa[i + bigoff] << off;
252                if (i + bigoff + 1 < BIGITS)
253                    tmp |= ret.m_mantissa[i + bigoff + 1] >> (16 - off);
254                ret.m_mantissa[i] = tmp;
255            }
256        }
257    }
258
259    return ret;
260}
261
262real real::operator *(real const &x) const
263{
264    real ret;
265
266    if (m_signexp << 1 == 0 || x.m_signexp << 1 == 0)
267    {
268        ret = (m_signexp << 1 == 0) ? *this : x;
269        ret.m_signexp ^= x.m_signexp & 0x80000000u;
270        return ret;
271    }
272
273    ret.m_signexp = (m_signexp ^ x.m_signexp) & 0x80000000u;
274    int e = (m_signexp & 0x7fffffffu) - (1 << 30) + 1
275          + (x.m_signexp & 0x7fffffffu) - (1 << 30) + 1;
276
277    /* Accumulate low order product; no need to store it, we just
278     * want the carry value */
279    uint64_t carry = 0;
280    for (int i = 0; i < BIGITS; i++)
281    {
282        for (int j = 0; j < i + 1; j++)
283            carry += (uint32_t)m_mantissa[BIGITS - 1 - j]
284                   * (uint32_t)x.m_mantissa[BIGITS - 1 + j - i];
285        carry >>= 16;
286    }
287
288    for (int i = 0; i < BIGITS; i++)
289    {
290        for (int j = i + 1; j < BIGITS; j++)
291            carry += (uint32_t)m_mantissa[BIGITS - 1 - j]
292                   * (uint32_t)x.m_mantissa[j - 1 - i];
293
294        carry += m_mantissa[BIGITS - 1 - i];
295        carry += x.m_mantissa[BIGITS - 1 - i];
296        ret.m_mantissa[BIGITS - 1 - i] = carry & 0xffffu;
297        carry >>= 16;
298    }
299
300    /* Renormalise in case we overflowed the mantissa */
301    if (carry)
302    {
303        carry--;
304        for (int i = 0; i < BIGITS; i++)
305        {
306            uint16_t tmp = ret.m_mantissa[i];
307            ret.m_mantissa[i] = (carry << 15) | (tmp >> 1);
308            carry = tmp & 0x0001u;
309        }
310        e++;
311    }
312
313    ret.m_signexp |= e + (1 << 30) - 1;
314
315    return ret;
316}
317
318real real::operator /(real const &x) const
319{
320    return *this * fres(x);
321}
322
323real &real::operator +=(real const &x)
324{
325    real tmp = *this;
326    return *this = tmp + x;
327}
328
329real &real::operator -=(real const &x)
330{
331    real tmp = *this;
332    return *this = tmp - x;
333}
334
335real &real::operator *=(real const &x)
336{
337    real tmp = *this;
338    return *this = tmp * x;
339}
340
341real &real::operator /=(real const &x)
342{
343    real tmp = *this;
344    return *this = tmp / x;
345}
346
347bool real::operator ==(real const &x) const
348{
349    if ((m_signexp << 1) == 0 && (x.m_signexp << 1) == 0)
350        return true;
351
352    if (m_signexp != x.m_signexp)
353        return false;
354
355    return memcmp(m_mantissa, x.m_mantissa, sizeof(m_mantissa)) == 0;
356}
357
358bool real::operator !=(real const &x) const
359{
360    return !(*this == x);
361}
362
363bool real::operator <(real const &x) const
364{
365    /* Ensure both numbers are positive */
366    if (m_signexp >> 31)
367        return (x.m_signexp >> 31) ? -*this > -x : true;
368
369    if (x.m_signexp >> 31)
370        return false;
371
372    /* Compare all relevant bits */
373    if (m_signexp != x.m_signexp)
374        return m_signexp < x.m_signexp;
375
376    for (int i = 0; i < BIGITS; i++)
377        if (m_mantissa[i] != x.m_mantissa[i])
378            return m_mantissa[i] < x.m_mantissa[i];
379
380    return false;
381}
382
383bool real::operator <=(real const &x) const
384{
385    return !(*this > x);
386}
387
388bool real::operator >(real const &x) const
389{
390    /* Ensure both numbers are positive */
391    if (m_signexp >> 31)
392        return (x.m_signexp >> 31) ? -*this < -x : false;
393
394    if (x.m_signexp >> 31)
395        return true;
396
397    /* Compare all relevant bits */
398    if (m_signexp != x.m_signexp)
399        return m_signexp > x.m_signexp;
400
401    for (int i = 0; i < BIGITS; i++)
402        if (m_mantissa[i] != x.m_mantissa[i])
403            return m_mantissa[i] > x.m_mantissa[i];
404
405    return false;
406}
407
408bool real::operator >=(real const &x) const
409{
410    return !(*this < x);
411}
412
413real fres(real const &x)
414{
415    if (!(x.m_signexp << 1))
416    {
417        real ret = x;
418        ret.m_signexp = x.m_signexp | 0x7fffffffu;
419        ret.m_mantissa[0] = 0;
420        return ret;
421    }
422
423    /* Use the system's float inversion to approximate 1/x */
424    union { float f; uint32_t x; } u = { 1.0f }, v = { 1.0f };
425    v.x |= (uint32_t)x.m_mantissa[0] << 7;
426    v.x |= (uint32_t)x.m_mantissa[1] >> 9;
427    v.f = 1.0 / v.f;
428
429    real ret;
430    ret.m_mantissa[0] = (v.x >> 7) & 0xffffu;
431    ret.m_mantissa[1] = (v.x << 9) & 0xffffu;
432
433    uint32_t sign = x.m_signexp & 0x80000000u;
434    ret.m_signexp = sign;
435
436    int exponent = (x.m_signexp & 0x7fffffffu) + 1;
437    exponent = -exponent + (v.x >> 23) - (u.x >> 23);
438    ret.m_signexp |= (exponent - 1) & 0x7fffffffu;
439
440    /* Five steps of Newton-Raphson seems enough for 32-bigit reals. */
441    real two = 2;
442    ret = ret * (two - ret * x);
443    ret = ret * (two - ret * x);
444    ret = ret * (two - ret * x);
445    ret = ret * (two - ret * x);
446    ret = ret * (two - ret * x);
447
448    return ret;
449}
450
451real sqrt(real const &x)
452{
453    /* if zero, return x */
454    if (!(x.m_signexp << 1))
455        return x;
456
457    /* if negative, return NaN */
458    if (x.m_signexp >> 31)
459    {
460        real ret;
461        ret.m_signexp = 0x7fffffffu;
462        ret.m_mantissa[0] = 0xffffu;
463        return ret;
464    }
465
466    /* Use the system's float inversion to approximate 1/sqrt(x). First
467     * we construct a float in the [1..4[ range that has roughly the same
468     * mantissa as our real. Its exponent is 0 or 1, depending on the
469     * partity of x. The final exponent is 0, -1 or -2. We use the final
470     * exponent and final mantissa to pre-fill the result. */
471    union { float f; uint32_t x; } u = { 1.0f }, v = { 2.0f };
472    v.x -= ((x.m_signexp & 1) << 23);
473    v.x |= (uint32_t)x.m_mantissa[0] << 7;
474    v.x |= (uint32_t)x.m_mantissa[1] >> 9;
475    v.f = 1.0 / sqrtf(v.f);
476
477    real ret;
478    ret.m_mantissa[0] = (v.x >> 7) & 0xffffu;
479    ret.m_mantissa[1] = (v.x << 9) & 0xffffu;
480
481    uint32_t sign = x.m_signexp & 0x80000000u;
482    ret.m_signexp = sign;
483
484    int exponent = (x.m_signexp & 0x7fffffffu) - ((1 << 30) - 1);
485    exponent = - (exponent / 2) + (v.x >> 23) - (u.x >> 23);
486    ret.m_signexp |= (exponent + ((1 << 30) - 1)) & 0x7fffffffu;
487
488    /* Five steps of Newton-Raphson seems enough for 32-bigit reals. */
489    real three = 3;
490    ret = ret * (three - ret * ret * x);
491    ret.m_signexp--;
492    ret = ret * (three - ret * ret * x);
493    ret.m_signexp--;
494    ret = ret * (three - ret * ret * x);
495    ret.m_signexp--;
496    ret = ret * (three - ret * ret * x);
497    ret.m_signexp--;
498    ret = ret * (three - ret * ret * x);
499    ret.m_signexp--;
500
501    return ret * x;
502}
503
504real exp(real const &x)
505{
506    int square = 0;
507
508    /* FIXME: this is slow. Find a better way to approximate exp(x) for
509     * large values. */
510    real tmp = x, one = 1.0;
511    while (tmp > one)
512    {
513        tmp.m_signexp--;
514        square++;
515    }
516
517    real ret = 1.0, fact = 1.0, xn = tmp;
518
519    for (int i = 1; i < 100; i++)
520    {
521        fact *= (real)i;
522        ret += xn / fact;
523        xn *= tmp;
524    }
525
526    for (int i = 0; i < square; i++)
527        ret = ret * ret;
528
529    return ret;
530}
531
532void real::print(int ndigits) const
533{
534    real const r1 = 1, r10 = 10;
535    real x = *this;
536
537    if (x.m_signexp >> 31)
538    {
539        printf("-");
540        x = -x;
541    }
542
543    /* Normalise x so that mantissa is in [1..9.999] */
544    int exponent = 0;
545    if (x.m_signexp)
546    {
547        for (real div = r1, newdiv; true; div = newdiv)
548        {
549            newdiv = div * r10;
550            if (x < newdiv)
551            {
552                x /= div;
553                break;
554            }
555            exponent++;
556        }
557        for (real mul = 1, newx; true; mul *= r10)
558        {
559            newx = x * mul;
560            if (newx >= r1)
561            {
562                x = newx;
563                break;
564            }
565            exponent--;
566        }
567    }
568
569    /* Print digits */
570    for (int i = 0; i < ndigits; i++)
571    {
572        int digit = (int)x;
573        printf("%i", digit);
574        if (i == 0)
575            printf(".");
576        x -= real(digit);
577        x *= r10;
578    }
579
580    /* Print exponent information */
581    if (exponent < 0)
582        printf("e-%i", -exponent);
583    else if (exponent > 0)
584        printf("e+%i", exponent);
585
586    printf("\n");
587}
588
589} /* namespace lol */
590
Note: See TracBrowser for help on using the repository browser.