Re: [boost] fixed-pt [improved]

From: Anthony Williams (anthony.ajw_at_[hidden])
Date: 2008-08-01 10:31:37

Neal Becker <ndbecker2_at_[hidden]> writes:

> Here is an improved version of fixed-pt, incorporating some of the
> suggestions given here:

>   self& operator*=(self const& x) { 
>     base_type tmp = (base_type)val * (base_type)x.val;
>     val = tmp >> (frac_bits);
>     return *this;
>   }

This is subject to overflow, even where the result isn't.

Suppose you multiply the max value by a fixed point representation of
1. val=INT_MAX, x.val=(1<<frac_bits), and the calculation of tmp
overflows.

>   self& operator/=(self const& x) { 
>     val = ((base_type)val << frac_bits) / (base_type)x.val;
>     return *this; }

This is also subject to overflow with the (val<<frac_bits).

Here's the code for multiplication and division from my DDJ article on
fixed-point arithmetic:

http://www.justsoftwaresolutions.co.uk/news/optimizing-applications-with-fixed-point-arithmetic.html

The internal rep is 64-bit. fixed_resolution_shift is your
frac_bits. For my code it is 28, but this doesn't actually matter for
multiplication and division. All of it could easily be extended to
allow for different numbers of fractional bits, apart from those parts
that need the lookup tables (sin, cos, atan, exp, log).

fixed& fixed::operator*=(fixed const& val)
{
    bool const val_negative=val.m_nVal<0;
    bool const this_negative=m_nVal<0;
    bool const negate=val_negative ^ this_negative;
    unsigned __int64 const other=val_negative?-val.m_nVal:val.m_nVal;
    unsigned __int64 const self=this_negative?-m_nVal:m_nVal;
    
    if(unsigned __int64 const self_upper=(self>>32))
    {
        m_nVal=(self_upper*other)<<(32-fixed_resolution_shift);
    }
    else
    {
        m_nVal=0;
    }
    if(unsigned __int64 const self_lower=(self&0xffffffff))
    {
        unsigned long const other_upper=static_cast<unsigned long>(other>>32);
        unsigned long const other_lower=static_cast<unsigned long>(other&0xffffffff);
        unsigned __int64 const lower_self_upper_other_res=self_lower*other_upper;
        unsigned __int64 const lower_self_lower_other_res=self_lower*other_lower;
        m_nVal+=(lower_self_upper_other_res<<(32-fixed_resolution_shift))
            + (lower_self_lower_other_res>>fixed_resolution_shift);
    }
    
    if(negate)
    {
        m_nVal=-m_nVal;
    }
    return *this;
}

fixed& fixed::operator/=(fixed const& divisor)
{
    if( !divisor.m_nVal)
    {
        m_nVal=fixed_max.m_nVal;
    }
    else
    {
        bool const negate_this=(m_nVal<0);
        bool const negate_divisor=(divisor.m_nVal<0);
        bool const negate=negate_this ^ negate_divisor;
        unsigned __int64 a=negate_this?-m_nVal:m_nVal;
        unsigned __int64 b=negate_divisor?-divisor.m_nVal:divisor.m_nVal;

        unsigned __int64 res=0;
    
        unsigned __int64 temp=b;
        bool const a_large=a>b;
        unsigned shift=fixed_resolution_shift;

        if(a_large)
        {
            unsigned __int64 const half_a=a>>1;
            while(temp<half_a)
            {
                temp<<=1;
                ++shift;
            }
        }
        unsigned __int64 d=1I64<<shift;
        if(a_large)
        {
            a-=temp;
            res+=d;
        }

        while(a && temp && shift)
        {
            unsigned right_shift=0;
            while(right_shift<shift && (temp>a))
            {
                temp>>=1;
                ++right_shift;
            }
            d>>=right_shift;
            shift-=right_shift;
            a-=temp;
            res+=d;
        }
        m_nVal=(negate?-(__int64)res:res);
    }
    
    return *this;
}

Anthony
-- 
Anthony Williams            | Just Software Solutions Ltd
Custom Software Development | http://www.justsoftwaresolutions.co.uk
Registered in England, Company Number 5478976.
Registered Office: 15 Carrallack Mews, St Just, Cornwall, TR19 7UL