source: trunk/src/real.cpp @ 1004

Last change on this file since 1004 was 1004, checked in by sam, 11 years ago

real: constrain sin() and cos() on real numbers so that they work properly
with large values. Until now they were evaluating the Taylor series even
for huge values.

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