bbgo/pkg/fixedpoint/div128.go

135 lines
3.2 KiB
Go

//go:build dnum
// +build dnum
// Copyright Suneido Software Corp. All rights reserved.
// Governed by the MIT license found in the LICENSE file.
package fixedpoint
import (
"math/bits"
)
const (
e16 = 1_0000_0000_0000_0000
longMask = 0xffffffff
divNumBase = 1 << 32
e16Hi = e16 >> 32
e16Lo = e16 & longMask
)
// returns (1e16 * dividend) / divisor
// Used by dnum divide
// Based on cSuneido code
// which is based on jSuneido code
// which is based on Java BigDecimal code
// which is based on Hacker's Delight and Knuth TAoCP Vol 2
// A bit simpler with unsigned types
func div128(dividend, divisor uint64) uint64 {
//check(dividend != 0)
//check(divisor != 0)
// multiply dividend * e16
d1Hi := dividend >> 32
d1Lo := dividend & longMask
product := uint64(e16Lo) * d1Lo
d0 := product & longMask
d1 := product >> 32
product = uint64(e16Hi)*d1Lo + d1
d1 = product & longMask
d2 := product >> 32
product = uint64(e16Lo)*d1Hi + d1
d1 = product & longMask
d2 += product >> 32
d3 := d2 >> 32
d2 &= longMask
product = e16Hi*d1Hi + d2
d2 = product & longMask
d3 = ((product >> 32) + d3) & longMask
dividendHi := make64(uint32(d3), uint32(d2))
dividendLo := make64(uint32(d1), uint32(d0))
// divide
return divide128(dividendHi, dividendLo, divisor)
}
func divide128(dividendHi, dividendLo, divisor uint64) uint64 {
// so we can shift dividend as much as divisor
// don't allow equals to avoid quotient overflow (by 1)
//check(dividendHi < divisor)
// maximize divisor (bit wise), since we're mostly using the top half
shift := uint(bits.LeadingZeros64(divisor))
divisor = divisor << shift
// split divisor
v1 := divisor >> 32
v0 := divisor & longMask
// matching shift
dls := dividendLo << shift
// split dividendLo
u1 := uint32(dls >> 32)
u0 := uint32(dls & longMask)
// tmp1 = top 64 of dividend << shift
tmp1 := (dividendHi << shift) | (dividendLo >> (64 - shift))
var q1, rtmp1 uint64
if v1 == 1 {
q1 = tmp1
rtmp1 = 0
} else {
//check(tmp1 >= 0)
q1 = tmp1 / v1 // DIVIDE top 64 / top 32
rtmp1 = tmp1 % v1 // remainder
}
// adjust if quotient estimate too large
//check(q1 < divNumBase)
for q1*v0 > make64(uint32(rtmp1), u1) {
// done about 5.5 per 10,000 divides
q1--
rtmp1 += v1
if rtmp1 >= divNumBase {
break
}
}
//check(q1 >= 0)
u2 := tmp1 & longMask // low half
// u2,u1 is the MIDDLE 64 bits of the dividend
tmp2 := mulsub(uint32(u2), uint32(u1), uint32(v1), uint32(v0), q1)
var q0, rtmp2 uint64
if v1 == 1 {
q0 = tmp2
rtmp2 = 0
} else {
q0 = tmp2 / v1 // DIVIDE dividend remainder 64 / divisor high 32
rtmp2 = tmp2 % v1
}
// adjust if quotient estimate too large
//check(q0 < divNumBase)
for q0*v0 > make64(uint32(rtmp2), u0) {
// done about .33 times per divide
q0--
rtmp2 += v1
if rtmp2 >= divNumBase {
break
}
//check(q0 < divNumBase)
}
//check(q1 <= math.MaxUint32)
//check(q0 <= math.MaxUint32)
return make64(uint32(q1), uint32(q0))
}
// mulsub returns u1,u0 - v1,v0 * q0
func mulsub(u1, u0, v1, v0 uint32, q0 uint64) uint64 {
tmp := uint64(u0) - q0*uint64(v0)
return make64(u1+uint32(tmp>>32)-uint32(q0*uint64(v1)), uint32(tmp&longMask))
}
func make64(hi, lo uint32) uint64 {
return uint64(hi)<<32 | uint64(lo)
}