ai_old/goja/ftoa/ftoa.go
2024-09-20 16:50:35 +08:00

700 lines
15 KiB
Go

package ftoa
import (
"math"
"math/big"
)
const (
exp_11 = 0x3ff00000
frac_mask1 = 0xfffff
bletch = 0x10
quick_max = 14
int_max = 14
)
var (
tens = [...]float64{
1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1e20, 1e21, 1e22,
}
bigtens = [...]float64{1e16, 1e32, 1e64, 1e128, 1e256}
big5 = big.NewInt(5)
big10 = big.NewInt(10)
p05 = []*big.Int{big5, big.NewInt(25), big.NewInt(125)}
pow5Cache [7]*big.Int
dtoaModes = []int{
ModeStandard: 0,
ModeStandardExponential: 0,
ModeFixed: 3,
ModeExponential: 2,
ModePrecision: 2,
}
)
/*
d must be > 0 and must not be Inf
mode:
0 ==> shortest string that yields d when read in
and rounded to nearest.
1 ==> like 0, but with Steele & White stopping rule;
e.g. with IEEE P754 arithmetic , mode 0 gives
1e23 whereas mode 1 gives 9.999999999999999e22.
2 ==> max(1,ndigits) significant digits. This gives a
return value similar to that of ecvt, except
that trailing zeros are suppressed.
3 ==> through ndigits past the decimal point. This
gives a return value similar to that from fcvt,
except that trailing zeros are suppressed, and
ndigits can be negative.
4,5 ==> similar to 2 and 3, respectively, but (in
round-nearest mode) with the tests of mode 0 to
possibly return a shorter string that rounds to d.
With IEEE arithmetic and compilation with
-DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
as modes 2 and 3 when FLT_ROUNDS != 1.
6-9 ==> Debugging modes similar to mode - 4: don't try
fast floating-point estimate (if applicable).
Values of mode other than 0-9 are treated as mode 0.
*/
func ftoa(d float64, mode int, biasUp bool, ndigits int, buf []byte) ([]byte, int) {
startPos := len(buf)
dblBits := make([]byte, 0, 8)
be, bbits, dblBits := d2b(d, dblBits)
dBits := math.Float64bits(d)
word0 := uint32(dBits >> 32)
word1 := uint32(dBits)
i := int((word0 >> exp_shift1) & (exp_mask >> exp_shift1))
var d2 float64
var denorm bool
if i != 0 {
d2 = setWord0(d, (word0&frac_mask1)|exp_11)
i -= bias
denorm = false
} else {
/* d is denormalized */
i = bbits + be + (bias + (p - 1) - 1)
var x uint64
if i > 32 {
x = uint64(word0)<<(64-i) | uint64(word1)>>(i-32)
} else {
x = uint64(word1) << (32 - i)
}
d2 = setWord0(float64(x), uint32((x>>32)-31*exp_mask))
i -= (bias + (p - 1) - 1) + 1
denorm = true
}
/* At this point d = f*2^i, where 1 <= f < 2. d2 is an approximation of f. */
ds := (d2-1.5)*0.289529654602168 + 0.1760912590558 + float64(i)*0.301029995663981
k := int(ds)
if ds < 0.0 && ds != float64(k) {
k-- /* want k = floor(ds) */
}
k_check := true
if k >= 0 && k < len(tens) {
if d < tens[k] {
k--
}
k_check = false
}
/* At this point floor(log10(d)) <= k <= floor(log10(d))+1.
If k_check is zero, we're guaranteed that k = floor(log10(d)). */
j := bbits - i - 1
var b2, s2, b5, s5 int
/* At this point d = b/2^j, where b is an odd integer. */
if j >= 0 {
b2 = 0
s2 = j
} else {
b2 = -j
s2 = 0
}
if k >= 0 {
b5 = 0
s5 = k
s2 += k
} else {
b2 -= k
b5 = -k
s5 = 0
}
/* At this point d/10^k = (b * 2^b2 * 5^b5) / (2^s2 * 5^s5), where b is an odd integer,
b2 >= 0, b5 >= 0, s2 >= 0, and s5 >= 0. */
if mode < 0 || mode > 9 {
mode = 0
}
try_quick := true
if mode > 5 {
mode -= 4
try_quick = false
}
leftright := true
var ilim, ilim1 int
switch mode {
case 0, 1:
ilim, ilim1 = -1, -1
ndigits = 0
case 2:
leftright = false
fallthrough
case 4:
if ndigits <= 0 {
ndigits = 1
}
ilim, ilim1 = ndigits, ndigits
case 3:
leftright = false
fallthrough
case 5:
i = ndigits + k + 1
ilim = i
ilim1 = i - 1
}
/* ilim is the maximum number of significant digits we want, based on k and ndigits. */
/* ilim1 is the maximum number of significant digits we want, based on k and ndigits,
when it turns out that k was computed too high by one. */
fast_failed := false
if ilim >= 0 && ilim <= quick_max && try_quick {
/* Try to get by with floating-point arithmetic. */
i = 0
d2 = d
k0 := k
ilim0 := ilim
ieps := 2 /* conservative */
/* Divide d by 10^k, keeping track of the roundoff error and avoiding overflows. */
if k > 0 {
ds = tens[k&0xf]
j = k >> 4
if (j & bletch) != 0 {
/* prevent overflows */
j &= bletch - 1
d /= bigtens[len(bigtens)-1]
ieps++
}
for ; j != 0; i++ {
if (j & 1) != 0 {
ieps++
ds *= bigtens[i]
}
j >>= 1
}
d /= ds
} else if j1 := -k; j1 != 0 {
d *= tens[j1&0xf]
for j = j1 >> 4; j != 0; i++ {
if (j & 1) != 0 {
ieps++
d *= bigtens[i]
}
j >>= 1
}
}
/* Check that k was computed correctly. */
if k_check && d < 1.0 && ilim > 0 {
if ilim1 <= 0 {
fast_failed = true
} else {
ilim = ilim1
k--
d *= 10.
ieps++
}
}
/* eps bounds the cumulative error. */
eps := float64(ieps)*d + 7.0
eps = setWord0(eps, _word0(eps)-(p-1)*exp_msk1)
if ilim == 0 {
d -= 5.0
if d > eps {
buf = append(buf, '1')
k++
return buf, k + 1
}
if d < -eps {
buf = append(buf, '0')
return buf, 1
}
fast_failed = true
}
if !fast_failed {
fast_failed = true
if leftright {
/* Use Steele & White method of only
* generating digits needed.
*/
eps = 0.5/tens[ilim-1] - eps
for i = 0; ; {
l := int64(d)
d -= float64(l)
buf = append(buf, byte('0'+l))
if d < eps {
return buf, k + 1
}
if 1.0-d < eps {
buf, k = bumpUp(buf, k)
return buf, k + 1
}
i++
if i >= ilim {
break
}
eps *= 10.0
d *= 10.0
}
} else {
/* Generate ilim digits, then fix them up. */
eps *= tens[ilim-1]
for i = 1; ; i++ {
l := int64(d)
d -= float64(l)
buf = append(buf, byte('0'+l))
if i == ilim {
if d > 0.5+eps {
buf, k = bumpUp(buf, k)
return buf, k + 1
} else if d < 0.5-eps {
buf = stripTrailingZeroes(buf, startPos)
return buf, k + 1
}
break
}
d *= 10.0
}
}
}
if fast_failed {
buf = buf[:startPos]
d = d2
k = k0
ilim = ilim0
}
}
/* Do we have a "small" integer? */
if be >= 0 && k <= int_max {
/* Yes. */
ds = tens[k]
if ndigits < 0 && ilim <= 0 {
if ilim < 0 || d < 5*ds || (!biasUp && d == 5*ds) {
buf = buf[:startPos]
buf = append(buf, '0')
return buf, 1
}
buf = append(buf, '1')
k++
return buf, k + 1
}
for i = 1; ; i++ {
l := int64(d / ds)
d -= float64(l) * ds
buf = append(buf, byte('0'+l))
if i == ilim {
d += d
if (d > ds) || (d == ds && (((l & 1) != 0) || biasUp)) {
buf, k = bumpUp(buf, k)
}
break
}
d *= 10.0
if d == 0 {
break
}
}
return buf, k + 1
}
m2 := b2
m5 := b5
var mhi, mlo *big.Int
if leftright {
if mode < 2 {
if denorm {
i = be + (bias + (p - 1) - 1 + 1)
} else {
i = 1 + p - bbits
}
/* i is 1 plus the number of trailing zero bits in d's significand. Thus,
(2^m2 * 5^m5) / (2^(s2+i) * 5^s5) = (1/2 lsb of d)/10^k. */
} else {
j = ilim - 1
if m5 >= j {
m5 -= j
} else {
j -= m5
s5 += j
b5 += j
m5 = 0
}
i = ilim
if i < 0 {
m2 -= i
i = 0
}
/* (2^m2 * 5^m5) / (2^(s2+i) * 5^s5) = (1/2 * 10^(1-ilim))/10^k. */
}
b2 += i
s2 += i
mhi = big.NewInt(1)
/* (mhi * 2^m2 * 5^m5) / (2^s2 * 5^s5) = one-half of last printed (when mode >= 2) or
input (when mode < 2) significant digit, divided by 10^k. */
}
/* We still have d/10^k = (b * 2^b2 * 5^b5) / (2^s2 * 5^s5). Reduce common factors in
b2, m2, and s2 without changing the equalities. */
if m2 > 0 && s2 > 0 {
if m2 < s2 {
i = m2
} else {
i = s2
}
b2 -= i
m2 -= i
s2 -= i
}
b := new(big.Int).SetBytes(dblBits)
/* Fold b5 into b and m5 into mhi. */
if b5 > 0 {
if leftright {
if m5 > 0 {
pow5mult(mhi, m5)
b.Mul(mhi, b)
}
j = b5 - m5
if j != 0 {
pow5mult(b, j)
}
} else {
pow5mult(b, b5)
}
}
/* Now we have d/10^k = (b * 2^b2) / (2^s2 * 5^s5) and
(mhi * 2^m2) / (2^s2 * 5^s5) = one-half of last printed or input significant digit, divided by 10^k. */
S := big.NewInt(1)
if s5 > 0 {
pow5mult(S, s5)
}
/* Now we have d/10^k = (b * 2^b2) / (S * 2^s2) and
(mhi * 2^m2) / (S * 2^s2) = one-half of last printed or input significant digit, divided by 10^k. */
/* Check for special case that d is a normalized power of 2. */
spec_case := false
if mode < 2 {
if (_word1(d) == 0) && ((_word0(d) & bndry_mask) == 0) &&
((_word0(d) & (exp_mask & (exp_mask << 1))) != 0) {
/* The special case. Here we want to be within a quarter of the last input
significant digit instead of one half of it when the decimal output string's value is less than d. */
b2 += log2P
s2 += log2P
spec_case = true
}
}
/* Arrange for convenient computation of quotients:
* shift left if necessary so divisor has 4 leading 0 bits.
*
* Perhaps we should just compute leading 28 bits of S once
* and for all and pass them and a shift to quorem, so it
* can do shifts and ors to compute the numerator for q.
*/
var zz int
if s5 != 0 {
S_bytes := S.Bytes()
var S_hiWord uint32
for idx := 0; idx < 4; idx++ {
S_hiWord = S_hiWord << 8
if idx < len(S_bytes) {
S_hiWord |= uint32(S_bytes[idx])
}
}
zz = 32 - hi0bits(S_hiWord)
} else {
zz = 1
}
i = (zz + s2) & 0x1f
if i != 0 {
i = 32 - i
}
/* i is the number of leading zero bits in the most significant word of S*2^s2. */
if i > 4 {
i -= 4
b2 += i
m2 += i
s2 += i
} else if i < 4 {
i += 28
b2 += i
m2 += i
s2 += i
}
/* Now S*2^s2 has exactly four leading zero bits in its most significant word. */
if b2 > 0 {
b = b.Lsh(b, uint(b2))
}
if s2 > 0 {
S.Lsh(S, uint(s2))
}
/* Now we have d/10^k = b/S and
(mhi * 2^m2) / S = maximum acceptable error, divided by 10^k. */
if k_check {
if b.Cmp(S) < 0 {
k--
b.Mul(b, big10) /* we botched the k estimate */
if leftright {
mhi.Mul(mhi, big10)
}
ilim = ilim1
}
}
/* At this point 1 <= d/10^k = b/S < 10. */
if ilim <= 0 && mode > 2 {
/* We're doing fixed-mode output and d is less than the minimum nonzero output in this mode.
Output either zero or the minimum nonzero output depending on which is closer to d. */
if ilim >= 0 {
i = b.Cmp(S.Mul(S, big5))
}
if ilim < 0 || i < 0 || i == 0 && !biasUp {
/* Always emit at least one digit. If the number appears to be zero
using the current mode, then emit one '0' digit and set decpt to 1. */
buf = buf[:startPos]
buf = append(buf, '0')
return buf, 1
}
buf = append(buf, '1')
k++
return buf, k + 1
}
var dig byte
if leftright {
if m2 > 0 {
mhi.Lsh(mhi, uint(m2))
}
/* Compute mlo -- check for special case
* that d is a normalized power of 2.
*/
mlo = mhi
if spec_case {
mhi = mlo
mhi = new(big.Int).Lsh(mhi, log2P)
}
/* mlo/S = maximum acceptable error, divided by 10^k, if the output is less than d. */
/* mhi/S = maximum acceptable error, divided by 10^k, if the output is greater than d. */
var z, delta big.Int
for i = 1; ; i++ {
z.DivMod(b, S, b)
dig = byte(z.Int64() + '0')
/* Do we yet have the shortest decimal string
* that will round to d?
*/
j = b.Cmp(mlo)
/* j is b/S compared with mlo/S. */
delta.Sub(S, mhi)
var j1 int
if delta.Sign() <= 0 {
j1 = 1
} else {
j1 = b.Cmp(&delta)
}
/* j1 is b/S compared with 1 - mhi/S. */
if (j1 == 0) && (mode == 0) && ((_word1(d) & 1) == 0) {
if dig == '9' {
var flag bool
buf = append(buf, '9')
if buf, flag = roundOff(buf, startPos); flag {
k++
buf = append(buf, '1')
}
return buf, k + 1
}
if j > 0 {
dig++
}
buf = append(buf, dig)
return buf, k + 1
}
if (j < 0) || ((j == 0) && (mode == 0) && ((_word1(d) & 1) == 0)) {
if j1 > 0 {
/* Either dig or dig+1 would work here as the least significant decimal digit.
Use whichever would produce a decimal value closer to d. */
b.Lsh(b, 1)
j1 = b.Cmp(S)
if (j1 > 0) || (j1 == 0 && (((dig & 1) == 1) || biasUp)) {
dig++
if dig == '9' {
buf = append(buf, '9')
buf, flag := roundOff(buf, startPos)
if flag {
k++
buf = append(buf, '1')
}
return buf, k + 1
}
}
}
buf = append(buf, dig)
return buf, k + 1
}
if j1 > 0 {
if dig == '9' { /* possible if i == 1 */
buf = append(buf, '9')
buf, flag := roundOff(buf, startPos)
if flag {
k++
buf = append(buf, '1')
}
return buf, k + 1
}
buf = append(buf, dig+1)
return buf, k + 1
}
buf = append(buf, dig)
if i == ilim {
break
}
b.Mul(b, big10)
if mlo == mhi {
mhi.Mul(mhi, big10)
} else {
mlo.Mul(mlo, big10)
mhi.Mul(mhi, big10)
}
}
} else {
var z big.Int
for i = 1; ; i++ {
z.DivMod(b, S, b)
dig = byte(z.Int64() + '0')
buf = append(buf, dig)
if i >= ilim {
break
}
b.Mul(b, big10)
}
}
/* Round off last digit */
b.Lsh(b, 1)
j = b.Cmp(S)
if (j > 0) || (j == 0 && (((dig & 1) == 1) || biasUp)) {
var flag bool
buf, flag = roundOff(buf, startPos)
if flag {
k++
buf = append(buf, '1')
return buf, k + 1
}
} else {
buf = stripTrailingZeroes(buf, startPos)
}
return buf, k + 1
}
func bumpUp(buf []byte, k int) ([]byte, int) {
var lastCh byte
stop := 0
if len(buf) > 0 && buf[0] == '-' {
stop = 1
}
for {
lastCh = buf[len(buf)-1]
buf = buf[:len(buf)-1]
if lastCh != '9' {
break
}
if len(buf) == stop {
k++
lastCh = '0'
break
}
}
buf = append(buf, lastCh+1)
return buf, k
}
func setWord0(d float64, w uint32) float64 {
dBits := math.Float64bits(d)
return math.Float64frombits(uint64(w)<<32 | dBits&0xffffffff)
}
func _word0(d float64) uint32 {
dBits := math.Float64bits(d)
return uint32(dBits >> 32)
}
func _word1(d float64) uint32 {
dBits := math.Float64bits(d)
return uint32(dBits)
}
func stripTrailingZeroes(buf []byte, startPos int) []byte {
bl := len(buf) - 1
for bl >= startPos && buf[bl] == '0' {
bl--
}
return buf[:bl+1]
}
/* Set b = b * 5^k. k must be nonnegative. */
func pow5mult(b *big.Int, k int) *big.Int {
if k < (1 << (len(pow5Cache) + 2)) {
i := k & 3
if i != 0 {
b.Mul(b, p05[i-1])
}
k >>= 2
i = 0
for {
if k&1 != 0 {
b.Mul(b, pow5Cache[i])
}
k >>= 1
if k == 0 {
break
}
i++
}
return b
}
return b.Mul(b, new(big.Int).Exp(big5, big.NewInt(int64(k)), nil))
}
func roundOff(buf []byte, startPos int) ([]byte, bool) {
i := len(buf)
for i != startPos {
i--
if buf[i] != '9' {
buf[i]++
return buf[:i+1], false
}
}
return buf[:startPos], true
}
func init() {
p := big.NewInt(625)
pow5Cache[0] = p
for i := 1; i < len(pow5Cache); i++ {
p = new(big.Int).Mul(p, p)
pow5Cache[i] = p
}
}