math/big: underflow panic in roundShortest

[bcmills]

	f := new(big.Float).SetPrec(big.MaxPrec).SetFloat64(1)

	b := f.Append(nil, 'g', -1)

panics with the following stack trace:

panic: underflow

goroutine 1 [running]:
math/big.nat.sub(0x0, 0x0, 0x0, 0x104402c0, 0x0, 0x6, 0x180010, 0x1, 0x1, 0x0, ...)
	/usr/local/go/src/math/big/nat.go:129 +0x3e0
math/big.roundShortest(0x10429f38, 0x10429f90)
	/usr/local/go/src/math/big/ftoa.go:195 +0x1c0
math/big.(*Float).Append(0x10429f90, 0x0, 0x0, 0x0, 0x10429f67, 0xffffffff, 0x104340b0, 0x371f, 0x0, 0x1040a130)
	/usr/local/go/src/math/big/ftoa.go:97 +0x680
main.main()
	/tmp/sandbox082194187/main.go:11 +0xc0

(https://play.golang.org/p/68rBsop-pW)

[ALTree]

Can reproduce on tip and on go1.8.3 so not a recent regression.

[gopherbot]

Change https://golang.org/cl/52011 mentions this issue: math/big: fix roundShortest overflow

[rsc]

Broken since Go 1.8 so not for Go 1.10 at this point.

[zairm]

I've been looking into this and have an idea on how I might approach this.
However, I'd like a review at this point before moving on incase anyone feels there is a better approach.

Reference:

// math/big/ftoa.go
func roundShortest(d *decimal, x *Float) {
  	// *** SOME CODE ***

  	// ---PART 1---
  	mant := nat(nil).set(x.mant)
  	exp := int(x.exp) - mant.bitLen()
  	s := mant.bitLen() - int(x.prec+1)
  	switch {
  	case s < 0:
  		mant = mant.shl(mant, uint(-s))
  	case s > 0:
  		mant = mant.shr(mant, uint(+s))
  	}

  	// ---PART 2---
  	exp += s

  	var lower decimal
  	var tmp nat
  	lower.init(tmp.sub(mant, natOne), exp)

  	var upper decimal
  	upper.init(tmp.add(mant, natOne), exp)

  	// *** SOME CODE ***
}

PART 1
0 <= mant.bitLen() <= maxInt Because mant.bitLen() returns a non-neg int.
0 <= x.prec <= maxUint32 because x.prec is a uint32
So s can range in [-maxUint32-1, maxInt]
We want the calls to mant.shl and mant.shr to make sense.

The switch statement in part 1 can be changed to something to the effect of

// s = mant.bitLen() - x.prec - 1
if mant.bitLen()-1 < x.prec { // s negative case
    if mant.bitLen() > 0 {
        // mant.bitLen() cancels out the -1 so that |mant.bitLen()-1 - x.prec| <= maxUint
        mant.shl(mant, x.prec - uint(mant.bitLen - 1))
    } else {
        mant.shl(mant, x.prec) // x.prec is possibly maxUint32
        mant.shl(mant, 1)
    }
} else if mant.bitLen()-1 > x.prec { // s pos case.
    //  bitLen() - (x.prec+1) is no larger than maxInt32 on a 32-bit machine.
    mant.shr(mant, x.prec - uint(mant.bitLen) - x.prec -1) 
}

In the above I have made the assumption that mant.shl can be split into 2 calls and retain the same result.

Part 2
Given definition of s and exp in PART 1 and the addition in PART 2, exp = x.exp - (x.prec+1) which can also overflow on 32-bit systems.
Similar to PART 1, we have to make sure that the solution works with lower.init(..., exp) and upper.init(..., exp), both of which require their second arg to be an int, which will be an int32 on 32-bit systems.

We can apply an approach similar to PART 1 here. Change decimal.init(m nat, shift int) to decimal.init(m nat, shift uint, sb sign). decimal.init takes the sign of shift to mean the direction of shift and (after possibly changing the sign) converts it into a uint to call nat.shl, nat.shr, or decimal.shr. We can use sb to determine direction of shift instead.

Next step is to deal with the fact that exp has range [-maxInt32-maxUint32-1, maxInt32-1].
To fix this, the decimal.init calls can be replaced with something to the effect of

if x.exp-1 < x.prec { // exp=x.exp - x.prec - 1 is neg
    if x.exp-1 >= 0 { // |exp| <= maxUint32
        lower.init(tmp.sub(mant, natOne), x.prec - uint(x.exp - 1), neg)
        upper.init(tmp.sub(mant, natOne), x.prec - uint(x.exp - 1), neg)
    } else { // |exp| can exceed maxUint32
        lower.init(tmp.sub(mant, natOne), x.prec, neg)
        lower.newShr(x.exp-1)
        upper.init(tmp.add(mant, natOne), x.prec, neg)
        upper.newShr(x.exp-1)
    }
} else { // exp=x.exp - x.prec - 1 is pos
    lower.init(tmp.sub(mant, natOne), uint(x.exp) - x.prec - 1, pos)
    upper.init(tmp.sub(mant, natOne), uint(x.exp) - x.prec - 1, pos)
}

decimal.newShr will be a new function.

[ALTree]

@zairm easiest way to do this is to actually upload a draft patch to gerrit and ask for feedback there. It's much easier to read it and it's also really easy for others to cherrypick it locally and test it. You can write

DO NOT SUBMIT

in the commit message to make sure the patch doesn't get submitted in case it's not ready yet.

[ALTree]

Anyway I feel that any patch that requires to change a substantial amount of code and/or slows down the common code path just to fix the corner case reported in this issue has a low chance to be accepted. That's the reason I gave up on the old (incorrect) patch I had, actually. I coudn't find a simple way to fix this and I didn't want to significantly change mat/big code in several place to fix what seems to be a corner case.

[zairm]

I coudn't find a simple way to fix this and I didn't want to significantly change mat/big code in several place to fix what seems to be a corner case.

This is why I was hesitant to start work on it myself. It appears to me there are no easy solutions to something that seems simple and unlikely enough to prohibit even a moderately significant change.

I severely hope I'm wrong but my hunch is a fix would require significant changes to mat/big/decimal.go or significant enough changes to mat/big/float.go so as to require a documentation change.