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

Last change on this file since 1513 was 1513, checked in by sam, 8 years ago

core: replace usage of sin() or std::sin() with lol::sin() where appropriate.

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