-
Notifications
You must be signed in to change notification settings - Fork 1
/
exp.go
72 lines (66 loc) · 3.11 KB
/
exp.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
package unsigned
// ----- ---- --- -- -
// Copyright 2019, 2020 The Axiom Foundation. All Rights Reserved.
//
// Licensed under the Apache License 2.0 (the "License"). You may not use
// this file except in compliance with the License. You can obtain a copy
// in the file LICENSE in the source distribution or at
// https://www.apache.org/licenses/LICENSE-2.0.txt
// - -- --- ---- -----
import (
"errors"
"math"
)
// This file contains an implementation of e^x (the exp function) that works for fractions
// between 0 and 1 in a 64-bit fixed point world.
// This frees us from the use of big math and it is also literally 25 times faster than the
// big package and has no memory allocation.
// ExpFrac calculates e^x, where x is a fraction numerator/denominator between
// 0 and 1. We use a Taylor Series expansion of e^x that converges well in the target range.
// This expansion is
// x^0/0! + x^1/1! + x^2/2! ...
// We can collapse the first two terms for convenience to 1+x.
// In addition, we make use of the fact that (numerator/denominator)^2 = numerator^2/denominator^2 so we can use muldiv
// and we require that denominator <= maxint32, and that numerator < denominator.
// Basically, we compute (denominator + numerator + numerator^2/2denominator + numerator^3/6denominator^2 ...) which is denominator times our desired result
// (so that we have the implied denominator).
//
// The return value is the numerator for the fraction; the denominator is unchanged.
// This fixed point calculation tends to produce values that are slightly off in the last digit
// (as compared to a floating point implementation) because of accumulated rounding errors.
// Therefore, what we do is scale the input fraction by multiplying both numerator and denominator by
// a scaling value and then divide by it again at the end.
// This means that the practical limit for denominator is maxint32 / 10, which is still larger than our
// napu multiplication factor of 100,000,000 (which is also the value we use for percentages).
func ExpFrac(numerator, denominator uint64) (uint64, error) {
rounder := uint64(10)
numerator *= rounder
denominator *= rounder
if denominator > (math.MaxUint64 / 2) {
return 0, errors.New("denominator too large")
}
if numerator > denominator {
return 0, errors.New("fraction must be between 0 and 1")
}
// start the sum at 1 + x, which is b/b + a/b, and we only care about the
// numerator, so it's just b+a
sum := denominator + numerator
// we accumulate a product by starting with the original numerator,
// then multiplying it by the fraction using muldiv; we don't square
// the denominator because there's an implied division by b in the result.
// In other words, to square 1200/10000, we muldiv(1200, 1200, 10000) and get 144;
// the result is 144/10000 which is correct.
// This converges at a rate of approximately one decimal digit per loop.
product := numerator
fact := uint64(1)
var err error
for i := uint64(2); product != 0; i++ {
product, err = MulDiv(product, numerator, denominator)
if err != nil {
return 0, err
}
fact *= i
sum += product / fact
}
return (sum + rounder/2) / rounder, nil
}