3 ms·

Interesting... I think we might be interpreting the statement differently? It seems like you're saying the statement claims the difficulty is in deciding whethe

Interesting... I think we might be interpreting the statement differently? It seems like you're saying the statement claims the difficulty is in deciding whether f(x) = 0 or f(x) = 1 when we see f(x) ≈ 0.5? If that is the only problem then the obvious solution would be "just use f(x) = 0.5 instead", right? (i.e. allow for an extra digit.) Somehow that's not how I've been interpreting the statement until now. Rather, the way I read it, we don't know how to predict the amount of work it takes to guarantee correct rounding to any number of digits (even if we allow for extra digits beyond those we really want)—not merely to a fixed number of digits. Or to put it another way, it sounded as if for pathological inputs even guaranteeing >= N digits of accuracy isn't possible with a predictable amount of effort. What led me to believe this is a couple of different statements from Kahan [1]: A desideratum generally acknowledged but not always achieved is: "In error by rather less than an ulp if by more than half." Reputable math libraries compute elementary transcendental functions mostly within slightly more than half an ulp and almost always well within one ulp. Both of these suggested to me that errors can occasionally exceed half an ulp by a non-negligible amount. So, if we were rounding to the nearest integer, this would imply occasionally not only rounding arbitrary overshoots like 0.5 + 1E-100 to zero (whose difficulty you're pointing out), but in fact even occasionally even producing 0 when the output is 0.7. That's far worse. Is this interpretation wrong? If it is, then why should errors sometimes not be "well within one ulp"? https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT

yorwba 6y ago
> If that is the only problem then the obvious solution would be "just use f(x) = 0.5 instead", right? (i.e. allow for an extra digit.) Yes, but in many cases that's not desirable. E.g. if you want to compute things in hardware where adding an extra digit means adding more transistors to your chip. (And then why not always use that extra digit instead of rounding to one digit less?) > mostly within slightly more than half an ulp The reason this doesn't say "at most half an ulp" but allows for "slightly more" is the table-makers dilemma. (If you keep reading after the part you quoted, you'll see that he's talking about precisely the case where a value is close to x.5) > almost always well within one ulp This is a different case, concerning functions that are expensive to compute to a high degree of precision (e.g. due to high sensitivity to rounding of intermediate results) but "reputable" libraries will still try really hard to give good results. In terms of the midpoint-radius representation, the first is about not knowing how much smaller the radius needs to be to achieve correct rounding to a predetermined number of digits, whereas the second is about needing to wait for a long (but usually predictable) time before the radius drops below a predetermined target.
- dataflow 6y ago
  > This is a different case, concerning functions that are expensive to compute to a high degree of precision (e.g. due to high sensitivity to rounding of intermediate results) but "reputable" libraries will still try really hard to give good results. This is what I was referring to. See where it says "Why can't y^w be rounded within half an ulp like SQRT? Because nobody knows how much computation it would cost."... that's why I was confused how a arbitrary-precision library would handle something like this. User asks to raise some long decimal to another long decimal power and it seems to me you might sometimes either have to return a very wrong result or spend so much time computing that it won't be worth it.
  - yorwba 6y ago
    > User asks to raise some long decimal to another long decimal power and it seems to me you might sometimes either have to return a very wrong result or spend so much time computing that it won't be worth it. y^w isn't that expensive. Usually it gets expressed roughly as exp(log(y^w)) = exp(w * log(y)), where exp and log are implemented using polynomial series. So if you want more precision, you just need to evaluate a few more polynomial terms.
    - dataflow 6y ago
      If just a few more terms would suffice then why does he say "nobody knows how much computation it would cost"? Are you saying he's just referring to the half-ulp + epsilon case here too? If so then I've definitely mixed these up!
      - yorwba 6y ago
        Yeah, the question he asks is "Why can't y^w be rounded within half an ulp like SQRT?" not "Why can't y^w be rounded within slightly more than half an ulp like SQRT?" It's really just that tiny epsilon difference that turns "just a few more terms" into "nobody knows how much it'd cost". The reason SQRT can actually be rounded within half an espilon in a predictable amount of time is that SQRT(x) >= y + 0.5 precisely when x >= (y + 0.5)^2 = y^2 + y + 0.25, which means you can determine whether to round up or not with a simple comparison.

yorwba 6y ago

> If that is the only problem then the obvious solution would be "just use f(x) = 0.5 instead", right? (i.e. allow for an extra digit.) Yes, but in many cases that's not desirable. E.g. if you want to compute things in hardware where adding an extra digit means adding more transistors to your chip. (And then why not always use that extra digit instead of rounding to one digit less?) > mostly within slightly more than half an ulp The reason this doesn't say "at most half an ulp" but allows for "slightly more" is the table-makers dilemma. (If you keep reading after the part you quoted, you'll see that he's talking about precisely the case where a value is close to x.5) > almost always well within one ulp This is a different case, concerning functions that are expensive to compute to a high degree of precision (e.g. due to high sensitivity to rounding of intermediate results) but "reputable" libraries will still try really hard to give good results. In terms of the midpoint-radius representation, the first is about not knowing how much smaller the radius needs to be to achieve correct rounding to a predetermined number of digits, whereas the second is about needing to wait for a long (but usually predictable) time before the radius drops below a predetermined target.

dataflow 6y ago
> This is a different case, concerning functions that are expensive to compute to a high degree of precision (e.g. due to high sensitivity to rounding of intermediate results) but "reputable" libraries will still try really hard to give good results. This is what I was referring to. See where it says "Why can't y^w be rounded within half an ulp like SQRT? Because nobody knows how much computation it would cost."... that's why I was confused how a arbitrary-precision library would handle something like this. User asks to raise some long decimal to another long decimal power and it seems to me you might sometimes either have to return a very wrong result or spend so much time computing that it won't be worth it.
- yorwba 6y ago
  > User asks to raise some long decimal to another long decimal power and it seems to me you might sometimes either have to return a very wrong result or spend so much time computing that it won't be worth it. y^w isn't that expensive. Usually it gets expressed roughly as exp(log(y^w)) = exp(w * log(y)), where exp and log are implemented using polynomial series. So if you want more precision, you just need to evaluate a few more polynomial terms.
  - dataflow 6y ago
    If just a few more terms would suffice then why does he say "nobody knows how much computation it would cost"? Are you saying he's just referring to the half-ulp + epsilon case here too? If so then I've definitely mixed these up!
    - yorwba 6y ago
      Yeah, the question he asks is "Why can't y^w be rounded within half an ulp like SQRT?" not "Why can't y^w be rounded within slightly more than half an ulp like SQRT?" It's really just that tiny epsilon difference that turns "just a few more terms" into "nobody knows how much it'd cost". The reason SQRT can actually be rounded within half an espilon in a predictable amount of time is that SQRT(x) >= y + 0.5 precisely when x >= (y + 0.5)^2 = y^2 + y + 0.25, which means you can determine whether to round up or not with a simple comparison.