source: trunk/src/real.cpp @ 974

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

core: implement real subtraction.

File size: 8.3 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)
26{
27    union { float f; uint32_t x; } u = { f };
28
29    uint32_t sign = u.x & 0x80000000u;
30    uint32_t exponent = (u.x >> 23) & 0xff;
31
32    switch (exponent)
33    {
34    case 0x00:
35    case 0xff:
36        m_signexp = sign | exponent;
37        break;
38    default:
39        m_signexp = sign | (exponent + (1 << 30) - (1 << 7));
40        break;
41    }
42
43    m_mantissa[0] = u.x >> 7;
44    m_mantissa[1] = u.x << 9;
45    memset(m_mantissa + 2, 0, sizeof(m_mantissa) - sizeof(m_mantissa[0]));
46}
47
48real::operator float() const
49{
50    union { float f; uint32_t x; } u;
51
52    uint32_t sign = m_signexp & 0x80000000u;
53    uint32_t exponent = m_signexp & 0x7fffffffu;
54    uint32_t mantissa = (m_mantissa[0] << 7) | (m_mantissa[1] >> 9);
55
56    int e = (int)exponent - (1 << 30) + (1 << 7);
57
58    if (e < 0)
59        u.x = sign;
60    else if (e >= 0xff)
61        u.x = sign | (0xff << 23);
62    else
63        u.x = sign | (e << 23) | mantissa;
64
65    return u.f;
66}
67
68real real::operator -() const
69{
70    real ret = *this;
71    ret.m_signexp ^= 0x80000000u;
72    return ret;
73}
74
75real real::operator +(real const &x) const
76{
77    if (x.m_signexp << 1 == 0)
78        return *this;
79
80    /* Ensure both arguments are positive. Otherwise, switch signs,
81     * or replace + with -). */
82    if (m_signexp >> 31)
83        return -(-*this + -x);
84
85    if (x.m_signexp >> 31)
86        return *this - x;
87
88    /* Ensure *this is the larger exponent (no need to be strictly larger,
89     * as in subtraction). Otherwise, switch. */
90    if ((m_signexp << 1) < (x.m_signexp << 1))
91        return x + *this;
92
93    real ret;
94
95    int e1 = m_signexp - (1 << 30) + 1;
96    int e2 = x.m_signexp - (1 << 30) + 1;
97
98    int bigoff = (e1 - e2) / (sizeof(uint16_t) * 8);
99    int off = e1 - e2 - bigoff * (sizeof(uint16_t) * 8);
100
101    ret.m_signexp = m_signexp;
102
103    uint32_t carry = 0;
104    for (int i = BIGITS; i--; )
105    {
106        carry += m_mantissa[i];
107        if (i - bigoff >= 0)
108            carry += x.m_mantissa[i - bigoff] >> off;
109        else if (i - bigoff == -1)
110            carry += 0x0001u >> off;
111
112        if (i - bigoff > 0)
113            carry += (x.m_mantissa[i - bigoff - 1] << (16 - off)) & 0xffffu;
114        else if (i - bigoff == 0)
115            carry += 0x0001u << (16 - off);
116
117        ret.m_mantissa[i] = carry;
118        carry >>= 16;
119    }
120
121    /* Renormalise in case we overflowed the mantissa */
122    if (carry)
123    {
124        carry--;
125        for (int i = 0; i < BIGITS; i++)
126        {
127            uint16_t tmp = ret.m_mantissa[i];
128            ret.m_mantissa[i] = (carry << 15) | (tmp >> 1);
129            carry = tmp & 0x0001u;
130        }
131        ret.m_signexp++;
132    }
133
134    return ret;
135}
136
137real real::operator -(real const &x) const
138{
139    if (x.m_signexp << 1 == 0)
140        return *this;
141
142    /* Ensure both arguments are positive. Otherwise, switch signs,
143     * or replace - with +). */
144    if (m_signexp >> 31)
145        return -(-*this + x);
146
147    if (x.m_signexp >> 31)
148        return (*this) + (-x);
149
150    /* Ensure *this is larger than x */
151    if (*this < x)
152        return -(x - *this);
153
154    real ret;
155
156    int e1 = m_signexp - (1 << 30) + 1;
157    int e2 = x.m_signexp - (1 << 30) + 1;
158
159    int bigoff = (e1 - e2) / (sizeof(uint16_t) * 8);
160    int off = e1 - e2 - bigoff * (sizeof(uint16_t) * 8);
161
162    ret.m_signexp = m_signexp;
163
164    int32_t carry = 0;
165    for (int i = 0; i < bigoff; i++)
166    {
167        carry -= x.m_mantissa[BIGITS - i];
168        carry = (carry & 0xffff0000u) | (carry >> 16);
169    }
170    carry -= x.m_mantissa[BIGITS - 1 - bigoff] & ((1 << off) - 1);
171    carry /= (1 << off);
172
173    for (int i = BIGITS; i--; )
174    {
175        carry += m_mantissa[i];
176        if (i - bigoff >= 0)
177            carry -= x.m_mantissa[i - bigoff] >> off;
178        else if (i - bigoff == -1)
179            carry -= 0x0001u >> off;
180
181        if (i - bigoff > 0)
182            carry -= (x.m_mantissa[i - bigoff - 1] << (16 - off)) & 0xffffu;
183        else if (i - bigoff == 0)
184            carry -= 0x0001u << (16 - off);
185
186        ret.m_mantissa[i] = carry;
187        carry = (carry & 0xffff0000u) | (carry >> 16);
188    }
189
190    carry += 1;
191
192    /* Renormalise if we underflowed the mantissa */
193    if (carry == 0)
194    {
195        /* How much do we need to shift the mantissa? FIXME: this could
196         * be computed above */
197        off = 0;
198        for (int i = 0; i < BIGITS; i++)
199        {
200            if (!ret.m_mantissa[i])
201            {
202                off += sizeof(uint16_t) * 8;
203                continue;
204            }
205
206            for (uint16_t tmp = ret.m_mantissa[i]; tmp < 0x8000u; tmp <<= 1)
207                off++;
208            break;
209        }
210        if (off == BIGITS * sizeof(uint16_t) * 8)
211            ret.m_signexp &= 0x80000000u;
212        else
213        {
214            off++; /* Shift one more to get rid of the leading one */
215            ret.m_signexp -= off;
216
217            bigoff = off / (sizeof(uint16_t) * 8);
218            off -= bigoff * sizeof(uint16_t) * 8;
219
220            for (int i = 0; i < BIGITS; i++)
221            {
222                uint16_t tmp = 0;
223                if (i + bigoff < BIGITS)
224                    tmp |= ret.m_mantissa[i + bigoff] << off;
225                if (i + bigoff + 1 < BIGITS)
226                    tmp |= ret.m_mantissa[i + bigoff + 1] >> (16 - off);
227                ret.m_mantissa[i] = tmp;
228            }
229        }
230    }
231
232    return ret;
233}
234
235real real::operator *(real const &x) const
236{
237    real ret;
238
239    ret.m_signexp = (m_signexp ^ x.m_signexp) & 0x80000000u;
240    int e = (m_signexp & 0x7fffffffu) - (1 << 30) + 1
241          + (x.m_signexp & 0x7fffffffu) - (1 << 30) + 1;
242
243    /* Accumulate low order product; no need to store it, we just
244     * want the carry value */
245    uint32_t carry = 0;
246    for (int i = 0; i < BIGITS; i++)
247    {
248        for (int j = 0; j < i + 1; j++)
249            carry += m_mantissa[BIGITS - 1 - j]
250                   * x.m_mantissa[BIGITS - 1 + j - i];
251        carry >>= 16;
252    }
253
254    for (int i = 0; i < BIGITS; i++)
255    {
256        for (int j = i + 1; j < BIGITS; j++)
257            carry += m_mantissa[BIGITS - 1 - j]
258                   * x.m_mantissa[j - 1 - i];
259
260        carry += m_mantissa[BIGITS - 1 - i];
261        carry += x.m_mantissa[BIGITS - 1 - i];
262        ret.m_mantissa[BIGITS - 1 - i] = carry & 0xffffu;
263        carry >>= 16;
264    }
265
266    /* Renormalise in case we overflowed the mantissa */
267    if (carry)
268    {
269        carry--;
270        for (int i = 0; i < BIGITS; i++)
271        {
272            uint16_t tmp = ret.m_mantissa[i];
273            ret.m_mantissa[i] = (carry << 15) | (tmp >> 1);
274            carry = tmp & 0x0001u;
275        }
276        e++;
277    }
278
279    ret.m_signexp |= e + (1 << 30) - 1;
280
281    return ret;
282}
283
284bool real::operator <(real const &x) const
285{
286    /* Ensure both numbers are positive */
287    if (m_signexp >> 31)
288        return (x.m_signexp >> 31) ? -*this > -x : true;
289
290    if (x.m_signexp >> 31)
291        return false;
292
293    /* Compare all relevant bits */
294    if (m_signexp != x.m_signexp)
295        return m_signexp < x.m_signexp;
296
297    for (int i = 0; i < BIGITS; i++)
298        if (m_mantissa[i] != x.m_mantissa[i])
299            return m_mantissa[i] < x.m_mantissa[i];
300
301    return false;
302}
303
304bool real::operator <=(real const &x) const
305{
306    return !(*this > x);
307}
308
309bool real::operator >(real const &x) const
310{
311    /* Ensure both numbers are positive */
312    if (m_signexp >> 31)
313        return (x.m_signexp >> 31) ? -*this < -x : false;
314
315    if (x.m_signexp >> 31)
316        return true;
317
318    /* Compare all relevant bits */
319    if (m_signexp != x.m_signexp)
320        return m_signexp > x.m_signexp;
321
322    for (int i = 0; i < BIGITS; i++)
323        if (m_mantissa[i] != x.m_mantissa[i])
324            return m_mantissa[i] > x.m_mantissa[i];
325
326    return false;
327}
328
329bool real::operator >=(real const &x) const
330{
331    return !(*this < x);
332}
333
334void real::print() const
335{
336    printf("%x  %08x  ", m_signexp >> 31, (m_signexp << 1) >> 1);
337    for (int i = 0; i < BIGITS; i++)
338        printf("%04x ", m_mantissa[i]);
339    printf("\n");
340}
341
342} /* namespace lol */
343
Note: See TracBrowser for help on using the repository browser.