Rounding Algorithms from Property Based Testing

Steffen Forkmann recently posted a blog about incorrect rounding in a twitter poll, and how coding a rounding strategy is a hard problem. This got me thinking about how many correct rounding algorithms there are.

In these type of problems it is important to look at what properties the algorithm should have. Property based testing is a great tool when doing this.

The key property required for a fair rounding algorithm is that rounded values increase with the weights. It doesn't make sense for a lower weight to have a greater rounded value. Symmetry in the results for negative weights and negative value to be distributed are also important. This can easily be achieved by mapping from the positive results, but a robust algorithm shouldn't need to do this.

In the blog it was proposed that adjusting the largest weight would work, but in general this can only work when the rounding needs a positive adjustment due to the increasing with weight property. For negative adjustments the smallest non-zero weight would need to be adjusted. It may also be unfair to adjust these values if they already have a large rounding error.

This algorithm minimises the absolute and relative rounding errors. I normally apply absolute then relative but the reverse order also works and may be more correct for certain problems.

/// Distribute integer n over an array of weights let n weights = let wt = Array . sum weights if Array . isEmpty weights || Seq . maxBy weights * ( n + 1 ) |> > = wt * Int32 . MaxValue then None else let inline f = ( if f < 0.0 then f - 0.5 else f + 0.5 ) let ns = Array . map ( (*) ( n / wt ) > > ) weights let inline wi ni d = let wc = 0.5 * d * wt / n let wni = ni * wt / n if wc > 0.0 then ( wni - wi ) 0.0 + wc else ( wi - wni ) 0.0 - wc let d = n - Array . sum ns for __ = 1 to d do let _ , _ , _ , i = Seq . mapi2 ( fun i wi ni -> let absErr = wi ni ( d ) let relErr = absErr / wi let weight = if wt > 0.0 then - wi else wi absErr , relErr , weight , i ) weights ns |> Seq . min ns . [ i ] <- ns . [ i ] + d Some ns

The twitter tricky test below is interesting. It is not clear which values should be adjusted down. Neither the largest or smallest weights look like good candidates. The error minimisation algorithm sensibly selects the second largest weight and keeps the correct order.

Testing also shows the algorithm is sensitive to weights that are close together and issues with machine precision. This directs how the error function needs to calculate.

let roundingTests = testList "rounding" [ test "empty" { let r = distribute 1 [| |] Expect . equal r None "empty" } test "n zero" { let r = distribute 0 [| 406.0 ; 348.0 ; 246.0 ; 0.0 |] Expect . equal r ( Some [| 0 ; 0 ; 0 ; 0 |] ) "zero" } test "twitter" { let r = distribute 100 [| 406.0 ; 348.0 ; 246.0 ; 0.0 |] Expect . equal r ( Some [| 40 ; 35 ; 25 ; 0 |] ) "40 etc" } test "twitter n negative" { let r = distribute - 100 [| 406.0 ; 348.0 ; 246.0 ; 0.0 |] Expect . equal r ( Some [| - 40 ; - 35 ; - 25 ; 0 |] ) "-40 etc" } test "twitter weights negative" { let r = distribute 100 [| - 406.0 ; - 348.0 ; - 246.0 ; - 0.0 |] Expect . equal r ( Some [| 40 ; 35 ; 25 ; 0 |] ) "40 etc" } test "twitter both negative" { let r = distribute - 100 [| - 406.0 ; - 348.0 ; - 246.0 ; - 0.0 |] Expect . equal r ( Some [| - 40 ; - 35 ; - 25 ; 0 |] ) "-40 etc" } test "twitter tricky" { let r = distribute 100 [| 404.0 ; 397.0 ; 57.0 ; 57.0 ; 57.0 ; 28.0 |] Expect . equal r ( Some [| 40 ; 39 ; 6 ; 6 ; 6 ; 3 |] ) "o no" } test "negative example" { let r1 = distribute 42 [| 1.5 ; 1.0 ; 39.5 ; - 1.0 ; 1.0 |] Expect . equal r1 ( Some [| 2 ; 1 ; 39 ; - 1 ; 1 |] ) "2 etc" let r2 = distribute - 42 [| 1.5 ; 1.0 ; 39.5 ; - 1.0 ; 1.0 |] Expect . equal r2 ( Some [| - 2 ; - 1 ; - 39 ; 1 ; - 1 |] ) "-2 etc" } testProp "ni total correctly" ( fun n ( Gen . RationalFloats ws ) -> distribute n ws |> Option . iter ( fun ns -> Expect . equal ( Array . sum ns ) n "sum ns = n" ) ) testProp "negative n returns opposite of positive n" ( fun n ( Gen . RationalFloats ws ) -> let r1 = distribute - n ws |> Option . map ( Array . map ( ~- ) ) let r2 = distribute n ws Expect . equal r1 r2 "r1 = r2" ) testProp "increase with weight" ( fun n ( Gen . RationalFloats ws ) -> let d = if Seq . sum ws > 0.0 <> ( n > 0 ) then - 1 else 1 distribute n ws |> Option . iter ( Seq . map ( (*) d ) > > Seq . zip ws > > Seq . sort > > Seq . pairwise > > Seq . iter ( fun ( ni1 , ni2 ) -> Expect . isLessThanOrEqual ni1 ni2 "ni1 <= ni2" ) ) ) testProp "smallest error" ( fun n ( Gen . RationalFloats ws ) randomChanges -> distribute n ws |> Option . iter ( fun ns -> let totalError ns = let wt = Seq . sum ws Seq . map2 ( fun ni wi -> float ni * wt / float n - wi |> abs ) ns ws |> Seq . sum let errorBefore = totalError ns let l = Array . length ns List . iter ( fun ( i , j ) -> ns . [ abs i % l ] <- ns . [ abs i % l ] - 1 ns . [ abs j % l ] <- ns . [ abs j % l ] + 1 ) randomChanges let errorAfter = totalError ns Expect . floatLessThanOrClose Accuracy . veryHigh errorBefore errorAfter "before <= after" ) ) ]

This is an example of how property based testing can actually help in algorithm design. It gives example failing cases that can steer you to a better solution.

I have seen this problem in order management systems where orders for a number of shares are to be allocated across a number of portfolios. The buy and sell orders have a number of partial fills, but ultimately everything needs to add up in a consistent and robust way.

Error minimisation is the best rounding algorithm I have found. By construction it is also the one with the smallest total rounding error.

I don't think there is any other general purpose algorithm that can satisfy these properties as well.

module Rounding

namespace System

val distribute : n:int -> weights:float [] -> int [] option





Distribute integer n over an array of weights

val n : int

val weights : float []

val wt : float

type Array =

member Clone : unit -> obj

member CopyTo : array:Array * index:int -> unit + 1 overload

member GetEnumerator : unit -> IEnumerator

member GetLength : dimension:int -> int

member GetLongLength : dimension:int -> int64

member GetLowerBound : dimension:int -> int

member GetUpperBound : dimension:int -> int

member GetValue : [<ParamArray>] indices:int[] -> obj + 7 overloads

member Initialize : unit -> unit

member IsFixedSize : bool

...

val sum : array:'T [] -> 'T (requires member ( + ) and member get_Zero)

val isEmpty : array:'T [] -> bool

module Seq



from Microsoft.FSharp.Collections

val maxBy : projection:('T -> 'U) -> source:seq<'T> -> 'T (requires comparison)

val abs : value:'T -> 'T (requires member Abs)

Multiple items

val float : value:'T -> float (requires member op_Explicit)



--------------------

type float = Double



--------------------

type float<'Measure> = float

type Int32 =

struct

member CompareTo : value:obj -> int + 1 overload

member Equals : obj:obj -> bool + 1 overload

member GetHashCode : unit -> int

member GetTypeCode : unit -> TypeCode

member ToString : unit -> string + 3 overloads

static val MaxValue : int

static val MinValue : int

static member Parse : s:string -> int + 3 overloads

static member TryParse : s:string * result:int -> bool + 1 overload

end

field int.MaxValue: int = 2147483647

union case Option.None: Option<'T>

val round : (float -> int)

val f : float

Multiple items

val int : value:'T -> int (requires member op_Explicit)



--------------------

type int = int32



--------------------

type int<'Measure> = int

val ns : int []

val map : mapping:('T -> 'U) -> array:'T [] -> 'U []

val absoluteError : (float -> 'a -> 'b -> float) (requires member op_Explicit and member op_Explicit)

val wi : float

val ni : 'a (requires member op_Explicit)

val d : 'b (requires member op_Explicit)

val wc : float

val wni : float

val min : e1:'T -> e2:'T -> 'T (requires comparison)

val d : int

val __ : int

val i : int

val mapi2 : mapping:(int -> 'T1 -> 'T2 -> 'U) -> source1:seq<'T1> -> source2:seq<'T2> -> seq<'U>

val ni : int

val absErr : float

val sign : value:'T -> int (requires member get_Sign)

val relErr : float

val weight : float

val min : source:seq<'T> -> 'T (requires comparison)

union case Option.Some: Value: 'T -> Option<'T>

Multiple items

union case RationalFloats.RationalFloats: float [] -> RationalFloats



--------------------

type RationalFloats = | RationalFloats of float []

val rationalFloats : obj

val fraction : (int -> int -> int -> float)

val a : int

val b : int

val c : int

val private config : obj

val typeof<'T> : Type

module Gen



from Rounding

Multiple items

union case Gen.RationalFloats.RationalFloats: float [] -> Gen.RationalFloats



--------------------

type RationalFloats = | RationalFloats of float []

val testProp : name:'a -> 'b

val name : 'a

val ptestProp : name:'a -> 'b

val ftestProp : stdgen:'a -> name:'b -> 'c

val stdgen : 'a

val name : 'b

val roundingTests : obj

module Option



from Microsoft.FSharp.Core

val iter : action:('T -> unit) -> option:'T option -> unit

val map : mapping:('T -> 'U) -> option:'T option -> 'U option

val sum : source:seq<'T> -> 'T (requires member ( + ) and member get_Zero)

val map : mapping:('T -> 'U) -> source:seq<'T> -> seq<'U>

val zip : source1:seq<'T1> -> source2:seq<'T2> -> seq<'T1 * 'T2>

val sort : source:seq<'T> -> seq<'T> (requires comparison)

val pairwise : source:seq<'T> -> seq<'T * 'T>

val iter : action:('T -> unit) -> source:seq<'T> -> unit

val map2 : mapping:('T1 -> 'T2 -> 'U) -> source1:seq<'T1> -> source2:seq<'T2> -> seq<'U>

val length : array:'T [] -> int

Multiple items

module List



from Microsoft.FSharp.Collections



--------------------

type List<'T> =

| ( [] )

| ( :: ) of Head: 'T * Tail: 'T list

interface IReadOnlyList<'T>

interface IReadOnlyCollection<'T>

interface IEnumerable

interface IEnumerable<'T>

member GetSlice : startIndex:int option * endIndex:int option -> 'T list

member Head : 'T

member IsEmpty : bool

member Item : index:int -> 'T with get

member Length : int

member Tail : 'T list

...

val iter : action:('T -> unit) -> list:'T list -> unit