From 3cf61177898fad4e6be85b52083e9c1b508d6568 Mon Sep 17 00:00:00 2001 From: Adrien Hopkins Date: Tue, 19 Sep 2023 19:44:35 -0500 Subject: factors.Score: Avoid overflow by using math/big When using really large numbers, factors.Score could overflow, which would cause an incorrect result. Using arbitrary-precision arithmetic fixes this. I only do so above 2^28, since below then factor sums are guaranteed to not overflow, and normal arithmetic is faster. --- factors/factors_test.go | 2 ++ factors/score.go | 27 ++++++++++++++++++++++++++- 2 files changed, 28 insertions(+), 1 deletion(-) diff --git a/factors/factors_test.go b/factors/factors_test.go index a18d30b..838186f 100644 --- a/factors/factors_test.go +++ b/factors/factors_test.go @@ -94,6 +94,8 @@ var factorScoreCases = map[uint]float64{ 10: 1.8, 12: 7.0 / 3.0, 120: 3.0, + // number that will use bigScore + 367567200: 62496.0 / 12155.0, } func TestFactorScore(t *testing.T) { diff --git a/factors/score.go b/factors/score.go index 0759a74..c777d11 100644 --- a/factors/score.go +++ b/factors/score.go @@ -1,6 +1,9 @@ package factors -import "math" +import ( + "math" + "math/big" +) // Score returns a "factor score" equal to the sum of the reciprocoals // of the number n's factors. @@ -20,6 +23,8 @@ import "math" func Score(n uint) float64 { if n == 0 { return math.NaN() + } else if n > maxSmall { + return bigScore(n) } factorSum := uint(0) @@ -29,6 +34,26 @@ func Score(n uint) float64 { return float64(factorSum) / float64(n) } +const maxSmall = 1 << 28 - 1 + +func bigScore(n uint) float64 { + factorSum := new(big.Int) + + // initialize bigFactor like this to avoid reallocating memory + bigFactor := new(big.Int).SetUint64(uint64(n)) + + for _, f := range Factors(n) { + bigFactor.SetUint64(uint64(f)) + factorSum.Add(factorSum, bigFactor) + } + + bigN := new(big.Int).SetUint64(uint64(n)) + score := new(big.Rat).SetFrac(factorSum, bigN) + + score64, _ := score.Float64() + return score64 +} + // BasicRank returns a rank describing how well a radix handles the simplest // fractions (1/2, 1/3, 1/4 and 1/5). Zero and one are not true radices, // but because this rank otherwise only depends on a radix's remainder -- cgit v1.2.3