Listing 2: Implementation of class fp

// fp.cpp
#include <math.h>
#include <stdlib.h>
#include <iostream>
#include <iomanip>
#include <sstream>
#include "fp.h"

const int BIAS = 127;
const int MIN_EXP = 0;
const int MAX_EXP = 255;
const int P = 15;
const long MIN_NORM = 1 << (P - 1);
const long ONE = 1 << P;

inline long fp::lf() const
    {    // sign-extend fraction and convert to long
    return long(is_neg ? -frac : frac);
    }

fp::fp()
    {    // construct
    convert(0.0);
    }

fp::fp(const char *str)
    {    // construct from string
    convert(strtod(str, NULL));
    }

fp::fp(double d)
    {    // construct from double
    convert(d);
    }

bool fp::eq(const fp& f) const
    {    // return this equals f
    return is_neg == f.is_neg && exp == f.exp && frac == f.frac;
    }

bool fp::lt(const fp& f) const
    {    // return this less than f
    return is_neg
        ? exp == f.exp
            ? f.frac < frac
            : f.exp < exp
        : exp == f.exp
            ? frac < f.frac
            : exp < f.exp;
    }

fp operator + (const fp& f1, const fp& f2)
    {    // add f1 to f2
    fp res = f1;
    res += f2;
    return res;
    }

fp& fp::operator += (const fp& f)
    {    // add f to this fp object
    return f.exp <= exp
        ? add(*this, f)
        : add(f, *this);    
    }

fp operator - (const fp& f1, const fp& f2)
    {    // subtract f2 from f1
    fp res = f1;
    return res -= f2;
    }

fp& fp::operator -= (const fp& f)
    {    // subtract f from this fp object
    return f.exp <= exp
        ? add(*this, -f)
        : add(-f, *this);    
    }

fp operator - (const fp& f)
    {    // negate this fp object
    fp res = f;
    res.is_neg = !res.is_neg;
    return res;
    }

fp operator * (const fp& f1, const fp& f2)
    {    // multiply f1 by f2
    fp res = f1;
    res *= f2;
    return res;
    }

fp& fp::operator *= (const fp& f)
    {    // multiply this fp object by f
    if (exp == 0 && frac == 0)
        is_neg = is_neg != f.is_neg;
    else if (f.exp == 0 && f.frac == 0)
        {    // handle zero second argument
        exp = 0;
        frac = 0;
        is_neg = is_neg != f.is_neg;
        }
    else
        normalize(exp + f.exp - BIAS, lf() * f.lf());
    return *this;    
    }

fp operator / (const fp& f1, const fp& f2)
    {    // divide f1 by f2
    fp res = f1;
    res /= f2;
    return res;
    }

fp& fp::operator /= (const fp& f)
    {    // divide this fp object by f
    if (exp == 0 && frac == 0)
        is_neg = is_neg != f.is_neg;
    else if (f.exp == 0 && f.frac == 0)
        exception(div_by_zero);
    else
        normalize(exp - f.exp + BIAS,
            ((lf() << P) / f.lf()) << P);
    return *this;    
    }

fp& fp::add(const fp& f1, const fp& f2)
    {    // set this object to sum of f1 and f2
    long fracr = f1.lf() << P;
    if (f1.exp - f2.exp < P + 2)
        fracr += ((f2.lf() << P) >> (f1.exp - f2.exp));
    return normalize(f1.exp, fracr);
    }

void fp::convert(double d)
    {    // convert double to fp object
    int e;
    double f = frexp(d, &e);
    if (f == 0.0)
        {    // handle zero specially
        exp = 0;
        frac = 0;
        is_neg = false;
        }
    else
        normalize(e + BIAS, long(f * (1L << P)) << P);
    }

fp& fp::normalize(int e, long fr)
    {    // set this object to normalized value of
        // exponent e and fraction fr
    if (fr == 0)
        {    // handle zero specially
        exp = 0;
        frac = 0;
        return *this;
        }
    if ((ONE << P) <= abs(fr))
        {    // handle fraction overflow
        fr >>= 1;
        ++e;
        }
    while (abs(fr) < (MIN_NORM << P))
        {    // scale left
        fr <<= 1;
        --e;
        }
    // check for round
    if (abs(fr) & (1 << (P - 1)))
        if (fr < 0)
            fr -= 1 << P;
        else    
            fr += 1 << P;
    if ((ONE << P) <= abs(fr))
        {    // handle fraction overflow
        fr >>= 1;
        ++e;
        }
    if (e < MIN_EXP)
        exception(exp_underflow);
    else if (MAX_EXP < e)
        exception(exp_overflow);
    exp = e;
    is_neg = fr < 0;
    frac = abs(fr) >> P;
    return *this;
    }

std::ostream& operator << (std::ostream& out, const fp& f)
    {    // insert into stream
    double d;
    if (f.exp == 0 && f.frac == 0)
        d = 0.0;
    else
        {    // convert non-zero value to double
        d = f.frac;
        d /= (1 << P);
        int e = f.exp - BIAS;
        if (e == 0)
            ;
        else if (e < 0)
            for (int i = 0; i < -e; ++i)
                d /= 2;
        else    
            for (int i = 0; i < e; ++i)
                d *= 2;
        if (f.is_neg)
            d = -d;
        }
    return out << d << ' ' << to_string(f);
    }

std::string to_string(const fp& f)
    {    // convert to text representation of internal data
    std::ostringstream str;
    str << '(' << (f.is_neg ? "-," : "+,")
        << (int)(f.exp - BIAS) << ','
        << std::hex << f.frac << std::dec << ')';
    return str.str();
    }

static char *except_desc[] =
    {    // descriptions of floating point exceptions
    "divide by zero",
    "exponent underflow",
    "exponent overflow"
    };

void fp::exception(fp_exception except)
    {    // handle floating point exceptions
    std::cerr << "floating point error: "
        << except_desc[except] << std::endl;
    std::exit(1);
    }

const fp fp::zero(0.0);
const fp fp::one(1.0);
const fp fp::two(2.0);