Experilous

In beginning the simulation code of my city builder over the past few weeks, I’ve run into some difficulties with integers that quite surprised me.  For a variety of reasons, I want to ensure that the simulation code is absolutely completely deterministic, regardless of platform, compiler, machine, or any other variable.  While I’d need to do further research to properly justify this statement in all its details, my impression is that this is much easier to guarantee if one sticks to integer computations, avoiding floating point representations of numbers.

Of course, if one uses just integers in all computations (or even the somewhat fancier fixed point variants), then loss of precision becomes more important to deal with, whereas with floating point numbers, one can often get by ignoring such details.  Division in particular is the prime culprit.

I noticed this issue when implementing a basic economic simulation (based on this article by Lars Doucet, which is in turn based on this research paper by Jonathon Doran).  It is designed to let supply and demand reasonably influence price beliefs of agents, and thus actual agreed upon prices within the simulation.  At an early stage of my simulation, though, I had very little which would meaningful dictate appropriate prices.  (Essentially, supply of everything was infinite.)  My hope was that prices would just fluctuate around their initial values, as there was no pressure for them to move in any specific direction or stabilize around a specific balanced value.  Unfortunately, they grew exponentially as the simulation ran. In fact, road construction never completed because the greater the price the city offered to pay to have the roads constructed, the more the construction company demanded for the work, so no agreement was ever reached between the two parties. Oops.

There were a variety of sloppy bugs in the code that were contributing to this effect, but even with those solved, I was noticing problems (or at least potential problems) with some of the price calculations.  I was just using basic integer division whenever I needed to find average values, which according to the C99 and C++11 specifications must round toward zero.  This consistent direction was introducing a bias into the calculations which was causing an inappropriate positive pressure of the price.

Well, there’s an easy fix for that:  Use a less biased rounding scheme.  The most common is to round 0.5 to the nearest even number.  Sounds easy enough.  I also had a few places in code where I wanted to always round up or round down (regardless of sign).  So I thought I’d hunt down some standard implementations of these operations for integer division and rounding.  Unfortunately I couldn’t find any.  I did find an attempt at a few of them on Stack Overflow, but it didn’t provide a solution for round-to-nearest-even, and the implementations all had branches, which my intuition suggested was unnecessary.

In the end, I implemented my own. They do work for all general combinations of positive and negative numbers (+/+, -/+, +/-, and -/-).  I have no idea if they word for unsigned types. Beyond some simple unit tests, I have not done any extensive testing, especially in the realm of performance or machine code analysis.  But they might serve as a starting point for someone who wants to take them further, or even as a satisfactory solution for anyone else who needs similar functionality.  I might come back to them in the future, but for now they do exactly what I need in my prototype code, so I’ve moved on.

Anyway, enough rambling.  Here’s some code.  (For what it’s worth, I declare this code to be public domain.) Credit for the Sign() function goes to user79758 on stackoverflow.

template <typename TInteger>
TInteger Sign(TInteger value)
{
    return (0 < value) - (value < 0);
}

template <typename TInteger>
TInteger DivisionHalfTowardNearestEven(TInteger numerator, TInteger denominator)
{
    TInteger numeratorDoubled = numerator * 2;
    TInteger quotientDoubled = numeratorDoubled / denominator;
    TInteger remainderDoubled = numeratorDoubled - quotientDoubled * denominator;
    //Round to nearest, with half always roundup up.
    TInteger quotient = (quotientDoubled + Sign(quotientDoubled)) / 2;
    TInteger remainder = numerator - quotient * denominator;
    //Given that half is always rounded up, we sometimes need to counter that and round down instead.
    //Subtract 1 whenever remainderDoubled was zero (implying that we were exactly at either n.0 or n.5),
    // remainder was not zero (implying that we were not exactly at n.0),
    // and our rounded up quotient was odd (we want even in this case).
    return quotient - (remainderDoubled == 0 && remainder != 0) * (quotient % 2);
}

template <typename TInteger>
TInteger CeilingDivision(TInteger numerator, TInteger denominator)
{
    TInteger quotient = numerator / denominator;
    TInteger remainder = numerator - quotient * denominator;
    //By default we round toward zero, which for negative numbers means rounding up.
    //If we have a remainder and a positive real-value quotient, we add 1 to compensate for rounding down.
    return quotient + (remainder * denominator > 0);
}

template <typename TInteger>
TInteger FloorDivision(TInteger numerator, TInteger denominator)
{
    TInteger quotient = numerator / denominator;
    TInteger remainder = numerator - quotient * denominator;
    //By default we round to zero, which for positive numbers means rounding down.
    //If we have a remainder and a negative real-value quotient, we subtract 1 to compensate for rounding up.
    return quotient - (remainder * denominator < 0);
}

template <typename TInteger>
TInteger DivisionTowardZero(TInteger numerator, TInteger denominator)
{
    return numerator / denominator;
}

template <typename TInteger>
TInteger DivisionTowardInfinity(TInteger numerator, TInteger denominator)
{
    TInteger quotient = numerator / denominator;
    TInteger remainder = numerator - quotient * denominator;
    //By default we round to zero, which is exactly the opposite that we want.
    //If we have a remainder, we compensate by adding positive or negative 1 to move away from zero.
    return quotient + Sign(remainder * denominator);
}

1 Comment

Comment by Brandon M. Petty — 2014/06/04 @ 12:36

Rounding based on a boolean expression is interesting. You could also do a sanity check for a denominator of Zero and return a number you optionally pass in representing “infinity”… like 0.

Leave a comment