views:

4372

answers:

4

Hi there, I was wondering if anyone here knows of any good resources for fixed point math in c#? I've seen things like this (http://2ddev.72dpiarmy.com/viewtopic.php?id=156) and this (http://stackoverflow.com/questions/79677/whats-the-best-way-to-do-fixed-point-math), and a number of discussions about whether decimal is really fixed point or actually floating point (update: responders have confirmed that it's definitely floating point), but I haven't seen a solid C# library for things like calculating cosine and sine.

My needs are simple -- I need the basic operators, plus cosine, sine, arctan2, PI... I think that's about it. Maybe sqrt. I'm programming a 2D RTS game, which I have largely working, but the unit movement when using floating-point math (doubles) has very small inaccuracies over time (10-30 minutes) across multiple machines, leading to desyncs. This is presently only between a 32 bit OS and a 64 bit OS, all the 32 bit machines seem to stay in sync without issue, which is what makes me think this is a floating point issue.

I was aware from this as a possible issue from the outset, and so have limited my use of non-integer position math as much as possible, but for smooth diagonal movement at varying speeds I'm calculating the angle between points in radians, then getting the x and y components of movement with sin and cos. That's the main issue. I'm also doing some calculations for line segment intersections, line-circle intersections, circle-rect intersections, etc, that also probably need to move from floating-point to fixed-point to avoid cross-machine issues.

If there's something open source in Java or VB or another comparable language, I could probably convert the code for my uses. The main priority for me is accuracy, although I'd like as little speed loss over present performance as possible. This whole fixed point math thing is very new to me, and I'm surprised by how little practical information on it there is on google -- most stuff seems to be either theory or dense C++ header files.

Anything you could do to point me in the right direction is much appreciated; if I can get this working, I plan to open-source the math functions I put together so that there will be a resource for other C# programmers out there.

UPDATE: I could definitely make a cosine/sine lookup table work for my purposes, but I don't think that would work for arctan2, since I'd need to generate a table with about 64,000x64,000 entries (yikes). If you know any programmatic explanations of efficient ways to calculate things like arctan2, that would be awesome. My math background is all right, but the advanced formulas and traditional math notation are very difficult for me to translate into code.

+3  A: 

Use 64bit integers in for example 1/1000 scale. You can add and subtract normally. When you need to multiply then multiply integers and then divide by 1000. When you need to sqrt, sin, cos etc. then convert to long double, divide by 1000, sqrt, multiply by 1000, convert to integer. The differences between machines should not matter then.

You can use another scale for faster divides, for example 1024 as x/1024 == x >> 10.

Tometzky
This is a good solution if it fits neatly into your code. You do want to use a power-of-2 shift if at all possible (256 tends to work exceptionally well). Also, you want to use custom fixed-point transcendental functions instead of a slow and error-prone conversion to floating-point and back.
kquinn
That's what I'm worried about, is the shift to floating point and back. I'd like to do sin, cos, etc, calculations in the 64bit integer format itself. Otherwise I like this solution, but I'm afraid the accuracy won't be what I need without that.
x4000
I guess I could create and use sine/cosine tables in that situation, but I actually need atan2, which I don't think I could make a table for.
x4000
I think you should try to benchmark it. Converting an integer to long double should be fast, divide by 1024 also, as it is just manipulating an exponent, sin/cos/tan/atan would not be faster anyway, multiplying by 1024 should be fast, converting to integer also. Differences would be rounded away.
Tometzky
how would I run atan2 on a long, though? The point is that I never want to have an intermediate double, because on 32bit versus 64bit machines there are differences in them, causing the desyncs. Speed is less my focus at this point.
x4000
There will be no differences, as you'd round them to the nearest 1/1000. At least as long as you would not use values larger than about 10^11.
Tometzky
I'm not sure I'd follow. In your statement, is "them" a double? I don't expect to have values too large (not 10^6 even, I don't think), but I've been getting inaccuracies between 64/32bit already with smaller values and doubles, so I'd really like to stay away from those if I can avoid it.
x4000
Your problem is that somefunction(double) on 64bit is different than somefunction(double) on 32bit computer. But round(somefunction(double/1024)*1024) will be the same on both architectures, as differences you see are much smaller than 1/1024.
Tometzky
Your differences, very small at first, accumulate with time. This is why you see them only after several minutes.
Tometzky
+1  A: 

As well as scaled integers, there are a few arbitrary precision numeric libraries which usually include a "BigRational" type, and fixed point is just a fixed power of ten denominator.

Richard
That word "BigRational" is helpful, that gives me something new to google on. I've found libraries like those from www.extremeoptimization.com, but that's way more than I need -- and $1,000, to boot. I'm looking for something no/low cost.
x4000
+13  A: 

Ok, here's what I've come up with for a fixed-point struct, based on the link in my original question but also including some fixes to how it was handling division and multiplication, and added logic for modules, comparisons, shifts, etc:

public struct FInt
{
    public long RawValue;
    public const int SHIFT_AMOUNT = 12; //12 is 4096

    public const long One = 1 << SHIFT_AMOUNT;
    public const int OneI = 1 << SHIFT_AMOUNT;
    public static FInt OneF = FInt.Create( 1, true );

    #region Constructors
    public static FInt Create( long StartingRawValue, bool UseMultiple )
    {
        FInt fInt;
        fInt.RawValue = StartingRawValue;
        if ( UseMultiple )
            fInt.RawValue = fInt.RawValue << SHIFT_AMOUNT;
        return fInt;
    }
    public static FInt Create( double DoubleValue )
    {
        FInt fInt;
        DoubleValue *= (double)One;
        fInt.RawValue = (int)Math.Round( DoubleValue );
        return fInt;
    }
    #endregion

    public int IntValue
    {
        get { return (int)( this.RawValue >> SHIFT_AMOUNT ); }
    }

    public int ToInt()
    {
        return (int)( this.RawValue >> SHIFT_AMOUNT );
    }

    public double ToDouble()
    {
        return (double)this.RawValue / (double)One;
    }

    public FInt Inverse
    {
        get { return FInt.Create( -this.RawValue, false ); }
    }

    #region FromParts
    /// <summary>
    /// Create a fixed-int number from parts.  For example, to create 1.5 pass in 1 and 500.
    /// </summary>
    /// <param name="PreDecimal">The number above the decimal.  For 1.5, this would be 1.</param>
    /// <param name="PostDecimal">The number below the decimal, to three digits.  
    /// For 1.5, this would be 500. For 1.005, this would be 5.</param>
    /// <returns>A fixed-int representation of the number parts</returns>
    public static FInt FromParts( int PreDecimal, int PostDecimal )
    {
        FInt f = FInt.Create( PreDecimal, true );
        if ( PostDecimal != 0 )
            f.RawValue += ( FInt.Create( PostDecimal ) / 1000 ).RawValue;

        return f;
    }
    #endregion

    #region *
    public static FInt operator *( FInt one, FInt other )
    {
        FInt fInt;
        fInt.RawValue = ( one.RawValue * other.RawValue ) >> SHIFT_AMOUNT;
        return fInt;
    }

    public static FInt operator *( FInt one, int multi )
    {
        return one * (FInt)multi;
    }

    public static FInt operator *( int multi, FInt one )
    {
        return one * (FInt)multi;
    }
    #endregion

    #region /
    public static FInt operator /( FInt one, FInt other )
    {
        FInt fInt;
        fInt.RawValue = ( one.RawValue << SHIFT_AMOUNT ) / ( other.RawValue );
        return fInt;
    }

    public static FInt operator /( FInt one, int divisor )
    {
        return one / (FInt)divisor;
    }

    public static FInt operator /( int divisor, FInt one )
    {
        return (FInt)divisor / one;
    }
    #endregion

    #region %
    public static FInt operator %( FInt one, FInt other )
    {
        FInt fInt;
        fInt.RawValue = ( one.RawValue ) % ( other.RawValue );
        return fInt;
    }

    public static FInt operator %( FInt one, int divisor )
    {
        return one % (FInt)divisor;
    }

    public static FInt operator %( int divisor, FInt one )
    {
        return (FInt)divisor % one;
    }
    #endregion

    #region +
    public static FInt operator +( FInt one, FInt other )
    {
        FInt fInt;
        fInt.RawValue = one.RawValue + other.RawValue;
        return fInt;
    }

    public static FInt operator +( FInt one, int other )
    {
        return one + (FInt)other;
    }

    public static FInt operator +( int other, FInt one )
    {
        return one + (FInt)other;
    }
    #endregion

    #region -
    public static FInt operator -( FInt one, FInt other )
    {
        FInt fInt;
        fInt.RawValue = one.RawValue - other.RawValue;
        return fInt;
    }

    public static FInt operator -( FInt one, int other )
    {
        return one - (FInt)other;
    }

    public static FInt operator -( int other, FInt one )
    {
        return (FInt)other - one;
    }
    #endregion

    #region ==
    public static bool operator ==( FInt one, FInt other )
    {
        return one.RawValue == other.RawValue;
    }

    public static bool operator ==( FInt one, int other )
    {
        return one == (FInt)other;
    }

    public static bool operator ==( int other, FInt one )
    {
        return (FInt)other == one;
    }
    #endregion

    #region !=
    public static bool operator !=( FInt one, FInt other )
    {
        return one.RawValue != other.RawValue;
    }

    public static bool operator !=( FInt one, int other )
    {
        return one != (FInt)other;
    }

    public static bool operator !=( int other, FInt one )
    {
        return (FInt)other != one;
    }
    #endregion

    #region >=
    public static bool operator >=( FInt one, FInt other )
    {
        return one.RawValue >= other.RawValue;
    }

    public static bool operator >=( FInt one, int other )
    {
        return one >= (FInt)other;
    }

    public static bool operator >=( int other, FInt one )
    {
        return (FInt)other >= one;
    }
    #endregion

    #region <=
    public static bool operator <=( FInt one, FInt other )
    {
        return one.RawValue <= other.RawValue;
    }

    public static bool operator <=( FInt one, int other )
    {
        return one <= (FInt)other;
    }

    public static bool operator <=( int other, FInt one )
    {
        return (FInt)other <= one;
    }
    #endregion

    #region >
    public static bool operator >( FInt one, FInt other )
    {
        return one.RawValue > other.RawValue;
    }

    public static bool operator >( FInt one, int other )
    {
        return one > (FInt)other;
    }

    public static bool operator >( int other, FInt one )
    {
        return (FInt)other > one;
    }
    #endregion

    #region <
    public static bool operator <( FInt one, FInt other )
    {
        return one.RawValue < other.RawValue;
    }

    public static bool operator <( FInt one, int other )
    {
        return one < (FInt)other;
    }

    public static bool operator <( int other, FInt one )
    {
        return (FInt)other < one;
    }
    #endregion

    public static explicit operator int( FInt src )
    {
        return (int)( src.RawValue >> SHIFT_AMOUNT );
    }

    public static explicit operator FInt( int src )
    {
        return FInt.Create( src, true );
    }

    public static explicit operator FInt( long src )
    {
        return FInt.Create( src, true );
    }

    public static explicit operator FInt( ulong src )
    {
        return FInt.Create( (long)src, true );
    }

    public static FInt operator <<( FInt one, int Amount )
    {
        return FInt.Create( one.RawValue << Amount, false );
    }

    public static FInt operator >>( FInt one, int Amount )
    {
        return FInt.Create( one.RawValue >> Amount, false );
    }

    public override bool Equals( object obj )
    {
        if ( obj is FInt )
            return ( (FInt)obj ).RawValue == this.RawValue;
        else
            return false;
    }

    public override int GetHashCode()
    {
        return RawValue.GetHashCode();
    }

    public override string ToString()
    {
        return this.RawValue.ToString();
    }
}

public struct FPoint
{
    public FInt X;
    public FInt Y;

    public static FPoint Create( FInt X, FInt Y )
    {
        FPoint fp;
        fp.X = X;
        fp.Y = Y;
        return fp;
    }

    public static FPoint FromPoint( Point p )
    {
        FPoint f;
        f.X = (FInt)p.X;
        f.Y = (FInt)p.Y;
        return f;
    }

    public static Point ToPoint( FPoint f )
    {
        return new Point( f.X.IntValue, f.Y.IntValue );
    }

    #region Vector Operations
    public static FPoint VectorAdd( FPoint F1, FPoint F2 )
    {
        FPoint result;
        result.X = F1.X + F2.X;
        result.Y = F1.Y + F2.Y;
        return result;
    }

    public static FPoint VectorSubtract( FPoint F1, FPoint F2 )
    {
        FPoint result;
        result.X = F1.X - F2.X;
        result.Y = F1.Y - F2.Y;
        return result;
    }

    public static FPoint VectorDivide( FPoint F1, int Divisor )
    {
        FPoint result;
        result.X = F1.X / Divisor;
        result.Y = F1.Y / Divisor;
        return result;
    }
    #endregion
}

Based on the comments from ShuggyCoUk, I see that this is in Q12 format. That's reasonably precise for my purposes. Of course, aside from the bugfixes, I already had this basic format before I asked my question. What I was looking for were ways to calculate Sqrt, Atan2, Sin, and Cos in C# using a structure like this. There aren't any other things that I know of in C# that will handle this, but in Java I managed to find the MathFP library by Onno Hommes. It's a liberal source software license, so I've converted some of his functions to my purposes in C# (with a fix to atan2, I think). Enjoy:

    #region PI, DoublePI
    public static FInt PI = FInt.Create( 12868, false ); //PI x 2^12
    public static FInt TwoPIF = PI * 2; //radian equivalent of 260 degrees
    public static FInt PIOver180F = PI / (FInt)180; //PI / 180
    #endregion

    #region Sqrt
    public static FInt Sqrt( FInt f, int NumberOfIterations )
    {
        if ( f.RawValue < 0 ) //NaN in Math.Sqrt
            throw new ArithmeticException( "Input Error" );
        if ( f.RawValue == 0 )
            return (FInt)0;
        FInt k = f + FInt.OneF >> 1;
        for ( int i = 0; i < NumberOfIterations; i++ )
            k = ( k + ( f / k ) ) >> 1;

        if ( k.RawValue < 0 )
            throw new ArithmeticException( "Overflow" );
        else
            return k;
    }

    public static FInt Sqrt( FInt f )
    {
        byte numberOfIterations = 8;
        if ( f.RawValue > 0x64000 )
            numberOfIterations = 12;
        if ( f.RawValue > 0x3e8000 )
            numberOfIterations = 16;
        return Sqrt( f, numberOfIterations );
    }
    #endregion

    #region Sin
    public static FInt Sin( FInt i )
    {
        FInt j = (FInt)0;
        for ( ; i < 0; i += FInt.Create( 25736, false ) ) ;
        if ( i > FInt.Create( 25736, false ) )
            i %= FInt.Create( 25736, false );
        FInt k = ( i * FInt.Create( 10, false ) ) / FInt.Create( 714, false );
        if ( i != 0 && i != FInt.Create( 6434, false ) && i != FInt.Create( 12868, false ) && 
            i != FInt.Create( 19302, false ) && i != FInt.Create( 25736, false ) )
            j = ( i * FInt.Create( 100, false ) ) / FInt.Create( 714, false ) - k * FInt.Create( 10, false );
        if ( k <= FInt.Create( 90, false ) )
            return sin_lookup( k, j );
        if ( k <= FInt.Create( 180, false ) )
            return sin_lookup( FInt.Create( 180, false ) - k, j );
        if ( k <= FInt.Create( 270, false ) )
            return sin_lookup( k - FInt.Create( 180, false ), j ).Inverse;
        else
            return sin_lookup( FInt.Create( 360, false ) - k, j ).Inverse;
    }

    private static FInt sin_lookup( FInt i, FInt j )
    {
        if ( j > 0 && j < FInt.Create( 10, false ) && i < FInt.Create( 90, false ) )
            return FInt.Create( SIN_TABLE[i.RawValue], false ) + 
                ( ( FInt.Create( SIN_TABLE[i.RawValue + 1], false ) - FInt.Create( SIN_TABLE[i.RawValue], false ) ) / 
                FInt.Create( 10, false ) ) * j;
        else
            return FInt.Create( SIN_TABLE[i.RawValue], false );
    }

    private static int[] SIN_TABLE = {
        0, 71, 142, 214, 285, 357, 428, 499, 570, 641, 
        711, 781, 851, 921, 990, 1060, 1128, 1197, 1265, 1333, 
        1400, 1468, 1534, 1600, 1665, 1730, 1795, 1859, 1922, 1985, 
        2048, 2109, 2170, 2230, 2290, 2349, 2407, 2464, 2521, 2577, 
        2632, 2686, 2740, 2793, 2845, 2896, 2946, 2995, 3043, 3091, 
        3137, 3183, 3227, 3271, 3313, 3355, 3395, 3434, 3473, 3510, 
        3547, 3582, 3616, 3649, 3681, 3712, 3741, 3770, 3797, 3823, 
        3849, 3872, 3895, 3917, 3937, 3956, 3974, 3991, 4006, 4020, 
        4033, 4045, 4056, 4065, 4073, 4080, 4086, 4090, 4093, 4095, 
        4096
    };
    #endregion

    private static FInt mul( FInt F1, FInt F2 )
    {
        return F1 * F2;
    }

    #region Cos, Tan, Asin
    public static FInt Cos( FInt i )
    {
        return Sin( i + FInt.Create( 6435, false ) );
    }

    public static FInt Tan( FInt i )
    {
        return Sin( i ) / Cos( i );
    }

    public static FInt Asin( FInt F )
    {
        bool isNegative = F < 0;
        F = Abs( F );

        if ( F > FInt.OneF )
            throw new ArithmeticException( "Bad Asin Input:" + F.ToDouble() );

        FInt f1 = mul( mul( mul( mul( FInt.Create( 145103 >> FInt.SHIFT_AMOUNT, false ), F ) -
            FInt.Create( 599880 >> FInt.SHIFT_AMOUNT, false ), F ) +
            FInt.Create( 1420468 >> FInt.SHIFT_AMOUNT, false ), F ) -
            FInt.Create( 3592413 >> FInt.SHIFT_AMOUNT, false ), F ) +
            FInt.Create( 26353447 >> FInt.SHIFT_AMOUNT, false );
        FInt f2 = PI / FInt.Create( 2, true ) - ( Sqrt( FInt.OneF - F ) * f1 );

        return isNegative ? f2.Inverse : f2;
    }
    #endregion

    #region ATan, ATan2
    public static FInt Atan( FInt F )
    {
        return Asin( F / Sqrt( FInt.OneF + ( F * F ) ) );
    }

    public static FInt Atan2( FInt F1, FInt F2 )
    {
        if ( F2.RawValue == 0 && F1.RawValue == 0 )
            return (FInt)0;

        FInt result = (FInt)0;
        if ( F2 > 0 )
            result = Atan( F1 / F2 );
        else if ( F2 < 0 )
        {
            if ( F1 >= 0 )
                result = ( PI - Atan( Abs( F1 / F2 ) ) );
            else
                result = ( PI - Atan( Abs( F1 / F2 ) ) ).Inverse;
        }
        else
            result = ( F1 >= 0 ? PI : PI.Inverse ) / FInt.Create( 2, true );

        return result;
    }
    #endregion

    #region Abs
    public static FInt Abs( FInt F )
    {
        if ( F < 0 )
            return F.Inverse;
        else
            return F;
    }
    #endregion

There are a number of other functions in Dr. Hommes' MathFP library, but they were beyond what I needed, and so I have not taken the time to translate them to C# (that process was made extra difficult by the fact that he was using a long, and I am using the FInt struct, which makes the conversion rules are a bit challenging to see immediately).

The accuracy of these functions as they are coded here is more than enough for my purposes, but if you need more you can increase the SHIFT AMOUNT on FInt. Just be aware that if you do so, the constants on Dr. Hommes' functions will then need to be divided by 4096 and then multiplied by whatever your new SHIFT AMOUNT requires. You're likely to run into some bugs if you do that and aren't careful, so be sure to run checks against the built-in Math functions to make sure that your results aren't being put off by incorrectly adjusting a constant.

So far, this FInt logic seems as fast, if not perhaps a bit faster, than the equivalent built in .net functions. That would obviously vary by machine, since the fp coprocessor would determine that, so I have not run specific benchmarks. But they are integrated into my game now, and I've seen a slight decrease in processor utilization compared to before (this is on a Q6600 quad core -- about a 1% drop in usage on average).

Thanks again to everyone who commented for your help. No one pointed me directly to what I was looking for, but you gave me some clues that helped me find it myself on google. I hope this code turns out to be useful for someone else, since there doesn't seem to be anything comparable in C# posted publicly.

x4000
It seems you calculated the asin value from the taylor series of the atan. And when you implement atan you calculated it from the asin function. I suggest calculating atan with taylor, because it is simpler and the calculate asin with atan.
Calmarius
Thank you! :) Very helpful library
Walt W
A: 

I created a similar fixedpoint struct. You get a performance hit using new() because it puts data onto the heap even though you are using a struct. See Google(C# Heap(ing) Vs Stack(ing) in .NET: Part I) the true power of using struct is the ablity to not use new and pass by value to the stack. My example below does the following on stack. 1. [result int ] on stack 2. [a int ] on stack 3. [b int ] on stack 4. [* ] operator on stack 5. value result returned no heap allocated costs.

    public static Num operator *(Num a, Num b)
    {
        Num result;
        result.NumValue = a.NumValue * b.NumValue;
        return result;
    }
Very interesting -- in profiling, I had noticed the performance hit on New, during the math operators in particular. I had not yet investigated it to find out precisely why. Very good point about this, I had forgotten that about new with structs. I have updated the code above to get rid of new, and the result is definitely a bit faster and less pointless load on the heap, too. Thanks!
x4000
-1: your code doesn't compile. You need `Num result = new Num();`. This is C#, not C++.
John Saunders