lolengine
/
lol
дзеркало https://github.com/lolengine/lol
1
0
Форк 0
Код Проблеми 0 Релізи 0 Вікі Активність
2082 Коміти
5 Гілки
1.3 GiB
 
 
 
Дерево: cdfb4e7abf
lol/src/math/real.cpp

1445 рядки
36 KiB
Неформатований Звинувачення Історія 

  1. //

  2. // Lol Engine

  3. //

  4. // Copyright: (c) 2010-2013 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://www.wtfpl.net/ for more details.

  9. //

  10. 

  11. #if defined HAVE_CONFIG_H

  12. #   include "config.h"

  13. #endif

  14. 

  15. #include <new>

  16. #include <cstring>

  17. #include <cstdio>

  18. #include <cstdlib>

  19. 

  20. #include "core.h"

  21. 

  22. using namespace std;

  23. 

  24. namespace lol

  25. {

  26. 

  27. /*

  28.  * First handle explicit specialisation of our templates.

  29.  *

  30.  * Initialisation order is not important because everything is

  31.  * done on demand, but here is the dependency list anyway:

  32.  *  - fast_log() requires R_1

  33.  *  - log() requires R_LN2

  34.  *  - re() require R_2

  35.  *  - exp() requires R_0, R_1, R_LN2

  36.  *  - sqrt() requires R_3

  37.  */

  38. 

  39. static real fast_log(real const &x);

  40. static real fast_pi();

  41. 

  42. #define LOL_CONSTANT_GETTER(name, value) \

  43.     template<> real const& real::name() \

  44.     { \

  45.         static real const ret = value; \

  46.         return ret; \

  47.     }

  48. 

  49. LOL_CONSTANT_GETTER(R_0,        (real)0.0);

  50. LOL_CONSTANT_GETTER(R_1,        (real)1.0);

  51. LOL_CONSTANT_GETTER(R_2,        (real)2.0);

  52. LOL_CONSTANT_GETTER(R_3,        (real)3.0);

  53. LOL_CONSTANT_GETTER(R_10,       (real)10.0);

  54. 

  55. LOL_CONSTANT_GETTER(R_LN2,      fast_log(R_2()));

  56. LOL_CONSTANT_GETTER(R_LN10,     log(R_10()));

  57. LOL_CONSTANT_GETTER(R_LOG2E,    re(R_LN2()));

  58. LOL_CONSTANT_GETTER(R_LOG10E,   re(R_LN10()));

  59. LOL_CONSTANT_GETTER(R_E,        exp(R_1()));

  60. LOL_CONSTANT_GETTER(R_PI,       fast_pi());

  61. LOL_CONSTANT_GETTER(R_PI_2,     R_PI() / 2);

  62. LOL_CONSTANT_GETTER(R_PI_3,     R_PI() / R_3());

  63. LOL_CONSTANT_GETTER(R_PI_4,     R_PI() / 4);

  64. LOL_CONSTANT_GETTER(R_1_PI,     re(R_PI()));

  65. LOL_CONSTANT_GETTER(R_2_PI,     R_1_PI() * 2);

  66. LOL_CONSTANT_GETTER(R_2_SQRTPI, re(sqrt(R_PI())) * 2);

  67. LOL_CONSTANT_GETTER(R_SQRT2,    sqrt(R_2()));

  68. LOL_CONSTANT_GETTER(R_SQRT3,    sqrt(R_3()));

  69. LOL_CONSTANT_GETTER(R_SQRT1_2,  R_SQRT2() / 2);

  70. 

  71. #undef LOL_CONSTANT_GETTER

  72. 

  73. /*

  74.  * Now carry on with the rest of the Real class.

  75.  */

  76. 

  77. template<> real::Real()

  78. {

  79.     m_mantissa = new uint32_t[BIGITS];

  80.     memset(m_mantissa, 0, BIGITS * sizeof(uint32_t));

  81.     m_signexp = 0;

  82. }

  83. 

  84. template<> real::Real(real const &x)

  85. {

  86.     m_mantissa = new uint32_t[BIGITS];

  87.     memcpy(m_mantissa, x.m_mantissa, BIGITS * sizeof(uint32_t));

  88.     m_signexp = x.m_signexp;

  89. }

  90. 

  91. template<> real const &real::operator =(real const &x)

  92. {

  93.     if (&x != this)

  94.     {

  95.         memcpy(m_mantissa, x.m_mantissa, BIGITS * sizeof(uint32_t));

  96.         m_signexp = x.m_signexp;

  97.     }

  98. 

  99.     return *this;

  100. }

  101. 

  102. template<> real::~Real()

  103. {

  104.     delete[] m_mantissa;

  105. }

  106. 

  107. template<> real::Real(float f) { new(this) real((double)f); }

  108. template<> real::Real(int i) { new(this) real((double)i); }

  109. template<> real::Real(unsigned int i) { new(this) real((double)i); }

  110. 

  111. template<> real::Real(double d)

  112. {

  113.     new(this) real();

  114. 

  115.     union { double d; uint64_t x; } u = { d };

  116. 

  117.     uint32_t sign = (u.x >> 63) << 31;

  118.     uint32_t exponent = (u.x << 1) >> 53;

  119. 

  120.     switch (exponent)

  121.     {

  122.     case 0x00:

  123.         m_signexp = sign;

  124.         break;

  125.     case 0x7ff:

  126.         m_signexp = sign | 0x7fffffffu;

  127.         break;

  128.     default:

  129.         m_signexp = sign | (exponent + (1 << 30) - (1 << 10));

  130.         break;

  131.     }

  132. 

  133.     m_mantissa[0] = (uint32_t)(u.x >> 20);

  134.     m_mantissa[1] = (uint32_t)(u.x << 12);

  135.     memset(m_mantissa + 2, 0, (BIGITS - 2) * sizeof(m_mantissa[0]));

  136. }

  137. 

  138. template<> real::operator float() const { return (float)(double)(*this); }

  139. template<> real::operator int() const { return (int)(double)(*this); }

  140. template<> real::operator unsigned() const { return (unsigned)(double)(*this); }

  141. 

  142. template<> real::operator double() const

  143. {

  144.     union { double d; uint64_t x; } u;

  145. 

  146.     /* Get sign */

  147.     u.x = m_signexp >> 31;

  148.     u.x <<= 11;

  149. 

  150.     /* Compute new exponent */

  151.     uint32_t exponent = (m_signexp << 1) >> 1;

  152.     int e = (int)exponent - (1 << 30) + (1 << 10);

  153. 

  154.     if (e < 0)

  155.         u.x <<= 52;

  156.     else if (e >= 0x7ff)

  157.     {

  158.         u.x |= 0x7ff;

  159.         u.x <<= 52;

  160.     }

  161.     else

  162.     {

  163.         u.x |= e;

  164. 

  165.         /* Store mantissa if necessary */

  166.         u.x <<= 32;

  167.         u.x |= m_mantissa[0];

  168.         u.x <<= 20;

  169.         u.x |= m_mantissa[1] >> 12;

  170.         /* Rounding */

  171.         u.x += (m_mantissa[1] >> 11) & 1;

  172.     }

  173. 

  174.     return u.d;

  175. }

  176. 

  177. /*

  178.  * Create a real number from an ASCII representation

  179.  */

  180. template<> real::Real(char const *str)

  181. {

  182.     real ret = 0;

  183.     int exponent = 0;

  184.     bool comma = false, nonzero = false, negative = false, finished = false;

  185. 

  186.     for (char const *p = str; *p && !finished; p++)

  187.     {

  188.         switch (*p)

  189.         {

  190.         case '-':

  191.         case '+':

  192.             if (p != str)

  193.                 break;

  194.             negative = (*p == '-');

  195.             break;

  196.         case '.':

  197.             if (comma)

  198.                 finished = true;

  199.             comma = true;

  200.             break;

  201.         case '0': case '1': case '2': case '3': case '4':

  202.         case '5': case '6': case '7': case '8': case '9':

  203.             if (nonzero)

  204.             {

  205.                 real x = ret + ret;

  206.                 x = x + x + ret;

  207.                 ret = x + x;

  208.             }

  209.             if (*p != '0')

  210.             {

  211.                 ret += (int)(*p - '0');

  212.                 nonzero = true;

  213.             }

  214.             if (comma)

  215.                 exponent--;

  216.             break;

  217.         case 'e':

  218.         case 'E':

  219.             exponent += atoi(p + 1);

  220.             finished = true;

  221.             break;

  222.         default:

  223.             finished = true;

  224.             break;

  225.         }

  226.     }

  227. 

  228.     if (exponent)

  229.         ret *= pow(R_10(), (real)exponent);

  230. 

  231.     if (negative)

  232.         ret = -ret;

  233. 

  234.     new(this) real(ret);

  235. }

  236. 

  237. template<> real real::operator +() const

  238. {

  239.     return *this;

  240. }

  241. 

  242. template<> real real::operator -() const

  243. {

  244.     real ret = *this;

  245.     ret.m_signexp ^= 0x80000000u;

  246.     return ret;

  247. }

  248. 

  249. template<> real real::operator +(real const &x) const

  250. {

  251.     if (x.m_signexp << 1 == 0)

  252.         return *this;

  253. 

  254.     /* Ensure both arguments are positive. Otherwise, switch signs,

  255.      * or replace + with -. */

  256.     if (m_signexp >> 31)

  257.         return -(-*this + -x);

  258. 

  259.     if (x.m_signexp >> 31)

  260.         return *this - (-x);

  261. 

  262.     /* Ensure *this has the larger exponent (no need for the mantissa to

  263.      * be larger, as in subtraction). Otherwise, switch. */

  264.     if ((m_signexp << 1) < (x.m_signexp << 1))

  265.         return x + *this;

  266. 

  267.     real ret;

  268. 

  269.     int e1 = m_signexp - (1 << 30) + 1;

  270.     int e2 = x.m_signexp - (1 << 30) + 1;

  271. 

  272.     int bigoff = (e1 - e2) / BIGIT_BITS;

  273.     int off = e1 - e2 - bigoff * BIGIT_BITS;

  274. 

  275.     if (bigoff > BIGITS)

  276.         return *this;

  277. 

  278.     ret.m_signexp = m_signexp;

  279. 

  280.     uint64_t carry = 0;

  281.     for (int i = BIGITS; i--; )

  282.     {

  283.         carry += m_mantissa[i];

  284.         if (i - bigoff >= 0)

  285.             carry += x.m_mantissa[i - bigoff] >> off;

  286. 

  287.         if (off && i - bigoff > 0)

  288.             carry += (x.m_mantissa[i - bigoff - 1] << (BIGIT_BITS - off)) & 0xffffffffu;

  289.         else if (i - bigoff == 0)

  290.             carry += (uint64_t)1 << (BIGIT_BITS - off);

  291. 

  292.         ret.m_mantissa[i] = (uint32_t)carry;

  293.         carry >>= BIGIT_BITS;

  294.     }

  295. 

  296.     /* Renormalise in case we overflowed the mantissa */

  297.     if (carry)

  298.     {

  299.         carry--;

  300.         for (int i = 0; i < BIGITS; i++)

  301.         {

  302.             uint32_t tmp = ret.m_mantissa[i];

  303.             ret.m_mantissa[i] = ((uint32_t)carry << (BIGIT_BITS - 1))

  304.                               | (tmp >> 1);

  305.             carry = tmp & 1u;

  306.         }

  307.         ret.m_signexp++;

  308.     }

  309. 

  310.     return ret;

  311. }

  312. 

  313. template<> real real::operator -(real const &x) const

  314. {

  315.     if (x.m_signexp << 1 == 0)

  316.         return *this;

  317. 

  318.     /* Ensure both arguments are positive. Otherwise, switch signs,

  319.      * or replace - with +. */

  320.     if (m_signexp >> 31)

  321.         return -(-*this + x);

  322. 

  323.     if (x.m_signexp >> 31)

  324.         return (*this) + (-x);

  325. 

  326.     /* Ensure *this is larger than x */

  327.     if (*this < x)

  328.         return -(x - *this);

  329. 

  330.     real ret;

  331. 

  332.     int e1 = m_signexp - (1 << 30) + 1;

  333.     int e2 = x.m_signexp - (1 << 30) + 1;

  334. 

  335.     int bigoff = (e1 - e2) / BIGIT_BITS;

  336.     int off = e1 - e2 - bigoff * BIGIT_BITS;

  337. 

  338.     if (bigoff > BIGITS)

  339.         return *this;

  340. 

  341.     ret.m_signexp = m_signexp;

  342. 

  343.     /* int64_t instead of uint64_t to preserve sign through shifts */

  344.     int64_t carry = 0;

  345.     for (int i = 0; i < bigoff; i++)

  346.     {

  347.         carry -= x.m_mantissa[BIGITS - 1 - i];

  348.         /* Emulates a signed shift */

  349.         carry >>= BIGIT_BITS;

  350.         carry |= carry << BIGIT_BITS;

  351.     }

  352.     if (bigoff < BIGITS)

  353.         carry -= x.m_mantissa[BIGITS - 1 - bigoff] & (((int64_t)1 << off) - 1);

  354.     carry /= (int64_t)1 << off;

  355. 

  356.     for (int i = BIGITS; i--; )

  357.     {

  358.         carry += m_mantissa[i];

  359.         if (i - bigoff >= 0)

  360.             carry -= x.m_mantissa[i - bigoff] >> off;

  361. 

  362.         if (off && i - bigoff > 0)

  363.             carry -= (x.m_mantissa[i - bigoff - 1] << (BIGIT_BITS - off)) & 0xffffffffu;

  364.         else if (i - bigoff == 0)

  365.             carry -= (uint64_t)1 << (BIGIT_BITS - off);

  366. 

  367.         ret.m_mantissa[i] = (uint32_t)carry;

  368.         carry >>= BIGIT_BITS;

  369.         carry |= carry << BIGIT_BITS;

  370.     }

  371. 

  372.     carry += 1;

  373. 

  374.     /* Renormalise if we underflowed the mantissa */

  375.     if (carry == 0)

  376.     {

  377.         /* How much do we need to shift the mantissa? FIXME: this could

  378.          * be computed above */

  379.         off = 0;

  380.         for (int i = 0; i < BIGITS; i++)

  381.         {

  382.             if (!ret.m_mantissa[i])

  383.             {

  384.                 off += BIGIT_BITS;

  385.                 continue;

  386.             }

  387. 

  388.             for (uint32_t tmp = ret.m_mantissa[i]; tmp < 0x80000000u; tmp <<= 1)

  389.                 off++;

  390.             break;

  391.         }

  392.         if (off == BIGITS * BIGIT_BITS)

  393.             ret.m_signexp &= 0x80000000u;

  394.         else

  395.         {

  396.             off++; /* Shift one more to get rid of the leading one */

  397.             ret.m_signexp -= off;

  398. 

  399.             bigoff = off / BIGIT_BITS;

  400.             off -= bigoff * BIGIT_BITS;

  401. 

  402.             for (int i = 0; i < BIGITS; i++)

  403.             {

  404.                 uint32_t tmp = 0;

  405.                 if (i + bigoff < BIGITS)

  406.                     tmp |= ret.m_mantissa[i + bigoff] << off;

  407.                 if (off && i + bigoff + 1 < BIGITS)

  408.                     tmp |= ret.m_mantissa[i + bigoff + 1] >> (BIGIT_BITS - off);

  409.                 ret.m_mantissa[i] = tmp;

  410.             }

  411.         }

  412.     }

  413. 

  414.     return ret;

  415. }

  416. 

  417. template<> real real::operator *(real const &x) const

  418. {

  419.     real ret;

  420. 

  421.     if (m_signexp << 1 == 0 || x.m_signexp << 1 == 0)

  422.     {

  423.         ret = (m_signexp << 1 == 0) ? *this : x;

  424.         ret.m_signexp ^= x.m_signexp & 0x80000000u;

  425.         return ret;

  426.     }

  427. 

  428.     ret.m_signexp = (m_signexp ^ x.m_signexp) & 0x80000000u;

  429.     int e = (m_signexp & 0x7fffffffu) - (1 << 30) + 1

  430.           + (x.m_signexp & 0x7fffffffu) - (1 << 30) + 1;

  431. 

  432.     /* Accumulate low order product; no need to store it, we just

  433.      * want the carry value */

  434.     uint64_t carry = 0, hicarry = 0, prev;

  435.     for (int i = 0; i < BIGITS; i++)

  436.     {

  437.         for (int j = 0; j < i + 1; j++)

  438.         {

  439.             prev = carry;

  440.             carry += (uint64_t)m_mantissa[BIGITS - 1 - j]

  441.                    * (uint64_t)x.m_mantissa[BIGITS - 1 + j - i];

  442.             if (carry < prev)

  443.                 hicarry++;

  444.         }

  445.         carry >>= BIGIT_BITS;

  446.         carry |= hicarry << BIGIT_BITS;

  447.         hicarry >>= BIGIT_BITS;

  448.     }

  449. 

  450.     for (int i = 0; i < BIGITS; i++)

  451.     {

  452.         for (int j = i + 1; j < BIGITS; j++)

  453.         {

  454.             prev = carry;

  455.             carry += (uint64_t)m_mantissa[BIGITS - 1 - j]

  456.                    * (uint64_t)x.m_mantissa[j - 1 - i];

  457.             if (carry < prev)

  458.                 hicarry++;

  459.         }

  460.         prev = carry;

  461.         carry += m_mantissa[BIGITS - 1 - i];

  462.         carry += x.m_mantissa[BIGITS - 1 - i];

  463.         if (carry < prev)

  464.             hicarry++;

  465.         ret.m_mantissa[BIGITS - 1 - i] = carry & 0xffffffffu;

  466.         carry >>= BIGIT_BITS;

  467.         carry |= hicarry << BIGIT_BITS;

  468.         hicarry >>= BIGIT_BITS;

  469.     }

  470. 

  471.     /* Renormalise in case we overflowed the mantissa */

  472.     if (carry)

  473.     {

  474.         carry--;

  475.         for (int i = 0; i < BIGITS; i++)

  476.         {

  477.             uint32_t tmp = (uint32_t)ret.m_mantissa[i];

  478.             ret.m_mantissa[i] = ((uint32_t)carry << (BIGIT_BITS - 1))

  479.                               | (tmp >> 1);

  480.             carry = tmp & 1u;

  481.         }

  482.         e++;

  483.     }

  484. 

  485.     ret.m_signexp |= e + (1 << 30) - 1;

  486. 

  487.     return ret;

  488. }

  489. 

  490. template<> real real::operator /(real const &x) const

  491. {

  492.     return *this * re(x);

  493. }

  494. 

  495. template<> real const &real::operator +=(real const &x)

  496. {

  497.     real tmp = *this;

  498.     return *this = tmp + x;

  499. }

  500. 

  501. template<> real const &real::operator -=(real const &x)

  502. {

  503.     real tmp = *this;

  504.     return *this = tmp - x;

  505. }

  506. 

  507. template<> real const &real::operator *=(real const &x)

  508. {

  509.     real tmp = *this;

  510.     return *this = tmp * x;

  511. }

  512. 

  513. template<> real const &real::operator /=(real const &x)

  514. {

  515.     real tmp = *this;

  516.     return *this = tmp / x;

  517. }

  518. 

  519. template<> bool real::operator ==(real const &x) const

  520. {

  521.     if ((m_signexp << 1) == 0 && (x.m_signexp << 1) == 0)

  522.         return true;

  523. 

  524.     if (m_signexp != x.m_signexp)

  525.         return false;

  526. 

  527.     return memcmp(m_mantissa, x.m_mantissa, BIGITS * sizeof(uint32_t)) == 0;

  528. }

  529. 

  530. template<> bool real::operator !=(real const &x) const

  531. {

  532.     return !(*this == x);

  533. }

  534. 

  535. template<> bool real::operator <(real const &x) const

  536. {

  537.     /* Ensure both numbers are positive */

  538.     if (m_signexp >> 31)

  539.         return (x.m_signexp >> 31) ? -*this > -x : true;

  540. 

  541.     if (x.m_signexp >> 31)

  542.         return false;

  543. 

  544.     /* Compare all relevant bits */

  545.     if (m_signexp != x.m_signexp)

  546.         return m_signexp < x.m_signexp;

  547. 

  548.     for (int i = 0; i < BIGITS; i++)

  549.         if (m_mantissa[i] != x.m_mantissa[i])

  550.             return m_mantissa[i] < x.m_mantissa[i];

  551. 

  552.     return false;

  553. }

  554. 

  555. template<> bool real::operator <=(real const &x) const

  556. {

  557.     return !(*this > x);

  558. }

  559. 

  560. template<> bool real::operator >(real const &x) const

  561. {

  562.     /* Ensure both numbers are positive */

  563.     if (m_signexp >> 31)

  564.         return (x.m_signexp >> 31) ? -*this < -x : false;

  565. 

  566.     if (x.m_signexp >> 31)

  567.         return true;

  568. 

  569.     /* Compare all relevant bits */

  570.     if (m_signexp != x.m_signexp)

  571.         return m_signexp > x.m_signexp;

  572. 

  573.     for (int i = 0; i < BIGITS; i++)

  574.         if (m_mantissa[i] != x.m_mantissa[i])

  575.             return m_mantissa[i] > x.m_mantissa[i];

  576. 

  577.     return false;

  578. }

  579. 

  580. template<> bool real::operator >=(real const &x) const

  581. {

  582.     return !(*this < x);

  583. }

  584. 

  585. template<> bool real::operator !() const

  586. {

  587.     return !(bool)*this;

  588. }

  589. 

  590. template<> real::operator bool() const

  591. {

  592.     /* A real is "true" if it is non-zero (exponent is non-zero) AND

  593.      * not NaN (exponent is not full bits OR higher order mantissa is zero) */

  594.     uint32_t exponent = m_signexp << 1;

  595.     return exponent && (~exponent || m_mantissa[0] == 0);

  596. }

  597. 

  598. template<> real min(real const &a, real const &b)

  599. {

  600.     return (a < b) ? a : b;

  601. }

  602. 

  603. template<> real max(real const &a, real const &b)

  604. {

  605.     return (a > b) ? a : b;

  606. }

  607. 

  608. template<> real clamp(real const &x, real const &a, real const &b)

  609. {

  610.     return (x < a) ? a : (x > b) ? b : x;

  611. }

  612. 

  613. template<> real re(real const &x)

  614. {

  615.     if (!(x.m_signexp << 1))

  616.     {

  617.         real ret = x;

  618.         ret.m_signexp = x.m_signexp | 0x7fffffffu;

  619.         ret.m_mantissa[0] = 0;

  620.         return ret;

  621.     }

  622. 

  623.     /* Use the system's float inversion to approximate 1/x */

  624.     union { float f; uint32_t x; } u = { 1.0f }, v = { 1.0f };

  625.     v.x |= x.m_mantissa[0] >> 9;

  626.     v.f = 1.0f / v.f;

  627. 

  628.     real ret;

  629.     ret.m_mantissa[0] = v.x << 9;

  630. 

  631.     uint32_t sign = x.m_signexp & 0x80000000u;

  632.     ret.m_signexp = sign;

  633. 

  634.     int exponent = (x.m_signexp & 0x7fffffffu) + 1;

  635.     exponent = -exponent + (v.x >> 23) - (u.x >> 23);

  636.     ret.m_signexp |= (exponent - 1) & 0x7fffffffu;

  637. 

  638.     /* FIXME: 1+log2(BIGITS) steps of Newton-Raphson seems to be enough for

  639.      * convergence, but this hasn't been checked seriously. */

  640.     for (int i = 1; i <= real::BIGITS; i *= 2)

  641.         ret = ret * (real::R_2() - ret * x);

  642. 

  643.     return ret;

  644. }

  645. 

  646. template<> real sqrt(real const &x)

  647. {

  648.     /* if zero, return x */

  649.     if (!(x.m_signexp << 1))

  650.         return x;

  651. 

  652.     /* if negative, return NaN */

  653.     if (x.m_signexp >> 31)

  654.     {

  655.         real ret;

  656.         ret.m_signexp = 0x7fffffffu;

  657.         ret.m_mantissa[0] = 0xffffu;

  658.         return ret;

  659.     }

  660. 

  661.     /* Use the system's float inversion to approximate 1/sqrt(x). First

  662.      * we construct a float in the [1..4[ range that has roughly the same

  663.      * mantissa as our real. Its exponent is 0 or 1, depending on the

  664.      * partity of x. The final exponent is 0, -1 or -2. We use the final

  665.      * exponent and final mantissa to pre-fill the result. */

  666.     union { float f; uint32_t x; } u = { 1.0f }, v = { 2.0f };

  667.     v.x -= ((x.m_signexp & 1) << 23);

  668.     v.x |= x.m_mantissa[0] >> 9;

  669.     v.f = 1.0f / sqrtf(v.f);

  670. 

  671.     real ret;

  672.     ret.m_mantissa[0] = v.x << 9;

  673. 

  674.     uint32_t sign = x.m_signexp & 0x80000000u;

  675.     ret.m_signexp = sign;

  676. 

  677.     uint32_t exponent = (x.m_signexp & 0x7fffffffu);

  678.     exponent = ((1 << 30) + (1 << 29) - 1) - (exponent + 1) / 2;

  679.     exponent = exponent + (v.x >> 23) - (u.x >> 23);

  680.     ret.m_signexp |= exponent & 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.         ret = ret * (real::R_3() - ret * ret * x);

  687.         ret.m_signexp--;

  688.     }

  689. 

  690.     return ret * x;

  691. }

  692. 

  693. template<> real cbrt(real const &x)

  694. {

  695.     /* if zero, return x */

  696.     if (!(x.m_signexp << 1))

  697.         return x;

  698. 

  699.     /* Use the system's float inversion to approximate cbrt(x). First

  700.      * we construct a float in the [1..8[ range that has roughly the same

  701.      * mantissa as our real. Its exponent is 0, 1 or 2, depending on the

  702.      * value of x. The final exponent is 0 or 1 (special case). We use

  703.      * the final exponent and final mantissa to pre-fill the result. */

  704.     union { float f; uint32_t x; } u = { 1.0f }, v = { 1.0f };

  705.     v.x += ((x.m_signexp % 3) << 23);

  706.     v.x |= x.m_mantissa[0] >> 9;

  707.     v.f = powf(v.f, 0.33333333333333333f);

  708. 

  709.     real ret;

  710.     ret.m_mantissa[0] = v.x << 9;

  711. 

  712.     uint32_t sign = x.m_signexp & 0x80000000u;

  713.     ret.m_signexp = sign;

  714. 

  715.     int exponent = (x.m_signexp & 0x7fffffffu) - (1 << 30) + 1;

  716.     exponent = exponent / 3 + (v.x >> 23) - (u.x >> 23);

  717.     ret.m_signexp |= (exponent + (1 << 30) - 1) & 0x7fffffffu;

  718. 

  719.     /* FIXME: 1+log2(BIGITS) steps of Newton-Raphson seems to be enough for

  720.      * convergence, but this hasn't been checked seriously. */

  721.     for (int i = 1; i <= real::BIGITS; i *= 2)

  722.     {

  723.         static real third = re(real::R_3());

  724.         ret = third * (x / (ret * ret) + (ret / 2));

  725.     }

  726. 

  727.     return ret;

  728. }

  729. 

  730. template<> real pow(real const &x, real const &y)

  731. {

  732.     if (!y)

  733.         return real::R_1();

  734.     if (!x)

  735.         return real::R_0();

  736.     if (x > real::R_0())

  737.         return exp(y * log(x));

  738.     else /* x < 0 */

  739.     {

  740.         /* Odd integer exponent */

  741.         if (y == (round(y / 2) * 2))

  742.             return exp(y * log(-x));

  743. 

  744.         /* Even integer exponent */

  745.         if (y == round(y))

  746.             return -exp(y * log(-x));

  747. 

  748.         /* FIXME: negative nth root */

  749.         return real::R_0();

  750.     }

  751. }

  752. 

  753. static real fast_fact(unsigned int x)

  754. {

  755.     real ret = real::R_1();

  756.     unsigned int i = 1, multiplier = 1, exponent = 0;

  757. 

  758.     for (;;)

  759.     {

  760.         if (i++ >= x)

  761.             /* Multiplication is a no-op if multiplier == 1 */

  762.             return ldexp(ret * multiplier, exponent);

  763. 

  764.         unsigned int tmp = i;

  765.         while ((tmp & 1) == 0)

  766.         {

  767.             tmp >>= 1;

  768.             exponent++;

  769.         }

  770.         if (multiplier * tmp / tmp != multiplier)

  771.         {

  772.             ret *= multiplier;

  773.             multiplier = 1;

  774.         }

  775.         multiplier *= tmp;

  776.     }

  777. }

  778. 

  779. template<> real gamma(real const &x)

  780. {

  781.     /* We use Spouge's formula. FIXME: precision is far from acceptable,

  782.      * especially with large values. We need to compute this with higher

  783.      * precision values in order to attain the desired accuracy. It might

  784.      * also be useful to sort the ck values by decreasing absolute value

  785.      * and do the addition in this order. */

  786.     int a = (int)ceilf(logf(2) / logf(2 * F_PI)

  787.                           * real::BIGITS * real::BIGIT_BITS);

  788. 

  789.     real ret = sqrt(real::R_PI() * 2);

  790.     real fact_k_1 = real::R_1();

  791. 

  792.     for (int k = 1; k < a; k++)

  793.     {

  794.         real a_k = (real)(a - k);

  795.         real ck = pow(a_k, (real)((float)k - 0.5)) * exp(a_k)

  796.                 / (fact_k_1 * (x + (real)(k - 1)));

  797.         ret += ck;

  798.         fact_k_1 *= (real)-k;

  799.     }

  800. 

  801.     ret *= pow(x + (real)(a - 1), x - (real::R_1() / 2));

  802.     ret *= exp(-x - (real)(a - 1));

  803. 

  804.     return ret;

  805. }

  806. 

  807. template<> real fabs(real const &x)

  808. {

  809.     real ret = x;

  810.     ret.m_signexp &= 0x7fffffffu;

  811.     return ret;

  812. }

  813. 

  814. template<> real abs(real const &x)

  815. {

  816.     return fabs(x);

  817. }

  818. 

  819. template<> real fract(real const &x)

  820. {

  821.     return x - floor(x);

  822. }

  823. 

  824. template<> real degrees(real const &x)

  825. {

  826.     static real mul = real(180) * real::R_1_PI();

  827. 

  828.     return x * mul;

  829. }

  830. 

  831. template<> real radians(real const &x)

  832. {

  833.     static real mul = real::R_PI() / real(180);

  834. 

  835.     return x * mul;

  836. }

  837. 

  838. static real fast_log(real const &x)

  839. {

  840.     /* This fast log method is tuned to work on the [1..2] range and

  841.      * no effort whatsoever was made to improve convergence outside this

  842.      * domain of validity. It can converge pretty fast, provided we use

  843.      * the following variable substitutions:

  844.      *    y = sqrt(x)

  845.      *    z = (y - 1) / (y + 1)

  846.      *

  847.      * And the following identities:

  848.      *    ln(x) = 2 ln(y)

  849.      *          = 2 ln((1 + z) / (1 - z))

  850.      *          = 4 z (1 + z^2 / 3 + z^4 / 5 + z^6 / 7...)

  851.      *

  852.      * Any additional sqrt() call would halve the convergence time, but

  853.      * would also impact the final precision. For now we stick with one

  854.      * sqrt() call. */

  855.     real y = sqrt(x);

  856.     real z = (y - real::R_1()) / (y + real::R_1()), z2 = z * z, zn = z2;

  857.     real sum = real::R_1();

  858. 

  859.     for (int i = 3; ; i += 2)

  860.     {

  861.         real newsum = sum + zn / (real)i;

  862.         if (newsum == sum)

  863.             break;

  864.         sum = newsum;

  865.         zn *= z2;

  866.     }

  867. 

  868.     return z * sum * 4;

  869. }

  870. 

  871. template<> real log(real const &x)

  872. {

  873.     /* Strategy for log(x): if x = 2^E*M then log(x) = E log(2) + log(M),

  874.      * with the property that M is in [1..2[, so fast_log() applies here. */

  875.     real tmp = x;

  876.     if (x.m_signexp >> 31 || x.m_signexp == 0)

  877.     {

  878.         tmp.m_signexp = 0xffffffffu;

  879.         tmp.m_mantissa[0] = 0xffffffffu;

  880.         return tmp;

  881.     }

  882.     tmp.m_signexp = (1 << 30) - 1;

  883.     return (real)(int)(x.m_signexp - (1 << 30) + 1) * real::R_LN2()

  884.            + fast_log(tmp);

  885. }

  886. 

  887. template<> real log2(real const &x)

  888. {

  889.     /* Strategy for log2(x): see log(x). */

  890.     real tmp = x;

  891.     if (x.m_signexp >> 31 || x.m_signexp == 0)

  892.     {

  893.         tmp.m_signexp = 0xffffffffu;

  894.         tmp.m_mantissa[0] = 0xffffffffu;

  895.         return tmp;

  896.     }

  897.     tmp.m_signexp = (1 << 30) - 1;

  898.     return (real)(int)(x.m_signexp - (1 << 30) + 1)

  899.            + fast_log(tmp) * real::R_LOG2E();

  900. }

  901. 

  902. template<> real log10(real const &x)

  903. {

  904.     return log(x) * real::R_LOG10E();

  905. }

  906. 

  907. static real fast_exp_sub(real const &x, real const &y)

  908. {

  909.     /* This fast exp method is tuned to work on the [-1..1] range and

  910.      * no effort whatsoever was made to improve convergence outside this

  911.      * domain of validity. The argument y is used for cases where we

  912.      * don't want the leading 1 in the Taylor series. */

  913.     real ret = real::R_1() - y, xn = x;

  914.     int i = 1;

  915. 

  916.     for (;;)

  917.     {

  918.         real newret = ret + xn;

  919.         if (newret == ret)

  920.             break;

  921.         ret = newret * ++i;

  922.         xn *= x;

  923.     }

  924. 

  925.     return ret / fast_fact(i);

  926. }

  927. 

  928. template<> real exp(real const &x)

  929. {

  930.     /* Strategy for exp(x): the Taylor series does not converge very fast

  931.      * with large positive or negative values.

  932.      *

  933.      * However, we know that the result is going to be in the form M*2^E,

  934.      * where M is the mantissa and E the exponent. We first try to predict

  935.      * a value for E, which is approximately log2(exp(x)) = x / log(2).

  936.      *

  937.      * Let E0 be an integer close to x / log(2). We need to find a value x0

  938.      * such that exp(x) = 2^E0 * exp(x0). We get x0 = x - E0 log(2).

  939.      *

  940.      * Thus the final algorithm:

  941.      *  int E0 = x / log(2)

  942.      *  real x0 = x - E0 log(2)

  943.      *  real x1 = exp(x0)

  944.      *  return x1 * 2^E0

  945.      */

  946.     int e0 = x / real::R_LN2();

  947.     real x0 = x - (real)e0 * real::R_LN2();

  948.     real x1 = fast_exp_sub(x0, real::R_0());

  949.     x1.m_signexp += e0;

  950.     return x1;

  951. }

  952. 

  953. template<> real exp2(real const &x)

  954. {

  955.     /* Strategy for exp2(x): see strategy in exp(). */

  956.     int e0 = x;

  957.     real x0 = x - (real)e0;

  958.     real x1 = fast_exp_sub(x0 * real::R_LN2(), real::R_0());

  959.     x1.m_signexp += e0;

  960.     return x1;

  961. }

  962. 

  963. template<> real sinh(real const &x)

  964. {

  965.     /* We cannot always use (exp(x)-exp(-x))/2 because we'll lose

  966.      * accuracy near zero. We only use this identity for |x|>0.5. If

  967.      * |x|<=0.5, we compute exp(x)-1 and exp(-x)-1 instead. */

  968.     bool near_zero = (fabs(x) < real::R_1() / 2);

  969.     real x1 = near_zero ? fast_exp_sub(x, real::R_1()) : exp(x);

  970.     real x2 = near_zero ? fast_exp_sub(-x, real::R_1()) : exp(-x);

  971.     return (x1 - x2) / 2;

  972. }

  973. 

  974. template<> real tanh(real const &x)

  975. {

  976.     /* See sinh() for the strategy here */

  977.     bool near_zero = (fabs(x) < real::R_1() / 2);

  978.     real x1 = near_zero ? fast_exp_sub(x, real::R_1()) : exp(x);

  979.     real x2 = near_zero ? fast_exp_sub(-x, real::R_1()) : exp(-x);

  980.     real x3 = near_zero ? x1 + x2 + real::R_2() : x1 + x2;

  981.     return (x1 - x2) / x3;

  982. }

  983. 

  984. template<> real cosh(real const &x)

  985. {

  986.     /* No need to worry about accuracy here; maybe the last bit is slightly

  987.      * off, but that's about it. */

  988.     return (exp(x) + exp(-x)) / 2;

  989. }

  990. 

  991. template<> real frexp(real const &x, int *exp)

  992. {

  993.     if (!x)

  994.     {

  995.         *exp = 0;

  996.         return x;

  997.     }

  998. 

  999.     real ret = x;

  1000.     int exponent = (ret.m_signexp & 0x7fffffffu) - (1 << 30) + 1;

  1001.     *exp = exponent + 1;

  1002.     ret.m_signexp -= exponent + 1;

  1003.     return ret;

  1004. }

  1005. 

  1006. template<> real ldexp(real const &x, int exp)

  1007. {

  1008.     real ret = x;

  1009.     if (ret)

  1010.         ret.m_signexp += exp;

  1011.     return ret;

  1012. }

  1013. 

  1014. template<> real modf(real const &x, real *iptr)

  1015. {

  1016.     real absx = fabs(x);

  1017.     real tmp = floor(absx);

  1018. 

  1019.     *iptr = copysign(tmp, x);

  1020.     return copysign(absx - tmp, x);

  1021. }

  1022. 

  1023. template<> real ulp(real const &x)

  1024. {

  1025.     real ret = real::R_1();

  1026.     if (x)

  1027.         ret.m_signexp = x.m_signexp + 1 - real::BIGITS * real::BIGIT_BITS;

  1028.     else

  1029.         ret.m_signexp = 0;

  1030.     return ret;

  1031. }

  1032. 

  1033. template<> real nextafter(real const &x, real const &y)

  1034. {

  1035.     if (x == y)

  1036.         return x;

  1037.     else if (x < y)

  1038.         return x + ulp(x);

  1039.     else

  1040.         return x - ulp(x);

  1041. }

  1042. 

  1043. template<> real copysign(real const &x, real const &y)

  1044. {

  1045.     real ret = x;

  1046.     ret.m_signexp &= 0x7fffffffu;

  1047.     ret.m_signexp |= y.m_signexp & 0x80000000u;

  1048.     return ret;

  1049. }

  1050. 

  1051. template<> real floor(real const &x)

  1052. {

  1053.     /* Strategy for floor(x):

  1054.      *  - if negative, return -ceil(-x)

  1055.      *  - if zero or negative zero, return x

  1056.      *  - if less than one, return zero

  1057.      *  - otherwise, if e is the exponent, clear all bits except the

  1058.      *    first e. */

  1059.     if (x < -real::R_0())

  1060.         return -ceil(-x);

  1061.     if (!x)

  1062.         return x;

  1063.     if (x < real::R_1())

  1064.         return real::R_0();

  1065. 

  1066.     real ret = x;

  1067.     int exponent = x.m_signexp - (1 << 30) + 1;

  1068. 

  1069.     for (int i = 0; i < real::BIGITS; i++)

  1070.     {

  1071.         if (exponent <= 0)

  1072.             ret.m_mantissa[i] = 0;

  1073.         else if (exponent < real::BIGIT_BITS)

  1074.             ret.m_mantissa[i] &= ~((1 << (real::BIGIT_BITS - exponent)) - 1);

  1075. 

  1076.         exponent -= real::BIGIT_BITS;

  1077.     }

  1078. 

  1079.     return ret;

  1080. }

  1081. 

  1082. template<> real ceil(real const &x)

  1083. {

  1084.     /* Strategy for ceil(x):

  1085.      *  - if negative, return -floor(-x)

  1086.      *  - if x == floor(x), return x

  1087.      *  - otherwise, return floor(x) + 1 */

  1088.     if (x < -real::R_0())

  1089.         return -floor(-x);

  1090.     real ret = floor(x);

  1091.     if (x == ret)

  1092.         return ret;

  1093.     else

  1094.         return ret + real::R_1();

  1095. }

  1096. 

  1097. template<> real round(real const &x)

  1098. {

  1099.     if (x < real::R_0())

  1100.         return -round(-x);

  1101. 

  1102.     return floor(x + (real::R_1() / 2));

  1103. }

  1104. 

  1105. template<> real fmod(real const &x, real const &y)

  1106. {

  1107.     if (!y)

  1108.         return real::R_0(); /* FIXME: return NaN */

  1109. 

  1110.     if (!x)

  1111.         return x;

  1112. 

  1113.     real tmp = round(x / y);

  1114.     return x - tmp * y;

  1115. }

  1116. 

  1117. template<> real sin(real const &x)

  1118. {

  1119.     int switch_sign = x.m_signexp & 0x80000000u;

  1120. 

  1121.     real absx = fmod(fabs(x), real::R_PI() * 2);

  1122.     if (absx > real::R_PI())

  1123.     {

  1124.         absx -= real::R_PI();

  1125.         switch_sign = !switch_sign;

  1126.     }

  1127. 

  1128.     if (absx > real::R_PI_2())

  1129.         absx = real::R_PI() - absx;

  1130. 

  1131.     real ret = real::R_0(), fact = real::R_1(), xn = absx, mx2 = -absx * absx;

  1132.     int i = 1;

  1133.     for (;;)

  1134.     {

  1135.         real newret = ret + xn;

  1136.         if (newret == ret)

  1137.             break;

  1138.         ret = newret * ((i + 1) * (i + 2));

  1139.         xn *= mx2;

  1140.         i += 2;

  1141.     }

  1142.     ret /= fast_fact(i);

  1143. 

  1144.     /* Propagate sign */

  1145.     if (switch_sign)

  1146.         ret.m_signexp ^= 0x80000000u;

  1147.     return ret;

  1148. }

  1149. 

  1150. template<> real cos(real const &x)

  1151. {

  1152.     return sin(real::R_PI_2() - x);

  1153. }

  1154. 

  1155. template<> real tan(real const &x)

  1156. {

  1157.     /* Constrain input to [-π,π] */

  1158.     real y = fmod(x, real::R_PI());

  1159. 

  1160.     /* Constrain input to [-π/2,π/2] */

  1161.     if (y < -real::R_PI_2())

  1162.         y += real::R_PI();

  1163.     else if (y > real::R_PI_2())

  1164.         y -= real::R_PI();

  1165. 

  1166.     /* In [-π/4,π/4] return sin/cos */

  1167.     if (fabs(y) <= real::R_PI_4())

  1168.         return sin(y) / cos(y);

  1169. 

  1170.     /* Otherwise, return cos/sin */

  1171.     if (y > real::R_0())

  1172.         y = real::R_PI_2() - y;

  1173.     else

  1174.         y = -real::R_PI_2() - y;

  1175. 

  1176.     return cos(y) / sin(y);

  1177. }

  1178. 

  1179. static inline real asinacos(real const &x, int is_asin, int is_negative)

  1180. {

  1181.     /* Strategy for asin(): in [-0.5..0.5], use a Taylor series around

  1182.      * zero. In [0.5..1], use asin(x) = π/2 - 2*asin(sqrt((1-x)/2)), and

  1183.      * in [-1..-0.5] just revert the sign.

  1184.      * Strategy for acos(): use acos(x) = π/2 - asin(x) and try not to

  1185.      * lose the precision around x=1. */

  1186.     real absx = fabs(x);

  1187.     int around_zero = (absx < (real::R_1() / 2));

  1188. 

  1189.     if (!around_zero)

  1190.         absx = sqrt((real::R_1() - absx) / 2);

  1191. 

  1192.     real ret = absx, xn = absx, x2 = absx * absx, fact1 = 2, fact2 = 1;

  1193.     for (int i = 1; ; i++)

  1194.     {

  1195.         xn *= x2;

  1196.         real mul = (real)(2 * i + 1);

  1197.         real newret = ret + ldexp(fact1 * xn / (mul * fact2), -2 * i);

  1198.         if (newret == ret)

  1199.             break;

  1200.         ret = newret;

  1201.         fact1 *= (real)((2 * i + 1) * (2 * i + 2));

  1202.         fact2 *= (real)((i + 1) * (i + 1));

  1203.     }

  1204. 

  1205.     if (is_negative)

  1206.         ret = -ret;

  1207. 

  1208.     if (around_zero)

  1209.         ret = is_asin ? ret : real::R_PI_2() - ret;

  1210.     else

  1211.     {

  1212.         real adjust = is_negative ? real::R_PI() : real::R_0();

  1213.         if (is_asin)

  1214.             ret = real::R_PI_2() - adjust - ret * 2;

  1215.         else

  1216.             ret = adjust + ret * 2;

  1217.     }

  1218. 

  1219.     return ret;

  1220. }

  1221. 

  1222. template<> real asin(real const &x)

  1223. {

  1224.     return asinacos(x, 1, x.m_signexp >> 31);

  1225. }

  1226. 

  1227. template<> real acos(real const &x)

  1228. {

  1229.     return asinacos(x, 0, x.m_signexp >> 31);

  1230. }

  1231. 

  1232. template<> real atan(real const &x)

  1233. {

  1234.     /* Computing atan(x): we choose a different Taylor series depending on

  1235.      * the value of x to help with convergence.

  1236.      *

  1237.      * If |x| < 0.5 we evaluate atan(y) near 0:

  1238.      *  atan(y) = y - y^3/3 + y^5/5 - y^7/7 + y^9/9 ...

  1239.      *

  1240.      * If 0.5 <= |x| < 1.5 we evaluate atan(1+y) near 0:

  1241.      *  atan(1+y) = π/4 + y/(1*2^1) - y^2/(2*2^1) + y^3/(3*2^2)

  1242.      *                  - y^5/(5*2^3) + y^6/(6*2^3) - y^7/(7*2^4)

  1243.      *                  + y^9/(9*2^5) - y^10/(10*2^5) + y^11/(11*2^6) ...

  1244.      *

  1245.      * If 1.5 <= |x| < 2 we evaluate atan(sqrt(3)+2y) near 0:

  1246.      *  atan(sqrt(3)+2y) = π/3 + 1/2 y - sqrt(3)/2 y^2/2

  1247.      *                         + y^3/3 - sqrt(3)/2 y^4/4 + 1/2 y^5/5

  1248.      *                         - 1/2 y^7/7 + sqrt(3)/2 y^8/8

  1249.      *                         - y^9/9 + sqrt(3)/2 y^10/10 - 1/2 y^11/11

  1250.      *                         + 1/2 y^13/13 - sqrt(3)/2 y^14/14

  1251.      *                         + y^15/15 - sqrt(3)/2 y^16/16 + 1/2 y^17/17 ...

  1252.      *

  1253.      * If |x| >= 2 we evaluate atan(y) near +∞:

  1254.      *  atan(y) = π/2 - y^-1 + y^-3/3 - y^-5/5 + y^-7/7 - y^-9/9 ...

  1255.      */

  1256.     real absx = fabs(x);

  1257. 

  1258.     if (absx < (real::R_1() / 2))

  1259.     {

  1260.         real ret = x, xn = x, mx2 = -x * x;

  1261.         for (int i = 3; ; i += 2)

  1262.         {

  1263.             xn *= mx2;

  1264.             real newret = ret + xn / (real)i;

  1265.             if (newret == ret)

  1266.                 break;

  1267.             ret = newret;

  1268.         }

  1269.         return ret;

  1270.     }

  1271. 

  1272.     real ret = 0;

  1273. 

  1274.     if (absx < (real::R_3() / 2))

  1275.     {

  1276.         real y = real::R_1() - absx;

  1277.         real yn = y, my2 = -y * y;

  1278.         for (int i = 0; ; i += 2)

  1279.         {

  1280.             real newret = ret + ldexp(yn / (real)(2 * i + 1), -i - 1);

  1281.             yn *= y;

  1282.             newret += ldexp(yn / (real)(2 * i + 2), -i - 1);

  1283.             yn *= y;

  1284.             newret += ldexp(yn / (real)(2 * i + 3), -i - 2);

  1285.             if (newret == ret)

  1286.                 break;

  1287.             ret = newret;

  1288.             yn *= my2;

  1289.         }

  1290.         ret = real::R_PI_4() - ret;

  1291.     }

  1292.     else if (absx < real::R_2())

  1293.     {

  1294.         real y = (absx - real::R_SQRT3()) / 2;

  1295.         real yn = y, my2 = -y * y;

  1296.         for (int i = 1; ; i += 6)

  1297.         {

  1298.             real newret = ret + ((yn / (real)i) / 2);

  1299.             yn *= y;

  1300.             newret -= (real::R_SQRT3() / 2) * yn / (real)(i + 1);

  1301.             yn *= y;

  1302.             newret += yn / (real)(i + 2);

  1303.             yn *= y;

  1304.             newret -= (real::R_SQRT3() / 2) * yn / (real)(i + 3);

  1305.             yn *= y;

  1306.             newret += (yn / (real)(i + 4)) / 2;

  1307.             if (newret == ret)

  1308.                 break;

  1309.             ret = newret;

  1310.             yn *= my2;

  1311.         }

  1312.         ret = real::R_PI_3() + ret;

  1313.     }

  1314.     else

  1315.     {

  1316.         real y = re(absx);

  1317.         real yn = y, my2 = -y * y;

  1318.         ret = y;

  1319.         for (int i = 3; ; i += 2)

  1320.         {

  1321.             yn *= my2;

  1322.             real newret = ret + yn / (real)i;

  1323.             if (newret == ret)

  1324.                 break;

  1325.             ret = newret;

  1326.         }

  1327.         ret = real::R_PI_2() - ret;

  1328.     }

  1329. 

  1330.     /* Propagate sign */

  1331.     ret.m_signexp |= (x.m_signexp & 0x80000000u);

  1332.     return ret;

  1333. }

  1334. 

  1335. template<> real atan2(real const &y, real const &x)

  1336. {

  1337.     if (!y)

  1338.     {

  1339.         if ((x.m_signexp >> 31) == 0)

  1340.             return y;

  1341.         if (y.m_signexp >> 31)

  1342.             return -real::R_PI();

  1343.         return real::R_PI();

  1344.     }

  1345. 

  1346.     if (!x)

  1347.     {

  1348.         if (y.m_signexp >> 31)

  1349.             return -real::R_PI();

  1350.         return real::R_PI();

  1351.     }

  1352. 

  1353.     /* FIXME: handle the Inf and NaN cases */

  1354.     real z = y / x;

  1355.     real ret = atan(z);

  1356.     if (x < real::R_0())

  1357.         ret += (y > real::R_0()) ? real::R_PI() : -real::R_PI();

  1358.     return ret;

  1359. }

  1360. 

  1361. template<> void real::hexprint() const

  1362. {

  1363.     std::printf("%08x", m_signexp);

  1364.     for (int i = 0; i < BIGITS; i++)

  1365.         std::printf(" %08x", m_mantissa[i]);

  1366. }

  1367. 

  1368. template<> void real::sprintf(char *str, int ndigits) const;

  1369. 

  1370. template<> void real::print(int ndigits) const

  1371. {

  1372.     char *buf = new char[ndigits + 32 + 10];

  1373.     real::sprintf(buf, ndigits);

  1374.     std::printf("%s", buf);

  1375.     delete[] buf;

  1376. }

  1377. 

  1378. template<> void real::sprintf(char *str, int ndigits) const

  1379. {

  1380.     real x = *this;

  1381. 

  1382.     if (x.m_signexp >> 31)

  1383.     {

  1384.         *str++ = '-';

  1385.         x = -x;

  1386.     }

  1387. 

  1388.     if (!x)

  1389.     {

  1390.         std::strcpy(str, "0.0");

  1391.         return;

  1392.     }

  1393. 

  1394.     /* Normalise x so that mantissa is in [1..9.999] */

  1395.     /* FIXME: better use int64_t when the cast is implemented */

  1396.     int exponent = ceil(log10(x));

  1397.     x /= pow(R_10(), (real)exponent);

  1398. 

  1399.     if (x < R_1())

  1400.     {

  1401.         x *= R_10();

  1402.         exponent--;

  1403.     }

  1404. 

  1405.     /* Print digits */

  1406.     for (int i = 0; i < ndigits; i++)

  1407.     {

  1408.         int digit = (int)floor(x);

  1409.         *str++ = '0' + digit;

  1410.         if (i == 0)

  1411.             *str++ = '.';

  1412.         x -= real(digit);

  1413.         x *= R_10();

  1414.     }

  1415. 

  1416.     /* Print exponent information */

  1417.     if (exponent)

  1418.         str += std::sprintf(str, "e%c%i",

  1419.                             exponent >= 0 ? '+' : '-', lol::abs(exponent));

  1420. 

  1421.     *str++ = '\0';

  1422. }

  1423. 

  1424. static real fast_pi()

  1425. {

  1426.     /* Approximate Pi using Machin's formula: 16*atan(1/5)-4*atan(1/239) */

  1427.     real ret = 0, x0 = 5, x1 = 239;

  1428.     real const m0 = -x0 * x0, m1 = -x1 * x1, r16 = 16, r4 = 4;

  1429. 

  1430.     for (int i = 1; ; i += 2)

  1431.     {

  1432.         real newret = ret + r16 / (x0 * (real)i) - r4 / (x1 * (real)i);

  1433.         if (newret == ret)

  1434.             break;

  1435.         ret = newret;

  1436.         x0 *= m0;

  1437.         x1 *= m1;

  1438.     }

  1439. 

  1440.     return ret;

  1441. }

  1442. 

  1443. } /* namespace lol */