source: trunk/src/real.cpp @ 1016

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

core: encode real numbers using uint32_t rather than uint16_t.

File size: 24.7 KB
Line 
1//
2// Lol Engine
3//
4// Copyright: (c) 2010-2011 Sam Hocevar <sam@hocevar.net>
5//   This program is free software; you can redistribute it and/or
6//   modify it under the terms of the Do What The Fuck You Want To
7//   Public License, Version 2, as published by Sam Hocevar. See
8//   http://sam.zoy.org/projects/COPYING.WTFPL for more details.
9//
10
11#if defined HAVE_CONFIG_H
12#   include "config.h"
13#endif
14
15#include <cstring>
16#include <cstdio>
17
18#include "core.h"
19
20using namespace std;
21
22namespace lol
23{
24
25static int const BIGIT_BITS = 32;
26
27real::real(float f) { *this = (double)f; }
28real::real(int i) { *this = (double)i; }
29real::real(unsigned int i) { *this = (double)i; }
30
31real::real(double d)
32{
33    union { double d; uint64_t x; } u = { d };
34
35    uint32_t sign = (u.x >> 63) << 31;
36    uint32_t exponent = (u.x << 1) >> 53;
37
38    switch (exponent)
39    {
40    case 0x00:
41        m_signexp = sign;
42        break;
43    case 0x7ff:
44        m_signexp = sign | 0x7fffffffu;
45        break;
46    default:
47        m_signexp = sign | (exponent + (1 << 30) - (1 << 10));
48        break;
49    }
50
51    m_mantissa[0] = u.x >> 20;
52    m_mantissa[1] = u.x << 12;
53    memset(m_mantissa + 2, 0, (BIGITS - 2) * sizeof(m_mantissa[0]));
54}
55
56real::operator float() const { return (float)(double)(*this); }
57real::operator int() const { return (int)(double)(*this); }
58real::operator unsigned int() const { return (unsigned int)(double)(*this); }
59
60real::operator double() const
61{
62    union { double d; uint64_t x; } u;
63
64    /* Get sign */
65    u.x = m_signexp >> 31;
66    u.x <<= 11;
67
68    /* Compute new exponent */
69    uint32_t exponent = (m_signexp << 1) >> 1;
70    int e = (int)exponent - (1 << 30) + (1 << 10);
71
72    if (e < 0)
73        u.x <<= 52;
74    else if (e >= 0x7ff)
75    {
76        u.x |= 0x7ff;
77        u.x <<= 52;
78    }
79    else
80    {
81        u.x |= e;
82
83        /* Store mantissa if necessary */
84        u.x <<= 32;
85        u.x |= m_mantissa[0];
86        u.x <<= 20;
87        u.x |= m_mantissa[1] >> 12;
88        /* Rounding */
89        u.x += (m_mantissa[1] >> 11) & 1;
90    }
91
92    return u.d;
93}
94
95real real::operator +() const
96{
97    return *this;
98}
99
100real real::operator -() const
101{
102    real ret = *this;
103    ret.m_signexp ^= 0x80000000u;
104    return ret;
105}
106
107real real::operator +(real const &x) const
108{
109    if (x.m_signexp << 1 == 0)
110        return *this;
111
112    /* Ensure both arguments are positive. Otherwise, switch signs,
113     * or replace + with -. */
114    if (m_signexp >> 31)
115        return -(-*this + -x);
116
117    if (x.m_signexp >> 31)
118        return *this - (-x);
119
120    /* Ensure *this has the larger exponent (no need for the mantissa to
121     * be larger, as in subtraction). Otherwise, switch. */
122    if ((m_signexp << 1) < (x.m_signexp << 1))
123        return x + *this;
124
125    real ret;
126
127    int e1 = m_signexp - (1 << 30) + 1;
128    int e2 = x.m_signexp - (1 << 30) + 1;
129
130    int bigoff = (e1 - e2) / BIGIT_BITS;
131    int off = e1 - e2 - bigoff * BIGIT_BITS;
132
133    if (bigoff > BIGITS)
134        return *this;
135
136    ret.m_signexp = m_signexp;
137
138    uint64_t carry = 0;
139    for (int i = BIGITS; i--; )
140    {
141        carry += m_mantissa[i];
142        if (i - bigoff >= 0)
143            carry += x.m_mantissa[i - bigoff] >> off;
144
145        if (off && i - bigoff > 0)
146            carry += (x.m_mantissa[i - bigoff - 1] << (BIGIT_BITS - off)) & 0xffffffffu;
147        else if (i - bigoff == 0)
148            carry += (uint64_t)1 << (BIGIT_BITS - off);
149
150        ret.m_mantissa[i] = carry;
151        carry >>= BIGIT_BITS;
152    }
153
154    /* Renormalise in case we overflowed the mantissa */
155    if (carry)
156    {
157        carry--;
158        for (int i = 0; i < BIGITS; i++)
159        {
160            uint32_t tmp = ret.m_mantissa[i];
161            ret.m_mantissa[i] = (carry << (BIGIT_BITS - 1)) | (tmp >> 1);
162            carry = tmp & 1u;
163        }
164        ret.m_signexp++;
165    }
166
167    return ret;
168}
169
170real real::operator -(real const &x) const
171{
172    if (x.m_signexp << 1 == 0)
173        return *this;
174
175    /* Ensure both arguments are positive. Otherwise, switch signs,
176     * or replace - with +. */
177    if (m_signexp >> 31)
178        return -(-*this + x);
179
180    if (x.m_signexp >> 31)
181        return (*this) + (-x);
182
183    /* Ensure *this is larger than x */
184    if (*this < x)
185        return -(x - *this);
186
187    real ret;
188
189    int e1 = m_signexp - (1 << 30) + 1;
190    int e2 = x.m_signexp - (1 << 30) + 1;
191
192    int bigoff = (e1 - e2) / BIGIT_BITS;
193    int off = e1 - e2 - bigoff * BIGIT_BITS;
194
195    if (bigoff > BIGITS)
196        return *this;
197
198    ret.m_signexp = m_signexp;
199
200    int64_t carry = 0;
201    for (int i = 0; i < bigoff; i++)
202    {
203        carry -= x.m_mantissa[BIGITS - i];
204        /* Emulates a signed shift */
205        carry >>= BIGIT_BITS;
206        carry |= carry << BIGIT_BITS;
207    }
208    carry -= x.m_mantissa[BIGITS - 1 - bigoff] & (((int64_t)1 << off) - 1);
209    carry /= (int64_t)1 << off;
210
211    for (int i = BIGITS; i--; )
212    {
213        carry += m_mantissa[i];
214        if (i - bigoff >= 0)
215            carry -= x.m_mantissa[i - bigoff] >> off;
216
217        if (off && i - bigoff > 0)
218            carry -= (x.m_mantissa[i - bigoff - 1] << (BIGIT_BITS - off)) & 0xffffffffu;
219        else if (i - bigoff == 0)
220            carry -= (uint64_t)1 << (BIGIT_BITS - off);
221
222        ret.m_mantissa[i] = carry;
223        carry >>= BIGIT_BITS;
224        carry |= carry << BIGIT_BITS;
225    }
226
227    carry += 1;
228
229    /* Renormalise if we underflowed the mantissa */
230    if (carry == 0)
231    {
232        /* How much do we need to shift the mantissa? FIXME: this could
233         * be computed above */
234        off = 0;
235        for (int i = 0; i < BIGITS; i++)
236        {
237            if (!ret.m_mantissa[i])
238            {
239                off += BIGIT_BITS;
240                continue;
241            }
242
243            for (uint32_t tmp = ret.m_mantissa[i]; tmp < 0x80000000u; tmp <<= 1)
244                off++;
245            break;
246        }
247        if (off == BIGITS * BIGIT_BITS)
248            ret.m_signexp &= 0x80000000u;
249        else
250        {
251            off++; /* Shift one more to get rid of the leading one */
252            ret.m_signexp -= off;
253
254            bigoff = off / BIGIT_BITS;
255            off -= bigoff * BIGIT_BITS;
256
257            for (int i = 0; i < BIGITS; i++)
258            {
259                uint32_t tmp = 0;
260                if (i + bigoff < BIGITS)
261                    tmp |= ret.m_mantissa[i + bigoff] << off;
262                if (off && i + bigoff + 1 < BIGITS)
263                    tmp |= ret.m_mantissa[i + bigoff + 1] >> (BIGIT_BITS - off);
264                ret.m_mantissa[i] = tmp;
265            }
266        }
267    }
268
269    return ret;
270}
271
272real real::operator *(real const &x) const
273{
274    real ret;
275
276    if (m_signexp << 1 == 0 || x.m_signexp << 1 == 0)
277    {
278        ret = (m_signexp << 1 == 0) ? *this : x;
279        ret.m_signexp ^= x.m_signexp & 0x80000000u;
280        return ret;
281    }
282
283    ret.m_signexp = (m_signexp ^ x.m_signexp) & 0x80000000u;
284    int e = (m_signexp & 0x7fffffffu) - (1 << 30) + 1
285          + (x.m_signexp & 0x7fffffffu) - (1 << 30) + 1;
286
287    /* Accumulate low order product; no need to store it, we just
288     * want the carry value */
289    uint64_t carry = 0, hicarry = 0, prev;
290    for (int i = 0; i < BIGITS; i++)
291    {
292        for (int j = 0; j < i + 1; j++)
293        {
294            prev = carry;
295            carry += (uint64_t)m_mantissa[BIGITS - 1 - j]
296                   * (uint64_t)x.m_mantissa[BIGITS - 1 + j - i];
297            if (carry < prev)
298                hicarry++;
299        }
300        carry >>= BIGIT_BITS;
301        carry |= hicarry << BIGIT_BITS;
302        hicarry >>= BIGIT_BITS;
303    }
304
305    for (int i = 0; i < BIGITS; i++)
306    {
307        for (int j = i + 1; j < BIGITS; j++)
308        {
309            prev = carry;
310            carry += (uint64_t)m_mantissa[BIGITS - 1 - j]
311                   * (uint64_t)x.m_mantissa[j - 1 - i];
312            if (carry < prev)
313                hicarry++;
314        }
315        prev = carry;
316        carry += m_mantissa[BIGITS - 1 - i];
317        carry += x.m_mantissa[BIGITS - 1 - i];
318        if (carry < prev)
319            hicarry++;
320        ret.m_mantissa[BIGITS - 1 - i] = carry & 0xffffffffu;
321        carry >>= BIGIT_BITS;
322        carry |= hicarry << BIGIT_BITS;
323        hicarry >>= BIGIT_BITS;
324    }
325
326    /* Renormalise in case we overflowed the mantissa */
327    if (carry)
328    {
329        carry--;
330        for (int i = 0; i < BIGITS; i++)
331        {
332            uint32_t tmp = ret.m_mantissa[i];
333            ret.m_mantissa[i] = (carry << (BIGIT_BITS - 1)) | (tmp >> 1);
334            carry = tmp & 1u;
335        }
336        e++;
337    }
338
339    ret.m_signexp |= e + (1 << 30) - 1;
340
341    return ret;
342}
343
344real real::operator /(real const &x) const
345{
346    return *this * re(x);
347}
348
349real &real::operator +=(real const &x)
350{
351    real tmp = *this;
352    return *this = tmp + x;
353}
354
355real &real::operator -=(real const &x)
356{
357    real tmp = *this;
358    return *this = tmp - x;
359}
360
361real &real::operator *=(real const &x)
362{
363    real tmp = *this;
364    return *this = tmp * x;
365}
366
367real &real::operator /=(real const &x)
368{
369    real tmp = *this;
370    return *this = tmp / x;
371}
372
373real real::operator <<(int x) const
374{
375    real tmp = *this;
376    return tmp <<= x;
377}
378
379real real::operator >>(int x) const
380{
381    real tmp = *this;
382    return tmp >>= x;
383}
384
385real &real::operator <<=(int x)
386{
387    if (m_signexp << 1)
388        m_signexp += x;
389    return *this;
390}
391
392real &real::operator >>=(int x)
393{
394    if (m_signexp << 1)
395        m_signexp -= x;
396    return *this;
397}
398
399bool real::operator ==(real const &x) const
400{
401    if ((m_signexp << 1) == 0 && (x.m_signexp << 1) == 0)
402        return true;
403
404    if (m_signexp != x.m_signexp)
405        return false;
406
407    return memcmp(m_mantissa, x.m_mantissa, sizeof(m_mantissa)) == 0;
408}
409
410bool real::operator !=(real const &x) const
411{
412    return !(*this == x);
413}
414
415bool real::operator <(real const &x) const
416{
417    /* Ensure both numbers are positive */
418    if (m_signexp >> 31)
419        return (x.m_signexp >> 31) ? -*this > -x : true;
420
421    if (x.m_signexp >> 31)
422        return false;
423
424    /* Compare all relevant bits */
425    if (m_signexp != x.m_signexp)
426        return m_signexp < x.m_signexp;
427
428    for (int i = 0; i < BIGITS; i++)
429        if (m_mantissa[i] != x.m_mantissa[i])
430            return m_mantissa[i] < x.m_mantissa[i];
431
432    return false;
433}
434
435bool real::operator <=(real const &x) const
436{
437    return !(*this > x);
438}
439
440bool real::operator >(real const &x) const
441{
442    /* Ensure both numbers are positive */
443    if (m_signexp >> 31)
444        return (x.m_signexp >> 31) ? -*this < -x : false;
445
446    if (x.m_signexp >> 31)
447        return true;
448
449    /* Compare all relevant bits */
450    if (m_signexp != x.m_signexp)
451        return m_signexp > x.m_signexp;
452
453    for (int i = 0; i < BIGITS; i++)
454        if (m_mantissa[i] != x.m_mantissa[i])
455            return m_mantissa[i] > x.m_mantissa[i];
456
457    return false;
458}
459
460bool real::operator >=(real const &x) const
461{
462    return !(*this < x);
463}
464
465bool real::operator !() const
466{
467    return !(bool)*this;
468}
469
470real::operator bool() const
471{
472    /* A real is "true" if it is non-zero (exponent is non-zero) AND
473     * not NaN (exponent is not full bits OR higher order mantissa is zero) */
474    uint32_t exponent = m_signexp << 1;
475    return exponent && (~exponent || m_mantissa[0] == 0);
476}
477
478real re(real const &x)
479{
480    if (!(x.m_signexp << 1))
481    {
482        real ret = x;
483        ret.m_signexp = x.m_signexp | 0x7fffffffu;
484        ret.m_mantissa[0] = 0;
485        return ret;
486    }
487
488    /* Use the system's float inversion to approximate 1/x */
489    union { float f; uint32_t x; } u = { 1.0f }, v = { 1.0f };
490    v.x |= x.m_mantissa[0] >> 9;
491    v.f = 1.0 / v.f;
492
493    real ret;
494    ret.m_mantissa[0] = v.x << 9;
495
496    uint32_t sign = x.m_signexp & 0x80000000u;
497    ret.m_signexp = sign;
498
499    int exponent = (x.m_signexp & 0x7fffffffu) + 1;
500    exponent = -exponent + (v.x >> 23) - (u.x >> 23);
501    ret.m_signexp |= (exponent - 1) & 0x7fffffffu;
502
503    /* FIXME: log2(BIGITS) steps of Newton-Raphson seems to be enough for
504     * convergence, but this hasn't been checked seriously. */
505    for (int i = 2; i < real::BIGITS; i *= 2)
506        ret = ret * (real::R_2 - ret * x);
507
508    return ret;
509}
510
511real sqrt(real const &x)
512{
513    /* if zero, return x */
514    if (!(x.m_signexp << 1))
515        return x;
516
517    /* if negative, return NaN */
518    if (x.m_signexp >> 31)
519    {
520        real ret;
521        ret.m_signexp = 0x7fffffffu;
522        ret.m_mantissa[0] = 0xffffu;
523        return ret;
524    }
525
526    /* Use the system's float inversion to approximate 1/sqrt(x). First
527     * we construct a float in the [1..4[ range that has roughly the same
528     * mantissa as our real. Its exponent is 0 or 1, depending on the
529     * partity of x. The final exponent is 0, -1 or -2. We use the final
530     * exponent and final mantissa to pre-fill the result. */
531    union { float f; uint32_t x; } u = { 1.0f }, v = { 2.0f };
532    v.x -= ((x.m_signexp & 1) << 23);
533    v.x |= x.m_mantissa[0] >> 9;
534    v.f = 1.0 / sqrtf(v.f);
535
536    real ret;
537    ret.m_mantissa[0] = v.x << 9;
538
539    uint32_t sign = x.m_signexp & 0x80000000u;
540    ret.m_signexp = sign;
541
542    uint32_t exponent = (x.m_signexp & 0x7fffffffu);
543    exponent = ((1 << 30) + (1 << 29) - 1) - (exponent + 1) / 2;
544    exponent = exponent + (v.x >> 23) - (u.x >> 23);
545    ret.m_signexp |= exponent & 0x7fffffffu;
546
547    /* FIXME: log2(BIGITS) steps of Newton-Raphson seems to be enough for
548     * convergence, but this hasn't been checked seriously. */
549    for (int i = 2; i < real::BIGITS; i *= 2)
550    {
551        ret = ret * (real::R_3 - ret * ret * x);
552        ret.m_signexp--;
553    }
554
555    return ret * x;
556}
557
558real fabs(real const &x)
559{
560    real ret = x;
561    ret.m_signexp &= 0x7fffffffu;
562    return ret;
563}
564
565static real fast_log(real const &x)
566{
567    /* This fast log method is tuned to work on the [1..2] range and
568     * no effort whatsoever was made to improve convergence outside this
569     * domain of validity. It can converge pretty fast, provided we use
570     * the following variable substitutions:
571     *    y = sqrt(x)
572     *    z = (y - 1) / (y + 1)
573     *
574     * And the following identities:
575     *    ln(x) = 2 ln(y)
576     *          = 2 ln((1 + z) / (1 - z))
577     *          = 4 z (1 + z^2 / 3 + z^4 / 5 + z^6 / 7...)
578     *
579     * Any additional sqrt() call would halve the convergence time, but
580     * would also impact the final precision. For now we stick with one
581     * sqrt() call. */
582    real y = sqrt(x);
583    real z = (y - real::R_1) / (y + real::R_1), z2 = z * z, zn = z2;
584    real sum = real::R_1;
585
586    for (int i = 3; ; i += 2)
587    {
588        real newsum = sum + zn / (real)i;
589        if (newsum == sum)
590            break;
591        sum = newsum;
592        zn *= z2;
593    }
594
595    return z * (sum << 2);
596}
597
598real log(real const &x)
599{
600    /* Strategy for log(x): if x = 2^E*M then log(x) = E log(2) + log(M),
601     * with the property that M is in [1..2[, so fast_log() applies here. */
602    real tmp = x;
603    if (x.m_signexp >> 31 || x.m_signexp == 0)
604    {
605        tmp.m_signexp = 0xffffffffu;
606        tmp.m_mantissa[0] = 0xffffffffu;
607        return tmp;
608    }
609    tmp.m_signexp = (1 << 30) - 1;
610    return (real)(x.m_signexp - (1 << 30) + 1) * real::R_LN2 + fast_log(tmp);
611}
612
613real exp(real const &x)
614{
615    /* Strategy for exp(x): the Taylor series does not converge very fast
616     * with large positive or negative values.
617     *
618     * However, we know that the result is going to be in the form M*2^E,
619     * where M is the mantissa and E the exponent. We first try to predict
620     * a value for E, which is approximately log2(exp(x)) = x / log(2).
621     *
622     * Let E0 be an integer close to x / log(2). We need to find a value x0
623     * such that exp(x) = 2^E0 * exp(x0). We get x0 = x - E0 log(2).
624     *
625     * Thus the final algorithm:
626     *  int E0 = x / log(2)
627     *  real x0 = x - E0 log(2)
628     *  real x1 = exp(x0)
629     *  return x1 * 2^E0
630     */
631    int e0 = x / real::R_LN2;
632    real x0 = x - (real)e0 * real::R_LN2;
633    real x1 = real::R_1, fact = real::R_1, xn = x0;
634
635    for (int i = 1; ; i++)
636    {
637        fact *= (real)i;
638        real newx1 = x1 + xn / fact;
639        if (newx1 == x1)
640            break;
641        x1 = newx1;
642        xn *= x0;
643    }
644
645    x1.m_signexp += e0;
646    return x1;
647}
648
649real floor(real const &x)
650{
651    /* Strategy for floor(x):
652     *  - if negative, return -ceil(-x)
653     *  - if zero or negative zero, return x
654     *  - if less than one, return zero
655     *  - otherwise, if e is the exponent, clear all bits except the
656     *    first e. */
657    if (x < -real::R_0)
658        return -ceil(-x);
659    if (!x)
660        return x;
661    if (x < real::R_1)
662        return real::R_0;
663
664    real ret = x;
665    int exponent = x.m_signexp - (1 << 30) + 1;
666
667    for (int i = 0; i < real::BIGITS; i++)
668    {
669        if (exponent <= 0)
670            ret.m_mantissa[i] = 0;
671        else if (exponent < BIGIT_BITS)
672            ret.m_mantissa[i] &= ~((1 << (BIGIT_BITS - exponent)) - 1);
673
674        exponent -= BIGIT_BITS;
675    }
676
677    return ret;
678}
679
680real ceil(real const &x)
681{
682    /* Strategy for ceil(x):
683     *  - if negative, return -floor(-x)
684     *  - if x == floor(x), return x
685     *  - otherwise, return floor(x) + 1 */
686    if (x < -real::R_0)
687        return -floor(-x);
688    real ret = floor(x);
689    if (x == ret)
690        return ret;
691    else
692        return ret + real::R_1;
693}
694
695real round(real const &x)
696{
697    if (x < real::R_0)
698        return -round(-x);
699
700    return floor(x + (real::R_1 >> 1));
701}
702
703real fmod(real const &x, real const &y)
704{
705    if (!y)
706        return real::R_0; /* FIXME: return NaN */
707
708    if (!x)
709        return x;
710
711    real tmp = round(x / y);
712    return x - tmp * y;
713}
714
715real sin(real const &x)
716{
717    bool switch_sign = x.m_signexp & 0x80000000u;
718
719    real absx = fmod(fabs(x), real::R_PI << 1);
720    if (absx > real::R_PI)
721    {
722        absx -= real::R_PI;
723        switch_sign = !switch_sign;
724    }
725
726    if (absx > real::R_PI_2)
727        absx = real::R_PI - absx;
728
729    real ret = real::R_0, fact = real::R_1, xn = absx, x2 = absx * absx;
730    for (int i = 1; ; i += 2)
731    {
732        real newret = ret + xn / fact;
733        if (newret == ret)
734            break;
735        ret = newret;
736        xn *= x2;
737        fact *= (real)(-(i + 1) * (i + 2));
738    }
739
740    /* Propagate sign */
741    if (switch_sign)
742        ret.m_signexp ^= 0x80000000u;
743    return ret;
744}
745
746real cos(real const &x)
747{
748    return sin(real::R_PI_2 - x);
749}
750
751static real asinacos(real const &x, bool is_asin, bool is_negative)
752{
753    /* Strategy for asin(): in [-0.5..0.5], use a Taylor series around
754     * zero. In [0.5..1], use asin(x) = π/2 - 2*asin(sqrt((1-x)/2)), and
755     * in [-1..-0.5] just revert the sign.
756     * Strategy for acos(): use acos(x) = π/2 - asin(x) and try not to
757     * lose the precision around x=1. */
758    real absx = fabs(x);
759    bool around_zero = (absx < (real::R_1 >> 1));
760
761    if (!around_zero)
762        absx = sqrt((real::R_1 - absx) >> 1);
763
764    real ret = absx, xn = absx, x2 = absx * absx, fact1 = 2, fact2 = 1;
765    for (int i = 1; ; i++)
766    {
767        xn *= x2;
768        real mul = (real)(2 * i + 1);
769        real newret = ret + ((fact1 * xn / (mul * fact2)) >> (i * 2));
770        if (newret == ret)
771            break;
772        ret = newret;
773        fact1 *= (real)((2 * i + 1) * (2 * i + 2));
774        fact2 *= (real)((i + 1) * (i + 1));
775    }
776
777    if (is_negative)
778        ret = -ret;
779
780    if (around_zero)
781        ret = is_asin ? ret : real::R_PI_2 - ret;
782    else
783    {
784        real adjust = is_negative ? real::R_PI : real::R_0;
785        if (is_asin)
786            ret = real::R_PI_2 - adjust - (ret << 1);
787        else
788            ret = adjust + (ret << 1);
789    }
790
791    return ret;
792}
793
794real asin(real const &x)
795{
796    return asinacos(x, true, x.m_signexp >> 31);
797}
798
799real acos(real const &x)
800{
801    return asinacos(x, false, x.m_signexp >> 31);
802}
803
804real atan(real const &x)
805{
806    /* Computing atan(x): we choose a different Taylor series depending on
807     * the value of x to help with convergence.
808     *
809     * If |x| < 0.5 we evaluate atan(y) near 0:
810     *  atan(y) = y - y^3/3 + y^5/5 - y^7/7 + y^9/9 ...
811     *
812     * If 0.5 <= |x| < 1.5 we evaluate atan(1+y) near 0:
813     *  atan(1+y) = π/4 + y/(1*2^1) - y^2/(2*2^1) + y^3/(3*2^2)
814     *                  - y^5/(5*2^3) + y^6/(6*2^3) - y^7/(7*2^4)
815     *                  + y^9/(9*2^5) - y^10/(10*2^5) + y^11/(11*2^6) ...
816     *
817     * If 1.5 <= |x| < 2 we evaluate atan(sqrt(3)+2y) near 0:
818     *  atan(sqrt(3)+2y) = π/3 + 1/2 y - sqrt(3)/2 y^2/2
819     *                         + y^3/3 - sqrt(3)/2 y^4/4 + 1/2 y^5/5
820     *                         - 1/2 y^7/7 + sqrt(3)/2 y^8/8
821     *                         - y^9/9 + sqrt(3)/2 y^10/10 - 1/2 y^11/11
822     *                         + 1/2 y^13/13 - sqrt(3)/2 y^14/14
823     *                         + y^15/15 - sqrt(3)/2 y^16/16 + 1/2 y^17/17 ...
824     *
825     * If |x| >= 2 we evaluate atan(y) near +∞:
826     *  atan(y) = π/2 - y^-1 + y^-3/3 - y^-5/5 + y^-7/7 - y^-9/9 ...
827     */
828    real absx = fabs(x);
829
830    if (absx < (real::R_1 >> 1))
831    {
832        real ret = x, xn = x, mx2 = -x * x;
833        for (int i = 3; ; i += 2)
834        {
835            xn *= mx2;
836            real newret = ret + xn / (real)i;
837            if (newret == ret)
838                break;
839            ret = newret;
840        }
841        return ret;
842    }
843
844    real ret = 0;
845
846    if (absx < (real::R_3 >> 1))
847    {
848        real y = real::R_1 - absx;
849        real yn = y, my2 = -y * y;
850        for (int i = 0; ; i += 2)
851        {
852            real newret = ret + ((yn / (real)(2 * i + 1)) >> (i + 1));
853            yn *= y;
854            newret += (yn / (real)(2 * i + 2)) >> (i + 1);
855            yn *= y;
856            newret += (yn / (real)(2 * i + 3)) >> (i + 2);
857            if (newret == ret)
858                break;
859            ret = newret;
860            yn *= my2;
861        }
862        ret = real::R_PI_4 - ret;
863    }
864    else if (absx < real::R_2)
865    {
866        real y = (absx - real::R_SQRT3) >> 1;
867        real yn = y, my2 = -y * y;
868        for (int i = 1; ; i += 6)
869        {
870            real newret = ret + ((yn / (real)i) >> 1);
871            yn *= y;
872            newret -= (real::R_SQRT3 >> 1) * yn / (real)(i + 1);
873            yn *= y;
874            newret += yn / (real)(i + 2);
875            yn *= y;
876            newret -= (real::R_SQRT3 >> 1) * yn / (real)(i + 3);
877            yn *= y;
878            newret += (yn / (real)(i + 4)) >> 1;
879            if (newret == ret)
880                break;
881            ret = newret;
882            yn *= my2;
883        }
884        ret = real::R_PI_3 + ret;
885    }
886    else
887    {
888        real y = re(absx);
889        real yn = y, my2 = -y * y;
890        ret = y;
891        for (int i = 3; ; i += 2)
892        {
893            yn *= my2;
894            real newret = ret + yn / (real)i;
895            if (newret == ret)
896                break;
897            ret = newret;
898        }
899        ret = real::R_PI_2 - ret;
900    }
901
902    /* Propagate sign */
903    ret.m_signexp |= (x.m_signexp & 0x80000000u);
904    return ret;
905}
906
907void real::print(int ndigits) const
908{
909    real const r1 = 1, r10 = 10;
910    real x = *this;
911
912    if (x.m_signexp >> 31)
913    {
914        printf("-");
915        x = -x;
916    }
917
918    /* Normalise x so that mantissa is in [1..9.999] */
919    int exponent = 0;
920    if (x.m_signexp)
921    {
922        for (real div = r1, newdiv; true; div = newdiv)
923        {
924            newdiv = div * r10;
925            if (x < newdiv)
926            {
927                x /= div;
928                break;
929            }
930            exponent++;
931        }
932        for (real mul = 1, newx; true; mul *= r10)
933        {
934            newx = x * mul;
935            if (newx >= r1)
936            {
937                x = newx;
938                break;
939            }
940            exponent--;
941        }
942    }
943
944    /* Print digits */
945    for (int i = 0; i < ndigits; i++)
946    {
947        int digit = (int)x;
948        printf("%i", digit);
949        if (i == 0)
950            printf(".");
951        x -= real(digit);
952        x *= r10;
953    }
954
955    /* Print exponent information */
956    if (exponent < 0)
957        printf("e-%i", -exponent);
958    else if (exponent > 0)
959        printf("e+%i", exponent);
960
961    printf("\n");
962}
963
964static real fast_pi()
965{
966    /* Approximate Pi using Machin's formula: 16*atan(1/5)-4*atan(1/239) */
967    real ret = 0.0, x0 = 5.0, x1 = 239.0;
968    real const m0 = -x0 * x0, m1 = -x1 * x1, r16 = 16.0, r4 = 4.0;
969
970    for (int i = 1; ; i += 2)
971    {
972        real newret = ret + r16 / (x0 * (real)i) - r4 / (x1 * (real)i);
973        if (newret == ret)
974            break;
975        ret = newret;
976        x0 *= m0;
977        x1 *= m1;
978    }
979
980    return ret;
981}
982
983real const real::R_0        = (real)0.0;
984real const real::R_1        = (real)1.0;
985real const real::R_2        = (real)2.0;
986real const real::R_3        = (real)3.0;
987real const real::R_10       = (real)10.0;
988
989real const real::R_LN2      = fast_log(R_2);
990real const real::R_LN10     = log(R_10);
991real const real::R_LOG2E    = re(R_LN2);
992real const real::R_LOG10E   = re(R_LN10);
993real const real::R_E        = exp(R_1);
994real const real::R_PI       = fast_pi();
995real const real::R_PI_2     = R_PI >> 1;
996real const real::R_PI_3     = R_PI / R_3;
997real const real::R_PI_4     = R_PI >> 2;
998real const real::R_1_PI     = re(R_PI);
999real const real::R_2_PI     = R_1_PI << 1;
1000real const real::R_2_SQRTPI = re(sqrt(R_PI)) << 1;
1001real const real::R_SQRT2    = sqrt(R_2);
1002real const real::R_SQRT3    = sqrt(R_3);
1003real const real::R_SQRT1_2  = R_SQRT2 >> 1;
1004
1005} /* namespace lol */
1006
Note: See TracBrowser for help on using the repository browser.