2

I would like to get to the code that makes floats and doubles work. The reason why I'm trying is that I need to introduce a new floating point type defined by the military that is not supported natively by Microsoft.

Has this been done before or is this even possible? If so, is there files that I can retrieve or steps to follow to do it myself?

2
  • 4
    A software emulation of an FPU should be a good starting point.
    – Yegor
    Commented Jun 11, 2015 at 5:07
  • 4
    Open Watcom is a compiler that dates back (it wasn't open then) to when Intel CPUs didn't (neccessarily) have FPUs, so it comes with FPU emulation source. However, unless you are comfortable with the mathematical concepts (approximating sqrt with Newton; approximating trig functions with taylor expansions), much of what's happening will seem like magic to you. Commented Jun 12, 2015 at 11:14

2 Answers 2

3
  1. For starters you can encode/decode the military format to IEEE 754

    this way you can use FPU and standard math stuff so at start everything would work. The code would look like this:

        military fabs(military x) 
         {
         double _x;
         _x=military2double(x);
         _x=fabs(_x);
         return double2military(_x);
         }
    
  2. convert functions to your format

    so start rewriting math function one by one to handle military format natively

        military fabs(military x) 
         {
         if (x<0) x=-x;
         return x;
         }
    

    the military should be a class from this point and you need to code the basic operators and functions into it so it handles as common variable data type. Do not forget to optimize what you can for example fabs can be done by single bitwise AND mask out the mantissa sign bit unless your military format does not use two's complement.

Here is how one of mine floating point number classes look like:

//---------------------------------------------------------------------------
const int _arbnum_max_siz=1<<18;            // mantisa limit [32bit words] 1<<18 = MByte
const int _arbnum_min_exp=-0x80000000;      // min exponent limit for base = 2
const int _arbnum_max_exp=+0x7FFFFFFD;      // max exponent limit for base = 2
const int _arbnum_nan_exp=+0x7FFFFFFE;      // not a number exponent constatn
const int _arbnum_inf_exp=+0x7FFFFFFF;      // infinity exponent constatn
//---------------------------------------------------------------------------
enum _arbnum_error_enum
    {
    _arbnum_error_AllOk=0,
    _arbnum_error_NotEnoughMemory,          // memory allocation
    _arbnum_error_CorruptedSignum,          // bad combination of sig,exp,siz,dat
    _arbnum_error_ZeroSize,                 // unallocated number usage
    _arbnum_error_UnexpectedCharInFormat,   // str()
    _arbnum_error_TooBigNumber,             // str() dec num is bigger than expected
    };
//---------------------------------------------------------------------------

// |-----|---------------------------|---------------|------|
// | sig | MSB      mantisa      LSB |   exponent    | bits |
// |-----|---------------------------|---------------|------|
// | +1  | 0.(0      ...          0) | 2^0           |   0  | +zero
// | -1  | 0.(0      ...          0) | 2^0           |   0  | -zero
// |-----|---------------------------|---------------|------|
// | +1  | 1.(dat[0] ... dat[siz-1]) | 2^exp         |   n  | +number
// | -1  | 1.(dat[0] ... dat[siz-1]) | 2^exp         |   n  | -number
// |-----|---------------------------|---------------|------|
// | +1  | 1.0                       | 2^+0x7FFFFFFE |   1  | +infinity
// | -1  | 1.0                       | 2^+0x7FFFFFFE |   1  | -infinity
// |-----|---------------------------|---------------|------|
//         ^                               ^             ^
//         ^                               ^          used bits count in mantisa
//         ^                          exponent of msb of dat[0]
//         1. is included in mantisa dat[0] !!!
//---------------------------------------------------------------------------
class arbnum
    {
public:
    ALU32 alu;
    // dat is LSDW first ... MSDW last
    DWORD *dat; int siz,exp,sig,bits;
    arbnum()                { dat=NULL; siz=0; exp=1; sig=+1; bits=0; }
    arbnum(int a)           { dat=NULL; siz=0; exp=1; sig=+1; bits=0; *this=a; }
    arbnum(DWORD a)         { dat=NULL; siz=0; exp=1; sig=+1; bits=0; *this=a; }
    arbnum(float a)         { dat=NULL; siz=0; exp=1; sig=+1; bits=0; *this=a; }
    arbnum(double a)        { dat=NULL; siz=0; exp=1; sig=+1; bits=0; *this=a; }
    arbnum(AnsiString a)    { dat=NULL; siz=0; exp=1; sig=+1; bits=0; *this=a; }
    arbnum(arbnum& a)       { dat=NULL; siz=0; exp=1; sig=+1; bits=0; *this=a; }
    ~arbnum()               { _free(); }
    arbnum* operator = (const arbnum *a) { *this=*a; return this; }
    arbnum* operator = (const arbnum &a) { _alloc(a.siz); for(int i=0;(i<siz)&&(i<a.siz);i++) dat[i]=a.dat[i]; exp=a.exp; sig=a.sig; bits=a.bits; return this; }
    // error handling
    void _error(int err)
        {
        err=1/0;
        }
    // memory management
    void _free();
    void _alloc(int _siz);
    void _realloc(int _siz);
    void _normalize(int _fbits=-1);
    // export functions
    AnsiString str(const char *cs="");          // cs: "5.3" ".3" "5" "s5.3" " 5.3" "M.M" b,d,h e
    int        geti32();
    DWORD      getu32();
    float      getf32();
    double     getf64();
    // assign operators
    arbnum operator = (int x);
    arbnum operator = (DWORD x);
    arbnum operator = (float x);
    arbnum operator = (double x);
    arbnum operator = (AnsiString x);           // nahra cislo zo stringu (+/-)dec/bin/hex(d/b/h) podla posledneho znaku
    // arithmetic operators
    arbnum operator + ();
    arbnum operator - ();
    arbnum operator ++ ();                      // ++x  return this+1
    arbnum operator -- ();                      // --x  return this-1
    arbnum operator ++ (int y);                 // x++  return this, and then this+1
    arbnum operator -- (int y);                 // x--  return this, and then this-1
    arbnum operator << (int n);
    arbnum operator >> (int n);
    arbnum operator <<=(int n)          { this[0]=this[0]<<n; return *this; };
    arbnum operator >>=(int n)          { this[0]=this[0]>>n; return *this; };
    arbnum operator & (const arbnum &x);
    arbnum operator | (const arbnum &x);
    arbnum operator ^ (const arbnum &x);
    arbnum operator + (const arbnum &x);
    arbnum operator - (const arbnum &x);
    arbnum operator * (const arbnum &x);
    arbnum operator / (const arbnum &x);
    arbnum operator % (const arbnum &x);
    arbnum operator &=(const arbnum &x) { this[0]=this[0]&x; return *this; };
    arbnum operator |=(const arbnum &x) { this[0]=this[0]|x; return *this; };
    arbnum operator ^=(const arbnum &x) { this[0]=this[0]^x; return *this; };
    arbnum operator +=(const arbnum &x) { this[0]=this[0]+x; return *this; };
    arbnum operator -=(const arbnum &x) { this[0]=this[0]-x; return *this; };
    arbnum operator *=(const arbnum &x) { this[0]=this[0]*x; return *this; };
    arbnum operator /=(const arbnum &x) { this[0]=this[0]/x; return *this; };
    arbnum operator %=(const arbnum &x) { this[0]=this[0]%x; return *this; };
    // bool operators
    int operator !  (               ) { return   iszero();                       }
    int operator && (const arbnum &x) { return ((isnonzero())&&(x.isnonzero())); }
    int operator || (const arbnum &x) { return ((isnonzero())||(x.isnonzero())); }
    int operator == (const arbnum &x);
    int operator != (const arbnum &x);
    int operator <  (const arbnum &x);
    int operator >  (const arbnum &x);
    int operator <= (const arbnum &x);
    int operator >= (const arbnum &x);
    int  isnan();                               // ( this == nan )
    int  isinf();                               // ( this == inf )
    int  isone();                               // ( this == +1 )
    int  iszero();                              // ( this == 0 )
    int  isnonzero();                           // ( this != 0 )
    int  isinteger();                           // ( this==integer(this))   ...-2,-1,0,+1,+2,...
    int  isreal();                              // ( this!=integer(this))   12.3456
    int  isfraction();                          // ( 1   >     abs(this))   0.555 -0.789
    // bit functions
    void bitset(int e);                         // this|=2^e
    void bitres(int e);                         // this&=2^e        [ normalize! ]
    void bitnot(int e);                         // this^=2^e        [ normalize! ]
    int  bitget(int e);                         //      =this & 2^e
    int  bits2size(int n);                      // = num of DWORDs needed for n>=0 bits
    int  nibits();                              // = num of used integer bits (from number)
    int  nfbits();                              // = num of used fractional bits (from number)
    int  mibits();                              // = num of used integer bits (from mantisa)
    int  mfbits();                              // = num of used fractional bits (from mantisa)
    // low level arithmetics
    void nan();                                 // this =+nan
    void inf();                                 // this =+inf
    void one();                                 // this =+1
    void zero();                                // this =+0
    void integer();                             // this = odreze desatinnu cast
    void fraction();                            // this = odreze celu cast
    void round();                               // this = zaokruhli podla 0.5
    void floor();                               // this = zaokruhli nadol
    void ceil ();                               // this = zaokruhli nahor
    void overflow();                            // this = zaokruhli ak zistil pretecenie (????????????.FFFFFFFFFFFFFF??h     [****]
    void shl(int n);                            // mantisa <<= n    [ normalize! ]
    void shr(int n);                            // mantisa >>= n    [ normalize! ]
    void inci();                                // this = sig (|this|+1)
    void deci();                                // this = sig (|this|-1)
    void incf();                                // this = sig (|this|+lsb)
    void decf();                                // this = sig (|this|-lsb)
    void add(const arbnum &x,const arbnum &y);  // this = |x|+|y|                    kazi sig !!! [ normalize! ]
    void sub(const arbnum &x,const arbnum &y);  // this = |x|-|y|  ; |x|>=|y|        kazi sig !!! [ normalize! ]
    void mul(const arbnum &x,const arbnum &y);  // this = |x|*|y|                    kazi sig !!! [ normalize! ]
    void _mul_karatsuba(DWORD *dst,DWORD *x,DWORD *y,int n);
    void mul_karatsuba(const arbnum &x,const arbnum &y);
    void mul_NTT(const arbnum &x,const arbnum &y);
    void sqr(const arbnum &x);                  // this = |x|*|x|                    kazi sig !!! [ normalize! ]
    void sqr_NTT(const arbnum &x);
    void div(const arbnum &x,const arbnum &y,int acc=-1);   // this = |x|/|y|  ; y = <0.5,1)     kazi sig !!! [ normalize! ]
    void div_binary(const arbnum &x,const arbnum &y);
    int  geq(const arbnum &x,const arbnum &y);  //      = |x|>=|y| ; 0 <, 1 >, 2==
    int int_log2(int x)
        {
        if (x<=0) return 0;
        int y=0,a=1;
        for (;a<=x;a<<=1,y++); y--;
        return y;
        }
    // friend arithmetics
    friend arbnum integer (const arbnum &x);                    // = integer(x)
    friend arbnum fraction(const arbnum &x);                    // = fraction(x)
    friend arbnum round   (const arbnum &x);                    // = round(x);
    friend arbnum floor   (const arbnum &x);                    // = floor(x)
    friend arbnum ceil    (const arbnum &x);                    // = ceil(x)
    friend arbnum abs     (const arbnum &x);                    // = |x|
    friend arbnum sqrt    (const arbnum &x);                    // = |x|^0.5
    friend arbnum exp2    (const arbnum &x);                    // = 2^x
    friend arbnum expe    (const arbnum &x);                    // = e^x
    friend arbnum exp10   (const arbnum &x);                    // = 10^x
    friend arbnum log2    (const arbnum &x);                    // = log2(x)
    friend arbnum loge    (const arbnum &x);                    // = ln(x)
    friend arbnum log10   (const arbnum &x);                    // = log10(x)
    friend arbnum powi    (const arbnum &x,const arbnum &y);    // = x^int(y)
    friend arbnum powi    (const arbnum &x,int y);              // = x^int(y)
    friend arbnum pow     (const arbnum &x,const arbnum &y);    // = x^y
    friend arbnum sin     (const arbnum &x);                    // = sin(x)
    friend arbnum cos     (const arbnum &x);                    // = cos(x)
    friend arbnum tan     (const arbnum &x);                    // = tan(x)
    friend arbnum cotan   (const arbnum &x);                    // = cotan(x)
    friend arbnum asin    (const arbnum &x);                    // = asin(x)
    friend arbnum acos    (const arbnum &x);                    // = acos(x)
    friend arbnum atan    (const arbnum &x);                    // = atan(x)
    friend arbnum acotan  (const arbnum &x);                    // = acotan(x)
    friend arbnum gcd     (arbnum x,arbnum y);                  // = gcd (int(|x|),int(|y|))  binary greatest common divisor
    friend arbnum fact    (const DWORD &x);                     // = x! = (x/2)!^2 * mul[ p(i)^e(i) ]
    friend arbnum fact_aprox(const DWORD &x);                   // ~ x!
    };
//---------------------------------------------------------------------------

[notes]

Without more info about your target platform or military format it is impossible to decide which approach is better in respect to speed,accuracy. We need info like mantissa bits/format exponent bits/format, do you have FPU at disposal (and for which format), what float types you have natively at disposal etc.

1
  • @CodeMonkeyV2 finished editing. sorry some comments are not in English but the code should be straightforward enough anyway ...
    – Spektre
    Commented Jun 15, 2015 at 7:51
1

The IEEE 754 standard describes how the floating point representation works under the hood.

2
  • Sure, as a representation. What I need is the code that makes it work, what is happening behind the scenes. Like how Microsoft's std::string can be traced to view the code that makes it function. Commented Jun 10, 2015 at 21:48
  • 2
    There is no such code in the modern OS, as it is implemented directly by the CPU.
    – Spl3en
    Commented Jun 10, 2015 at 22:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.