centflip

Find the exact inputs where f64 interest, fee or FX code rounds to the wrong cent, with a proof that the list is complete.

The problem

Money code written as amount * rate * days / 360 in f64 and then rounded to cents passes its unit tests, yet books the wrong cent on inputs whose exact value is a half-cent (or within a few ulps of one) and whose f64 result lands on the other side of it. Which inputs those are depends on the rate, the day count, the order of operations and whether the compiler fused a*b + c into one FMA instruction, so neither code review nor random testing finds them reliably. The usual responses are to rewrite everything with a decimal type without knowing whether the float code was wrong, or to ship it. centflip takes the formula and the input domain and lists every input that rounds wrongly, or proves that none does.

How it works

Each formula is parsed into a tree whose shape is the f64 evaluation order (P*r*d/360 is ((P*r)*d)/360). Every input can be evaluated three ways:

• exactly, in rationals on num-bigint, with decimal literals taken at their decimal value;
• in f64 without contraction: each literal is its nearest double and each operation rounds in tree order;
• in f64 with FMA contraction: a*b + c, c + a*b, a*b - c, c - a*b become one fused multiply-add.

An input is mis-rounded when the f64 result rounds to a different cent from the exact result. The double is rounded from its exact binary value, so no second float error comes in at that step. Half-up (ties away from zero, as in f64::round) and half-even are both checked.

The search never guesses. It relies on one fact: the f64 cent can differ from the exact cent only if a half-cent boundary lies between the two values. So:

1. Error bound. An interval pass over the tree (outward-rounded, using next_up/next_down) gives, for a whole block of inputs, an enclosure of the exact values and a bound E on |f64 - exact| that holds for both FMA settings. It uses the standard model |fl(z) - z| <= 2^-53 |z| per operation and carries the propagated error through + - * / and fused multiply-add. Blocks are split at powers of two in the scanned variable, so E scales with each block's own magnitude.

2. Affine formulas: continued-fraction search. With every variable except one integer x fixed, P*r*d/360, fees, FX chains and P*(1+r/12)^n all have an exact value of the form a*x + b with rational a, b. An input needs checking only if 100(a*x + b) lies within 100E of k + 1/2, which is a Diophantine closeness condition on frac(alpha*x + beta).

   • alpha and beta are replaced by their last continued-fraction convergents p/q (denominator at most 2^36) and r/s (at most 2^24). The approximation error, computed exactly, is added to the tolerance. For decimal rates the convergent is the rate itself, so this step is exact. For compounded rates with denominators like 240^60 it costs about 2^-45 of a cent.
   • The condition becomes (A*x + B) mod M in [-t, t] with M = q*s <= 2^60, all in i128.
   • The smallest solution of L <= (A*j) mod M <= R comes from a Euclid-style recursion (A, M) -> (M mod A, A). It follows the continued-fraction quotients of A/M, so each hit costs O(log M).
   • Hits repeat with period M / gcd(A, M). One period is solved and then replicated along the range.

   The result is the complete list of candidates. Every other input is at least E from a half-cent, so its f64 value rounds to the same cent. Each candidate is then checked exactly.

3. Other formulas: interval pruning. The range is bisected. Any sub-range whose enclosure, widened by E, contains no half-cent is discarded with a proof. Sub-ranges of 16 inputs or fewer are checked one input at a time.

4. Precision limit. At |exact| >= 2^46 (about 7.04e13) consecutive doubles are more than a cent apart. Those inputs are reported as a warning region and are not analysed.

Worked example: P*0.07/100 at P = 250. The exact value is 0.175, a half-cent. 0.07 is stored as 0.07000000000000000666…, and 250 times that is 17.5000000000000016…, less than half an ulp above 17.5, so the product rounds to exactly 17.5. Then 17.5/100 rounds to the nearest double to 0.175, which is 0.17499999999999998889…. Half-up books 0.17 where the exact answer is 0.18. The continued-fraction step reaches this input without visiting the others: the slope 100 * 0.0007 = 7/100 gives M = 200, and the hits are exactly the P ≡ 50 (mod 100).

Install and usage

Requires Rust 1.86 or newer.

git clone <this repository> centflip && cd centflip
cargo build --release
cargo test

Scan simple interest on 3.66 billion inputs for half-even errors:

$ cargo run --release -q -- scan 'P*r*d/360' --P 1..10000000 --r 0.0525 --d 1..366 --mode half-even --limit 6
formula   P*r*d/360
domain    P = 1..=10000000, r = 0.0525, d = 1..=366
rounding  half-even; FMA contraction off and on
search    continued-fraction search along P
inputs    3660000000 in domain, 58354150 checked exactly, 4.27 s

27608437 mis-rounded results (half-up: 0, half-even: 27608437; evaluation-order: 27608437, fma-contraction: 0)

P     d  mode       exact value  f64 value                  exact  f64   diff  cause
240   1  half-even  0.035        0.03499999999999999639...  0.04   0.03  -1    evaluation-order
1200  1  half-even  0.175        0.17499999999999998889...  0.18   0.17  -1    evaluation-order
2640  1  half-even  0.385        0.38500000000000000888...  0.38   0.39  +1    evaluation-order
3120  1  half-even  0.455        0.45499999999999996003...  0.46   0.45  -1    evaluation-order
3600  1  half-even  0.525        0.52500000000000002220...  0.52   0.53  +1    evaluation-order
4080  1  half-even  0.595        0.59499999999999997335...  0.60   0.59  -1    evaluation-order
... 27608431 more (raise --limit or use --format csv)

f64 value is the exact decimal expansion of the double, cut after 20 digits. --format csv prints every row in full.

Separate the errors that FMA contraction introduces:

$ cargo run --release -q -- scan 'P*r + fee' --P 1..1000000 --r 0.0725 --fee 0.35 --mode half-up --fma on --limit 4
formula   P*r + fee
domain    P = 1..=1000000, r = 0.0725, fee = 0.35
rounding  half-up; FMA contraction on
search    continued-fraction search along P
inputs    1000000 in domain, 250000 checked exactly, 0.04 s

229334 mis-rounded results (half-up: 229334, half-even: 0; evaluation-order: 208384, fma-contraction: 20950)

P   mode     exact value  f64 value                  exact  f64   diff  cause
2   half-up  0.495        0.49499999999999999555...  0.50   0.49  -1    evaluation-order
6   half-up  0.785        0.78499999999999992006...  0.79   0.78  -1    evaluation-order
10  half-up  1.075        1.07499999999999995559...  1.08   1.07  -1    evaluation-order
14  half-up  1.365        1.36499999999999976907...  1.37   1.36  -1    evaluation-order
... 229330 more (raise --limit or use --format csv)

With --fma on, evaluation-order means the input is also wrong without contraction, and fma-contraction means only the fused evaluation is wrong.

A formula with no mis-rounded input. The exact cents are P*d/73, which is never a half-cent, and the search checks no input at all:

$ cargo run --release -q -- scan 'P*r*d/365' --P 1..100000000 --r 0.05 --d 1..366
formula   P*r*d/365
domain    P = 1..=100000000, r = 0.05, d = 1..=366
rounding  half-up, half-even; FMA contraction off and on
search    continued-fraction search along P
inputs    36600000000 in domain, 0 checked exactly, 0.06 s

No mis-rounded input. Each input was either checked exactly or shown to lie farther from a half-cent boundary than the proven f64 error bound, so its f64 result rounds to the same cent as the exact result.

Formula syntax: numbers, variables, + - * / ^, parentheses. ^ takes an integer exponent (constant, or a variable bound to a range that is iterated in the outer loop) and is evaluated as repeated multiplication. Every variable gets a flag: --P 1..10000000 is an inclusive integer range, and --r 0.0525 or --d 30 is a constant. Options: --mode half-up|half-even|both, --fma off|on|both, --strategy auto|interval|brute, --limit N, --format text|csv. The exit status is 0 when there is no mis-rounded input, 1 when some were found, and 2 on usage or evaluation errors.

The library API is the same: Problem::new(formula, bindings), then scan(&problem, &options, |counterexample| ...).

Results

Reproduce with cargo run --release --example bench, which takes about five minutes. -- --quick uses domains 100 times smaller.

Measured on an Apple Silicon (arm64) Mac under macOS, single-threaded, release build:

case	domain (inputs)	counterexamples	search time	inputs checked exactly	brute force (est.)	speedup	fuzzer samples in same time	fuzzer recall
P*r*d/360, r=0.0525, P 1..1e8, d 1..366	3.66e10	649,554,505	56.7 s	5.84e8	2394 s	42x	8.10e8	2.19%
P*(1+r/12)^n, r=0.05, P 1..1e7, n 1..60	6.00e8	749,241	0.16 s	8.57e5	2650 s	16,937x	3.69e4	0.00%

How the columns were measured:

• Search is scan with the default strategy. Counterexamples are counted per rounding mode, and one input can appear twice, once for half-up and once for half-even.
• Brute force runs Strategy::Brute, which applies the same per-input check to every input. Running it on the full domains would take 40 minutes or more, so it runs on a slice: P up to 1e6 for the first row and 1e5 for the second, with the outer variable over its full range. On that slice it must produce the same counterexample set as the search; the benchmark compares an order-independent hash and the count, and exits non-zero if they differ. Both slices matched. The slice time is then scaled by domain size.
• Fuzzer draws inputs uniformly at random for the same wall-clock time the search took, using the same per-input check. Recall is the number of distinct counterexamples it found divided by the total.

The two rows show the two regimes:

• Decimal rates. The exact cents are P*d*7/480, so exact half-cents recur periodically. The search had to check 1.6% of inputs, and it reports 649.6 million mode-labelled counterexamples, about 1.8% of the input count. The output is large, and the search's cost tracks its size. The speedup is 42x, and a same-budget fuzzer sees 2.2% of the set.
• Compounded rates. Mis-roundings are rare (0.12%), but each exact check needs big-integer arithmetic (about 4 µs per input). The search checks fewer than a million inputs out of 600 million. In the same 0.16 s the fuzzer checked 37,000 random inputs and found none of the 749,241.

Correctness evidence, all in cargo test:

• Exhaustive comparison on the headline formula. For P*0.0525*d/360 over every P up to 1e5 and every d from 1 to 366 (36.6 million inputs), the search must match an oracle that uses only u128 arithmetic and the f64 expression written directly in Rust.
• Brute-force comparison on more formulas. Fees with refunds, FX chains, compound interest, contractible formulas, nonlinear formulas, exact ties and multi-variable formulas are each compared against brute force under all three strategies.
• Property tests. Randomized affine formulas and domains check four claims: the search equals brute force; every pruned interval, enumerated, has no mis-rounding; the candidate generator covers every close input; and the enclosure bounds every f64 error.
• Independent re-check of every reported counterexample. Each one is re-checked with big-integer rounding from the double's bits.
• FMA attribution. Inputs labelled fma-contraction must round correctly without FMA, and inputs labelled evaluation-order must round wrongly without it. --fma off and --fma on must report exactly the inputs that are wrong under that evaluation.
• Hand-derived cases and edge cases. 1.005, 0.1 + 0.2 + 0.005, and P*0.07/100 at P=250, all derived in the test comments; negative amounts; exact half-cents; amounts past the 2^46 precision limit; and error reporting.

Design notes

Filter with a proven bound, then decide exactly. The core of the project is the tolerance E, not the arithmetic. With a sound bound on |f64 - exact|, finding candidates becomes a pure number-theory question with an exact answer, and each candidate is then settled by exact arithmetic. The bound can be loose: a loose bound only adds candidates, and each candidate is checked exactly. So the analysis uses simple running-error rules plus relative inflation instead of a tight but fragile model. The completeness tests hold the whole pipeline to account: search and brute force must produce the same set, and they are compared on 36.6 million inputs of the headline formula and on randomized affine formulas. The cost is that the search does work proportional to its output. When a formula mis-rounds a large share of its domain, it cannot be much faster than that share.

Approximate the slope, not the arithmetic. Compounded rates have exact slopes with thousands of digits. Running the modular recursion in big integers would be simple but slow, so the slope and offset are replaced by bounded continued-fraction convergents and the exact approximation error goes into the tolerance. Decimal rates lose nothing, since the convergent equals the rate. Huge denominators cost a tolerance of about 2^-45 cent. All arithmetic stays in i128, and the proof stays intact because the error is added, not ignored.

The tree is the evaluation order. Rust, C and Java compilers do not reorder floating-point operations under default flags, so the parser keeps the written order. ^ is expanded into repeated multiplication instead of calling powi, whose results are not specified across platforms. The only transformation modelled is FMA contraction, and it is modelled explicitly. That is why a cause can be attributed: fma-contraction means the non-fused evaluation of the same tree rounds correctly.

Exact rounding of the double. The f64 result is rounded to cents from its exact binary value (mantissa times a power of two, rounded in integers). A second float step such as (x * 100.0).round() would add error of its own. Mixing it in would blur which formula a counterexample belongs to.

The precision limit belongs to each input. Inputs are "beyond f64 precision" when their exact result is at least 2^46. Because this is a property of the exact value, brute force and the search agree on exactly which inputs are excluded, and the tests compare those counts too.

Limitations

• Variables are integers. Amounts with cents must be written as integer cents with an explicit /100, and that division is part of the evaluation order being tested. Rates and other decimal constants are fixed per scan, not ranged.
• Only one variable gets the smart search. The widest range that does not appear in an exponent is searched. Every other ranged variable is iterated in full, so --d 1..366 costs 366 searches.
• The float model is IEEE binary64 with round-to-nearest. x87 80-bit intermediates, -ffast-math reassociation and contraction across statements are not modelled. FMA contraction applies to a product that is a direct operand of + or -, with the left product taking precedence, which is one common choice among several.
• The rounding step is modelled as exact. Code that rounds with (x*100.0).round()/100.0, or rounds an intermediate result such as a daily rate, needs that step written into the formula, and there is no syntax for rounding inside a formula.
• No pow with fractional exponents, exp, ln, min, max or floor.
• Interval pruning is weak for steep formulas. If the value moves more than about a cent per unit step, most sub-ranges contain a half-cent and the search checks nearly every input. It is still correct, but no faster than brute force.
• Decimal rates produce dense output. For P*0.0525*d/360, exact half-cents recur periodically, and the f64 result often lands on the wrong side of them. About 1.6% of inputs need an exact check. The search spends its time checking and reporting those, and the speedup comes from skipping the rest.
• The search is single-threaded. It is not parallelized, and neither are the baselines.
• Near the precision limit, from about 2^40 to 2^46, results are analysed and reported, but doubles are only a fraction of a cent apart there, so almost every half-cent input is a counterexample.

License

MIT. See LICENSE.