Even powers

[andreasnoack]

I don't have the book handy, but this seems wrong wider than necessary

julia> x = Ubound34(-1,1)             
UnumX.Interval{UnumX.Unum{3,4,UInt64}}
[-1.0,1.0]                            

julia> x*x                            
UnumX.Interval{UnumX.Unum{3,4,UInt64}}
[-1.0,1.0]                            

julia> abs(x)*abs(x)                  
UnumX.Interval{UnumX.Unum{3,4,UInt64}}
[0.0,1.0]    

Maybe it is worth checking if x === y in *.

I think @alanedelman mentioned that the intervals were growing faster than expected and this might be a reason why. I wanted to run an SVD on Ubound34s, but because this issue, the norm, used in the scaling of the elementary reflectors, included negative numbers.

[simonbyrne]

I don't think checking for === would be correct (in fact, === is the same as ==). This should probably be done through ^ (which is currently unimplemented).

[JohnLGustafson]

The book deals extensively with this problem, both the specific cases of x*x (page 141) and any other routine where x appears more than once ("single-use expression" or SUE, to the interval arithmetic community). Chapter 16. If you want to use an algorithm crafted for floats, such as SVD or just about any LINPACK-type routine, you need to implement the guess function and apply it after every operation to mimic the way floats replace ranges with specific floats (rounding).

In particular, please read Chapter 17 on solving equations. Unum solutions to linear algebra problems are far superior to the ones you can obtain with floats, but don't expect to simply change a classic routines data type from float to unum and get a nice result unless you also have the unums set up to mimic rounding error (Chapter 18).

[simonbyrne]

It might be possible to do something like that via macros?

[JohnLGustafson]

Macros or some kind of environment setting, not unlike a compiler flag. My preference would be that it be shown explicitly, like "guessu(timesu(x,y))" to show a rounded value for the unum product of x and y. A simple preprocessor could apply the guess function to every unum operation, and it would make the code look messy but it would then expose all the place where there is rounding.

Then, you can try taking some of the guess calls out and see what happens. If you only guess, say, after every ten operations, it will pick the midpoint of the ubound range and I suspect that will often be a better guess than the one made by floats. Like a delayed rounding.

It also should become obvious where some loops should be replaced by fused operations. In LINPACK, for example, any call to the BLAS routines SDOT or DDOT should be replaced with the exact unum dot product. That will make a massive improvement in any routine that uses dot products as its innermost loop. Some routines use a different loop nesting that places SAXPY or DAXPY as the innermost loop, and those could use fused multiply-adds (followed by the guess function) but will not show as much benefit from unum accuracy control. If we use Level 2 BLAS and higher, it should always be possible to make an exact dot product the kernel of an array-times-array operation.