var tests = [...]struct {
	name string
	in   float64 // used for error messages, not an input
	got  float64
	want float64
}{
	{"sqrt0", 0, math.Sqrt(0), 0},
	{"sqrt1", 1, math.Sqrt(1), 1},
	{"sqrt2", 2, math.Sqrt(2), math.Sqrt2},
	{"sqrt4", 4, math.Sqrt(4), 2},
	{"sqrt100", 100, math.Sqrt(100), 10},
	{"sqrt101", 101, math.Sqrt(101), 10.04987562112089},
}

var nanTests = [...]struct {
	name string

	Pi  = 3.14159265358979323846264338327950288419716939937510582097494459 // https://oeis.org/A000796
	Phi = 1.61803398874989484820458683436563811772030917980576286213544862 // https://oeis.org/A001622

	Sqrt2   = 1.41421356237309504880168872420969807856967187537694807317667974 // https://oeis.org/A002193
	SqrtE   = 1.64872127070012814684865078781416357165377610071014801157507931 // https://oeis.org/A019774

	switch {
	case math.IsNaN(x) || math.IsInf(x, 1):
		return x
	case x < 0:
		return math.NaN()
	case x == 0:
		return math.Inf(-1)
	}

	// reduce
	f1, ki := math.Frexp(x)
	if f1 < math.Sqrt2/2 {
		f1 *= 2
		ki--
	}
	f := f1 - 1
	k := float64(ki)

	// compute
	s := float64(f / (2 + f))
	s2 := float64(s * s)
	s4 := float64(s2 * s2)

	)

	// special cases
	switch {
	case IsNaN(x) || IsInf(x, 1):
		return x
	case x < 0:
		return NaN()
	case x == 0:
		return Inf(-1)
	}

	// reduce
	f1, ki := Frexp(x)
	if f1 < Sqrt2/2 {
		f1 *= 2
		ki--
	}
	f := f1 - 1
	k := float64(ki)

	// compute
	s := f / (2 + f)
	s2 := s * s
	s4 := s2 * s2
	t1 := s2 * (L1 + s4*(L3+s4*(L5+s4*L7)))

	ORPD    X0, X2 // X2= f1
	SHRQ    $52, BX
	ANDL    $0x7FF, BX
	SUBL    $0x3FE, BX
	XORPS   X1, X1 // break dependency for CVTSL2SD
	CVTSL2SD BX, X1 // x1= k, x2= f1
	// if f1 < math.Sqrt2/2 { k -= 1; f1 *= 2 }
	MOVSD   $HSqrt2, X0 // x0= 0.7071, x1= k, x2= f1
	CMPSD   X2, X0, 5 // cmpnlt; x0= 0 or ^0, x1= k, x2 = f1
	MOVSD   $1.0, X3 // x0= 0 or ^0, x1= k, x2 = f1, x3= 1

			return true
		}
	}
	return false
}

// alpha returns 1 if i is 0 and returns √2 otherwise.
func alpha(i int) float64 {
	if i == 0 {
		return 1
	}
	return math.Sqrt2
}

var cosines [32]float64 // cosines[k] = cos(π/2 * k/8)

func init() {
	for k := range cosines {
		cosines[k] = math.Cos(math.Pi * float64(k) / 16)
	}
}