1 | // Copyright (c) 2016 The mathutil Authors. All rights reserved.
|
---|
2 | // Use of this source code is governed by a BSD-style
|
---|
3 | // license that can be found in the LICENSE file.
|
---|
4 |
|
---|
5 | package mathutil // import "modernc.org/mathutil"
|
---|
6 |
|
---|
7 | import (
|
---|
8 | "math/big"
|
---|
9 | )
|
---|
10 |
|
---|
11 | type float struct {
|
---|
12 | n *big.Int
|
---|
13 | fracBits int
|
---|
14 | maxFracBits int
|
---|
15 | }
|
---|
16 |
|
---|
17 | func newFloat(n *big.Int, fracBits, maxFracBits int) float {
|
---|
18 | f := float{n: n, fracBits: fracBits, maxFracBits: maxFracBits}
|
---|
19 | f.normalize()
|
---|
20 | return f
|
---|
21 | }
|
---|
22 |
|
---|
23 | func (f *float) normalize() {
|
---|
24 | n := f.n.BitLen()
|
---|
25 | if n == 0 {
|
---|
26 | return
|
---|
27 | }
|
---|
28 |
|
---|
29 | if n := f.fracBits - f.maxFracBits; n > 0 {
|
---|
30 | bit := f.n.Bit(n - 1)
|
---|
31 | f.n.Rsh(f.n, uint(n))
|
---|
32 | if bit != 0 {
|
---|
33 | f.n.Add(f.n, _1)
|
---|
34 | }
|
---|
35 | f.fracBits -= n
|
---|
36 | }
|
---|
37 |
|
---|
38 | var i int
|
---|
39 | for ; f.fracBits > 0 && i <= f.fracBits && f.n.Bit(i) == 0; i++ {
|
---|
40 | f.fracBits--
|
---|
41 | }
|
---|
42 |
|
---|
43 | if i != 0 {
|
---|
44 | f.n.Rsh(f.n, uint(i))
|
---|
45 | }
|
---|
46 | }
|
---|
47 |
|
---|
48 | func (f *float) eq1() bool { return f.fracBits == 0 && f.n.BitLen() == 1 }
|
---|
49 | func (f *float) ge2() bool { return f.n.BitLen() > f.fracBits+1 }
|
---|
50 |
|
---|
51 | func (f *float) div2() {
|
---|
52 | f.fracBits++
|
---|
53 | f.normalize()
|
---|
54 | }
|
---|
55 |
|
---|
56 | // BinaryLog computes the binary logarithm of n. The result consists of a
|
---|
57 | // characteristic and a mantissa having precision mantissaBits. The value of
|
---|
58 | // the binary logarithm is
|
---|
59 | //
|
---|
60 | // characteristic + mantissa*(2^-mantissaBits)
|
---|
61 | //
|
---|
62 | // BinaryLog panics for n <= 0 or mantissaBits < 0.
|
---|
63 | func BinaryLog(n *big.Int, mantissaBits int) (characteristic int, mantissa *big.Int) {
|
---|
64 | if n.Sign() <= 0 || mantissaBits < 0 {
|
---|
65 | panic("invalid argument of BinaryLog")
|
---|
66 | }
|
---|
67 |
|
---|
68 | characteristic = n.BitLen() - 1
|
---|
69 | mantissa = big.NewInt(0)
|
---|
70 | x := newFloat(n, characteristic, mantissaBits)
|
---|
71 | for ; mantissaBits != 0 && !x.eq1(); mantissaBits-- {
|
---|
72 | x.sqr()
|
---|
73 | mantissa.Lsh(mantissa, 1)
|
---|
74 | if x.ge2() {
|
---|
75 | mantissa.SetBit(mantissa, 0, 1)
|
---|
76 | x.div2()
|
---|
77 | }
|
---|
78 | }
|
---|
79 | return characteristic, mantissa
|
---|
80 | }
|
---|