sqrt.gno

  1package physics
  2
  3import "math"
  4
  5// fsqrt is math.Sqrt, bit for bit, at a fraction of its gas.
  6//
  7// gno's math.Sqrt computes the root one bit at a time in interpreted code
  8// (~170k gas); the engine takes several roots every millisecond step. Here a
  9// float Newton iteration gets within a unit in the last place, and an exact
 10// integer test settles the last bit: for x = M * 2^e (e even) the correctly
 11// rounded root is R * 2^q, q = e/2 - 26, where R is the integer nearest
 12// sqrt(N), N = M * 2^52, and R is that integer exactly when
 13// R*R - R < N <= R*R + R (no square root of an integer lies halfway between
 14// two integers). N - R*R is small, so it is computed exactly mod 2^64. A
 15// Newton result one ulp off is moved by that ulp; anything else (rare) goes
 16// to sqrtSlow. Zero, subnormals, infinities and NaN go to math.Sqrt.
 17func fsqrt(x float64) float64 {
 18	b := math.Float64bits(x)
 19	if b-1<<52 >= 0x7fe0000000000000 { // negative, zero, subnormal, infinite or NaN
 20		return math.Sqrt(x)
 21	}
 22	y := math.Float64frombits((b >> 1) + 0x1ff8000000000000)
 23	y = 0.5 * (y + x/y)
 24	y = 0.5 * (y + x/y)
 25	y = 0.5 * (y + x/y)
 26	y = 0.5 * (y + x/y)
 27	yb := math.Float64bits(y)
 28	// x = m * 2^(ex), y = R * 2^q: y is the root when R*R - R < N <= R*R + R
 29	// for N = m * 2^(ex-2q), which is m << 52 or m << 53 when y is in the
 30	// root's binade.
 31	s := b>>52 - yb>>52<<1 + 1075
 32	if s-52 > 1 {
 33		return sqrtSlow(x)
 34	}
 35	R := yb&(1<<52-1) | 1<<52
 36	N := (b&(1<<52-1) | 1<<52) << s
 37	d := int64(N - R*R)
 38	if d > int64(R) {
 39		// One ulp low: R+1 is the root when N - (R+1)^2 = d - 2R - 1 is in range.
 40		if d2 := d - int64(R<<1) - 1; d2 <= int64(R+1) && d2 > -int64(R+1) && R+1 < 1<<53 {
 41			return math.Float64frombits(yb + 1)
 42		}
 43		return sqrtSlow(x)
 44	}
 45	if d <= -int64(R) {
 46		if d2 := d + int64(R<<1) - 1; d2 <= int64(R-1) && d2 > -int64(R-1) && R > 1<<52 {
 47			return math.Float64frombits(yb - 1)
 48		}
 49		return sqrtSlow(x)
 50	}
 51	return y
 52}
 53
 54// sqrtSlow settles the roots the fast path cannot: the general case.
 55func sqrtSlow(x float64) float64 {
 56	b := math.Float64bits(x)
 57	exp := int(b>>52) & 0x7ff
 58	if b>>63 != 0 || exp == 0 || exp == 0x7ff {
 59		return math.Sqrt(x)
 60	}
 61	// Halve the exponent, then four Newton steps: within an ulp.
 62	y := math.Float64frombits((b >> 1) + 0x1ff8000000000000)
 63	y = 0.5 * (y + x/y)
 64	y = 0.5 * (y + x/y)
 65	y = 0.5 * (y + x/y)
 66	y = 0.5 * (y + x/y)
 67
 68	m := b&(1<<52-1) | 1<<52
 69	e := exp - 1075 // x = m * 2^e
 70	if e&1 != 0 {
 71		m <<= 1
 72		e--
 73	}
 74	N := m << 52
 75	yb := math.Float64bits(y)
 76	if int(yb>>52) != e/2+1049 { // y is not in the root's binade: settle it the long way
 77		return rootOf(N, e/2-26, y)
 78	}
 79	R := yb&(1<<52-1) | 1<<52 // y = R * 2^q
 80	d := int64(N - R*R)
 81	if d > int64(R) {
 82		if d2 := int64(N - (R+1)*(R+1)); d2 > int64(R+1) || d2 <= -int64(R+1) {
 83			return rootOf(N, e/2-26, y)
 84		}
 85		return math.Float64frombits(yb + 1)
 86	}
 87	if d <= -int64(R) {
 88		if R == 1<<52 {
 89			return rootOf(N, e/2-26, y)
 90		}
 91		if d2 := int64(N - (R-1)*(R-1)); d2 > int64(R-1) || d2 <= -int64(R-1) {
 92			return rootOf(N, e/2-26, y)
 93		}
 94		return math.Float64frombits(yb - 1)
 95	}
 96	return y
 97}
 98
 99// rootOf is the general case: R from the estimate y, an integer Newton step,
100// then whole units until R is the nearest integer to sqrt(N).
101func rootOf(N uint64, q int, y float64) float64 {
102	R := uint64(y * math.Float64frombits(uint64(1023-q)<<52))
103	d := int64(N - R*R)
104	R = uint64(int64(R) + d/int64(2*R))
105	for {
106		d = int64(N - R*R)
107		if d > int64(R) {
108			R++
109		} else if d <= -int64(R) {
110			R--
111		} else {
112			break
113		}
114	}
115	return float64(R) * math.Float64frombits(uint64(q+1023)<<52)
116}