[This message was posted by Walter Mascarenhas of GeoCAD <[EMAIL PROTECTED]> to
the "FAST Protocol" discussion forum at http://fixprotocol.org/discuss/46. You
can reply to it on-line at http://fixprotocol.org/discuss/read/758d78cb -
PLEASE DO NOT REPLY BY MAIL.]
Dimitri,
I agree that Daniel´s approach is simpler
and is certainly faster. It is clearly the best
option if you do not worry about high precision:
who cares about the 5th digit in a million
dollar amount?
(it is not bad in this aspect too, it is
just no optimal)
Last night I played a little bit with this
stuff and come up with the code pasted below.
It is specific for the 10^-4 precision and
certainly not perfect. I believe that, on
the average, if gives full precision using
a few shifts and multiplies of 64bit integers
(which is quite adequate for the 64bit platforms)
The "average" aspect comes from the normalization
step, which is not considered by Daniel and can
be quite time consuming. In the algorithm below
I belive that in 90% of the time there will
be no need to normalize the result, in 9% it
will normalized with one division, in 0.09%
normalization will take 2 divisions and so on.
The same approach works for precision 0,1,2 and 3
but higher precision would require 96bits arithmetic
or even higher for very high precisions.
Anyway, here is the code:
Walter
//
// Given a double x >= 0 and and integer precision
// DoubleToDecimal4 returns a mantissa m and
// assigns an exponent to e exp so that
//
// a) |x - m * 10^(-e)| <= 0.5 * 10^(-4)
// b) either x = m = 0 or (m % 10) != 0
//
// The function may trhow an exception if
// |x| >= 10^250 because for some x in this
// range case no such e and m exist.
//
#define FORMAT_SENSITIVE 1
__int64 DoubleToDecimal4(double x, int& exp)
{
int e2;
int sign;
__int64 mLong;
if( x < 0 )
{
x = -x;
sign = -1;
}
else
{
sign = 1;
}
double md = frexp(x,&e2);
mLong = (__int64) _scalb(md,53);
// At this point we have x = 2^(e2 - 53) * mLong (exactly, no rounding)
int e;
if( e2 < 49) // mLong * (2^(e2 - 49) * 5^4 is not an integer
{
int shift = 49 - e2;
if( shift >= 64 ) // for instance in the case of denormals
{
exp = 0;
return 0;
}
mLong *= 625; // 625 = 5^4 < 2^10 and since mLong <= 2^53 there is no
overflow here
__int64 mask = (LongOne << shift) - 1;
__int64 remainder = mLong & mask;
if( (remainder << 1) >= mask )
{
mLong = (mLong >> shift) + 1;
}
else
{
mLong >>= shift;
}
if( mLong == 0 )
{
exp = 0;
return 0;
}
e = 4;
}
else
{
// here e2 >= 49 e mLong > 0
// Since x = 2^(e2 - 49) * 5^4 * mLong * 10^-4 there will
// be no rounding in this case.
switch(e2)
{
case 49:
{
if( mLong & 1 )
{
exp = 4;
return 625 * mLong;
}
mLong >>= 1;
// there is no break here
}
case 50:
{
if( mLong & 1 )
{
exp = 3;
return 125 * mLong;
}
mLong >>= 1;
// there is no break here
}
case 51:
{
if( mLong & 1 )
{
exp = 2;
return 25 * mLong;
}
mLong >>= 1;
// there is no break here
}
case 52:
{
if( mLong & 1 )
{
exp = 1;
return 5 * mLong;
}
mLong >>= 1;
}
case 53:
{
e = 0;
break;
}
default:
{
e = 0; // in the next while we try to divide mLong by 5 instead
// of multiply it by 2 in order to avoid overflow in mLong
do
{
__int64 remainder = (mLong % 5);
if( remainder )
{
int shift = e2 - 53;
if( mLong >= (LongOne << (63 - shift)) )
{
throw std::exception("Overflow in function DoubleToDecimal4");
}
exp = e;
return sign * (mLong << shift);
}
--e;
mLong /= 5;
}
while( --e2 > 53 );
}
}
}
// normalizing mLong
while( (mLong % 10) == 0 )
{
--e;
mLong /= 10;
}
exp = e;
return sign * mLong;
}
> Thanks, Walter,
>
> you are right - using byte vector for mantissa essentially erases most
> of FAST compression benefits. As for the algorithm, I think Daniel's
> approach is simpler and is probably more efficient, especially in the
> first case (there is only one multiplication operation). In the current
> implementation, we perform the conversion as double->string->FAST
> scaled number, so I am curious how much faster is double->FAST decimal
> approach is.
>
> Dimitry
>
> > Daniel and Dimitri,
> >
> > These floating point issues can get very tricky. I,
> > for instance, have been mislead by them many, many times.
> >
> > One point I am trying to make is that you don´t have the
> > 12.34500 value in the double world to start with. There is no such a
> > thing as a nice double with value x = 12.34500. The closest you
> > will ever get to x is xr equal to
> > 123450000000000006394884621840901672840118408203125 x 10^(-49) with
> > an error of about 6.4 x 10^(-16).
> >
> > As a result, in most cases the double produced by rounding will have
> > its least signficant decimal different from zero (for instance
> > 12.34500 would look like 13.3450000000001) and you would have no
> > compression: the full 64 bits of the number would need to be sent down
> > the wire, because you would never know if that last one was there on
> > purpose or it was caused by rounding. Or, in the bit vector
> > alternative suggest a few messages earlier, the bit vector would have
> > its maximum length most of the time.
> >
> > as for Daniel´s suggestion,
> >
> > >> float x = 123.4500 int64_t mantissa = (int64_t)(x * 10000); int32_t
> > >> exponent = 4;
> >
> > I would do something different, but his suggestion is much simpler and
> > he is right about the details: it can get very messy in general cases.
> >
> > I would go roughly like this (assuming x > 0, that we have a
> > unsigned int with 128 bits and disregarding some details, so that
> > you get the idea):
> >
> > int e; double xm = frexp(x,&e); UINT64 m = (UINT64) _scalb(xm,53); //
> > now x = 2^(e - 53) m exactly
> >
> > UINT128 result = Mult(m,5^4); // now x = result * 2^(e - 49) * 10^(-
> > 4), exactly
> >
> > if( e>= 49 ) { result <<= (e - 49); // now x = result * 10^(-4)
> > exactly } else { UINT64 one = 1; UINT64 mask = ((one << (49 - e)) -
> > 1); UINT128 reminder = result & mask; if( 2* reminder > mask ) result
> > = (result >> (49 - e)) + 1; else result >>= 49 - e; } Hope this helps
> > (but don´t trust the details...)
> >
> > Walter.
> >
> >
> >
> >
> > > Dimity, There are two ways to look at this problem. If you know
> > > ahead of time that you are going to always use a fixed number of
> > > decimal precision (again, let's say your business rules require a
> > > precision of 4 decimal points), and you have already either rounded
> > > your float to those 4 points of precision, or truncation is okay for
> > > your application, then the conversion is simple:
> > >
> > > float x = 123.4500 int64_t mantissa = (int64_t)(x * 10000); int32_t
> > > exponent = 4;
> > >
> > >
> > > The more complicated case is when you are trying to generically
> > > convert a floating point number to a FAST scaled decimal with
> > > a.) not losing any precision, and
> > > b.) optimizing the exponent as to keep the mantissa as small as
> > > possible.
> > >
> > > For this, take a look at modf() in the standard C library. It breaks
> > > a float into the whole and fractional parts. You can then cast or
> > > convert the float whole and fractional parts to integers. Next, you
> > > would remove any unnecessary precision from the fractional integer
> > > by using a mod and divide by 10 while there are trailing zeros left.
> > >
> > > Finally, to create the FAST scaled decimal mantissa, you must
> > > determine the FAST exponent by inspecting the size of the fractional
> > > integer, and then adjust the whole integer by that factor, and add
> > > back the fractional part.
> > >
> > > I am working on C/C++ an example for you...
> > >
> > > /Daniel
[You can unsubscribe from this discussion group by sending a message to
mailto:[EMAIL PROTECTED]
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"Financial Information eXchange" group.
To post to this group, send email to FIX-Protocol@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at
http://groups.google.com/group/FIX-Protocol?hl=en
-~----------~----~----~----~------~----~------~--~---
