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

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

math: add an sprintf() method to real numbers, and ensure they are always
fully initialised.

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