Long Double

longdouble.coffee

two_sum = (a, b) ->
	s = a + b; d = s - a
	
	return [s, (a - (s - d)) + (b - d)]

two_diff = (a, b) ->
	s = a - b; d = s - a
	
	return [s, (a - (s - d)) - (b + d)]

SPLITTER     =  134217729.0               # = 2^27 + 1
SPLIT_THRESH =  6.69692879491417e+299     # = 2^996


split = (a) ->
	if (a > SPLIT_THRESH || a < -SPLIT_THRESH)
		a *= 3.7252902984619140625e-09    # 2^-28
		
		temp = SPLITTER * a
		
		hi = temp - (temp - a); lo = a - hi
		
		hi *= 268435456.0                 # 2^28
		lo *= 268435456.0                 # 2^28
	else
		temp = SPLITTER * a
		hi = temp - (temp - a); lo = a - hi
	
	return [hi, lo]

if (false) # Fused Multiply-Add is never implemented in JS
	two_prod = (a, b) ->
		p = a * b
		
		return [p, Math.fma(a, b, -p)]
	
else
	two_prod = (a, b) ->
		p = a * b
		
		[a_hi, a_lo] = split(a)
		[b_hi, b_lo] = split(b)
		
		error = ((a_hi * b_hi - p) + a_hi * b_lo + a_lo * b_hi) + a_lo * b_lo
		
		return [p, error]
	

quick_two_sum = (a, b) ->
	s = a + b
	
	return [s, b - (s - a)]

quick_two_diff = (a, b) ->
	s = a - b
	
	return [s, (a - s) - b]

class CSLongDouble
	constructor: (@sum, @error) ->
	
	negate: () ->
		return new CSLongDouble(-@sum, -@error)
	
	add: (other) ->
		[s1, s2] = two_sum(@sum, other.sum)
		[t1, t2] = two_sum(@error, other.error)
		
		s2 += t1
		
		[s1, s2] = quick_two_sum(s1, s2)
		
		s2 += t2
		
		[s1, s2] = quick_two_sum(s1, s2)
		
		return new CSLongDouble(s1, s2);
	
	addFast: (other) ->
		[s, e] = two_sum(sum, other.sum)
		
		e += (@error + other.error)
		
		[s, e] = quick_two_sum(s, e)
		
		return new CSLongDouble(s, e)
	
	addDouble: (other) ->
		[s1, s2] = two_sum(@sum, other)
		
		s2 += @error
		
		[s1, s2] = quick_two_sum(s1, s2)
		
		return new CSLongDouble(s1, s2)
	
	sub: (other) ->
		[s1, s2] = two_diff(@sum, other.sum)
		[t1, t2] = two_diff(@error, other.error)
		
		s2 += t1
		
		[s1, s2] = quick_two_sum(s1, s2)
		
		s2 += t2
		
		[s1, s2] = quick_two_sum(s1, s2)
		
		return new CSLongDouble(s1, s2)
	
	subFast: (other) ->
		[s, e] = two_diff(@sum, other.sum)
		
		e += @error; e -= other.error
		
		[s, e] = quick_two_sum(s, e)
		
		return new CSLongDouble(s, e)
	
	subDouble: (other) ->
		[s1, s2] = two_diff(@sum, other)
		
		s2 += @error
		
		[s1, s2] = quick_two_sum(s1, s2)
		
		return new CSLongDouble(s1, s2)
	
	mul: (other) ->
		[p1, p2] = two_prod(@sum, other.sum)
		
		p2 += (@sum * other.error + @error * other.sum);
		
		[p1, p2] = quick_two_sum(p1, p2)
		
		return new CSLongDouble(p1, p2)
	
	mulDouble: (other) ->
		[p1, p2] = two_prod(@sum, other)
		
		p2 += (@error * other)
		
		[p1, p2] = quick_two_sum(p1, p2)
		
		return new CSLongDouble(p1, p2)
	
	div: (other) ->
		q1 = @sum / other.sum
		
		r = this.sub(other.mulDouble(q1))
		
		q2 = r.sum / other.sum
		
		r.decrement(other.mulDouble(q2))
		
		q3 = r.sum / other.sum
		
		[q1, q2] = quick_two_sum(q1, q2)
		
		r = new CSLongDouble(q1, q2)
		
		return r.addDouble(q3)
	
	divFast: (other) ->
		q1 = @sum / other.sum
		
		r = other.mulDouble(q1)
		
		[s1, s2] = two_diff(@sum, r.sum)
		
		s2 -= r.error
		s2 += @error
		
		q2 = (s1 + s2) / other.sum
		
		[r.sum, r.error] = quick_two_sum(q1, q2)
		
		return r
	
	divDouble: (other) ->
		q1 = @sum / other
		
		[p1, p2] = two_prod(q1, other)
		
		[s1, s2] = two_diff(@sum, p1)
		
		s2 += @error
		s2 -= p2
		
		q2 = (s1 + s2) / other
		
		[q1, q2] = quick_two_sum(q1, q2)
		
		return new CSLongDouble(q1, q2)
	
	increment: (other) ->
		[s1, s2] = two_sum(@sum, other.sum)
		[t1, t2] = two_sum(@error, other.error)
		
		s2 += t1
		
		[s1, s2] = quick_two_sum(s1, s2)
		
		s2 += t2;
		
		[@sum, @error] = quick_two_sum(s1, s2)
		
		return this
	
	incrementFast: (other) ->
		[s, e] = two_sum(@sum, other.sum)
		
		e += @error; e += other.error
		
		[@sum, @error] = quick_two_sum(s, e)
		
		return this
	
	incrementDouble: (other) ->
		[s1, s2] = two_sum(@sum, other)
		
		s2 += @error
		
		[@sum, @error] = quick_two_sum(s1, s2)
		
		return this
	
	decrement: (other) ->
		[s1, s2] = two_diff(@sum, other.sum)
		[t1, t2] = two_diff(@error, other.error)
		
		s2 += t1
		
		[s1, s2] = quick_two_sum(s1, s2)
		
		s2 += t2
		
		[@sum, @error] = quick_two_sum(s1, s2)
		
		return this
	
	decrementFast: (other) ->
		[s, e] = two_diff(@sum, other.sum)
		
		e += @error; e -= other.error
		
		[@sum, @error] = quick_two_sum(s, e)
		
		return this
	
	decrementDouble: (other) ->
		[s1, s2] = two_diff(@sum, other)
		
		s2 += @error
		
		[@sum, @error] = quick_two_sum(s1, s2)
		
		return this
	

window.CSLongDouble = CSLongDouble

test.html

<html>
  <head>
    <title>Long Double</title>
    
    <script type='application/javascript' src='longdouble.js'></script>
    
    <script type='text/javascript'>
      a = new CSLongDouble(1.0, 0.0)
      b = new CSLongDouble(2.0, 0.0)
      c = new CSLongDouble(3.0, 0.0)
      
      console.log("a + b", a.add(b))
      console.log("b + c", b.add(c))
      
      console.log("a + 3", a.addDouble(3))
      console.log("b + 7", b.addDouble(7))
      
      console.log("a - b", a.sub(b))
      console.log("b - c", b.sub(c))
      
      console.log("a - 3", a.subDouble(3))
      console.log("b - 7", b.subDouble(7))
      
      console.log("a * b", a.mul(b))
      console.log("b * c", b.mul(c))
      
      console.log("a * 3", a.mulDouble(3))
      console.log("b * 7", b.mulDouble(7))
      
      console.log("a / b", a.div(b))
      console.log("b / c", b.div(c))
      
      console.log("a / 3", a.divDouble(3))
      console.log("b / 7", b.divDouble(7))
    </script>
  </head>
</html>