0

For which input does the following function reduce to the identity (modulo 1 ulp)?

public static double invertTwice(double value) {
  double inverse = 1.0 / value;
  return 1.0 / inverse;
}

Also, what is the largest interval [min, max] for which the function is the identity (modulo 1 ulp)?

How come values in the range [1E-320...1E-310] are mapped to +Infinity when inverted? (A subsequent inversion of which will result in 0.0!)

Motivation: the general relation

a / b == a * (1.0 / b)

does not hold for all double pairs, even when b != 0.0, for instance a = 0.0 and b = 4.9E-324.

This has to be taken into account when normalizing vectors, for instance

[0.0 b 0.0] * (1.0 / b) == [NaN Infinity NaN]

but

[0.0/b b/b 0.0/b] == [0.0 1.0 0.0]

Remark: My question is not about normalizing vectors, but about what I wrote above. Thank you!

datahaki
  • 600
  • 7
  • 23
  • sounds similar to [this question](https://stackoverflow.com/questions/44623331/does-a-floating-point-reciprocal-always-round-trip/44626041#44626041). The `invertTwice` function is the identity on the interval `[1.0, sqrt(2))`, and on the result of scaling that interval by any power of two, provided underflow and overflow are avoided. On `[sqrt(2), 2.0)`, it's somewhat random, with values close to `sqrt(2)` _usually_ giving the identity, and values close to `2.0` giving the identity with "probability" around 0.5. – Mark Dickinson Jul 31 '17 at 10:12
  • thanks for the pointer and the partial answer. How come values in the range [1e-323...1e-310] all map to Infinity when inverted? (i am using Java) – datahaki Jul 31 '17 at 10:33
  • Because the max representable IEEE 754 double precision float is approximately 1.8e308, so for anything smaller than around 5.56e-309, its reciprocal is larger than the largest representable float, so "rounds" to infinity. – Mark Dickinson Jul 31 '17 at 11:46
  • Ah, I see you edited to include the 1 ulp condition. That's a different question altogether. Note that very large floats (near the 1.8e308 boundary) can lose more than 1 ulp on a double inversion, thanks to the loss of precision in the subnormal range. Do you need "up to 1 ulp", or is "up to a handful of ulps" good enough? – Mark Dickinson Jul 31 '17 at 11:48
  • thanks for the insights! you are right: "up to a handful of ulps" would be good enough. i mainly care about closeness under inversion. in my case, repeated inversion should not diverge to Infinity. – datahaki Jul 31 '17 at 17:31

1 Answers1

0

using bisection search, I found the interval of double values that is invariant under 2x inversion

value == 1.0 / (1.0 / value)

to be

[5.562684646268010E-309 ... 1.7976931348623151E308]
datahaki
  • 600
  • 7
  • 23