sqrt(negative)

[Tomek Czajka]

I think the fallback value for sqrt(negative) should be 0, not ∞. Much more useful.

This makes the function continuous at 0.

Imagine you're making a calculation such as sqrt(some formula - some other formula). Due to rounding errors you might end up with sqrt(-tiny), while the perfectly accurate answer is sqrt(+tiny). You would want 0 and to continue with the computation, which is perfectly fine.

Of course the programmer could just do sqrt(max(x, 0)) instead, so there is a workaround.

[MitchAlsup]

IEEE 754-2018 has a statement that::

SQRT(+0.0) = +0.0

and

SQRT(-0.0) = -0.0

and

and SQRT(-value) = NaN

I do not take a position on SQRT(-value) wrt unums; but definitional compatibility seems a worthy goal.

[Tomek Czajka]

My understanding is that with posits, you get a NaN **and** a fallback value of infinity after handling of NaN. Whereas with IEEE 754 it's just a NaN.

I'm saying it's better to make it a NaN with a fallback value is 0 rather than falling back to infinity.

Frankly I would prefer if there were no NaNs at all. E.g. consider this computation:

sqrt(1000.3 - 1000.0 - 0.3)

With 64-bit IEEE 754 floats it returns NaN. I'm not sure if it's a NaN with posits, but some similar formula is going to be a NaN.

But it's not really a bug, it's just a normal rounding error. If you're computing square roots at 0 with some rounding errors earlier in the computation, you expect some of them to involve negative numbers because of rounding. The proper way to make this computation numerically stable is to compute sqrt(max(x, 0)).

In other worlds, the function sqrt(max(x, 0)) is more useful than the function "if(x>=0) then sqrt(x) else NaN/infinity".

[MitchAlsup]

In an algebraic calculation, SQRT(-v) = COMPLEX( 0.0, SQRT(+value) )

And probably not something you want a function unit to do to your register file.

But in a real (sic) sense I agree with the notion of 0.0 rather than Infinity.

However, I bet Kahan has 17+ reasons that something other than 0.0 is a good choice, here, too.

Infinity has the property that "most SW won't deal with Infinity as gently as it deals with 0.0."

quote: " But it's not really a bug, it's just a normal rounding error. If you're computing square roots at 0 with some rounding errors earlier in the computation, you expect some of them to involve negative numbers because of rounding. "

The bug, here, is not in the definition of SQRT, but in the calculations leading to SQRT.

Do not try to fix SQRT to handle bugs elsewhere in your application.

[John Gustafson]

This has gotten to the point where I really have to jump in. Some quick points.

Defining the square root of "negative zero" as "negative zero" is one of the most hilarious and indefensible clauses in the IEEE Std 754™ document. As I said in my book, the only explanation was that it was the 1980s, it was Berkeley, and the use of controlled substances may have been involved. The rules for "negative zero" sometimes treat it as "negative epsilon", sometimes "tests equal positive zero" but if you take the reciprocal of negative zero you get –∞, infinitely distinguishable from the reciprocal of positive zero, +∞.

Unum Type I arithmetic made an attempt to be upward-compatible with IEEE Std 754, but drew the line at having trillions of ways to express "NaN" and supporting the nonsensical attempts to impart meaning to "negative zero" and "positive zero." Unum Type II and Type III arithmetic discard upward compatibility with IEEE Std 754 and are a clean-slate design. I think I'll attach the latest version of the Posit Standard so people can start taking a look at how much simpler it is (12 pages) compared to IEEE Std 754 (~100 pages).

There is widespread misunderstanding of how to construct math libraries that are correctly rounded for every input argument, that you have to use the method of Ziv and it can be unbelievably and unpredictably expensive. You don't. All you have to do is find an approximate that rounds the same way the mathematical function rounds, and that is a much, much less demanding goal. That's the Minefield Method. Santosh Nagarakotte and Jay Lim and others at Rutgers have automated the method brilliantly, and except for the problem of x^y for 32-bit x and y, I'm comfortable insisting that transcendental functions (and every other math library function) up to 32-bit inputs must return a correctly-rounded output.

Intel's Math Kernel Library (MKL) has thousands of cases that are incorrectly rounded, leading to irreproducibility since the gcc math library also has thousands of cases that are incorrectly rounded, but not the same cases. Now, get this: Nagarakotte and Lim applied the Minefield Method to the usual 32-bit trig, exponential, and transcendental functions. producing perfect rounding for every input, no exceptions. And it's faster than MKL. You heard me. Perfect rounding is faster than Intel's MKL or gcc, if you use the insight that correct rounding need not require a hyper-accurate approximation for the rare hard-to-round point. I just has to pass through the ULP-wide gap that leads it to round correctly.

John

P.S. I'm sending the attached version out for ratification by the Posit Working Group. We should have this finalized in a few days, unless someone is a stubborn holdout!

[John Gustafson]

More importantly, do not ask hardware to protect you from software bugs involving the square root. Fix the software in the rare cases where a negative zero might arise, and let the hardware be as simple and power-efficient and small chip area as possible.

Too much of IEEE Std 754 was based on building "code babysitters" into hardware. Most vendors have stripped those out for efficiency, and only support the sillier rules in 754 with software traps. People discover that those rules cause lower performance, and soon find the compiler flag that switches them off.

I meant to explain yesterday… the reason you cannot define 0^0 as 1 is that you have no idea what real number was each input, since 754 floats round small magnitude values to 0 (creating a 100% relative error). The base could have been positive, zero, or negative. The exponent could have been positive, zero, or negative. The possible set of real numbers that could be produced by an exact calculation is huge; therefore the output is indeterminate.

The general principle is simple: Imagine an open neighborhood in the complex plane containing the input values. Treat those neighborhoods as the inputs, and look at the output set. Can it be made arbitrarily close to a single real value with sufficiently small input neighborhoods? If it can, then you can take that real number as a well-defined limit and round it to the nearest representable quantity. That is the rule used in the Posit Standard. If you cannot get an arbitrarily small neighborhood containing a single real value no matter how small the neighborhoods of the inputs are, then you must declare the result indeterminate, or Not-a-Real (NaR).

Similarly, atan2(0, 0) is indeterminate. Every argument from –π to π is possible in any neighborhood of the inputs.

There is a human urge to want an exact answer and reject indeterminacy. That's an urge best resisted. If you willfully kid yourself that a result is deterministic when it isn't, the truth is likely to cause catastrophic surprises in program execution.

John

-- 
You received this message because you are subscribed to the Google Groups "Unum Computing" group.
To unsubscribe from this group and stop receiving emails from it, send an email to unum-computin...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/unum-computing/6d55b203-1667-4901-8e47-f1e8d0d47d2dn%40googlegroups.com.

[Tomek Czajka]

I feel there is an inconsistency.

In the case of sqrt you reject "babysitting" for rounding errors and expect the programmer to protect him/herself, while in the 0^0 you reject the calculation for "babysitting" reasons, to protect the programmer from rounding errors.

Perhaps the exponent is *not* rounded in some use cases, in which case the babysitting rule prevents the programmer from getting useful work done? This may actually often be the case: people often use integral exponents, which are not rounded but perfectly accurate.

And of course this rule is inconsistent with something like ceiling(0.0) which you could similarly argue should return NaR, because you never know if the exact value was slightly above or slightly below 0.

-Tomek

[Tomek Czajka]

As a constructive alternative, it may make sense to leave pow(0, 0) and pow(-1, 1), etc, as NaR, but in that case it seems useful to provide pown(x, n) that works for any integer n, when either n >= 0 or x > 0.

[MitchAlsup]

I like the clarity and brevity of this rule. Well done.

[MitchAlsup]

On Monday, January 10, 2022 at 3:52:39 PM UTC-6 tcz...@gmail.com wrote:

I feel there is an inconsistency.

In the case of sqrt you reject "babysitting" for rounding errors and expect the programmer to protect him/herself, while in the 0^0 you reject the calculation for "babysitting" reasons, to protect the programmer from rounding errors.

 

I do not see the inconsistency. In both cases one needs to fix the problem before getting to the F(~0) or F(~0,~0) 

{Where tilde ~ means rounded to actual}

Perhaps the exponent is *not* rounded in some use cases, in which case the babysitting rule prevents the programmer from getting useful work done?

The work the programmer should be doing is finding the sources of ~0s getting to the arguments of functions for which actual 0.0 is out of domain.

 

This may actually often be the case: people often use integral exponents, which are not rounded but perfectly accurate.

And of course this rule is inconsistent with something like ceiling(0.0) which you could similarly argue should return NaR, because you never know if the exact value was slightly above or slightly below 0.

But ceil() works across the entire Real-domain.

Calculations paths leading to ceil which get improperly rounded can change the result of ceil even though ceil is obeying all the numeric properties it was designed to obey.

 

-Tomek

[John Gustafson]

The neighborhood argument does not apply to the following functions in the Posit Standard:

negate

abs

sign

nearestInt

ceil

floor

next

prior

See Section 5.1.

If you want (0.0)^0 to evaluate to exact 1.0, call

compound(–1.0, 0)

and it will return 1 since the second input argument is an integer (Section 5.8).

Also, the examples I saw where 0^0 occurs seem to be as part of a series. If you initialize a value x = 1.0 and then multiply it by itself in a loop, you certainly will get 1.0 if you execute the loop zero times. Then x, x², x³, etc. That accumulates rounding errors after each multiplication, however.

John

--
You received this message because you are subscribed to the Google Groups "Unum Computing" group.
To unsubscribe from this group and stop receiving emails from it, send an email to unum-computin...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/unum-computing/3d465e26-cee8-45b5-8569-75a036d1c611n%40googlegroups.com.