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)

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.

The two parts. They are protected so the operators can read another value's parts without exposing them publicly.

Returns:

  • (Float) 


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

def high
  @high
end

#low ⇒ Float (readonly)

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.

The two parts. They are protected so the operators can read another value's parts without exposing them publicly.

Returns:

  • (Float) 


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

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



246
247
248
249
 
# File 'lib/horologium/numeric/two_part_float.rb', line 246

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:



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

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 fraction stays in the low part

Horologium::Numeric::TwoPartFloat.normalize(2.0, 0.3) ==
  Horologium::Numeric::TwoPartFloat.new(2.0, 0.3)
# => true

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:



174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
 
# File 'lib/horologium/numeric/two_part_float.rb', line 174

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



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

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



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

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



260
261
262
263
264
265
266
267
268
269
270
 
# File 'lib/horologium/numeric/two_part_float.rb', line 260

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



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

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

Multiplies by a plain number, keeping the precision a single Float would lose. The number is a scalar, not another two-part value.

Examples:

Horologium::Numeric::TwoPartFloat.new(1.5) * 2 ==
  Horologium::Numeric::TwoPartFloat.new(3.0)
# => true

Parameters:

  • scalar (Numeric) —

    the number to multiply by

Returns:



86
87
88
89
90
91
92
93
 
# File 'lib/horologium/numeric/two_part_float.rb', line 86

def *(scalar) # rubocop:disable Naming/BinaryOperatorParameterName
  factor = scalar.to_f
  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

Adds another value and keeps the extra precision from both parts.

Examples:

Horologium::Numeric::TwoPartFloat.new(1.5) +
  Horologium::Numeric::TwoPartFloat.new(2.25) ==
  Horologium::Numeric::TwoPartFloat.new(3.75)
# => true

Parameters:

Returns:



48
49
50
51
52
53
54
55
56
 
# File 'lib/horologium/numeric/two_part_float.rb', line 48

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

Subtracts another value and keeps the extra precision from both parts.

Examples:

Horologium::Numeric::TwoPartFloat.new(1.5) -
  Horologium::Numeric::TwoPartFloat.new(0.25) ==
  Horologium::Numeric::TwoPartFloat.new(1.25)
# => true

Parameters:

Returns:



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

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

Divides by a plain number, keeping the precision a single Float would lose. The number is a scalar, not another two-part value.

Examples:

Horologium::Numeric::TwoPartFloat.new(3.0) / 2 ==
  Horologium::Numeric::TwoPartFloat.new(1.5)
# => true

Parameters:

  • scalar (Numeric) —

    the number to divide by

Returns:

Raises:

  • (ZeroDivisionError) 


104
105
106
107
108
109
110
111
112
113
114
115
116
 
# File 'lib/horologium/numeric/two_part_float.rb', line 104

def /(scalar) # rubocop:disable Naming/BinaryOperatorParameterName
  divisor = scalar.to_f
  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

Two values are equal when their high and low parts are both equal. This compares the stored parts, not the number they add up to.

Examples:

obj = Horologium::Numeric::TwoPartFloat.new(1.0, 0.5)
obj == Horologium::Numeric::TwoPartFloat.new(1.0, 0.5) # => true

Parameters:

  • other (Object) —

    the object to compare

Returns:

  • (Boolean) —

    true when other is a TwoPartFloat with equal parts



139
140
141
 
# File 'lib/horologium/numeric/two_part_float.rb', line 139

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

#eql?(other) ⇒ Boolean

The stricter equality used for Hash keys and Sets. It matches when the high and low parts are eql?, so an integer part and a float part differ.

Parameters:

  • other (Object) —

    the object to compare

Returns:

  • (Boolean) —

    true when other is a TwoPartFloat with eql? parts



148
149
150
 
# File 'lib/horologium/numeric/two_part_float.rb', line 148

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 ]) 


305
306
307
 
# File 'lib/horologium/numeric/two_part_float.rb', line 305

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

#hash ⇒ Integer

Returns a hash built from the two parts.

Returns:

  • (Integer) —

    a hash built from the two parts



153
154
155
 
# File 'lib/horologium/numeric/two_part_float.rb', line 153

def hash
  [high, low].hash
end

#to_r ⇒ Rational

The value as a Rational: the two parts added with no loss. Each Float is an exact rational, so their sum is exact and keeps the low part's extra precision.

Examples:

Horologium::Numeric::TwoPartFloat.new(1.0, 1e-16).to_r ==
  1.to_r + Rational(1e-16)
# => true

Returns:

  • (Rational) —

    the sum of the high and low parts



127
128
129
 
# File 'lib/horologium/numeric/two_part_float.rb', line 127

def to_r
  high.to_r + low.to_r
end

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

Parameters:

  • left (Float)
  • right (Float)

Returns:

  • ([ Float, Float ]) 


301
302
303
 
# File 'lib/horologium/numeric/two_part_float.rb', line 301

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 ]) 


309
310
311
 
# File 'lib/horologium/numeric/two_part_float.rb', line 309

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 ]) 


297
298
299
 
# File 'lib/horologium/numeric/two_part_float.rb', line 297

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