Neal Becker <ndbecker2@gmail.com> 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-fix... 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

Neal, Please excuse me if I'm widening the scope of the discussion too far. I use an internally developed fixed_point_10 class for financial work to combat base2 vs base10 representation issues (e.g. it is impossible to represent 0.10 exactly as a float, double or base2 fixed_pt). The implementation is fundamentally very similar to youf fixed_pt class, but instead of specifying frac_bits you specify frac_digits. Instead of using left/right shifts for scaling it multiplies/divides by a compile time calculated denominator based on frac_digits. It seems like it might be possible to extend your implementation to support both base2/base10, either by specifying an integer base to go along with frac_bits (mabye better named frac_places given the generic base) or by specifying the denominator directly (e.g. 10000 or 32768) intead of the bits/digits. Shift is probably faster than multiply/divide when the denominator is a power of 2. However, given a denominator known at compile time, the optimizer could (no guarantees) choose shifts. Allowing specification of base (or numerator) opens up oddball cases like base 3, etc. This could be constrained by a STATIC_ASSERT, or just allowed to be. Does fixed_pt_10 deserve to be its own class to allow fixed_pt to concentrate on the DSP crowd? I don't know, but I thought I'd throw it out while the discussion is still young. Also, I assume input/output stream operators will come when/if this library becomes a reality? Thanks, Paul Rose

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