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);
}
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.