source: trunk/src/real.cpp @ 1023

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

core: implement log10, sinh and cosh for real numbers.

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