ieee754 and newton-raphson

sqrt.go

package main

import (
	"fmt"
	"math"
	"unsafe"
)

// IEEE-754 floating-point representation.
func ieee754BinRep(x float64) uint64 {
	return *(*uint64)(unsafe.Pointer(&x))
}

// Turn IEEE-754 floating-point representation
// into true floating-point form.
func ieee754ToFloatingPoint(x uint64) float64 {
	return *(*float64)(unsafe.Pointer(&x))
}

// Precisely calculating square root
// using Newton-Raphson method.
func Sqrt(x float64) float64 {
	switch {
	case x == 0 || math.IsNaN(x) || math.IsInf(x, 1):
		return x
	case x < 0:
		return math.NaN()
	}
	
	ix  := ieee754BinRep(x)
	exp := int((ix >> 32) & 0xffffffff)
	
	if exp == 0 {
		for ix & (1 << 32) == 0 {
			ix <<= 1
			exp--
		}
		exp++
	}
	
	exp -= 1023
	ix &^= 0xffffffff << 32
	ix |= 1 << 32
	
	if exp & 1 == 1 {
		ix <<= 1
	}
	
	exp >>= 1
	ix  <<= 1
	
	var q, s uint64
	
	r := uint64(1 << 33)
	
	for r != 0 {
		t := s + r
		if t <= ix {
			s = t + r
			ix -= t
			q += r
		}
		ix <<= 1
		r  >>= 1
	}
	
	if ix != 0 {
		q += q & 1
	}
	
	ix = q >> 1 + uint64(exp - 1 + 1023) << 32
	return ieee754ToFloatingPoint(ix)
}

func main() {
	// see :)
	fmt.Println(Sqrt(2), math.Sqrt(2))
}