github.com/cznic/mathutil@v0.0.0-20181122101859-297441e03548/binarylog.go (about)

     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
     6  
     7  import (
     8  	"math/big"
     9  
    10  	"github.com/remyoudompheng/bigfft"
    11  )
    12  
    13  type float struct {
    14  	n           *big.Int
    15  	fracBits    int
    16  	maxFracBits int
    17  }
    18  
    19  func newFloat(n *big.Int, fracBits, maxFracBits int) float {
    20  	f := float{n: n, fracBits: fracBits, maxFracBits: maxFracBits}
    21  	f.normalize()
    22  	return f
    23  }
    24  
    25  func (f *float) normalize() {
    26  	n := f.n.BitLen()
    27  	if n == 0 {
    28  		return
    29  	}
    30  
    31  	if n := f.fracBits - f.maxFracBits; n > 0 {
    32  		bit := f.n.Bit(n - 1)
    33  		f.n.Rsh(f.n, uint(n))
    34  		if bit != 0 {
    35  			f.n.Add(f.n, _1)
    36  		}
    37  		f.fracBits -= n
    38  	}
    39  
    40  	var i int
    41  	for ; f.fracBits > 0 && i <= f.fracBits && f.n.Bit(i) == 0; i++ {
    42  		f.fracBits--
    43  	}
    44  
    45  	if i != 0 {
    46  		f.n.Rsh(f.n, uint(i))
    47  	}
    48  }
    49  
    50  func (f *float) eq1() bool { return f.fracBits == 0 && f.n.BitLen() == 1 }
    51  func (f *float) ge2() bool { return f.n.BitLen() > f.fracBits+1 }
    52  
    53  func (f *float) div2() {
    54  	f.fracBits++
    55  	f.normalize()
    56  }
    57  
    58  func (f *float) sqr() {
    59  	f.n = bigfft.Mul(f.n, f.n)
    60  	f.fracBits *= 2
    61  	f.normalize()
    62  }
    63  
    64  // BinaryLog computes the binary logarithm of n. The result consists of a
    65  // characteristic and a mantissa having precision mantissaBits. The value of
    66  // the binary logarithm is
    67  //
    68  //	characteristic + mantissa*(2^-mantissaBits)
    69  //
    70  // BinaryLog panics for n <= 0 or mantissaBits < 0.
    71  func BinaryLog(n *big.Int, mantissaBits int) (characteristic int, mantissa *big.Int) {
    72  	if n.Sign() <= 0 || mantissaBits < 0 {
    73  		panic("invalid argument of BinaryLog")
    74  	}
    75  
    76  	characteristic = n.BitLen() - 1
    77  	mantissa = big.NewInt(0)
    78  	x := newFloat(n, characteristic, mantissaBits)
    79  	for ; mantissaBits != 0 && !x.eq1(); mantissaBits-- {
    80  		x.sqr()
    81  		mantissa.Lsh(mantissa, 1)
    82  		if x.ge2() {
    83  			mantissa.SetBit(mantissa, 0, 1)
    84  			x.div2()
    85  		}
    86  	}
    87  	return characteristic, mantissa
    88  }