# 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.

## Error minimisation algorithm

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 distribute n weights =
let wt = Array.sum weights
if Array.isEmpty weights
|| Seq.maxBy abs weights * float(abs n+1) |> abs
>= abs wt * float Int32.MaxValue then None
else
let inline round f = int(if f<0.0 then f-0.5 else f+0.5)
let ns = Array.map ((*)(float n / wt) >> round) weights

let inline absoluteError wi ni d =
let wc = 0.5 * float d * wt / float n
let wni = float ni * wt / float n
if wc > 0.0 then min (wni-wi) 0.0 + wc
else min (wi-wni) 0.0 - wc

let d = n - Array.sum ns
for __ = 1 to abs d do
let _,_,_,i =
Seq.mapi2 (fun i wi ni ->
let absErr = absoluteError wi ni (sign d)
let relErr = absErr / abs wi
let weight = if wt > 0.0 then -wi else wi
absErr, relErr, weight, i
) weights ns |> Seq.min
ns.[i] <- ns.[i] + sign d
Some ns


## Testing

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"
}

let r = distribute 100 [|406.0;348.0;246.0;0.0|]
Expect.equal r (Some [|40;35;25;0|]) "40 etc"
}

let r = distribute -100 [|406.0;348.0;246.0;0.0|]
Expect.equal r (Some [|-40;-35;-25;0|]) "-40 etc"
}

let r = distribute 100 [|-406.0;-348.0;-246.0;-0.0|]
Expect.equal r (Some [|40;35;25;0|]) "40 etc"
}

let r = distribute -100 [|-406.0;-348.0;-246.0;-0.0|]
Expect.equal r (Some [|-40;-35;-25;0|]) "-40 etc"
}

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"
)
)

]


## Conclusion

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