Ignore:
Timestamp:
Sep 27, 2011, 7:09:28 PM (11 years ago)
Author:
sam
Message:

core: add exp() for real numbers and fix the == operator.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/real.cpp

    r987 r988  
    347347bool real::operator ==(real const &x) const
    348348{
     349    if ((m_signexp << 1) == 0 && (x.m_signexp << 1) == 0)
     350        return true;
     351
    349352    if (m_signexp != x.m_signexp)
    350353        return false;
     
    499502}
    500503
     504real exp(real const &x)
     505{
     506    int square = 0;
     507
     508    /* FIXME: this is slow. Find a better way to approximate exp(x) for
     509     * large values. */
     510    real tmp = x, one = 1.0;
     511    while (tmp > one)
     512    {
     513        tmp.m_signexp--;
     514        square++;
     515    }
     516
     517    real ret = 1.0, fact = 1.0, xn = tmp;
     518
     519    for (int i = 1; i < 100; i++)
     520    {
     521        fact *= (real)i;
     522        ret += xn / fact;
     523        xn *= tmp;
     524    }
     525
     526    for (int i = 0; i < square; i++)
     527        ret = ret * ret;
     528
     529    return ret;
     530}
     531
    501532void real::print(int ndigits) const
    502533{
     
    512543    /* Normalise x so that mantissa is in [1..9.999] */
    513544    int exponent = 0;
    514     for (real div = r1, newdiv; true; div = newdiv)
    515     {
    516         newdiv = div * r10;
    517         if (x < newdiv)
     545    if (x.m_signexp)
     546    {
     547        for (real div = r1, newdiv; true; div = newdiv)
    518548        {
    519             x /= div;
    520             break;
     549            newdiv = div * r10;
     550            if (x < newdiv)
     551            {
     552                x /= div;
     553                break;
     554            }
     555            exponent++;
    521556        }
    522         exponent++;
    523     }
    524     for (real mul = 1, newx; true; mul *= r10)
    525     {
    526         newx = x * mul;
    527         if (newx >= r1)
     557        for (real mul = 1, newx; true; mul *= r10)
    528558        {
    529             x = newx;
    530             break;
     559            newx = x * mul;
     560            if (newx >= r1)
     561            {
     562                x = newx;
     563                break;
     564            }
     565            exponent--;
    531566        }
    532         exponent--;
    533567    }
    534568
Note: See TracChangeset for help on using the changeset viewer.