135 lines
3.2 KiB
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)
|
|
}
|