Source file
src/runtime/mkfastlog2table.go
1
2
3
4
5
6
7
8
9
10 package main
11
12 import (
13 "bytes"
14 "fmt"
15 "log"
16 "math"
17 "os"
18 )
19
20 func main() {
21 var buf bytes.Buffer
22
23 fmt.Fprintln(&buf, "// Code generated by mkfastlog2table.go; DO NOT EDIT.")
24 fmt.Fprintln(&buf, "// Run go generate from src/runtime to update.")
25 fmt.Fprintln(&buf, "// See mkfastlog2table.go for comments.")
26 fmt.Fprintln(&buf)
27 fmt.Fprintln(&buf, "package runtime")
28 fmt.Fprintln(&buf)
29 fmt.Fprintln(&buf, "const fastlogNumBits =", fastlogNumBits)
30 fmt.Fprintln(&buf)
31
32 fmt.Fprintln(&buf, "var fastlog2Table = [1<<fastlogNumBits + 1]float64{")
33 table := computeTable()
34 for _, t := range table {
35 fmt.Fprintf(&buf, "\t%v,\n", t)
36 }
37 fmt.Fprintln(&buf, "}")
38
39 if err := os.WriteFile("fastlog2table.go", buf.Bytes(), 0644); err != nil {
40 log.Fatalln(err)
41 }
42 }
43
44 const fastlogNumBits = 5
45
46 func computeTable() []float64 {
47 fastlog2Table := make([]float64, 1<<fastlogNumBits+1)
48 for i := 0; i <= (1 << fastlogNumBits); i++ {
49 fastlog2Table[i] = log2(1.0 + float64(i)/(1<<fastlogNumBits))
50 }
51 return fastlog2Table
52 }
53
54
55
56 func log2(x float64) float64 {
57 frac, exp := math.Frexp(x)
58
59
60 if frac == 0.5 {
61 return float64(exp - 1)
62 }
63 return float64(nlog(frac)*(1/math.Ln2)) + float64(exp)
64 }
65
66
67
68 func nlog(x float64) float64 {
69 const (
70 Ln2Hi = 6.93147180369123816490e-01
71 Ln2Lo = 1.90821492927058770002e-10
72 L1 = 6.666666666666735130e-01
73 L2 = 3.999999999940941908e-01
74 L3 = 2.857142874366239149e-01
75 L4 = 2.222219843214978396e-01
76 L5 = 1.818357216161805012e-01
77 L6 = 1.531383769920937332e-01
78 L7 = 1.479819860511658591e-01
79 )
80
81
82 switch {
83 case math.IsNaN(x) || math.IsInf(x, 1):
84 return x
85 case x < 0:
86 return math.NaN()
87 case x == 0:
88 return math.Inf(-1)
89 }
90
91
92 f1, ki := math.Frexp(x)
93 if f1 < math.Sqrt2/2 {
94 f1 *= 2
95 ki--
96 }
97 f := f1 - 1
98 k := float64(ki)
99
100
101 s := float64(f / (2 + f))
102 s2 := float64(s * s)
103 s4 := float64(s2 * s2)
104 t1 := s2 * float64(L1+float64(s4*float64(L3+float64(s4*float64(L5+float64(s4*L7))))))
105 t2 := s4 * float64(L2+float64(s4*float64(L4+float64(s4*L6))))
106 R := float64(t1 + t2)
107 hfsq := float64(0.5 * f * f)
108 return float64(k*Ln2Hi) - ((hfsq - (float64(s*float64(hfsq+R)) + float64(k*Ln2Lo))) - f)
109 }
110
View as plain text