Class: Horologium::Numeric::TwoPartFloat

Inherits:
Object
  • Object
show all
Defined in:
lib/horologium/numeric/two_part_float.rb,
sig/horologium/numeric/two_part_float.rbs

Overview

A number stored as the sum of two Floats, a high part and a low part. The low part carries the bits that do not fit in the high part, so the pair holds about twice the precision of a single Float.

The value is frozen on creation, and every operation returns a new one.

The representation and the way it is normalized follow the convention used by ERFA (a relicensed version of the SOFA library). The error-free transformations come from Shewchuk's work on robust floating-point arithmetic, and go back to Knuth and Dekker.

Examples:

Keeping precision a single Float would lose

a = Horologium::Numeric::TwoPartFloat.new(1.0)
b = Horologium::Numeric::TwoPartFloat.new(1e-16)
a + b == Horologium::Numeric::TwoPartFloat.new(1.0, 1e-16)
# => true
1.0 + 1e-16 == 1.0
# => true (the correction is lost)

See Also:

Constant Summary collapse

SPLIT_FACTOR =

The constant Dekker's split multiplies by, 2**27 + 1. It cuts a Float's 53-bit mantissa in two, so the halves multiply with no rounding.

Returns:

  • (Float) 
134_217_729.0

Instance Attribute Summary collapse

Class Method Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(high, low = 0.0) ⇒ TwoPartFloat

Returns a new instance of TwoPartFloat.

Parameters:

  • high (Float) —

    the high part

  • low (Float) (defaults to: 0.0) —

    the low part



33
34
35
36
37
 
# File 'lib/horologium/numeric/two_part_float.rb', line 33

def initialize(high, low = 0.0)
  @high = high
  @low = low
  freeze
end

Instance Attribute Details

#high ⇒ Float (readonly)

The two parts. Read them to pass the value to a foreign kernel that takes a day and a fraction of a day, such as an ERFA binding or a Chebyshev ephemeris segment. To add them together, use #to_f.

Returns:

  • (Float) 


259
260
261
 
# File 'lib/horologium/numeric/two_part_float.rb', line 259

def high
  @high
end

#low ⇒ Float (readonly)

The two parts. Read them to pass the value to a foreign kernel that takes a day and a fraction of a day, such as an ERFA binding or a Chebyshev ephemeris segment. To add them together, use #to_f.

Returns:

  • (Float) 


259
260
261
 
# File 'lib/horologium/numeric/two_part_float.rb', line 259

def low
  @low
end

Class Method Details

.fast_two_sum(larger, smaller) ⇒ Array(Float, Float)

This method is part of a private API. You should avoid using this method if possible, as it may be removed or be changed in the future.

Dekker's fast-two-sum, a quicker version of two_sum. It is only correct when larger is at least as large as smaller in magnitude, so use it when you already know which part is bigger.

Parameters:

  • larger (Float) —

    the part with the larger magnitude

  • smaller (Float) —

    the part with the smaller magnitude

Returns:

  • (Array(Float, Float)) —

    the sum and its rounding error



214
215
216
217
 
# File 'lib/horologium/numeric/two_part_float.rb', line 214

def self.fast_two_sum(larger, smaller)
  sum = larger + smaller
  [sum, smaller - (sum - larger)]
end

.from_real(value) ⇒ Horologium::Numeric::TwoPartFloat

Builds a two-part float from a single real number, keeping the precision a single Float would lose. The high part is the nearest Float and the low part carries the remainder.

Examples:

An integer past a Float's reach keeps its last digit

Horologium::Numeric::TwoPartFloat.from_real(2**53 + 1).to_r ==
  2**53 + 1
# => true

Parameters:

  • value (Numeric) —

    the number to represent

Returns:



169
170
171
172
 
# File 'lib/horologium/numeric/two_part_float.rb', line 169

def self.from_real(value)
  high = value.to_f
  new(high, (value.to_r - high.to_r).to_f)
end

.normalize(high, low = 0.0) ⇒ Horologium::Numeric::TwoPartFloat

Rebuilds a (high, low) pair into a canonical form: high is the nearest integer and low is the leftover fraction, between -0.5 and 0.5. Use it when the parts do not already follow that form, for example when low is larger than one half.

Examples:

A low part above one half carries into the high part

Horologium::Numeric::TwoPartFloat.normalize(2.0, 0.75) ==
  Horologium::Numeric::TwoPartFloat.new(3.0, -0.25)
# => true

Parameters:

  • high (Float) —

    the high part

  • low (Float) (defaults to: 0.0) —

    the low part

Returns:



142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
 
# File 'lib/horologium/numeric/two_part_float.rb', line 142

def self.normalize(high, low = 0.0)
  rounded = high.round.to_f
  remainder = high - rounded
  sum, error = two_sum(remainder, low)
  if sum.abs > 0.5
    carry = sum.round.to_f
    carry_sum, carry_error = two_sum(sum - carry, error)
    new_high = rounded + carry
    new_low = carry_sum + carry_error
  else
    new_high = rounded
    new_low = sum + error
  end

  new(new_high, new_low)
end

.split(value) ⇒ Array(Float, Float)

This method is part of a private API. You should avoid using this method if possible, as it may be removed or be changed in the future.

Splits a Float into a high and a low half that add back to the value and each multiply with no rounding. It is Dekker's method, cutting the mantissa in two with SPLIT_FACTOR.

Parameters:

  • value (Float)

Returns:

  • (Array(Float, Float)) —

    the high half and the low half



247
248
249
250
251
252
 
# File 'lib/horologium/numeric/two_part_float.rb', line 247

def self.split(value)
  scaled = SPLIT_FACTOR * value
  high = scaled - (scaled - value)
  low = value - high
  [high, low]
end

.two_diff(left, right) ⇒ Array(Float, Float)

This method is part of a private API. You should avoid using this method if possible, as it may be removed or be changed in the future.

Subtracts right from left and also returns the rounding error. Adding the two results back together gives left - right with no loss. It is the two-difference companion of two_sum and works for any two Floats.

Parameters:

  • left (Float)
  • right (Float)

Returns:

  • (Array(Float, Float)) —

    the difference and its rounding error



197
198
199
200
201
202
203
204
 
# File 'lib/horologium/numeric/two_part_float.rb', line 197

def self.two_diff(left, right)
  difference = left - right
  right_virtual = difference - left
  [
    difference,
    (left - (difference - right_virtual)) - (right + right_virtual)
  ]
end

.two_product(left, right) ⇒ Array(Float, Float)

This method is part of a private API. You should avoid using this method if possible, as it may be removed or be changed in the future.

Multiplies left and right and also returns the rounding error of the product. Adding the two results back together gives left * right with no loss. Ruby has no fused multiply-add, so each operand is split in two and the partial products stay exact.

Parameters:

  • left (Float)
  • right (Float)

Returns:

  • (Array(Float, Float)) —

    the product and its rounding error



228
229
230
231
232
233
234
235
236
237
238
 
# File 'lib/horologium/numeric/two_part_float.rb', line 228

def self.two_product(left, right)
  product = left * right
  left_high, left_low = split(left)
  right_high, right_low = split(right)
  error =
    ((left_high * right_high - product) +
      left_high * right_low +
      left_low * right_high) +
    left_low * right_low
  [product, error]
end

.two_sum(left, right) ⇒ Array(Float, Float)

This method is part of a private API. You should avoid using this method if possible, as it may be removed or be changed in the future.

Adds left and right and also returns the rounding error of the addition. Adding the two results back together gives left + right with no loss. This is the two-sum algorithm, due to Knuth. It works for any two Floats.

Parameters:

  • left (Float)
  • right (Float)

Returns:

  • (Array(Float, Float)) —

    the sum and its rounding error



183
184
185
186
187
 
# File 'lib/horologium/numeric/two_part_float.rb', line 183

def self.two_sum(left, right)
  sum = left + right
  right_virtual = sum - left
  [sum, (left - (sum - right_virtual)) + (right - right_virtual)]
end

Instance Method Details

#*(scalar) ⇒ Horologium::Numeric::TwoPartFloat

Parameters:

  • scalar (Integer, Float, Rational)

Returns:

Raises:

  • (ArgumentError) —

    when given anything but a plain number



66
67
68
69
70
71
72
73
 
# File 'lib/horologium/numeric/two_part_float.rb', line 66

def *(scalar) # rubocop:disable Naming/BinaryOperatorParameterName
  factor = scalar_float(scalar)
  high, low = two_sum(@high, @low)
  product, product_error = two_product(high, factor)
  result_high, result_low =
    fast_two_sum(product, product_error + low * factor)
  self.class.new(result_high, result_low)
end

#+(other) ⇒ Horologium::Numeric::TwoPartFloat

Parameters:

Returns:



41
42
43
44
45
46
47
48
49
 
# File 'lib/horologium/numeric/two_part_float.rb', line 41

def +(other)
  high_sum, high_error = two_sum(@high, other.high)
  low_sum, low_error = two_sum(@low, other.low)
  result_high, result_low = fast_two_sum(
    high_sum,
    high_error + low_error + low_sum
  )
  self.class.new(result_high, result_low)
end

#-(other) ⇒ Horologium::Numeric::TwoPartFloat

Parameters:

Returns:



53
54
55
56
57
58
59
60
61
 
# File 'lib/horologium/numeric/two_part_float.rb', line 53

def -(other)
  high_diff, high_error = two_diff(@high, other.high)
  low_diff, low_error = two_diff(@low, other.low)
  result_high, result_low = fast_two_sum(
    high_diff,
    high_error + low_error + low_diff
  )
  self.class.new(result_high, result_low)
end

#/(scalar) ⇒ Horologium::Numeric::TwoPartFloat

Parameters:

  • scalar (Integer, Float, Rational)

Returns:

Raises:

  • (ArgumentError) —

    when given anything but a plain number

  • (ZeroDivisionError) —

    when dividing by zero



79
80
81
82
83
84
85
86
87
88
89
90
91
 
# File 'lib/horologium/numeric/two_part_float.rb', line 79

def /(scalar) # rubocop:disable Naming/BinaryOperatorParameterName
  divisor = scalar_float(scalar)
  raise ZeroDivisionError, "divided by 0" if divisor.zero?

  high, low = two_sum(@high, @low)
  quotient = high / divisor
  product, product_error = two_product(quotient, divisor)
  remainder, remainder_error = two_diff(high, product)
  correction = (remainder_error + low) - product_error
  next_quotient = (remainder + correction) / divisor
  result_high, result_low = fast_two_sum(quotient, next_quotient)
  self.class.new(result_high, result_low)
end

#==(other) ⇒ Boolean

Compares the stored parts, not the number they add up to.

Parameters:

  • other (Object)

Returns:

  • (Boolean) 


114
115
116
 
# File 'lib/horologium/numeric/two_part_float.rb', line 114

def ==(other)
  other.is_a?(self.class) && high == other.high && low == other.low
end

#eql?(other) ⇒ Boolean

Parameters:

  • other (Object)

Returns:

  • (Boolean) 


120
121
122
 
# File 'lib/horologium/numeric/two_part_float.rb', line 120

def eql?(other)
  other.is_a?(self.class) && high.eql?(other.high) && low.eql?(other.low)
end

#fast_two_sum(larger, smaller) ⇒ [ Float, Float ]

Parameters:

  • larger (Float)
  • smaller (Float)

Returns:

  • ([ Float, Float ]) 


289
290
291
 
# File 'lib/horologium/numeric/two_part_float.rb', line 289

def fast_two_sum(larger, smaller)
  self.class.fast_two_sum(larger, smaller)
end

#hash ⇒ Integer

Returns:

  • (Integer) 


125
126
127
 
# File 'lib/horologium/numeric/two_part_float.rb', line 125

def hash
  [high, low].hash
end

#scalar_float(scalar) ⇒ Float

A plain number, as a Float. Multiplying or dividing by another two-part value would collapse it to one Float and lose the low part, so it is refused rather than silently rounded.

Parameters:

  • scalar (Integer, Float, Rational) —

    the number to check

Returns:

  • (Float) —

    the number as a Float

Raises:

  • (ArgumentError) —

    when it is not a plain number



270
271
272
273
274
275
276
277
278
279
 
# File 'lib/horologium/numeric/two_part_float.rb', line 270

def scalar_float(scalar)
  case scalar
  when Integer, Float, Rational
    scalar.to_f
  else
    raise ArgumentError,
      "a TwoPartFloat multiplies and divides by a plain number, " \
      "got a #{scalar.class}"
  end
end

#to_f ⇒ Float

The value as a single Float. One Float cannot hold what two hold, so the extra precision the split carries is dropped here. Do it at the end, once the arithmetic is done.

Returns:

  • (Float) 


106
107
108
 
# File 'lib/horologium/numeric/two_part_float.rb', line 106

def to_f
  high + low
end

#to_r ⇒ Rational

The two parts added with no loss. Each Float is an exact rational, so their sum is exact and keeps the low part.

Returns:

  • (Rational) 


97
98
99
 
# File 'lib/horologium/numeric/two_part_float.rb', line 97

def to_r
  high.to_r + low.to_r
end

#two_diff(left, right) ⇒ [ Float, Float ]

Parameters:

  • left (Float)
  • right (Float)

Returns:

  • ([ Float, Float ]) 


285
286
287
 
# File 'lib/horologium/numeric/two_part_float.rb', line 285

def two_diff(left, right)
  self.class.two_diff(left, right)
end

#two_product(left, right) ⇒ [ Float, Float ]

Parameters:

  • left (Float)
  • right (Float)

Returns:

  • ([ Float, Float ]) 


293
294
295
 
# File 'lib/horologium/numeric/two_part_float.rb', line 293

def two_product(left, right)
  self.class.two_product(left, right)
end

#two_sum(left, right) ⇒ [ Float, Float ]

Parameters:

  • left (Float)
  • right (Float)

Returns:

  • ([ Float, Float ]) 


281
282
283
 
# File 'lib/horologium/numeric/two_part_float.rb', line 281

def two_sum(left, right)
  self.class.two_sum(left, right)
end