source: trunk/src/math/real.cpp @ 1893

Last change on this file since 1893 was 1893, checked in by sam, 7 years ago

math: refactor real number constant declarations so that they are only
computed on demand with static initialisation.

  • Property svn:keywords set to Id
File size: 35.9 KB
Line 
1//
2// Lol Engine
3//
4// Copyright: (c) 2010-2012 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#if defined _XBOX
16#   define _USE_MATH_DEFINES /* for M_PI */
17#   include <xtl.h>
18#   undef near /* Fuck Microsoft */
19#   undef far /* Fuck Microsoft again */
20#elif defined WIN32
21#   define _USE_MATH_DEFINES /* for M_PI */
22#   define WIN32_LEAN_AND_MEAN
23#   include <windows.h>
24#   undef near /* Fuck Microsoft */
25#   undef far /* Fuck Microsoft again */
26#endif
27
28#include <new>
29#include <cstring>
30#include <cstdio>
31#include <cstdlib>
32
33#include "core.h"
34
35using namespace std;
36
37namespace lol
38{
39
40/*
41 * First handle explicit specialisation of our templates.
42 *
43 * Initialisation order is not important because everything is
44 * done on demand, but here is the dependency list anyway:
45 *  - fast_log() requires R_1
46 *  - log() requires R_LN2
47 *  - re() require R_2
48 *  - exp() requires R_0, R_1, R_LN2
49 *  - sqrt() requires R_3
50 */
51
52static real fast_log(real const &x);
53static real fast_pi();
54
55#define LOL_CONSTANT_GETTER(name, value) \
56    template<> real const& real::name() \
57    { \
58        static real const ret = value; \
59        return ret; \
60    }
61
62LOL_CONSTANT_GETTER(R_0,        (real)0.0);
63LOL_CONSTANT_GETTER(R_1,        (real)1.0);
64LOL_CONSTANT_GETTER(R_2,        (real)2.0);
65LOL_CONSTANT_GETTER(R_3,        (real)3.0);
66LOL_CONSTANT_GETTER(R_10,       (real)10.0);
67
68LOL_CONSTANT_GETTER(R_LN2,      fast_log(R_2()));
69LOL_CONSTANT_GETTER(R_LN10,     log(R_10()));
70LOL_CONSTANT_GETTER(R_LOG2E,    re(R_LN2()));
71LOL_CONSTANT_GETTER(R_LOG10E,   re(R_LN10()));
72LOL_CONSTANT_GETTER(R_E,        exp(R_1()));
73LOL_CONSTANT_GETTER(R_PI,       fast_pi());
74LOL_CONSTANT_GETTER(R_PI_2,     R_PI() / 2);
75LOL_CONSTANT_GETTER(R_PI_3,     R_PI() / R_3());
76LOL_CONSTANT_GETTER(R_PI_4,     R_PI() / 4);
77LOL_CONSTANT_GETTER(R_1_PI,     re(R_PI()));
78LOL_CONSTANT_GETTER(R_2_PI,     R_1_PI() * 2);
79LOL_CONSTANT_GETTER(R_2_SQRTPI, re(sqrt(R_PI())) * 2);
80LOL_CONSTANT_GETTER(R_SQRT2,    sqrt(R_2()));
81LOL_CONSTANT_GETTER(R_SQRT3,    sqrt(R_3()));
82LOL_CONSTANT_GETTER(R_SQRT1_2,  R_SQRT2() / 2);
83
84#undef LOL_CONSTANT_GETTER
85
86/*
87 * Now carry on with the rest of the Real class.
88 */
89
90template<> real::Real()
91{
92    m_mantissa = new uint32_t[BIGITS];
93    memset(m_mantissa, 0, BIGITS * sizeof(uint32_t));
94    m_signexp = 0;
95}
96
97template<> real::Real(real const &x)
98{
99    m_mantissa = new uint32_t[BIGITS];
100    memcpy(m_mantissa, x.m_mantissa, BIGITS * sizeof(uint32_t));
101    m_signexp = x.m_signexp;
102}
103
104template<> real const &real::operator =(real const &x)
105{
106    if (&x != this)
107    {
108        memcpy(m_mantissa, x.m_mantissa, BIGITS * sizeof(uint32_t));
109        m_signexp = x.m_signexp;
110    }
111
112    return *this;
113}
114
115template<> real::~Real()
116{
117    delete[] m_mantissa;
118}
119
120template<> real::Real(float f) { new(this) real((double)f); }
121template<> real::Real(int i) { new(this) real((double)i); }
122template<> real::Real(unsigned int i) { new(this) real((double)i); }
123
124template<> real::Real(double d)
125{
126    new(this) real();
127
128    union { double d; uint64_t x; } u = { d };
129
130    uint32_t sign = (u.x >> 63) << 31;
131    uint32_t exponent = (u.x << 1) >> 53;
132
133    switch (exponent)
134    {
135    case 0x00:
136        m_signexp = sign;
137        break;
138    case 0x7ff:
139        m_signexp = sign | 0x7fffffffu;
140        break;
141    default:
142        m_signexp = sign | (exponent + (1 << 30) - (1 << 10));
143        break;
144    }
145
146    m_mantissa[0] = (uint32_t)(u.x >> 20);
147    m_mantissa[1] = (uint32_t)(u.x << 12);
148    memset(m_mantissa + 2, 0, (BIGITS - 2) * sizeof(m_mantissa[0]));
149}
150
151template<> real::operator float() const { return (float)(double)(*this); }
152template<> real::operator int() const { return (int)(double)(*this); }
153template<> real::operator unsigned() const { return (unsigned)(double)(*this); }
154
155template<> real::operator double() const
156{
157    union { double d; uint64_t x; } u;
158
159    /* Get sign */
160    u.x = m_signexp >> 31;
161    u.x <<= 11;
162
163    /* Compute new exponent */
164    uint32_t exponent = (m_signexp << 1) >> 1;
165    int e = (int)exponent - (1 << 30) + (1 << 10);
166
167    if (e < 0)
168        u.x <<= 52;
169    else if (e >= 0x7ff)
170    {
171        u.x |= 0x7ff;
172        u.x <<= 52;
173    }
174    else
175    {
176        u.x |= e;
177
178        /* Store mantissa if necessary */
179        u.x <<= 32;
180        u.x |= m_mantissa[0];
181        u.x <<= 20;
182        u.x |= m_mantissa[1] >> 12;
183        /* Rounding */
184        u.x += (m_mantissa[1] >> 11) & 1;
185    }
186
187    return u.d;
188}
189
190/*
191 * Create a real number from an ASCII representation
192 */
193template<> real::Real(char const *str)
194{
195    real ret = 0;
196    int exponent = 0;
197    bool comma = false, nonzero = false, negative = false, finished = false;
198
199    for (char const *p = str; *p && !finished; p++)
200    {
201        switch (*p)
202        {
203        case '-':
204        case '+':
205            if (p != str)
206                break;
207            negative = (*p == '-');
208            break;
209        case '.':
210            if (comma)
211                finished = true;
212            comma = true;
213            break;
214        case '0': case '1': case '2': case '3': case '4':
215        case '5': case '6': case '7': case '8': case '9':
216            if (nonzero)
217            {
218                real x = ret + ret;
219                x = x + x + ret;
220                ret = x + x;
221            }
222            if (*p != '0')
223            {
224                ret += (int)(*p - '0');
225                nonzero = true;
226            }
227            if (comma)
228                exponent--;
229            break;
230        case 'e':
231        case 'E':
232            exponent += atoi(p + 1);
233            finished = true;
234            break;
235        default:
236            finished = true;
237            break;
238        }
239    }
240
241    if (exponent)
242        ret *= pow(R_10(), (real)exponent);
243
244    if (negative)
245        ret = -ret;
246
247    new(this) real(ret);
248}
249
250template<> real real::operator +() const
251{
252    return *this;
253}
254
255template<> real real::operator -() const
256{
257    real ret = *this;
258    ret.m_signexp ^= 0x80000000u;
259    return ret;
260}
261
262template<> real real::operator +(real const &x) const
263{
264    if (x.m_signexp << 1 == 0)
265        return *this;
266
267    /* Ensure both arguments are positive. Otherwise, switch signs,
268     * or replace + with -. */
269    if (m_signexp >> 31)
270        return -(-*this + -x);
271
272    if (x.m_signexp >> 31)
273        return *this - (-x);
274
275    /* Ensure *this has the larger exponent (no need for the mantissa to
276     * be larger, as in subtraction). Otherwise, switch. */
277    if ((m_signexp << 1) < (x.m_signexp << 1))
278        return x + *this;
279
280    real ret;
281
282    int e1 = m_signexp - (1 << 30) + 1;
283    int e2 = x.m_signexp - (1 << 30) + 1;
284
285    int bigoff = (e1 - e2) / BIGIT_BITS;
286    int off = e1 - e2 - bigoff * BIGIT_BITS;
287
288    if (bigoff > BIGITS)
289        return *this;
290
291    ret.m_signexp = m_signexp;
292
293    uint64_t carry = 0;
294    for (int i = BIGITS; i--; )
295    {
296        carry += m_mantissa[i];
297        if (i - bigoff >= 0)
298            carry += x.m_mantissa[i - bigoff] >> off;
299
300        if (off && i - bigoff > 0)
301            carry += (x.m_mantissa[i - bigoff - 1] << (BIGIT_BITS - off)) & 0xffffffffu;
302        else if (i - bigoff == 0)
303            carry += (uint64_t)1 << (BIGIT_BITS - off);
304
305        ret.m_mantissa[i] = (uint32_t)carry;
306        carry >>= BIGIT_BITS;
307    }
308
309    /* Renormalise in case we overflowed the mantissa */
310    if (carry)
311    {
312        carry--;
313        for (int i = 0; i < BIGITS; i++)
314        {
315            uint32_t tmp = ret.m_mantissa[i];
316            ret.m_mantissa[i] = ((uint32_t)carry << (BIGIT_BITS - 1))
317                              | (tmp >> 1);
318            carry = tmp & 1u;
319        }
320        ret.m_signexp++;
321    }
322
323    return ret;
324}
325
326template<> real real::operator -(real const &x) const
327{
328    if (x.m_signexp << 1 == 0)
329        return *this;
330
331    /* Ensure both arguments are positive. Otherwise, switch signs,
332     * or replace - with +. */
333    if (m_signexp >> 31)
334        return -(-*this + x);
335
336    if (x.m_signexp >> 31)
337        return (*this) + (-x);
338
339    /* Ensure *this is larger than x */
340    if (*this < x)
341        return -(x - *this);
342
343    real ret;
344
345    int e1 = m_signexp - (1 << 30) + 1;
346    int e2 = x.m_signexp - (1 << 30) + 1;
347
348    int bigoff = (e1 - e2) / BIGIT_BITS;
349    int off = e1 - e2 - bigoff * BIGIT_BITS;
350
351    if (bigoff > BIGITS)
352        return *this;
353
354    ret.m_signexp = m_signexp;
355
356    /* int64_t instead of uint64_t to preserve sign through shifts */
357    int64_t carry = 0;
358    for (int i = 0; i < bigoff; i++)
359    {
360        carry -= x.m_mantissa[BIGITS - 1 - i];
361        /* Emulates a signed shift */
362        carry >>= BIGIT_BITS;
363        carry |= carry << BIGIT_BITS;
364    }
365    if (bigoff < BIGITS)
366        carry -= x.m_mantissa[BIGITS - 1 - bigoff] & (((int64_t)1 << off) - 1);
367    carry /= (int64_t)1 << off;
368
369    for (int i = BIGITS; i--; )
370    {
371        carry += m_mantissa[i];
372        if (i - bigoff >= 0)
373            carry -= x.m_mantissa[i - bigoff] >> off;
374
375        if (off && i - bigoff > 0)
376            carry -= (x.m_mantissa[i - bigoff - 1] << (BIGIT_BITS - off)) & 0xffffffffu;
377        else if (i - bigoff == 0)
378            carry -= (uint64_t)1 << (BIGIT_BITS - off);
379
380        ret.m_mantissa[i] = (uint32_t)carry;
381        carry >>= BIGIT_BITS;
382        carry |= carry << BIGIT_BITS;
383    }
384
385    carry += 1;
386
387    /* Renormalise if we underflowed the mantissa */
388    if (carry == 0)
389    {
390        /* How much do we need to shift the mantissa? FIXME: this could
391         * be computed above */
392        off = 0;
393        for (int i = 0; i < BIGITS; i++)
394        {
395            if (!ret.m_mantissa[i])
396            {
397                off += BIGIT_BITS;
398                continue;
399            }
400
401            for (uint32_t tmp = ret.m_mantissa[i]; tmp < 0x80000000u; tmp <<= 1)
402                off++;
403            break;
404        }
405        if (off == BIGITS * BIGIT_BITS)
406            ret.m_signexp &= 0x80000000u;
407        else
408        {
409            off++; /* Shift one more to get rid of the leading one */
410            ret.m_signexp -= off;
411
412            bigoff = off / BIGIT_BITS;
413            off -= bigoff * BIGIT_BITS;
414
415            for (int i = 0; i < BIGITS; i++)
416            {
417                uint32_t tmp = 0;
418                if (i + bigoff < BIGITS)
419                    tmp |= ret.m_mantissa[i + bigoff] << off;
420                if (off && i + bigoff + 1 < BIGITS)
421                    tmp |= ret.m_mantissa[i + bigoff + 1] >> (BIGIT_BITS - off);
422                ret.m_mantissa[i] = tmp;
423            }
424        }
425    }
426
427    return ret;
428}
429
430template<> real real::operator *(real const &x) const
431{
432    real ret;
433
434    if (m_signexp << 1 == 0 || x.m_signexp << 1 == 0)
435    {
436        ret = (m_signexp << 1 == 0) ? *this : x;
437        ret.m_signexp ^= x.m_signexp & 0x80000000u;
438        return ret;
439    }
440
441    ret.m_signexp = (m_signexp ^ x.m_signexp) & 0x80000000u;
442    int e = (m_signexp & 0x7fffffffu) - (1 << 30) + 1
443          + (x.m_signexp & 0x7fffffffu) - (1 << 30) + 1;
444
445    /* Accumulate low order product; no need to store it, we just
446     * want the carry value */
447    uint64_t carry = 0, hicarry = 0, prev;
448    for (int i = 0; i < BIGITS; i++)
449    {
450        for (int j = 0; j < i + 1; j++)
451        {
452            prev = carry;
453            carry += (uint64_t)m_mantissa[BIGITS - 1 - j]
454                   * (uint64_t)x.m_mantissa[BIGITS - 1 + j - i];
455            if (carry < prev)
456                hicarry++;
457        }
458        carry >>= BIGIT_BITS;
459        carry |= hicarry << BIGIT_BITS;
460        hicarry >>= BIGIT_BITS;
461    }
462
463    for (int i = 0; i < BIGITS; i++)
464    {
465        for (int j = i + 1; j < BIGITS; j++)
466        {
467            prev = carry;
468            carry += (uint64_t)m_mantissa[BIGITS - 1 - j]
469                   * (uint64_t)x.m_mantissa[j - 1 - i];
470            if (carry < prev)
471                hicarry++;
472        }
473        prev = carry;
474        carry += m_mantissa[BIGITS - 1 - i];
475        carry += x.m_mantissa[BIGITS - 1 - i];
476        if (carry < prev)
477            hicarry++;
478        ret.m_mantissa[BIGITS - 1 - i] = carry & 0xffffffffu;
479        carry >>= BIGIT_BITS;
480        carry |= hicarry << BIGIT_BITS;
481        hicarry >>= BIGIT_BITS;
482    }
483
484    /* Renormalise in case we overflowed the mantissa */
485    if (carry)
486    {
487        carry--;
488        for (int i = 0; i < BIGITS; i++)
489        {
490            uint32_t tmp = (uint32_t)ret.m_mantissa[i];
491            ret.m_mantissa[i] = ((uint32_t)carry << (BIGIT_BITS - 1))
492                              | (tmp >> 1);
493            carry = tmp & 1u;
494        }
495        e++;
496    }
497
498    ret.m_signexp |= e + (1 << 30) - 1;
499
500    return ret;
501}
502
503template<> real real::operator /(real const &x) const
504{
505    return *this * re(x);
506}
507
508template<> real const &real::operator +=(real const &x)
509{
510    real tmp = *this;
511    return *this = tmp + x;
512}
513
514template<> real const &real::operator -=(real const &x)
515{
516    real tmp = *this;
517    return *this = tmp - x;
518}
519
520template<> real const &real::operator *=(real const &x)
521{
522    real tmp = *this;
523    return *this = tmp * x;
524}
525
526template<> real const &real::operator /=(real const &x)
527{
528    real tmp = *this;
529    return *this = tmp / x;
530}
531
532template<> bool real::operator ==(real const &x) const
533{
534    if ((m_signexp << 1) == 0 && (x.m_signexp << 1) == 0)
535        return true;
536
537    if (m_signexp != x.m_signexp)
538        return false;
539
540    return memcmp(m_mantissa, x.m_mantissa, BIGITS * sizeof(uint32_t)) == 0;
541}
542
543template<> bool real::operator !=(real const &x) const
544{
545    return !(*this == x);
546}
547
548template<> bool real::operator <(real const &x) const
549{
550    /* Ensure both numbers are positive */
551    if (m_signexp >> 31)
552        return (x.m_signexp >> 31) ? -*this > -x : true;
553
554    if (x.m_signexp >> 31)
555        return false;
556
557    /* Compare all relevant bits */
558    if (m_signexp != x.m_signexp)
559        return m_signexp < x.m_signexp;
560
561    for (int i = 0; i < BIGITS; i++)
562        if (m_mantissa[i] != x.m_mantissa[i])
563            return m_mantissa[i] < x.m_mantissa[i];
564
565    return false;
566}
567
568template<> bool real::operator <=(real const &x) const
569{
570    return !(*this > x);
571}
572
573template<> bool real::operator >(real const &x) const
574{
575    /* Ensure both numbers are positive */
576    if (m_signexp >> 31)
577        return (x.m_signexp >> 31) ? -*this < -x : false;
578
579    if (x.m_signexp >> 31)
580        return true;
581
582    /* Compare all relevant bits */
583    if (m_signexp != x.m_signexp)
584        return m_signexp > x.m_signexp;
585
586    for (int i = 0; i < BIGITS; i++)
587        if (m_mantissa[i] != x.m_mantissa[i])
588            return m_mantissa[i] > x.m_mantissa[i];
589
590    return false;
591}
592
593template<> bool real::operator >=(real const &x) const
594{
595    return !(*this < x);
596}
597
598template<> bool real::operator !() const
599{
600    return !(bool)*this;
601}
602
603template<> real::operator bool() const
604{
605    /* A real is "true" if it is non-zero (exponent is non-zero) AND
606     * not NaN (exponent is not full bits OR higher order mantissa is zero) */
607    uint32_t exponent = m_signexp << 1;
608    return exponent && (~exponent || m_mantissa[0] == 0);
609}
610
611template<> real min(real const &a, real const &b)
612{
613    return (a < b) ? a : b;
614}
615
616template<> real max(real const &a, real const &b)
617{
618    return (a > b) ? a : b;
619}
620
621template<> real clamp(real const &x, real const &a, real const &b)
622{
623    return (x < a) ? a : (x > b) ? b : x;
624}
625
626template<> real re(real const &x)
627{
628    if (!(x.m_signexp << 1))
629    {
630        real ret = x;
631        ret.m_signexp = x.m_signexp | 0x7fffffffu;
632        ret.m_mantissa[0] = 0;
633        return ret;
634    }
635
636    /* Use the system's float inversion to approximate 1/x */
637    union { float f; uint32_t x; } u = { 1.0f }, v = { 1.0f };
638    v.x |= x.m_mantissa[0] >> 9;
639    v.f = 1.0f / v.f;
640
641    real ret;
642    ret.m_mantissa[0] = v.x << 9;
643
644    uint32_t sign = x.m_signexp & 0x80000000u;
645    ret.m_signexp = sign;
646
647    int exponent = (x.m_signexp & 0x7fffffffu) + 1;
648    exponent = -exponent + (v.x >> 23) - (u.x >> 23);
649    ret.m_signexp |= (exponent - 1) & 0x7fffffffu;
650
651    /* FIXME: 1+log2(BIGITS) steps of Newton-Raphson seems to be enough for
652     * convergence, but this hasn't been checked seriously. */
653    for (int i = 1; i <= real::BIGITS; i *= 2)
654        ret = ret * (real::R_2() - ret * x);
655
656    return ret;
657}
658
659template<> real sqrt(real const &x)
660{
661    /* if zero, return x */
662    if (!(x.m_signexp << 1))
663        return x;
664
665    /* if negative, return NaN */
666    if (x.m_signexp >> 31)
667    {
668        real ret;
669        ret.m_signexp = 0x7fffffffu;
670        ret.m_mantissa[0] = 0xffffu;
671        return ret;
672    }
673
674    /* Use the system's float inversion to approximate 1/sqrt(x). First
675     * we construct a float in the [1..4[ range that has roughly the same
676     * mantissa as our real. Its exponent is 0 or 1, depending on the
677     * partity of x. The final exponent is 0, -1 or -2. We use the final
678     * exponent and final mantissa to pre-fill the result. */
679    union { float f; uint32_t x; } u = { 1.0f }, v = { 2.0f };
680    v.x -= ((x.m_signexp & 1) << 23);
681    v.x |= x.m_mantissa[0] >> 9;
682    v.f = 1.0f / sqrtf(v.f);
683
684    real ret;
685    ret.m_mantissa[0] = v.x << 9;
686
687    uint32_t sign = x.m_signexp & 0x80000000u;
688    ret.m_signexp = sign;
689
690    uint32_t exponent = (x.m_signexp & 0x7fffffffu);
691    exponent = ((1 << 30) + (1 << 29) - 1) - (exponent + 1) / 2;
692    exponent = exponent + (v.x >> 23) - (u.x >> 23);
693    ret.m_signexp |= exponent & 0x7fffffffu;
694
695    /* FIXME: 1+log2(BIGITS) steps of Newton-Raphson seems to be enough for
696     * convergence, but this hasn't been checked seriously. */
697    for (int i = 1; i <= real::BIGITS; i *= 2)
698    {
699        ret = ret * (real::R_3() - ret * ret * x);
700        ret.m_signexp--;
701    }
702
703    return ret * x;
704}
705
706template<> real cbrt(real const &x)
707{
708    /* if zero, return x */
709    if (!(x.m_signexp << 1))
710        return x;
711
712    /* Use the system's float inversion to approximate cbrt(x). First
713     * we construct a float in the [1..8[ range that has roughly the same
714     * mantissa as our real. Its exponent is 0, 1 or 2, depending on the
715     * value of x. The final exponent is 0 or 1 (special case). We use
716     * the final exponent and final mantissa to pre-fill the result. */
717    union { float f; uint32_t x; } u = { 1.0f }, v = { 1.0f };
718    v.x += ((x.m_signexp % 3) << 23);
719    v.x |= x.m_mantissa[0] >> 9;
720    v.f = powf(v.f, 0.33333333333333333f);
721
722    real ret;
723    ret.m_mantissa[0] = v.x << 9;
724
725    uint32_t sign = x.m_signexp & 0x80000000u;
726    ret.m_signexp = sign;
727
728    int exponent = (x.m_signexp & 0x7fffffffu) - (1 << 30) + 1;
729    exponent = exponent / 3 + (v.x >> 23) - (u.x >> 23);
730    ret.m_signexp |= (exponent + (1 << 30) - 1) & 0x7fffffffu;
731
732    /* FIXME: 1+log2(BIGITS) steps of Newton-Raphson seems to be enough for
733     * convergence, but this hasn't been checked seriously. */
734    for (int i = 1; i <= real::BIGITS; i *= 2)
735    {
736        static real third = re(real::R_3());
737        ret = third * (x / (ret * ret) + (ret / 2));
738    }
739
740    return ret;
741}
742
743template<> real pow(real const &x, real const &y)
744{
745    if (!y)
746        return real::R_1();
747    if (!x)
748        return real::R_0();
749    if (x > real::R_0())
750        return exp(y * log(x));
751    else /* x < 0 */
752    {
753        /* Odd integer exponent */
754        if (y == (round(y / 2) * 2))
755            return exp(y * log(-x));
756
757        /* Even integer exponent */
758        if (y == round(y))
759            return -exp(y * log(-x));
760
761        /* FIXME: negative nth root */
762        return real::R_0();
763    }
764}
765
766static real fast_fact(int x)
767{
768    real ret = real::R_1();
769    int i = 1, multiplier = 1, exponent = 0;
770
771    for (;;)
772    {
773        if (i++ >= x)
774            /* Multiplication is a no-op if multiplier == 1 */
775            return ldexp(ret * multiplier, exponent);
776
777        int tmp = i;
778        while ((tmp & 1) == 0)
779        {
780            tmp >>= 1;
781            exponent++;
782        }
783        if (multiplier * tmp / tmp != multiplier)
784        {
785            ret *= multiplier;
786            multiplier = 1;
787        }
788        multiplier *= tmp;
789    }
790}
791
792template<> real gamma(real const &x)
793{
794    /* We use Spouge's formula. FIXME: precision is far from acceptable,
795     * especially with large values. We need to compute this with higher
796     * precision values in order to attain the desired accuracy. It might
797     * also be useful to sort the ck values by decreasing absolute value
798     * and do the addition in this order. */
799    int a = ceilf(logf(2) / logf(2 * M_PI) * real::BIGITS * real::BIGIT_BITS);
800
801    real ret = sqrt(real::R_PI() * 2);
802    real fact_k_1 = real::R_1();
803
804    for (int k = 1; k < a; k++)
805    {
806        real a_k = (real)(a - k);
807        real ck = pow(a_k, (real)((float)k - 0.5)) * exp(a_k)
808                / (fact_k_1 * (x + (real)(k - 1)));
809        ret += ck;
810        fact_k_1 *= (real)-k;
811    }
812
813    ret *= pow(x + (real)(a - 1), x - (real::R_1() / 2));
814    ret *= exp(-x - (real)(a - 1));
815
816    return ret;
817}
818
819template<> real fabs(real const &x)
820{
821    real ret = x;
822    ret.m_signexp &= 0x7fffffffu;
823    return ret;
824}
825
826static real fast_log(real const &x)
827{
828    /* This fast log method is tuned to work on the [1..2] range and
829     * no effort whatsoever was made to improve convergence outside this
830     * domain of validity. It can converge pretty fast, provided we use
831     * the following variable substitutions:
832     *    y = sqrt(x)
833     *    z = (y - 1) / (y + 1)
834     *
835     * And the following identities:
836     *    ln(x) = 2 ln(y)
837     *          = 2 ln((1 + z) / (1 - z))
838     *          = 4 z (1 + z^2 / 3 + z^4 / 5 + z^6 / 7...)
839     *
840     * Any additional sqrt() call would halve the convergence time, but
841     * would also impact the final precision. For now we stick with one
842     * sqrt() call. */
843    real y = sqrt(x);
844    real z = (y - real::R_1()) / (y + real::R_1()), z2 = z * z, zn = z2;
845    real sum = real::R_1();
846
847    for (int i = 3; ; i += 2)
848    {
849        real newsum = sum + zn / (real)i;
850        if (newsum == sum)
851            break;
852        sum = newsum;
853        zn *= z2;
854    }
855
856    return z * sum * 4;
857}
858
859template<> real log(real const &x)
860{
861    /* Strategy for log(x): if x = 2^E*M then log(x) = E log(2) + log(M),
862     * with the property that M is in [1..2[, so fast_log() applies here. */
863    real tmp = x;
864    if (x.m_signexp >> 31 || x.m_signexp == 0)
865    {
866        tmp.m_signexp = 0xffffffffu;
867        tmp.m_mantissa[0] = 0xffffffffu;
868        return tmp;
869    }
870    tmp.m_signexp = (1 << 30) - 1;
871    return (real)(int)(x.m_signexp - (1 << 30) + 1) * real::R_LN2()
872           + fast_log(tmp);
873}
874
875template<> real log2(real const &x)
876{
877    /* Strategy for log2(x): see log(x). */
878    real tmp = x;
879    if (x.m_signexp >> 31 || x.m_signexp == 0)
880    {
881        tmp.m_signexp = 0xffffffffu;
882        tmp.m_mantissa[0] = 0xffffffffu;
883        return tmp;
884    }
885    tmp.m_signexp = (1 << 30) - 1;
886    return (real)(int)(x.m_signexp - (1 << 30) + 1)
887           + fast_log(tmp) * real::R_LOG2E();
888}
889
890template<> real log10(real const &x)
891{
892    return log(x) * real::R_LOG10E();
893}
894
895static real fast_exp_sub(real const &x, real const &y)
896{
897    /* This fast exp method is tuned to work on the [-1..1] range and
898     * no effort whatsoever was made to improve convergence outside this
899     * domain of validity. The argument y is used for cases where we
900     * don't want the leading 1 in the Taylor series. */
901    real ret = real::R_1() - y, xn = x;
902    int i = 1;
903
904    for (;;)
905    {
906        real newret = ret + xn;
907        if (newret == ret)
908            break;
909        ret = newret * ++i;
910        xn *= x;
911    }
912
913    return ret / fast_fact(i);
914}
915
916template<> real exp(real const &x)
917{
918    /* Strategy for exp(x): the Taylor series does not converge very fast
919     * with large positive or negative values.
920     *
921     * However, we know that the result is going to be in the form M*2^E,
922     * where M is the mantissa and E the exponent. We first try to predict
923     * a value for E, which is approximately log2(exp(x)) = x / log(2).
924     *
925     * Let E0 be an integer close to x / log(2). We need to find a value x0
926     * such that exp(x) = 2^E0 * exp(x0). We get x0 = x - E0 log(2).
927     *
928     * Thus the final algorithm:
929     *  int E0 = x / log(2)
930     *  real x0 = x - E0 log(2)
931     *  real x1 = exp(x0)
932     *  return x1 * 2^E0
933     */
934    int e0 = x / real::R_LN2();
935    real x0 = x - (real)e0 * real::R_LN2();
936    real x1 = fast_exp_sub(x0, real::R_0());
937    x1.m_signexp += e0;
938    return x1;
939}
940
941template<> real exp2(real const &x)
942{
943    /* Strategy for exp2(x): see strategy in exp(). */
944    int e0 = x;
945    real x0 = x - (real)e0;
946    real x1 = fast_exp_sub(x0 * real::R_LN2(), real::R_0());
947    x1.m_signexp += e0;
948    return x1;
949}
950
951template<> real sinh(real const &x)
952{
953    /* We cannot always use (exp(x)-exp(-x))/2 because we'll lose
954     * accuracy near zero. We only use this identity for |x|>0.5. If
955     * |x|<=0.5, we compute exp(x)-1 and exp(-x)-1 instead. */
956    bool near_zero = (fabs(x) < real::R_1() / 2);
957    real x1 = near_zero ? fast_exp_sub(x, real::R_1()) : exp(x);
958    real x2 = near_zero ? fast_exp_sub(-x, real::R_1()) : exp(-x);
959    return (x1 - x2) / 2;
960}
961
962template<> real tanh(real const &x)
963{
964    /* See sinh() for the strategy here */
965    bool near_zero = (fabs(x) < real::R_1() / 2);
966    real x1 = near_zero ? fast_exp_sub(x, real::R_1()) : exp(x);
967    real x2 = near_zero ? fast_exp_sub(-x, real::R_1()) : exp(-x);
968    real x3 = near_zero ? x1 + x2 + real::R_2() : x1 + x2;
969    return (x1 - x2) / x3;
970}
971
972template<> real cosh(real const &x)
973{
974    /* No need to worry about accuracy here; maybe the last bit is slightly
975     * off, but that's about it. */
976    return (exp(x) + exp(-x)) / 2;
977}
978
979template<> real frexp(real const &x, int *exp)
980{
981    if (!x)
982    {
983        *exp = 0;
984        return x;
985    }
986
987    real ret = x;
988    int exponent = (ret.m_signexp & 0x7fffffffu) - (1 << 30) + 1;
989    *exp = exponent + 1;
990    ret.m_signexp -= exponent + 1;
991    return ret;
992}
993
994template<> real ldexp(real const &x, int exp)
995{
996    real ret = x;
997    if (ret)
998        ret.m_signexp += exp;
999    return ret;
1000}
1001
1002template<> real modf(real const &x, real *iptr)
1003{
1004    real absx = fabs(x);
1005    real tmp = floor(absx);
1006
1007    *iptr = copysign(tmp, x);
1008    return copysign(absx - tmp, x);
1009}
1010
1011template<> real ulp(real const &x)
1012{
1013    real ret = real::R_1();
1014    if (x)
1015        ret.m_signexp = x.m_signexp + 1 - real::BIGITS * real::BIGIT_BITS;
1016    else
1017        ret.m_signexp = 0;
1018    return ret;
1019}
1020
1021template<> real nextafter(real const &x, real const &y)
1022{
1023    if (x == y)
1024        return x;
1025    else if (x < y)
1026        return x + ulp(x);
1027    else
1028        return x - ulp(x);
1029}
1030
1031template<> real copysign(real const &x, real const &y)
1032{
1033    real ret = x;
1034    ret.m_signexp &= 0x7fffffffu;
1035    ret.m_signexp |= y.m_signexp & 0x80000000u;
1036    return ret;
1037}
1038
1039template<> real floor(real const &x)
1040{
1041    /* Strategy for floor(x):
1042     *  - if negative, return -ceil(-x)
1043     *  - if zero or negative zero, return x
1044     *  - if less than one, return zero
1045     *  - otherwise, if e is the exponent, clear all bits except the
1046     *    first e. */
1047    if (x < -real::R_0())
1048        return -ceil(-x);
1049    if (!x)
1050        return x;
1051    if (x < real::R_1())
1052        return real::R_0();
1053
1054    real ret = x;
1055    int exponent = x.m_signexp - (1 << 30) + 1;
1056
1057    for (int i = 0; i < real::BIGITS; i++)
1058    {
1059        if (exponent <= 0)
1060            ret.m_mantissa[i] = 0;
1061        else if (exponent < real::BIGIT_BITS)
1062            ret.m_mantissa[i] &= ~((1 << (real::BIGIT_BITS - exponent)) - 1);
1063
1064        exponent -= real::BIGIT_BITS;
1065    }
1066
1067    return ret;
1068}
1069
1070template<> real ceil(real const &x)
1071{
1072    /* Strategy for ceil(x):
1073     *  - if negative, return -floor(-x)
1074     *  - if x == floor(x), return x
1075     *  - otherwise, return floor(x) + 1 */
1076    if (x < -real::R_0())
1077        return -floor(-x);
1078    real ret = floor(x);
1079    if (x == ret)
1080        return ret;
1081    else
1082        return ret + real::R_1();
1083}
1084
1085template<> real round(real const &x)
1086{
1087    if (x < real::R_0())
1088        return -round(-x);
1089
1090    return floor(x + (real::R_1() / 2));
1091}
1092
1093template<> real fmod(real const &x, real const &y)
1094{
1095    if (!y)
1096        return real::R_0(); /* FIXME: return NaN */
1097
1098    if (!x)
1099        return x;
1100
1101    real tmp = round(x / y);
1102    return x - tmp * y;
1103}
1104
1105template<> real sin(real const &x)
1106{
1107    int switch_sign = x.m_signexp & 0x80000000u;
1108
1109    real absx = fmod(fabs(x), real::R_PI() * 2);
1110    if (absx > real::R_PI())
1111    {
1112        absx -= real::R_PI();
1113        switch_sign = !switch_sign;
1114    }
1115
1116    if (absx > real::R_PI_2())
1117        absx = real::R_PI() - absx;
1118
1119    real ret = real::R_0(), fact = real::R_1(), xn = absx, mx2 = -absx * absx;
1120    int i = 1;
1121    for (;;)
1122    {
1123        real newret = ret + xn;
1124        if (newret == ret)
1125            break;
1126        ret = newret * ((i + 1) * (i + 2));
1127        xn *= mx2;
1128        i += 2;
1129    }
1130    ret /= fast_fact(i);
1131
1132    /* Propagate sign */
1133    if (switch_sign)
1134        ret.m_signexp ^= 0x80000000u;
1135    return ret;
1136}
1137
1138template<> real cos(real const &x)
1139{
1140    return sin(real::R_PI_2() - x);
1141}
1142
1143template<> real tan(real const &x)
1144{
1145    /* Constrain input to [-π,π] */
1146    real y = fmod(x, real::R_PI());
1147
1148    /* Constrain input to [-π/2,π/2] */
1149    if (y < -real::R_PI_2())
1150        y += real::R_PI();
1151    else if (y > real::R_PI_2())
1152        y -= real::R_PI();
1153
1154    /* In [-π/4,π/4] return sin/cos */
1155    if (fabs(y) <= real::R_PI_4())
1156        return sin(y) / cos(y);
1157
1158    /* Otherwise, return cos/sin */
1159    if (y > real::R_0())
1160        y = real::R_PI_2() - y;
1161    else
1162        y = -real::R_PI_2() - y;
1163
1164    return cos(y) / sin(y);
1165}
1166
1167static inline real asinacos(real const &x, int is_asin, int is_negative)
1168{
1169    /* Strategy for asin(): in [-0.5..0.5], use a Taylor series around
1170     * zero. In [0.5..1], use asin(x) = π/2 - 2*asin(sqrt((1-x)/2)), and
1171     * in [-1..-0.5] just revert the sign.
1172     * Strategy for acos(): use acos(x) = π/2 - asin(x) and try not to
1173     * lose the precision around x=1. */
1174    real absx = fabs(x);
1175    int around_zero = (absx < (real::R_1() / 2));
1176
1177    if (!around_zero)
1178        absx = sqrt((real::R_1() - absx) / 2);
1179
1180    real ret = absx, xn = absx, x2 = absx * absx, fact1 = 2, fact2 = 1;
1181    for (int i = 1; ; i++)
1182    {
1183        xn *= x2;
1184        real mul = (real)(2 * i + 1);
1185        real newret = ret + ldexp(fact1 * xn / (mul * fact2), -2 * i);
1186        if (newret == ret)
1187            break;
1188        ret = newret;
1189        fact1 *= (real)((2 * i + 1) * (2 * i + 2));
1190        fact2 *= (real)((i + 1) * (i + 1));
1191    }
1192
1193    if (is_negative)
1194        ret = -ret;
1195
1196    if (around_zero)
1197        ret = is_asin ? ret : real::R_PI_2() - ret;
1198    else
1199    {
1200        real adjust = is_negative ? real::R_PI() : real::R_0();
1201        if (is_asin)
1202            ret = real::R_PI_2() - adjust - ret * 2;
1203        else
1204            ret = adjust + ret * 2;
1205    }
1206
1207    return ret;
1208}
1209
1210template<> real asin(real const &x)
1211{
1212    return asinacos(x, 1, x.m_signexp >> 31);
1213}
1214
1215template<> real acos(real const &x)
1216{
1217    return asinacos(x, 0, x.m_signexp >> 31);
1218}
1219
1220template<> real atan(real const &x)
1221{
1222    /* Computing atan(x): we choose a different Taylor series depending on
1223     * the value of x to help with convergence.
1224     *
1225     * If |x| < 0.5 we evaluate atan(y) near 0:
1226     *  atan(y) = y - y^3/3 + y^5/5 - y^7/7 + y^9/9 ...
1227     *
1228     * If 0.5 <= |x| < 1.5 we evaluate atan(1+y) near 0:
1229     *  atan(1+y) = π/4 + y/(1*2^1) - y^2/(2*2^1) + y^3/(3*2^2)
1230     *                  - y^5/(5*2^3) + y^6/(6*2^3) - y^7/(7*2^4)
1231     *                  + y^9/(9*2^5) - y^10/(10*2^5) + y^11/(11*2^6) ...
1232     *
1233     * If 1.5 <= |x| < 2 we evaluate atan(sqrt(3)+2y) near 0:
1234     *  atan(sqrt(3)+2y) = π/3 + 1/2 y - sqrt(3)/2 y^2/2
1235     *                         + y^3/3 - sqrt(3)/2 y^4/4 + 1/2 y^5/5
1236     *                         - 1/2 y^7/7 + sqrt(3)/2 y^8/8
1237     *                         - y^9/9 + sqrt(3)/2 y^10/10 - 1/2 y^11/11
1238     *                         + 1/2 y^13/13 - sqrt(3)/2 y^14/14
1239     *                         + y^15/15 - sqrt(3)/2 y^16/16 + 1/2 y^17/17 ...
1240     *
1241     * If |x| >= 2 we evaluate atan(y) near +∞:
1242     *  atan(y) = π/2 - y^-1 + y^-3/3 - y^-5/5 + y^-7/7 - y^-9/9 ...
1243     */
1244    real absx = fabs(x);
1245
1246    if (absx < (real::R_1() / 2))
1247    {
1248        real ret = x, xn = x, mx2 = -x * x;
1249        for (int i = 3; ; i += 2)
1250        {
1251            xn *= mx2;
1252            real newret = ret + xn / (real)i;
1253            if (newret == ret)
1254                break;
1255            ret = newret;
1256        }
1257        return ret;
1258    }
1259
1260    real ret = 0;
1261
1262    if (absx < (real::R_3() / 2))
1263    {
1264        real y = real::R_1() - absx;
1265        real yn = y, my2 = -y * y;
1266        for (int i = 0; ; i += 2)
1267        {
1268            real newret = ret + ldexp(yn / (real)(2 * i + 1), -i - 1);
1269            yn *= y;
1270            newret += ldexp(yn / (real)(2 * i + 2), -i - 1);
1271            yn *= y;
1272            newret += ldexp(yn / (real)(2 * i + 3), -i - 2);
1273            if (newret == ret)
1274                break;
1275            ret = newret;
1276            yn *= my2;
1277        }
1278        ret = real::R_PI_4() - ret;
1279    }
1280    else if (absx < real::R_2())
1281    {
1282        real y = (absx - real::R_SQRT3()) / 2;
1283        real yn = y, my2 = -y * y;
1284        for (int i = 1; ; i += 6)
1285        {
1286            real newret = ret + ((yn / (real)i) / 2);
1287            yn *= y;
1288            newret -= (real::R_SQRT3() / 2) * yn / (real)(i + 1);
1289            yn *= y;
1290            newret += yn / (real)(i + 2);
1291            yn *= y;
1292            newret -= (real::R_SQRT3() / 2) * yn / (real)(i + 3);
1293            yn *= y;
1294            newret += (yn / (real)(i + 4)) / 2;
1295            if (newret == ret)
1296                break;
1297            ret = newret;
1298            yn *= my2;
1299        }
1300        ret = real::R_PI_3() + ret;
1301    }
1302    else
1303    {
1304        real y = re(absx);
1305        real yn = y, my2 = -y * y;
1306        ret = y;
1307        for (int i = 3; ; i += 2)
1308        {
1309            yn *= my2;
1310            real newret = ret + yn / (real)i;
1311            if (newret == ret)
1312                break;
1313            ret = newret;
1314        }
1315        ret = real::R_PI_2() - ret;
1316    }
1317
1318    /* Propagate sign */
1319    ret.m_signexp |= (x.m_signexp & 0x80000000u);
1320    return ret;
1321}
1322
1323template<> real atan2(real const &y, real const &x)
1324{
1325    if (!y)
1326    {
1327        if ((x.m_signexp >> 31) == 0)
1328            return y;
1329        if (y.m_signexp >> 31)
1330            return -real::R_PI();
1331        return real::R_PI();
1332    }
1333
1334    if (!x)
1335    {
1336        if (y.m_signexp >> 31)
1337            return -real::R_PI();
1338        return real::R_PI();
1339    }
1340
1341    /* FIXME: handle the Inf and NaN cases */
1342    real z = y / x;
1343    real ret = atan(z);
1344    if (x < real::R_0())
1345        ret += (y > real::R_0()) ? real::R_PI() : -real::R_PI();
1346    return ret;
1347}
1348
1349template<> void real::hexprint() const
1350{
1351    std::printf("%08x", m_signexp);
1352    for (int i = 0; i < BIGITS; i++)
1353        std::printf(" %08x", m_mantissa[i]);
1354    std::printf("\n");
1355}
1356
1357template<> void real::sprintf(char *str, int ndigits) const;
1358
1359template<> void real::print(int ndigits) const
1360{
1361    char *buf = new char[ndigits + 32 + 10];
1362    real::sprintf(buf, ndigits);
1363    std::printf("%s\n", buf);
1364    delete[] buf;
1365}
1366
1367template<> void real::sprintf(char *str, int ndigits) const
1368{
1369    real x = *this;
1370
1371    if (x.m_signexp >> 31)
1372    {
1373        *str++ = '-';
1374        x = -x;
1375    }
1376
1377    if (!x)
1378    {
1379        std::strcpy(str, "0.0\n");
1380        return;
1381    }
1382
1383    /* Normalise x so that mantissa is in [1..9.999] */
1384    /* FIXME: better use int64_t when the cast is implemented */
1385    int exponent = ceil(log10(x));
1386    x /= pow(R_10(), (real)exponent);
1387
1388    if (x < R_1())
1389    {
1390        x *= R_10();
1391        exponent--;
1392    }
1393
1394    /* Print digits */
1395    for (int i = 0; i < ndigits; i++)
1396    {
1397        int digit = (int)floor(x);
1398        *str++ = '0' + digit;
1399        if (i == 0)
1400            *str++ = '.';
1401        x -= real(digit);
1402        x *= R_10();
1403    }
1404
1405    /* Print exponent information */
1406    if (exponent)
1407        str += std::sprintf(str, "e%c%i", exponent > 0 ? '+' : '-', -exponent);
1408
1409    *str++ = '\0';
1410}
1411
1412static real fast_pi()
1413{
1414    /* Approximate Pi using Machin's formula: 16*atan(1/5)-4*atan(1/239) */
1415    real ret = 0, x0 = 5, x1 = 239;
1416    real const m0 = -x0 * x0, m1 = -x1 * x1, r16 = 16, r4 = 4;
1417
1418    for (int i = 1; ; i += 2)
1419    {
1420        real newret = ret + r16 / (x0 * (real)i) - r4 / (x1 * (real)i);
1421        if (newret == ret)
1422            break;
1423        ret = newret;
1424        x0 *= m0;
1425        x1 *= m1;
1426    }
1427
1428    return ret;
1429}
1430
1431} /* namespace lol */
1432
Note: See TracBrowser for help on using the repository browser.