# How would I use Supposition.jl on computer arithmetic corner cases?

**URL:** https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418
**Category:** General Usage
**Tags:** question
**Created:** [September 15, 2024, 1:42am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418 "2024-09-15T01:42:14Z")
**Posts on this page:** 19
**Page:** 1

<div class="post-metadata">

### Author: ![kapple](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kapple/32/218915_2.png) [@kapple](https://discourse.julialang.org/u/kapple)
#### Post date: [September 15, 2024, 1:42am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/1 "2024-09-15T01:42:14Z")

</div>

How would one implement a Supposition.jl test for something like the following MWE?

```julia
@check function phased(x = Data.Floats())
  cosd(x) == sind(x + 90)
end

```

which fails due to floating point stuff

```julia
x = 1.8014398509481896e16
(cosd(x), sind(x + 90))

```

which returns two different numbers.

Just like how `@test_broken` is an admission that a functionality is currently not behaving as desired,  
I’d imagine maybe splitting into two tests?

#1. An expected-passing `@check` test on a range of `Supposition.Data.Floats` (using `Supposition.Data.Satisfying`?) - different range for each concrete floating point numeric type.  
#2. An admission of broken functionality, somewhat a `@check_broken` on the floating point range complement of #1.

It’d be interesting to run such tests on Julia Base, and packages like SpecialFunctions.jl which can build awareness for people to know the specific limitations of the numeric types they’re implementing. And for their use cases, they might consider doing things like [what Milan did](https://www.youtube.com/watch?v=qba-7kdXiJg) for SpeedyWeather.jl (in his case, mapping the important dynamics to a custom numeric type containing a finer collection of values at the important ranges for his modelling).

I’m going to be annoying and tag @Sukera.

Edit: The test also fails with `isapprox` instead of `==`, thanks for bringing that up @jar1.

---

<div class="post-metadata">

### Author: ![jar1](https://avatars.discourse-cdn.com/v4/letter/j/c0e974/32.png) [@jar1](https://discourse.julialang.org/u/jar1)
#### Post date: [September 15, 2024, 1:51am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/2 "2024-09-15T01:51:58Z")

</div>

You might check out [`isapprox`](https://docs.julialang.org/en/v1/base/math/#Base.isapprox) and use that instead of `==` for floating-point comparisons.

---

<div class="post-metadata">

### Author: ![kapple](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kapple/32/218915_2.png) [@kapple](https://discourse.julialang.org/u/kapple)
#### Post date: [September 15, 2024, 1:54am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/3 "2024-09-15T01:54:58Z")

</div>

I forgot to say that it actually fails `isapprox` as well.

```julia
@check function phased(x = Data.Floats())
  isapprox(cosd(x), sind(x + 90))
end

```

```julia
Found counterexample
  Context: phased
  Arguments:
      x::Float64 = 1.8014398509481896e16

```

My original post before editing had `isapprox` in it because I tested with `==` first, then with `isapprox`, both of which failed, then I copy-pasted to here.

---

<div class="post-metadata">

### Author: ![kapple](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kapple/32/218915_2.png) [@kapple](https://discourse.julialang.org/u/kapple)
#### Post date: [September 15, 2024, 2:20am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/5 "2024-09-15T02:20:40Z")

</div>

So, more useful would be error plots I suppose.

Though from the context of, I’m a package developer, I have implemented some math functions, and I’d like to verify them, so I’ll use Supposition.jl as the closest thing to formal verification that we have, and implement them as part of my tests. In _that_ context, I’m wondering how to go about it.

And in the documentation, the error plots can be provided, with reference to the Supposition.jl testing results.

I would think this important for say, e.g. inclusion as uncertainty bounds in implemented numerical models.

---

<div class="post-metadata">

### Author: ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)
#### Post date: [September 15, 2024, 6:24am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/6 "2024-09-15T06:24:08Z")

</div>

This is an interesting question, thank you for bringing it up! I think the most relevant section of the docs here is [Alignment of Documentation · Supposition.jl Documentation](https://seelengrab.github.io/Supposition.jl/stable/Examples/docalignment.html), since it talks about these kinds of issues a bit in the context of floating point values in particular. Documenting these kinds of behaviors is IMO the best course of action.

Unlike with the precision of `sincosd`, I don’t think there’s any guarantee about the precision of the basic `sincosd` relative to `cosd` and how it should behave when shifting the phase of the input by 90°. So whether this inexactness is considered “broken” or not is really not that easy to answer! What I would do is basically the same as what @jar1 suggested - use `isapprox` with an absolute tolerance and check that the output is always within that tolerance. For a single operation like `sind` on a `Float64`, this is typically [1 ULP](https://en.wikipedia.org/wiki/Unit_in_the_last_place), and the result is compared to `sind` on a `BigFloat`. For your example though, you’d still have to compare the `Float64` results of `sind` & `cosd` with `isapprox`.

Note also that it’s going to be _extremely_ difficult to guarantee a large number of trigonometric identities to hold, even within some chosen precision. See e.g. [this floating point favorite of mine](https://discourse.julialang.org/t/array-ordering-and-naive-summation/1929) for an example where it’s hard.

* * *

All that being said - your intuition of using error plots as a guiding principle is spot on! In this example, I’d try to plot the error over the range you’re interested in and just take a look at how bad it really is. Since this is about trigonometry, I’d forego the `*d` versions and use `sin`/`cos` directly, producing a plot of all `Float16` numbers in [0, 1].

Maybe the overall error is already acceptable, maybe the worst case is MUCH worse than the example Supposition.jl found. Remember, it tries to find the minimal _input_ that produces the same test result, which may not necessarily be the minimal _error_ that you’re actually interested in! You could try to find that with `target!(abs(sind(x+90) - cosd(x)))` (see [this doc section](https://seelengrab.github.io/Supposition.jl/stable/Examples/target.html)), which looks like it could be simple enough to work out.

---

<div class="post-metadata">

### Author: ![kapple](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kapple/32/218915_2.png) [@kapple](https://discourse.julialang.org/u/kapple)
#### Post date: [September 15, 2024, 6:36am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/7 "2024-09-15T06:36:29Z")

</div>

Speaking hypothetically here, suppose there exists a package, let’s call it ComputationalErrors.jl (or maybe even just include the following in Supposition.jl). It has a macro, say `@uncertainty`. Analogous to `@btime` which outputs a measure of performance, `@uncertainty` would output a measure of uncertainty, and would be used with similar syntax as `@check`. It could then be used in package authors’ documentations, being calculated autonomously on each doc build. I feel like that would be in the spirit of what Supposition.jl is pointing toward. What is everyone’s thoughts on this idea?

It could be used in Julia’s documentation for math functions, on SpecialFunctions.jl, on automating the difference of accuracy calculations for different architectures/parallelisations/etc. It can inform on what uncertainty bounds to use as inputs when using a package’s functions in one’s own models. Is that something that people would find useful? I would personally find this more acceptable than manually trying to find an `isapprox` tolerance for things that should theoretically ideally be `isequal`.

Edit: Stefan’s post on that summation behaviour is crazy!

---

<div class="post-metadata">

### Author: ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)
#### Post date: [September 15, 2024, 7:19am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/8 "2024-09-15T07:19:27Z")

</div>

It’s not a bad idea, but I don’t think it’s a good fit for Supposition.jl. You can probably achieve this exact thing by using IntervalArithmetic.jl and checking that the true upper & lower bounds returned are within your margin of error, without having to add anything to Supposition.jl.

---

<div class="post-metadata">

### Author: ![mleseach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleseach/32/202840_2.png) [@mleseach](https://discourse.julialang.org/u/mleseach)
#### Post date: [September 15, 2024, 9:05am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/9 "2024-09-15T09:05:47Z")

</div>

I think in this case, you might want to either restrict the domain or check for relative precision with respect to x.

```julia
@check function phased(x = Data.Floats())
    if !isfinite(x) || !isfinite(x + 90)
        return true
    end

    isapprox(cosd(x), sind(x + 90); atol=eps(x)/2)
end

```

Here, `eps(x)/2` represents the rounding error. For functions with larger variations, you might want to use Taylor expansions to approximate the error. For example:

```julia
# First order approximation
atol = C * abs(df(x) * eps(x)/2)

# Second order approximation (useful if f' can be zero)
atol = C * abs(df(x) * eps(x)/2 + ddf(x) * (big(eps(x))/2)^2 / 2)

# You can extend this further for higher orders

```

Where `C` is a constant that you choose to adjust the tolerance level. In the second-order case, I use `BigFloat` because `eps(x)^2` can result in infinity (`Inf`) in floating-point arithmetic for large values of x.

---

<div class="post-metadata">

### Author: ![kapple](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kapple/32/218915_2.png) [@kapple](https://discourse.julialang.org/u/kapple)
#### Post date: [September 15, 2024, 9:15am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/10 "2024-09-15T09:15:43Z")

</div>

Oh it didn’t cross my mind to use a Taylor expansion in defining a tolerance, that’s a good idea. Would work well with autodiff utilities as well.

And I guess that would make sense if I wanted a more flexible tolerance system for my property tests. But I’m thinking of the scenario: I pick up Julia for the first time, I choose a package which implements some mathematics, and I would use said package with the expectation that the inplementations would hold to mathematical properties. So a package author could get in the mindset of giving their users the awareness: “this is what my package provides, and these are its specific limitations based on how it holds/doesn’t hold mathematical relationships.” For this perspective, I’m not looking to pass tests. I’m looking to measure deviation from objective truth.

Regarding using IntervalArithmetic, that would produce a range of outputs based on a range of inputs. But my thinking is of evaluating an incorrectness of expected properties based on single inputs. So we’re back to the error graph.

(On that note, sorry @Sukera for the premature “Solution” labelling and for un-solutioning your response, it took me a moment to process my thoughts.)

---

<div class="post-metadata">

### Author: ![mleseach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleseach/32/202840_2.png) [@mleseach](https://discourse.julialang.org/u/mleseach)
#### Post date: [September 15, 2024, 9:16am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/11 "2024-09-15T09:16:46Z")

</div>

I was going to suggest implementing something like that by hand (I wasn’t aware of this package), but I was unsure how to handle non-monotonic functions.

---

<div class="post-metadata">

### Author: ![kapple](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kapple/32/218915_2.png) [@kapple](https://discourse.julialang.org/u/kapple)
#### Post date: [September 15, 2024, 9:19am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/12 "2024-09-15T09:19:54Z")

</div>

@mleseach are you talking about the Taylor expansion method for deciding a tolerance? I would guess that the Taylor expansion would be centered at the value the test is running on, so it should work for non-monotonic functions.

What I would understand it wouldn’t work on is functions that are not C1, i.e. continuous and differentiable.

Edit: I just noticed you were replying to Sukera who was talking about using IntervalArithmetic, nevermind.

---

<div class="post-metadata">

### Author: ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)
#### Post date: [September 16, 2024, 6:55am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/13 "2024-09-16T06:55:03Z")

</div>

I think in this case the problem lies not with Supposition.jl but rather floating point arithmetic tbh. The property you test, i.e. `cosd(x) == sind(x + 90)` indeed does not hold for _all_ floats. However, this is due to finite precision, because for very large `x` you have `x+90 == x` (which Supposition.jl correctly identifies!). So I think you should rather try to fix your test to be more meaningful.  
I suggest splitting it into two parts:

- You test that `cosd(x) == sind(x + 90)` for some reasonable range, like `-720 < x < 720`
- You test that it correctly wraps the input, i.e. `cosd(x) = cosd(x % 360)`  
Together these tests essentially cover the same the ground as before but you avoid the trouble with floating point additions.

Perhaps you could combine these two cases like so:

```julia
@check function phased(x = Data.Floats())
  x = x % 720 # or some other reduction to a reasonable interval
  cosd(x) == sind(x + 90)
end

```

But I think leaving the cases separate makes is clearer and actually tests a bit more thorough.

---

<div class="post-metadata">

### Author: ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)
#### Post date: [September 16, 2024, 7:15am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/14 "2024-09-16T07:15:28Z")

</div>

Specifically the `x` value in the original post is an integer valued `Float64` which is large enough that `x + 90` cannot be represented exactly as a `Float64`.

```julia
julia> x = 1.8014398509481896e16
1.8014398509481896e16

julia> Int(x + 90) - Int(x)
88

julia> (cosd(Int(x)), sind(Int(x) + 90))
(0.9135454576426009, 0.9135454576426009)

```

---

<div class="post-metadata">

### Author: ![kapple](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kapple/32/218915_2.png) [@kapple](https://discourse.julialang.org/u/kapple)
#### Post date: [September 17, 2024, 3:38am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/15 "2024-09-17T03:38:20Z")

</div>

I like the point and demonstration of how my test can be improved. Splitting into two parts is actually a good idea, thank you!

I think the point I’m trying to make in this thread is the fact that Julia’s `cosd` doesn’t `mod` the input, and therefore, what’s the appropriate caveat of usage informing users on the potential errors? How large can `x` be before we expect inequality? And beyond that, how large can `x` be before, say, `acosd(cosd(x))` starts behaving in quantifiably undesirable ways?

That said, I am still starting to reach the same conclusion that Supposition.jl isn’t scoped for the presentation I am looking for. It is, however, very good for bringing awareness of the corner cases/corner domain regions, and for verifying that the presentation is correct.

Splitting into cases will be applied to, as I’ve listed above, the domain of inputs where a property exactly holds, the domain of inputs where a property approximately holds and present the error (plot or otherwise), and the domain of inputs where you get mangled garbage returned.

---

<div class="post-metadata">

### Author: ![kapple](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kapple/32/218915_2.png) [@kapple](https://discourse.julialang.org/u/kapple)
#### Post date: [September 17, 2024, 3:48am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/16 "2024-09-17T03:48:26Z")

</div>

Oh I see, that’s a good point I hadn’t realised.

So applying your observation to the point I’m trying to convey and ask about, in order for a user to get accurate output, they would need to:

1. Know about floating point results and arithmetic specific to the function(s) they’re using.
2. Implement appropriate workarounds to achieve the accuracy they want. In the example of `cosd` as you demonstrated: the internal wrapping with `Int` to their input to the trig function.

The `cosd` example isn’t the best example because I’d imagine pretty much no one works with inputs with magnitudes that large.

But there are other cases where such gotchas have been asked and discussed.

We give the generalised caveat on floating point arithmetic, but it’d be nice to have a library of demonstrations and numerical quantifications on the errors of floating point arithmetic in the context of the many ubiquitous math functions that get called throughout the ecosystem.

---

<div class="post-metadata">

### Author: ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [September 17, 2024, 8:25am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/17 "2024-09-17T08:25:42Z")

</div>

> [@kapple](#):
>
> the fact that Julia’s `cosd` doesn’t `mod` the input

It does, though.

> <https://github.com/JuliaLang/julia/blob/02549d5c54ec8ac8c7e62ea470803350bb3d1899/base/special/trig.jl#L1247-L1263>

---

<div class="post-metadata">

### Author: ![kapple](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kapple/32/218915_2.png) [@kapple](https://discourse.julialang.org/u/kapple)
#### Post date: [September 17, 2024, 8:37am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/18 "2024-09-17T08:37:13Z")

</div>

Oh so this issue is attributed more to the floating point sum in the test as pointed out by others above, rather than the way cosd parses the input. Okay cool, thank you.

---

<div class="post-metadata">

### Author: ![jar1](https://avatars.discourse-cdn.com/v4/letter/j/c0e974/32.png) [@jar1](https://discourse.julialang.org/u/jar1)
#### Post date: [September 19, 2024, 5:15am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/19 "2024-09-19T05:15:15Z")

</div>

Relevant discussion in [History of `isapprox` in languages - #26 by Richard\_Fateman](https://discourse.julialang.org/t/history-of-isapprox-in-languages/119347/26) and following.

---

<div class="post-metadata">

### Author: ![kapple](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kapple/32/218915_2.png) [@kapple](https://discourse.julialang.org/u/kapple)
#### Post date: [September 19, 2024, 5:23am UTC](https://discourse.julialang.org/t/how-would-i-use-supposition-jl-on-computer-arithmetic-corner-cases/119418/20 "2024-09-19T05:23:13Z")

</div>

Thanks heaps! I’ve actually been following that very interesting discussion and thinking about its relevance to what I’m thinking about in this thread. Thank you for the cross reference 😁
