Rounding Algorithms

[Duthomhas]

There are a zillion different ways to round floating point values to integers. C and C++ provide a couple basic ones in <math.h> or <cmath>.

There are two general categories of rounding algorithm: those that are symmetric about zero and those that are biased in some way.

Biased Rounding
A bias is a mathematical notion used in statistics to indicate that samples are not exactly representative of their true values; they are skewed in one way or another.

For example, the <cmath> floor() function is biased towards negative infinity, because it always chooses the lower integer number --that is, it always chooses the number closer to negative infinity.

floor( 7.5 ) --> 7
floor( -7.5 ) --> -8


Suppose your local big city and wanted to know how fast people are driving on a particular freeway. The first step is to gather the exact speeds that individual drivers are going, and the second would be to convert all the individual values into a single value that would represent the normal rate of speed. For simplicity here, we will just use the average value.

Let us say that the equipment which samples a motorist's speed is far more accurate than the computer's ability to store it. (This is actually not uncommon.) Again, for simplicity, we will say that the computer stores the speed as integers.

With a sampling of twelve motorists, we get the following speeds (in miles per hour)

49.087 57.901 28.500 46.738 51.270 53.096
44.795 47.218 46.347 45.989 47.582 50.563

A quick calculation shows that the average speed is 47.424 mph.

If the city were to simply use the floor() function to convert these to integers, it would get

49 57 28 46 51 53
44 47 46 45 47 50

which averages to 46.916 --> 46 mph (remember, integer arithmetic!)

Either way, the sampling is off by about a whole mile per hour. I don't think that the city would actually care about a single mile per hour, but this does illustrate the bias, or tendancy of the floor() function, to make the numbers closer to negative infinity, thuse skewing the data to an inaccurate number.

This was just a simple example that came off the top of my head, but in many sciences and statistical surveys, that difference can mean quite a lot. Suppose that the Apollo missed the moon by 1%? Suppose a pharmaceutical company put too much iron in a daily vitamin pill by 1%? Suppose a construction firm miscalculated the stresses a bridge can take by 1%? In all these scenarios the results would prove deadly. One percent is a lot.


Symmetric Rounding
A special case of bias is centered about zero. Let us fix the floor() function to tend towards zero.

double floor0( double value )
  {
  if (value < 0.0)
    return ceil( value );
  else
    return floor( value );
  }

Now, the absolute value of the result will always be the same:

floor0( -7.7 ) --> -7 floor0( 7.7 ) --> 7
floor0( -7.5 ) --> -7 floor0( 7.5 ) --> 7
floor0( -7.3 ) --> -7 floor0( 7.3 ) --> 7

Enough about that.


Unbiased Rounding
So, how do we handle these biases? By adding some rule that accounts for it.

Let us apply knowledge we all learned in gradeschool: In arithmetic rounding we round up if the next digit is 5 or more, and if it less than 5 we round down. We write ourselves a little function to do just that:

double round( double value )
  {
  return floor( value + 0.5 );
  }

The problem is that this is still biased. We have actually reversed the bias of floor() from negative infinity to positive infinity, because we always choose to round up when exactly halfway between two values:

round( 10.3 ) --> 10
round( 10.5 ) --> 11
round( 10.7 ) --> 11

You can actually see the bias in the above table: the result tends towards 11 and away from 10.

This brings us to the trick: which way do we round when exactly halfway between two values?

One very popular method is variously called "bankers rounding", "round to even", "convergent rounding", and even "unbiased rounding", to name a few. It works by skewing the bias itself.

Given a number exactly halfway between two values, round to the even value (zero is considered even here).

round( 1.7 ) --> 2 round( 2.7 ) --> 3
round( 1.5 ) --> 2 round( 2.5 ) --> 2
round( 1.3 ) --> 1 round( 2.3 ) --> 2

For random data this is very convenient. Bankers like it because money deposited and withdrawn is random. (There are trends, mind you, but you cannot predict exactly how much will be deposited and withdrawn.) The important point is that the bankers rounding is still biased if the data is biased. It is only unbiased with random data.

One solution is called "alternate rounding". It works by simply choosing to bias up or down every other time.

round( 1.5 ) --> 2
round( 1.5 ) --> 1
round( 1.5 ) --> 2
round( 1.5 ) --> 1
etc

This is not always useful though.

The only way to eliminate all bias is to use a random bias... This, of course, is impossible to generate in your typical PC, but it still goes toward solving the problem quite nicely.

If the sample is exactly halfway between two integers, it chooses one or the other randomly.

Of course, the Achilles heel of this method is the random number generator you use. The default pseudorandom generator for C and C++ is not that great. The Mersenne Twister is by far the most popular high-quality pseudorandom number generator, but it is non-trivial to implement so I will not include it below.

Anyway, what follows is a convenient, simple library you are free to use. I'll even permit you to cannabalize it at will (since the algorithms are so obvious...)

I may update it sometime in the future if I can figure out a way around the default epsilon issue. Feel free to make suggestions for improvements.

:-)

[Duthomhas]

rounding-algorithms.hpp

// rounding-algorithms.hpp
//
// General Rounding Algorithms
// Copyright (c) 2008 Michael Thomas Greer
// Boost Software License - Version 1.0 - August 17th, 2003
//
// Permission is hereby granted, free of charge, to any person or organization
// obtaining a copy of the software and accompanying documentation covered by
// this license (the "Software") to use, reproduce, display, distribute,
// execute, and transmit the Software, and to prepare derivative works of the
// Software, and to permit third-parties to whom the Software is furnished to
// do so, all subject to the following:
//
// The copyright notices in the Software and this entire statement, including
// the above license grant, this restriction and the following disclaimer,
// must be included in all copies of the Software, in whole or in part, and
// all derivative works of the Software, unless such copies or derivative
// works are solely in the form of machine-executable object code generated by
// a source language processor.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
// SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
// FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
//----------------------------------------------------------------------------
// Reference
// <http://www.pldesignline.com/howto/showArticle.jhtml;?articleID=175801189>
//
//----------------------------------------------------------------------------
// In this library, symmetric functions are indicated by a zero at the end
// of the function name.
//
// If you want a different default epsilon make sure to change
//
//   #define ROUNDING_EPSILON 0.001
//
// to whatever you want it to be. (I wanted to make it so that you could 
// define a different default epsilon each time you #included the file, but
// I haven't figured out how to get around the template restrictions yet.)
//

#ifndef ROUNDING_ALGORITHMS_HPP
#define ROUNDING_ALGORITHMS_HPP

#ifndef ROUNDING_EPSILON
#define ROUNDING_EPSILON 0.0000001
#endif

#include <cmath>
#include <cstdlib>
#include <ciso646>

namespace rounding
  {

  //--------------------------------------------------------------------------
  // round down
  // Bias: -Infinity
  using std::floor;

  //--------------------------------------------------------------------------
  // round up
  // Bias: +Infinity
  using std::ceil;

  //--------------------------------------------------------------------------
  // symmetric round down
  // Bias: towards zero
  template <typename FloatType>
  FloatType floor0( const FloatType& value )
    {
    FloatType result = std::floor( std::fabs( value ) );
    return (value < 0.0) ? -result : result;
    }

  //--------------------------------------------------------------------------
  // A common alias for floor0()
  // (notwithstanding hardware quirks)
  template <typename FloatType>
  inline
  FloatType trunc( const FloatType& value )
    {
    return floor0( value );
    }

  //--------------------------------------------------------------------------
  // symmetric round up
  // Bias: away from zero
  template <typename FloatType>
  FloatType ceil0( const FloatType& value )
    {
    FloatType result = std::ceil( std::fabs( value ) );
    return (value < 0.0) ? -result : result;
    }

  //--------------------------------------------------------------------------
  // Common rounding: round half up
  // Bias: +Infinity
  template <typename FloatType>
  FloatType roundhalfup( const FloatType& value )
    {
    return std::floor( value +0.5 );
    }

  //--------------------------------------------------------------------------
  // Round half down
  // Bias: -Infinity
  template <typename FloatType>
  FloatType roundhalfdown( const FloatType& value )
    {
    return std::ceil( value -0.5 );
    }

  //--------------------------------------------------------------------------
  // symmetric round half down
  // Bias: towards zero
  template <typename FloatType>
  FloatType roundhalfdown0( const FloatType& value )
    {
    FloatType result = roundhalfdown( std::fabs( value ) );
    return (value < 0.0) ? -result : result;
    }

  //--------------------------------------------------------------------------
  // symmetric round half up
  // Bias: away from zero
  template <typename FloatType>
  FloatType roundhalfup0( const FloatType& value )
    {
    FloatType result = roundhalfup( std::fabs( value ) );
    return (value < 0.0) ? -result : result;
    }

  //--------------------------------------------------------------------------
  // round half even (banker's rounding)
  // Bias: none
  template <typename FloatType>
  FloatType roundhalfeven(
    const FloatType& value,
    const FloatType& epsilon = ROUNDING_EPSILON
    ) {
    if (value < 0.0) return -roundhalfeven <FloatType> ( -value, epsilon );

    FloatType ipart;
    std::modf( value, &ipart );

    // If 'value' is exctly halfway between two integers
    if ((value -(ipart +0.5)) < epsilon)
      {
      // If 'ipart' is even then return 'ipart'
      if (std::fmod( ipart, 2.0 ) < epsilon)
        return ipart;

      // Else return the nearest even integer
      return ceil0( ipart +0.5 );
      }

    // Otherwise use the usual round to closest
    // (Either symmetric half-up or half-down will do0
    return roundhalfup0( value );
    }

  //--------------------------------------------------------------------------
  // round alternate
  // Bias: none for sequential calls
  bool _is_up = false;
  template <typename FloatType>
  FloatType roundalternate( const FloatType& value, int& is_up = _is_up )
    {
    if ((is_up != is_up))
      return roundhalfup( value );
    return roundhalfdown( value );
    }

  //--------------------------------------------------------------------------
  // symmetric round alternate
  // Bias: none for sequential calls
  template <typename FloatType>
  FloatType roundalternate0( const FloatType& value, int& is_up = _is_up )
    {
    if ((is_up != is_up))
      return roundhalfup0( value );
    return roundhalfdown0( value );
    }

  //--------------------------------------------------------------------------
  // round random
  // Bias: generator's bias
  template <typename FloatType, typename RandValue, typename RandomGenerator>
  FloatType roundrandom(
    const FloatType& value,
    const RandValue& mid,
    RandomGenerator& g
    ) {
    if (g() < mid)
      return roundhalfup0( value );
    return roundhalfdown0( value );
    }

  //--------------------------------------------------------------------------
  // default round random
  // Bias: rand()
  template <typename FloatType>
  FloatType roundrandom( const FloatType& value )
    {
    return roundrandom <FloatType, int, int(*)()> ( value, RAND_MAX /2, &rand );
    }
  }

#endif 

If you find a bug let me know!

[Mitsakos]

The roundhalfdown function doesn't work correctly. If you pass the number 12.5 to the function then you will have as a result 12, but using math it should be 13.
I think that this comes with the ceil() function. It rounds up only if the number after the decimal is greater than 0. When you have exactly 0 as decimal it doesn't.

Here is how I tried it:

#include "rounding-algorithms.h"
#include <iostream>
using namespace std;
using namespace rounding;

int main(){
   double num = 12.5;
   double num1 = 12.500001;
   cout << "12.5 after rounded: " << roundhalfdown<double>( num ) << '\n';
   cout << "12.500001 after rounded: " << roundhalfdown<double>( num1 ) << '\n';
}

If I have something wrong let me know. thanks

[helios]

The one true way to round a number:

//Named after an equivalent BASIC function.
long cint(double x){
	if (modf(x,0)>=.5)
		return x>=0?ceil(x):floor(x);
	else
		return x<0?ceil(x):floor(x);
}

[Duthomhas]

The round half down did exactly what it should have. Given a number exactly halfway between two integers (12.5) it rounds down to 12.

As an unrelated point of interest, if you are passing typed arguments to a template function you don't actually have to give template arguments:

double num = 12.5;
cout << roundhalfdown( num ) << endl;

num = 12.500001;
cout << roundhalfdown( num ) << endl;

One thing to remember about floating point numbers is that they are not exact. In fact, a lot more can go wrong with them than people realize. So, depending on your hardware and your clib's temperament that second one may or may not be considered the same as the first. IEEE single precision should handle it fine though.