Changeset 1016


Ignore:
Timestamp:
Oct 11, 2011, 12:51:26 AM (11 years ago)
Author:
sam
Message:

core: encode real numbers using uint32_t rather than uint16_t.

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/real.cpp

    r1013 r1016  
    2323{
    2424
     25static int const BIGIT_BITS = 32;
     26
    2527real::real(float f) { *this = (double)f; }
    2628real::real(int i) { *this = (double)i; }
    … …  
    4749    }
    4850
    49     m_mantissa[0] = u.x >> 36;
    50     m_mantissa[1] = u.x >> 20;
    51     m_mantissa[2] = u.x >> 4;
    52     m_mantissa[3] = u.x << 12;
    53     memset(m_mantissa + 4, 0, sizeof(m_mantissa) - 4 * sizeof(m_mantissa[0]));
     51    m_mantissa[0] = u.x >> 20;
     52    m_mantissa[1] = u.x << 12;
     53    memset(m_mantissa + 2, 0, (BIGITS - 2) * sizeof(m_mantissa[0]));
    5454}
    5555
    … …  
    8282
    8383        /* Store mantissa if necessary */
    84         u.x <<= 16;
     84        u.x <<= 32;
    8585        u.x |= m_mantissa[0];
    86         u.x <<= 16;
    87         u.x |= m_mantissa[1];
    88         u.x <<= 16;
    89         u.x |= m_mantissa[2];
    90         u.x <<= 4;
    91         u.x |= m_mantissa[3] >> 12;
     86        u.x <<= 20;
     87        u.x |= m_mantissa[1] >> 12;
    9288        /* Rounding */
    93         u.x += (m_mantissa[3] >> 11) & 1;
     89        u.x += (m_mantissa[1] >> 11) & 1;
    9490    }
    9591
    … …  
    132128    int e2 = x.m_signexp - (1 << 30) + 1;
    133129
    134     int bigoff = (e1 - e2) / (sizeof(uint16_t) * 8);
    135     int off = e1 - e2 - bigoff * (sizeof(uint16_t) * 8);
     130    int bigoff = (e1 - e2) / BIGIT_BITS;
     131    int off = e1 - e2 - bigoff * BIGIT_BITS;
    136132
    137133    if (bigoff > BIGITS)
    … …  
    140136    ret.m_signexp = m_signexp;
    141137
    142     uint32_t carry = 0;
     138    uint64_t carry = 0;
    143139    for (int i = BIGITS; i--; )
    144140    {
    … …  
    147143            carry += x.m_mantissa[i - bigoff] >> off;
    148144
    149         if (i - bigoff > 0)
    150             carry += (x.m_mantissa[i - bigoff - 1] << (16 - off)) & 0xffffu;
     145        if (off && i - bigoff > 0)
     146            carry += (x.m_mantissa[i - bigoff - 1] << (BIGIT_BITS - off)) & 0xffffffffu;
    151147        else if (i - bigoff == 0)
    152             carry += 0x0001u << (16 - off);
     148            carry += (uint64_t)1 << (BIGIT_BITS - off);
    153149
    154150        ret.m_mantissa[i] = carry;
    155         carry >>= 16;
     151        carry >>= BIGIT_BITS;
    156152    }
    157153
    … …  
    162158        for (int i = 0; i < BIGITS; i++)
    163159        {
    164             uint16_t tmp = ret.m_mantissa[i];
    165             ret.m_mantissa[i] = (carry << 15) | (tmp >> 1);
    166             carry = tmp & 0x0001u;
     160            uint32_t tmp = ret.m_mantissa[i];
     161            ret.m_mantissa[i] = (carry << (BIGIT_BITS - 1)) | (tmp >> 1);
     162            carry = tmp & 1u;
    167163        }
    168164        ret.m_signexp++;
    … …  
    194190    int e2 = x.m_signexp - (1 << 30) + 1;
    195191
    196     int bigoff = (e1 - e2) / (sizeof(uint16_t) * 8);
    197     int off = e1 - e2 - bigoff * (sizeof(uint16_t) * 8);
     192    int bigoff = (e1 - e2) / BIGIT_BITS;
     193    int off = e1 - e2 - bigoff * BIGIT_BITS;
    198194
    199195    if (bigoff > BIGITS)
    … …  
    202198    ret.m_signexp = m_signexp;
    203199
    204     int32_t carry = 0;
     200    int64_t carry = 0;
    205201    for (int i = 0; i < bigoff; i++)
    206202    {
    207203        carry -= x.m_mantissa[BIGITS - i];
    208         carry = (carry & 0xffff0000u) | (carry >> 16);
    209     }
    210     carry -= x.m_mantissa[BIGITS - 1 - bigoff] & ((1 << off) - 1);
    211     carry /= (1 << off);
     204        /* Emulates a signed shift */
     205        carry >>= BIGIT_BITS;
     206        carry |= carry << BIGIT_BITS;
     207    }
     208    carry -= x.m_mantissa[BIGITS - 1 - bigoff] & (((int64_t)1 << off) - 1);
     209    carry /= (int64_t)1 << off;
    212210
    213211    for (int i = BIGITS; i--; )
    … …  
    217215            carry -= x.m_mantissa[i - bigoff] >> off;
    218216
    219         if (i - bigoff > 0)
    220             carry -= (x.m_mantissa[i - bigoff - 1] << (16 - off)) & 0xffffu;
     217        if (off && i - bigoff > 0)
     218            carry -= (x.m_mantissa[i - bigoff - 1] << (BIGIT_BITS - off)) & 0xffffffffu;
    221219        else if (i - bigoff == 0)
    222             carry -= 0x0001u << (16 - off);
     220            carry -= (uint64_t)1 << (BIGIT_BITS - off);
    223221
    224222        ret.m_mantissa[i] = carry;
    225         carry = (carry & 0xffff0000u) | (carry >> 16);
     223        carry >>= BIGIT_BITS;
     224        carry |= carry << BIGIT_BITS;
    226225    }
    227226
    … …  
    238237            if (!ret.m_mantissa[i])
    239238            {
    240                 off += sizeof(uint16_t) * 8;
     239                off += BIGIT_BITS;
    241240                continue;
    242241            }
    243242
    244             for (uint16_t tmp = ret.m_mantissa[i]; tmp < 0x8000u; tmp <<= 1)
     243            for (uint32_t tmp = ret.m_mantissa[i]; tmp < 0x80000000u; tmp <<= 1)
    245244                off++;
    246245            break;
    247246        }
    248         if (off == BIGITS * sizeof(uint16_t) * 8)
     247        if (off == BIGITS * BIGIT_BITS)
    249248            ret.m_signexp &= 0x80000000u;
    250249        else
    … …  
    253252            ret.m_signexp -= off;
    254253
    255             bigoff = off / (sizeof(uint16_t) * 8);
    256             off -= bigoff * sizeof(uint16_t) * 8;
     254            bigoff = off / BIGIT_BITS;
     255            off -= bigoff * BIGIT_BITS;
    257256
    258257            for (int i = 0; i < BIGITS; i++)
    259258            {
    260                 uint16_t tmp = 0;
     259                uint32_t tmp = 0;
    261260                if (i + bigoff < BIGITS)
    262261                    tmp |= ret.m_mantissa[i + bigoff] << off;
    263                 if (i + bigoff + 1 < BIGITS)
    264                     tmp |= ret.m_mantissa[i + bigoff + 1] >> (16 - off);
     262                if (off && i + bigoff + 1 < BIGITS)
     263                    tmp |= ret.m_mantissa[i + bigoff + 1] >> (BIGIT_BITS - off);
    265264                ret.m_mantissa[i] = tmp;
    266265            }
    … …  
    288287    /* Accumulate low order product; no need to store it, we just
    289288     * want the carry value */
    290     uint64_t carry = 0;
     289    uint64_t carry = 0, hicarry = 0, prev;
    291290    for (int i = 0; i < BIGITS; i++)
    292291    {
    293292        for (int j = 0; j < i + 1; j++)
    294             carry += (uint32_t)m_mantissa[BIGITS - 1 - j]
    295                    * (uint32_t)x.m_mantissa[BIGITS - 1 + j - i];
    296         carry >>= 16;
     293        {
     294            prev = carry;
     295            carry += (uint64_t)m_mantissa[BIGITS - 1 - j]
     296                   * (uint64_t)x.m_mantissa[BIGITS - 1 + j - i];
     297            if (carry < prev)
     298                hicarry++;
     299        }
     300        carry >>= BIGIT_BITS;
     301        carry |= hicarry << BIGIT_BITS;
     302        hicarry >>= BIGIT_BITS;
    297303    }
    298304
    … …  
    300306    {
    301307        for (int j = i + 1; j < BIGITS; j++)
    302             carry += (uint32_t)m_mantissa[BIGITS - 1 - j]
    303                    * (uint32_t)x.m_mantissa[j - 1 - i];
    304 
     308        {
     309            prev = carry;
     310            carry += (uint64_t)m_mantissa[BIGITS - 1 - j]
     311                   * (uint64_t)x.m_mantissa[j - 1 - i];
     312            if (carry < prev)
     313                hicarry++;
     314        }
     315        prev = carry;
    305316        carry += m_mantissa[BIGITS - 1 - i];
    306317        carry += x.m_mantissa[BIGITS - 1 - i];
    307         ret.m_mantissa[BIGITS - 1 - i] = carry & 0xffffu;
    308         carry >>= 16;
     318        if (carry < prev)
     319            hicarry++;
     320        ret.m_mantissa[BIGITS - 1 - i] = carry & 0xffffffffu;
     321        carry >>= BIGIT_BITS;
     322        carry |= hicarry << BIGIT_BITS;
     323        hicarry >>= BIGIT_BITS;
    309324    }
    310325
    … …  
    315330        for (int i = 0; i < BIGITS; i++)
    316331        {
    317             uint16_t tmp = ret.m_mantissa[i];
    318             ret.m_mantissa[i] = (carry << 15) | (tmp >> 1);
    319             carry = tmp & 0x0001u;
     332            uint32_t tmp = ret.m_mantissa[i];
     333            ret.m_mantissa[i] = (carry << (BIGIT_BITS - 1)) | (tmp >> 1);
     334            carry = tmp & 1u;
    320335        }
    321336        e++;
    … …  
    473488    /* Use the system's float inversion to approximate 1/x */
    474489    union { float f; uint32_t x; } u = { 1.0f }, v = { 1.0f };
    475     v.x |= (uint32_t)x.m_mantissa[0] << 7;
    476     v.x |= (uint32_t)x.m_mantissa[1] >> 9;
     490    v.x |= x.m_mantissa[0] >> 9;
    477491    v.f = 1.0 / v.f;
    478492
    479493    real ret;
    480     ret.m_mantissa[0] = (v.x >> 7) & 0xffffu;
    481     ret.m_mantissa[1] = (v.x << 9) & 0xffffu;
     494    ret.m_mantissa[0] = v.x << 9;
    482495
    483496    uint32_t sign = x.m_signexp & 0x80000000u;
    … …  
    490503    /* FIXME: log2(BIGITS) steps of Newton-Raphson seems to be enough for
    491504     * convergence, but this hasn't been checked seriously. */
    492     for (int i = 1; i < real::BIGITS; i *= 2)
     505    for (int i = 2; i < real::BIGITS; i *= 2)
    493506        ret = ret * (real::R_2 - ret * x);
    494507
    … …  
    518531    union { float f; uint32_t x; } u = { 1.0f }, v = { 2.0f };
    519532    v.x -= ((x.m_signexp & 1) << 23);
    520     v.x |= (uint32_t)x.m_mantissa[0] << 7;
    521     v.x |= (uint32_t)x.m_mantissa[1] >> 9;
     533    v.x |= x.m_mantissa[0] >> 9;
    522534    v.f = 1.0 / sqrtf(v.f);
    523535
    524536    real ret;
    525     ret.m_mantissa[0] = (v.x >> 7) & 0xffffu;
    526     ret.m_mantissa[1] = (v.x << 9) & 0xffffu;
     537    ret.m_mantissa[0] = v.x << 9;
    527538
    528539    uint32_t sign = x.m_signexp & 0x80000000u;
    … …  
    536547    /* FIXME: log2(BIGITS) steps of Newton-Raphson seems to be enough for
    537548     * convergence, but this hasn't been checked seriously. */
    538     for (int i = 1; i < real::BIGITS; i *= 2)
     549    for (int i = 2; i < real::BIGITS; i *= 2)
    539550    {
    540551        ret = ret * (real::R_3 - ret * ret * x);
    … …  
    593604    {
    594605        tmp.m_signexp = 0xffffffffu;
    595         tmp.m_mantissa[0] = 0xffffu;
     606        tmp.m_mantissa[0] = 0xffffffffu;
    596607        return tmp;
    597608    }
    … …  
    658669        if (exponent <= 0)
    659670            ret.m_mantissa[i] = 0;
    660         else if (exponent < 8 * (int)sizeof(uint16_t))
    661             ret.m_mantissa[i] &= ~((1 << (16 - exponent)) - 1);
    662 
    663         exponent -= 8 * sizeof(uint16_t);
     671        else if (exponent < BIGIT_BITS)
     672            ret.m_mantissa[i] &= ~((1 << (BIGIT_BITS - exponent)) - 1);
     673
     674        exponent -= BIGIT_BITS;
    664675    }
    665676

  • trunk/src/real.h

    r1003 r1016  
    109109     * Newton-Raphson iterations) and real::print (the number of printed
    110110     * digits) */
    111     static int const BIGITS = 32;
     111    static int const BIGITS = 16;
    112112
    113113    uint32_t m_size;
    114114    uint32_t m_signexp;
    115     uint16_t m_mantissa[BIGITS];
     115    uint32_t m_mantissa[BIGITS];
    116116};
    117117
Note: See TracChangeset for help on using the changeset viewer.