source: trunk/src/real.cpp @ 998

Last change on this file since 998 was 998, checked in by sam, 9 years ago

core: add real constants for 0, 1, 2 and 10, and reorder static constants
so that their initialisation order ensures they have the proper values. Add
unit tests for these constants.

File size: 18.3 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 has the larger exponent (no need for the mantissa to
120     * be larger, 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    if (bigoff > BIGITS)
133        return *this;
134
135    ret.m_signexp = m_signexp;
136
137    uint32_t carry = 0;
138    for (int i = BIGITS; i--; )
139    {
140        carry += m_mantissa[i];
141        if (i - bigoff >= 0)
142            carry += x.m_mantissa[i - bigoff] >> off;
143
144        if (i - bigoff > 0)
145            carry += (x.m_mantissa[i - bigoff - 1] << (16 - off)) & 0xffffu;
146        else if (i - bigoff == 0)
147            carry += 0x0001u << (16 - off);
148
149        ret.m_mantissa[i] = carry;
150        carry >>= 16;
151    }
152
153    /* Renormalise in case we overflowed the mantissa */
154    if (carry)
155    {
156        carry--;
157        for (int i = 0; i < BIGITS; i++)
158        {
159            uint16_t tmp = ret.m_mantissa[i];
160            ret.m_mantissa[i] = (carry << 15) | (tmp >> 1);
161            carry = tmp & 0x0001u;
162        }
163        ret.m_signexp++;
164    }
165
166    return ret;
167}
168
169real real::operator -(real const &x) const
170{
171    if (x.m_signexp << 1 == 0)
172        return *this;
173
174    /* Ensure both arguments are positive. Otherwise, switch signs,
175     * or replace - with +. */
176    if (m_signexp >> 31)
177        return -(-*this + x);
178
179    if (x.m_signexp >> 31)
180        return (*this) + (-x);
181
182    /* Ensure *this is larger than x */
183    if (*this < x)
184        return -(x - *this);
185
186    real ret;
187
188    int e1 = m_signexp - (1 << 30) + 1;
189    int e2 = x.m_signexp - (1 << 30) + 1;
190
191    int bigoff = (e1 - e2) / (sizeof(uint16_t) * 8);
192    int off = e1 - e2 - bigoff * (sizeof(uint16_t) * 8);
193
194    if (bigoff > BIGITS)
195        return *this;
196
197    ret.m_signexp = m_signexp;
198
199    int32_t carry = 0;
200    for (int i = 0; i < bigoff; i++)
201    {
202        carry -= x.m_mantissa[BIGITS - i];
203        carry = (carry & 0xffff0000u) | (carry >> 16);
204    }
205    carry -= x.m_mantissa[BIGITS - 1 - bigoff] & ((1 << off) - 1);
206    carry /= (1 << off);
207
208    for (int i = BIGITS; i--; )
209    {
210        carry += m_mantissa[i];
211        if (i - bigoff >= 0)
212            carry -= x.m_mantissa[i - bigoff] >> off;
213
214        if (i - bigoff > 0)
215            carry -= (x.m_mantissa[i - bigoff - 1] << (16 - off)) & 0xffffu;
216        else if (i - bigoff == 0)
217            carry -= 0x0001u << (16 - off);
218
219        ret.m_mantissa[i] = carry;
220        carry = (carry & 0xffff0000u) | (carry >> 16);
221    }
222
223    carry += 1;
224
225    /* Renormalise if we underflowed the mantissa */
226    if (carry == 0)
227    {
228        /* How much do we need to shift the mantissa? FIXME: this could
229         * be computed above */
230        off = 0;
231        for (int i = 0; i < BIGITS; i++)
232        {
233            if (!ret.m_mantissa[i])
234            {
235                off += sizeof(uint16_t) * 8;
236                continue;
237            }
238
239            for (uint16_t tmp = ret.m_mantissa[i]; tmp < 0x8000u; tmp <<= 1)
240                off++;
241            break;
242        }
243        if (off == BIGITS * sizeof(uint16_t) * 8)
244            ret.m_signexp &= 0x80000000u;
245        else
246        {
247            off++; /* Shift one more to get rid of the leading one */
248            ret.m_signexp -= off;
249
250            bigoff = off / (sizeof(uint16_t) * 8);
251            off -= bigoff * sizeof(uint16_t) * 8;
252
253            for (int i = 0; i < BIGITS; i++)
254            {
255                uint16_t tmp = 0;
256                if (i + bigoff < BIGITS)
257                    tmp |= ret.m_mantissa[i + bigoff] << off;
258                if (i + bigoff + 1 < BIGITS)
259                    tmp |= ret.m_mantissa[i + bigoff + 1] >> (16 - off);
260                ret.m_mantissa[i] = tmp;
261            }
262        }
263    }
264
265    return ret;
266}
267
268real real::operator *(real const &x) const
269{
270    real ret;
271
272    if (m_signexp << 1 == 0 || x.m_signexp << 1 == 0)
273    {
274        ret = (m_signexp << 1 == 0) ? *this : x;
275        ret.m_signexp ^= x.m_signexp & 0x80000000u;
276        return ret;
277    }
278
279    ret.m_signexp = (m_signexp ^ x.m_signexp) & 0x80000000u;
280    int e = (m_signexp & 0x7fffffffu) - (1 << 30) + 1
281          + (x.m_signexp & 0x7fffffffu) - (1 << 30) + 1;
282
283    /* Accumulate low order product; no need to store it, we just
284     * want the carry value */
285    uint64_t carry = 0;
286    for (int i = 0; i < BIGITS; i++)
287    {
288        for (int j = 0; j < i + 1; j++)
289            carry += (uint32_t)m_mantissa[BIGITS - 1 - j]
290                   * (uint32_t)x.m_mantissa[BIGITS - 1 + j - i];
291        carry >>= 16;
292    }
293
294    for (int i = 0; i < BIGITS; i++)
295    {
296        for (int j = i + 1; j < BIGITS; j++)
297            carry += (uint32_t)m_mantissa[BIGITS - 1 - j]
298                   * (uint32_t)x.m_mantissa[j - 1 - i];
299
300        carry += m_mantissa[BIGITS - 1 - i];
301        carry += x.m_mantissa[BIGITS - 1 - i];
302        ret.m_mantissa[BIGITS - 1 - i] = carry & 0xffffu;
303        carry >>= 16;
304    }
305
306    /* Renormalise in case we overflowed the mantissa */
307    if (carry)
308    {
309        carry--;
310        for (int i = 0; i < BIGITS; i++)
311        {
312            uint16_t tmp = ret.m_mantissa[i];
313            ret.m_mantissa[i] = (carry << 15) | (tmp >> 1);
314            carry = tmp & 0x0001u;
315        }
316        e++;
317    }
318
319    ret.m_signexp |= e + (1 << 30) - 1;
320
321    return ret;
322}
323
324real real::operator /(real const &x) const
325{
326    return *this * re(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
347real &real::operator /=(real const &x)
348{
349    real tmp = *this;
350    return *this = tmp / x;
351}
352
353real real::operator <<(int x) const
354{
355    real tmp = *this;
356    return tmp <<= x;
357}
358
359real real::operator >>(int x) const
360{
361    real tmp = *this;
362    return tmp >>= x;
363}
364
365real &real::operator <<=(int x)
366{
367    if (m_signexp << 1)
368        m_signexp += x;
369    return *this;
370}
371
372real &real::operator >>=(int x)
373{
374    if (m_signexp << 1)
375        m_signexp -= x;
376    return *this;
377}
378
379bool real::operator ==(real const &x) const
380{
381    if ((m_signexp << 1) == 0 && (x.m_signexp << 1) == 0)
382        return true;
383
384    if (m_signexp != x.m_signexp)
385        return false;
386
387    return memcmp(m_mantissa, x.m_mantissa, sizeof(m_mantissa)) == 0;
388}
389
390bool real::operator !=(real const &x) const
391{
392    return !(*this == x);
393}
394
395bool real::operator <(real const &x) const
396{
397    /* Ensure both numbers are positive */
398    if (m_signexp >> 31)
399        return (x.m_signexp >> 31) ? -*this > -x : true;
400
401    if (x.m_signexp >> 31)
402        return false;
403
404    /* Compare all relevant bits */
405    if (m_signexp != x.m_signexp)
406        return m_signexp < x.m_signexp;
407
408    for (int i = 0; i < BIGITS; i++)
409        if (m_mantissa[i] != x.m_mantissa[i])
410            return m_mantissa[i] < x.m_mantissa[i];
411
412    return false;
413}
414
415bool real::operator <=(real const &x) const
416{
417    return !(*this > x);
418}
419
420bool real::operator >(real const &x) const
421{
422    /* Ensure both numbers are positive */
423    if (m_signexp >> 31)
424        return (x.m_signexp >> 31) ? -*this < -x : false;
425
426    if (x.m_signexp >> 31)
427        return true;
428
429    /* Compare all relevant bits */
430    if (m_signexp != x.m_signexp)
431        return m_signexp > x.m_signexp;
432
433    for (int i = 0; i < BIGITS; i++)
434        if (m_mantissa[i] != x.m_mantissa[i])
435            return m_mantissa[i] > x.m_mantissa[i];
436
437    return false;
438}
439
440bool real::operator >=(real const &x) const
441{
442    return !(*this < x);
443}
444
445bool real::operator !() const
446{
447    return !(bool)*this;
448}
449
450real::operator bool() const
451{
452    /* A real is "true" if it is non-zero (exponent is non-zero) AND
453     * not NaN (exponent is not full bits OR higher order mantissa is zero) */
454    uint32_t exponent = m_signexp << 1;
455    return exponent && (~exponent || m_mantissa[0] == 0);
456}
457
458real re(real const &x)
459{
460    if (!(x.m_signexp << 1))
461    {
462        real ret = x;
463        ret.m_signexp = x.m_signexp | 0x7fffffffu;
464        ret.m_mantissa[0] = 0;
465        return ret;
466    }
467
468    /* Use the system's float inversion to approximate 1/x */
469    union { float f; uint32_t x; } u = { 1.0f }, v = { 1.0f };
470    v.x |= (uint32_t)x.m_mantissa[0] << 7;
471    v.x |= (uint32_t)x.m_mantissa[1] >> 9;
472    v.f = 1.0 / v.f;
473
474    real ret;
475    ret.m_mantissa[0] = (v.x >> 7) & 0xffffu;
476    ret.m_mantissa[1] = (v.x << 9) & 0xffffu;
477
478    uint32_t sign = x.m_signexp & 0x80000000u;
479    ret.m_signexp = sign;
480
481    int exponent = (x.m_signexp & 0x7fffffffu) + 1;
482    exponent = -exponent + (v.x >> 23) - (u.x >> 23);
483    ret.m_signexp |= (exponent - 1) & 0x7fffffffu;
484
485    /* Five steps of Newton-Raphson seems enough for 32-bigit reals. */
486    real two = 2;
487    ret = ret * (two - ret * x);
488    ret = ret * (two - ret * x);
489    ret = ret * (two - ret * x);
490    ret = ret * (two - ret * x);
491    ret = ret * (two - ret * x);
492
493    return ret;
494}
495
496real sqrt(real const &x)
497{
498    /* if zero, return x */
499    if (!(x.m_signexp << 1))
500        return x;
501
502    /* if negative, return NaN */
503    if (x.m_signexp >> 31)
504    {
505        real ret;
506        ret.m_signexp = 0x7fffffffu;
507        ret.m_mantissa[0] = 0xffffu;
508        return ret;
509    }
510
511    /* Use the system's float inversion to approximate 1/sqrt(x). First
512     * we construct a float in the [1..4[ range that has roughly the same
513     * mantissa as our real. Its exponent is 0 or 1, depending on the
514     * partity of x. The final exponent is 0, -1 or -2. We use the final
515     * exponent and final mantissa to pre-fill the result. */
516    union { float f; uint32_t x; } u = { 1.0f }, v = { 2.0f };
517    v.x -= ((x.m_signexp & 1) << 23);
518    v.x |= (uint32_t)x.m_mantissa[0] << 7;
519    v.x |= (uint32_t)x.m_mantissa[1] >> 9;
520    v.f = 1.0 / sqrtf(v.f);
521
522    real ret;
523    ret.m_mantissa[0] = (v.x >> 7) & 0xffffu;
524    ret.m_mantissa[1] = (v.x << 9) & 0xffffu;
525
526    uint32_t sign = x.m_signexp & 0x80000000u;
527    ret.m_signexp = sign;
528
529    int exponent = (x.m_signexp & 0x7fffffffu) - ((1 << 30) - 1);
530    exponent = - (exponent / 2) + (v.x >> 23) - (u.x >> 23);
531    ret.m_signexp |= (exponent + ((1 << 30) - 1)) & 0x7fffffffu;
532
533    /* Five steps of Newton-Raphson seems enough for 32-bigit reals. */
534    real three = 3;
535    ret = ret * (three - ret * ret * x);
536    ret.m_signexp--;
537    ret = ret * (three - ret * ret * x);
538    ret.m_signexp--;
539    ret = ret * (three - ret * ret * x);
540    ret.m_signexp--;
541    ret = ret * (three - ret * ret * x);
542    ret.m_signexp--;
543    ret = ret * (three - ret * ret * x);
544    ret.m_signexp--;
545
546    return ret * x;
547}
548
549real fabs(real const &x)
550{
551    real ret = x;
552    ret.m_signexp &= 0x7fffffffu;
553    return ret;
554}
555
556static real fastlog(real const &x)
557{
558    /* This fast log method is tuned to work on the [1..2] range and
559     * no effort whatsoever was made to improve convergence outside this
560     * domain of validity. It can converge pretty fast, provided we use
561     * the following variable substitutions:
562     *    y = sqrt(x)
563     *    z = (y - 1) / (y + 1)
564     *
565     * And the following identities:
566     *    ln(x) = 2 ln(y)
567     *          = 2 ln((1 + z) / (1 - z))
568     *          = 4 z (1 + z^2 / 3 + z^4 / 5 + z^6 / 7...)
569     *
570     * Any additional sqrt() call would halve the convergence time, but
571     * would also impact the final precision. For now we stick with one
572     * sqrt() call. */
573    real y = sqrt(x);
574    real z = (y - (real)1) / (y + (real)1), z2 = z * z, zn = z2;
575    real sum = 1.0;
576
577    for (int i = 3; i < 200; i += 2)
578    {
579        sum += zn / (real)i;
580        zn *= z2;
581    }
582
583    return z * (sum << 2);
584}
585
586static real LOG_2 = fastlog((real)2);
587
588real log(real const &x)
589{
590    /* Strategy for log(x): if x = 2^E*M then log(x) = E log(2) + log(M),
591     * with the property that M is in [1..2[, so fastlog() applies here. */
592    real tmp = x;
593    if (x.m_signexp >> 31 || x.m_signexp == 0)
594    {
595        tmp.m_signexp = 0xffffffffu;
596        tmp.m_mantissa[0] = 0xffffu;
597        return tmp;
598    }
599    tmp.m_signexp = (1 << 30) - 1;
600    return (real)(x.m_signexp - (1 << 30) + 1) * LOG_2 + fastlog(tmp);
601}
602
603real exp(real const &x)
604{
605    /* Strategy for exp(x): the Taylor series does not converge very fast
606     * with large positive or negative values.
607     *
608     * However, we know that the result is going to be in the form M*2^E,
609     * where M is the mantissa and E the exponent. We first try to predict
610     * a value for E, which is approximately log2(exp(x)) = x / log(2).
611     *
612     * Let E0 be an integer close to x / log(2). We need to find a value x0
613     * such that exp(x) = 2^E0 * exp(x0). We get x0 = x - E0 log(2).
614     *
615     * Thus the final algorithm:
616     *  int E0 = x / log(2)
617     *  real x0 = x - E0 log(2)
618     *  real x1 = exp(x0)
619     *  return x1 * 2^E0
620     */
621    int e0 = x / LOG_2;
622    real x0 = x - (real)e0 * LOG_2;
623    real x1 = 1.0, fact = 1.0, xn = x0;
624
625    for (int i = 1; i < 100; i++)
626    {
627        fact *= (real)i;
628        x1 += xn / fact;
629        xn *= x0;
630    }
631
632    x1.m_signexp += e0;
633    return x1;
634}
635
636real sin(real const &x)
637{
638    real ret = 0.0, fact = 1.0, xn = x, x2 = x * x;
639
640    for (int i = 1; ; i += 2)
641    {
642        real newret = ret + xn / fact;
643        if (ret == newret)
644            break;
645        ret = newret;
646        xn *= x2;
647        fact *= (real)(-(i + 1) * (i + 2));
648    }
649
650    return ret;
651}
652
653real cos(real const &x)
654{
655    real ret = 0.0, fact = 1.0, xn = 1.0, x2 = x * x;
656
657    for (int i = 1; ; i += 2)
658    {
659        real newret = ret + xn / fact;
660        if (ret == newret)
661            break;
662        ret = newret;
663        xn *= x2;
664        fact *= (real)(-i * (i + 1));
665    }
666
667    return ret;
668}
669
670void real::print(int ndigits) const
671{
672    real const r1 = 1, r10 = 10;
673    real x = *this;
674
675    if (x.m_signexp >> 31)
676    {
677        printf("-");
678        x = -x;
679    }
680
681    /* Normalise x so that mantissa is in [1..9.999] */
682    int exponent = 0;
683    if (x.m_signexp)
684    {
685        for (real div = r1, newdiv; true; div = newdiv)
686        {
687            newdiv = div * r10;
688            if (x < newdiv)
689            {
690                x /= div;
691                break;
692            }
693            exponent++;
694        }
695        for (real mul = 1, newx; true; mul *= r10)
696        {
697            newx = x * mul;
698            if (newx >= r1)
699            {
700                x = newx;
701                break;
702            }
703            exponent--;
704        }
705    }
706
707    /* Print digits */
708    for (int i = 0; i < ndigits; i++)
709    {
710        int digit = (int)x;
711        printf("%i", digit);
712        if (i == 0)
713            printf(".");
714        x -= real(digit);
715        x *= r10;
716    }
717
718    /* Print exponent information */
719    if (exponent < 0)
720        printf("e-%i", -exponent);
721    else if (exponent > 0)
722        printf("e+%i", exponent);
723
724    printf("\n");
725}
726
727static real fast_pi()
728{
729    /* Approximate Pi using Machin's formula: 16*atan(1/5)-4*atan(1/239) */
730    real ret = 0.0, x0 = 5.0, x1 = 239.0;
731    real const m0 = -x0 * x0, m1 = -x1 * x1, r16 = 16.0, r4 = 4.0;
732
733    /* Degree 240 is required for 512-bit mantissa precision */
734    for (int i = 1; i < 240; i += 2)
735    {
736        ret += r16 / (x0 * (real)i) - r4 / (x1 * (real)i);
737        x0 *= m0;
738        x1 *= m1;
739    }
740
741    return ret;
742}
743
744real const real::R_0        = (real)0.0;
745real const real::R_1        = (real)1.0;
746real const real::R_2        = (real)2.0;
747real const real::R_10       = (real)10.0;
748
749real const real::R_E        = exp(R_1);
750real const real::R_LN2      = log(R_2);
751real const real::R_LN10     = log(R_10);
752real const real::R_LOG2E    = re(R_LN2);
753real const real::R_LOG10E   = re(R_LN10);
754real const real::R_PI       = fast_pi();
755real const real::R_PI_2     = R_PI >> 1;
756real const real::R_PI_4     = R_PI >> 2;
757real const real::R_1_PI     = re(R_PI);
758real const real::R_2_PI     = R_1_PI << 1;
759real const real::R_2_SQRTPI = re(sqrt(R_PI)) << 1;
760real const real::R_SQRT2    = sqrt(R_2);
761real const real::R_SQRT1_2  = R_SQRT2 >> 1;
762
763} /* namespace lol */
764
Note: See TracBrowser for help on using the repository browser.