Add exact rules for irrational numbers

[hyrodium]

This PR fixes https://github.com/JuliaMath/IrrationalConstants.jl/issues/20, and is an alternative to JuliaMath/Tau.jl#7.

On current master

julia> sqrt2^2
2.0000000000000004

julia> sin(quartπ)
0.7071067811865475

julia> inv(quartπ)
1.2732395447351628

With this PR

julia> sqrt2^2
2.0

julia> sin(quartπ)
invsqrt2 = 0.7071067811865...

julia> inv(quartπ)
fourinvπ = 1.2732395447351...

[coveralls]

Pull Request Test Coverage Report for Build 2239524464

Warning: This coverage report may be inaccurate.

This pull request's base commit is no longer the HEAD commit of its target branch. This means it includes changes from outside the original pull request, including, potentially, unrelated coverage changes.

• For more information on this, see Tracking coverage changes with pull request builds.
• To avoid this issue with future PRs, see these Recommended CI Configurations.
• For a quick fix, rebase this PR at GitHub. Your next report should be accurate.

Details

• 26 of 27 (96.3%) changed or added relevant lines in 1 file are covered.
• No unchanged relevant lines lost coverage.
• Overall coverage decreased (-3.7%) to 96.296%

---

Changes Missing Coverage	Covered Lines	Changed/Added Lines	%
src/rules.jl	26	27	96.3%
<!--	Total:	26	27

Totals
Change from base Build 1333249280:	-3.7%
Covered Lines:	26
Relevant Lines:	27

---

💛 - Coveralls

[codecov[bot]]

Codecov Report

Merging #14 (f1e6c86) into main (d226b95) will increase coverage by 100.00%. The diff coverage is 100.00%.

@@             Coverage Diff             @@
##           main       JuliaMath/Tau.jl#14        +/-   ##
===========================================
+ Coverage      0   100.00%   +100.00%     
===========================================
  Files         0         1         +1     
  Lines         0        28        +28     
===========================================
+ Hits          0        28        +28     

Impacted Files	Coverage Δ
src/rules.jl	100.00% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update d226b95...f1e6c86. Read the comment docs.

[devmotion]

I'm in favour of adding such improvements in principle. However, I think we should be consistent with Julia base and non-breaking, ie all shortcuts should return values of the same type as currently. Thus generally they should return Float64 but no irrationals.

[hyrodium]

Thus generally they should return Float64 but no irrationals.

For example, log(ℯ) is not Float64 but Int.

julia> log(ℯ)
1

non-breaking, ie all shortcuts should return values of the same type as currently.

This package is before v1.0, so I thought it would be okay to add breaking changes.

[devmotion]

It would be OK to add breaking changes here but IMO it's not good to be inconsistent with Julia base. In the sin(pi) PR triage made it clear that the trigonometric shortcuts should return Float64, and hence we should do the same. The definition for log(e) is somewhat special and older, and deviates from the common pattern in Julia base it seems.

[hyrodium]

In the sin(pi) PR triage made it clear that the trigonometric shortcuts should return Float64, and hence we should do the same.

In that PR(https://github.com/JuliaLang/julia/pull/42595), the discussion 1.0 vs 1 is not about accuracy, but we need to care about it in this PR.

On the current master

julia> sin(quartπ)^2
0.4999999999999999

With this PR

julia> sin(quartπ)^2
0.5

I have read the following comment in https://github.com/JuliaLang/julia/pull/35823#issuecomment-630367743.

The main benefit I see in returning integers is that it avoids unnecessary widening types. Returning Float64 for all integer input types doesn't seem that useful as a feature, but in generic code this could lead to Float32s being converted to Float64 accidentally. (...) I don't think it's unreasonable for users to expect <:Integer as output type if the output of a function can only ever be a whole number. My main worry about this change is that it would be a breaking change in some cases.

We don't need to care about breaking changes here, so I still think the return values can be Irrational. And also, there was no discussion about the return types in JuliaMath/Tau.jl#7.

@jishnub @JeffreySarnoff Do you have any thoughts on this?

[devmotion]

In that PR(JuliaLang/julia#42595), the discussion 1.0 vs 1 is not about accuracy, but we need to care about it in this PR.

I don't see why we can't improve accuracy while keeping the same return types as base (and as currently in this package). Using

sin(::typeof(quartpi)) = Float64(invsqrt2)

would keep consistent return types while speeding up calculations and potentially improving accuracy (at least it's definitely not worse).

Such changes would be much less controversial and easily justifiable. Hence even if we would want to change the behaviour at some point, they would be a simple non-breaking first step.

IMO it would be very unfortunate to introduce inconsistencies such as

julia> sin(pi) isa Float64
true

julia> sin(quartpi) isa Irrational
true

Apart from the inconsistencies that we would introduce, I also wonder more generally wether it's worth returning Irrationals. As @JeffreySarnoff pointed out in the other PR it is unclear where to stop (and as the commented out line in this PR shows some definitions are not even possible without introducing type piracy). In most downstream calculations at some point irrationals are promoted to floating point numbers anyway, and hence downstream packages already now have to take care of these promotions (and e.g. usually avoid multiplying them with integers) to ensure that the output is of the desired type and no undesired promotions occur. Hence loss of accuracy usually should not be an issue. And if your computation contains a term such as x * sin(quartpi)^2 it seems straightforward enough to just replace it with x/2 instead of relying on Julia and/or compiler optimizations.

[JeffreySarnoff]

Mathematically, elementary functions almost never return an integer value. To introduce a second return type would force other people's code to work with Unions as returned types (whether they know it or not), and this is disadvantageous.

[hyrodium]

Thanks for the detailed description! I still have some questions.

sin(::typeof(quartpi)) = Float64(invsqrt2) would keep consistent return types while speeding up calculations and potentially improving accuracy

As you may know, if we use this definition, then the sin(quartπ)^2 will not be precise.

julia> Float64(invsqrt2)^2
0.5000000000000001

And the speed is not relevant here:

julia> f(x) = sin(x)^2  # with this PR
f (generic function with 1 method)

julia> @code_llvm f(quartπ)
;  @ REPL[27]:1 within `f`
define double @julia_f_292() #0 {
top:
  ret double 5.000000e-01

julia> f(x) = sin(x)^2  # On the current master
f (generic function with 1 method)

julia> @code_llvm f(quartπ)
;  @ REPL[2]:1 within `f`
define double @julia_f_147() #0 {
top:
  ret double 0x3FDFFFFFFFFFFFFE
}

IMO it would be very unfortunate to introduce inconsistencies such as

Both sin(pi) and sin(quartπ) are Real, and there is no type-instability, so I think this is okay. Note that both π and quartπ are Irrational, but they have different concrete types.

julia> typeof(π) == typeof(quartπ)
false

As a different example, the power operator ^ returns different types.

julia> 2^2
4

julia> 2^(-2)
0.25

Both inputs are integer, but the output types are different. This is not problematic because the type-instability doesn't exist after compiling. Irrational numbers can be handled in compile-time as well.

I also wonder more generally wether it's worth returning Irrationals. As @JeffreySarnoff pointed out in the other PR it is unclear where to stop

I was also having the same concerns. I thought "return Irrational if we have the irrational, otherwise, return Float64" is the best here.

And if your computation contains a term such as x * sin(quartpi)^2 it seems straightforward enough to just replace it with x/2 instead of relying on Julia and/or compiler optimizations.

I agree with that, but I think the objection can be adapted to sin(π).

To introduce a second return type would force other people's code to work with Unions as returned types (whether they know it or not), and this is disadvantageous.

Is there any practical case to produce Union? For example, inv.([invπ, 3]) is not a Vector{Union{Irrational{:π}, Float64}} but a Vector{Float64}.

[JeffreySarnoff]

I meant, if we evaluated sinpi(1.0) === 0 (as an Int) while other arg values returned a Float64, then code like this would have a return type of Union{Int, Float64} for compilation purposes.

function a(x)
   sinpi(x)
end

It is much the same case as this

julia> function fn(x::Int)
            iseven(x) ? 1 : 1.0
          end

julia> Base.return_type(fn, (Int64,))
1-element Vector{Any}:
 Union{Float64, Int64}

[hyrodium]

It is much the same case as this

I see, but the function fn(x::Int) is not type-stable. My proposal such as sin(quartπ) isa Irrational is type-stable, so it doesn't seem disadvantageous.

[jishnub]

IIUC, the issue is with something like

julia> sin.(Irrational[π, quartπ])
2-element Vector{Real}:
        1.2246467991473532e-16
 invsqrt2 = 0.7071067811865...

vs. something like

ulia> sin.(Irrational[π, ℯ])
2-element Vector{Float64}:
 1.2246467991473532e-16
 0.41078129050290885

[hyrodium]

I think the example is similar to the following, so sin.(Irrational[π, quartπ]) isa Vector{Real} seems fine for me.

julia> abs2.(Real[1.,1//2])
2-element Vector{Real}:
  1.0
 1//4

julia> abs2.([1.,1//2])
2-element Vector{Float64}:
 1.0
 0.25

If we use Irrational[...] or Real[...], then it implies the output may produce type instability.

julia> isconcretetype(eltype(Irrational[π, quartπ]))
false

julia> isconcretetype(eltype(Real[1.,1//2]))
false

julia> isconcretetype(eltype([π, quartπ]))
true

julia> isconcretetype(eltype([1.,1//2]))
true

Note that we already have the following property because of log(ℯ) isa Int.

julia> log.(Irrational[π, ℯ])
2-element Vector{Real}:
 1.1447298858494002
 1

(I still could not agree that the definition of log(ℯ) is wrong.)

[JeffreySarnoff]

@hyrodium we may be looking at slightly different effects, my concern with providing Integers when the result is exact and integral while providing floats for most results is that the following does not allow the same downstream optimizations as would be available were only floats provided. Some specializations and simplifications with the flow within other users' functions which utilize this and pass along the results to other methods certainly could become unavailable.

function sinpi_example(x)
	  result = sinpi(x)
	  isinteger(result) ? Int(result) : result
end

julia> Base.return_types(sinpi_example, (Irrational,))
1-element Vector{Any}:
 Union{Float64, Int64}

[devmotion]

Both sin(pi) and sin(quartπ) are Real, and there is no type-instability, so I think this is okay.

Surely it is type stable but my point was that it is inconsistent - with the behaviour in base but also within this PR. As pointed out by others, this is problematic since it makes it more difficult to reason about the return type, both for the compiler and downstream developers.

To reiterate, I think it would be good to avoid these problems and just consistently ensure sin(::Irrational) isa Float64. This is also non-breaking and generally I think would be much less controversial due to the existing definitions of e.g sin(pi) = 0.0 in base. And, of course, this does not rule out that the behaviour might be changed at a later point, e.g., if base decides to change its definition.

[hyrodium]

we may be looking at slightly different effects, my concern with providing Integers when the result is exact and integral while providing floats for most results is that the following does not allow the same downstream optimizations as would be available were only floats provided.

"Exact result for exact input" doesn't imply type-instable, and the sinpi_example function doesn't seem eligible here. Note that the discussion in https://github.com/JuliaLang/julia/issues/35820 and https://github.com/JuliaLang/julia/pull/35823 is about whether sinpi(1) isa Int or sinpi(1) isa Float64, but not about sinpi(1.0) isa Int or sinpi(1.0) isa Float64.

[hyrodium]

this is problematic since it makes it more difficult to reason about the return type, both for the compiler and downstream developers.

Is this comment on the compiler true? I'm not much familiar with the inside of the Julia compiler, but I think if we have type-stable code, that will not be problematic. For the downstream developers, we can just add more documentation. If someone needs Float64, one can just write float(sin(quartπ)).

And, of course, this does not rule out that the behaviour might be changed at a later point, e.g., if base decides to change its definition.

I'm okay to separate this PR into [add rules] and [fix return type], but, I couldn't agree that the Julia base enforces sin(::Irrational) isa Float64.

• log(ℯ) seems correct for me.
• I think we just need sin(::Real) isa Real, not sin(::Irrational) isa Float64.
• My proposal will be sometimes useful to avoid undesired type promotions.

julia> sin(quartπ)
invsqrt2 = 0.7071067811865...

julia> f(x,a) = sin(x) + a
f (generic function with 1 method)

julia> f(3.0f0, 4.0f0)
4.14112f0

julia> f(quartπ, 4.0f0)
4.7071066f0

In the above example, f(π,4.0f0) will be still Float64. I think this should be fixed in Base to return sin(π) === false, just like zero(Irrational). I'm understanding the reason for sinpi(0) isa Float64 and sin(π) isa Float64 is just to avoid breaking. (cf. https://github.com/JuliaLang/julia/pull/35823#issuecomment-630367743)

[JeffreySarnoff]

@hyrodium At this point in the discussion it is safe to say that your position and reasoning are well understood by the thread participants. You ask "Is this comment on the compiler true?" It is true and well understood by the people who work with performance sensitive coding. In addition, it may increase precompilation time, as either each case needs to be followed into calling code or optimization paths dropped. I understand your preference for mathematical "honesty", and respect it. Nonetheless, your proposal would lead people to avoid using our work -- including me.

[hyrodium]

@JeffreySarnoff Thanks for summarizing the threads so far!

It is true and well understood by the people who work with performance sensitive coding. In addition, it may increase precompilation time, as either each case needs to be followed into calling code or optimization paths dropped.

it makes it more difficult to reason about the return type, both for the compiler and downstream developers.

Sorry to bother you, but I have one last question. Is there any performance benchmark for that? In my local environment, the compilation time with this PR is much less than the current main branch.

on the current main branch

julia> using IrrationalConstants

julia> t=time(); inv(inv(sin(quartπ)))^2; println(time()-t)
0.04355502128601074

with this PR

julia> using IrrationalConstants

julia> t=time(); inv(inv(sin(quartπ)))^2; println(time()-t)
0.000640869140625

Note that I could not compare the performance with @time macro; both results were 0.000001 seconds.

on the current main branch

julia> using IrrationalConstants

julia> @time inv(inv(sin(quartπ)))^2
  0.000001 seconds
0.4999999999999999

with this PR

julia> using IrrationalConstants

julia> @time inv(inv(sin(quartπ)))^2
  0.000001 seconds
0.5

I'm using the following environment.

julia> versioninfo()
Julia Version 1.7.2
Commit bf53498635 (2022-02-06 15:21 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: AMD Ryzen 7 2700X Eight-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, znver1)

Sorry if the method for measuring compile time is not correct.

Since Julia uses a JIT compilation, I thought it was normal not to worry about compilation time unless the time is significantly longer.

Nonetheless, your proposal would lead people to avoid using our work -- including me.

I'm not a maintainer of this repository, so it's okay to close this PR. Then, I will create a package IrrationalConstantRules.jl to provide exact rules. (The package will pirate types just like ColorTypes.jl and ColorVectorSpace.jl.)

[devmotion]

In my local environment, the compilation time with this PR is much less than the current main branch.

As @giordano mentioned, there are some problems with the implementation in the PR, regardless of the discussion above. So it would be good to fix these before comparing run or compilation times.

One could use @time using IrrationalConstants to benchmark how long it takes to load the model (one could also start Julia with --compiled-modules=no to disable precompilation). @time includes both the runtime and compilation time when running a function for the first time. If you want to benchmark run time I strongly recommend using BenchmarkTools. Generally, I don't think one should use something like t=time(); ....

I'm not a maintainer of this repository, so it's okay to close this PR.

I would be happy if you could fix the implementation issue and change it to a non-breaking PR. But I understand if you don't want to spend time on this if you have a different goal in mind. In this case I might open a PR when I find the time for it. The main issue with a package such as IrrationalConstantRules.jl is that it will contain most/only type piracy (as you already mention) which affects even packages that only use it as a indirect dependency. This will severely limit its adoption in the Julia ecosystem, I assume. (Out of curiosity, what's the type piracy in ColorTypes?)

[hyrodium]

One could use @time using IrrationalConstants to benchmark how long it takes to load the model

Thanks! I measured the time.

On commit 1db0ac (enforce the output type with Float64)

julia> @time using IrrationalConstants
  0.029152 seconds (25.73 k allocations: 2.263 MiB, 34.32% compilation time)

On commit b0a9e4 (update testset names)

julia> @time using IrrationalConstants
  0.027626 seconds (25.73 k allocations: 2.264 MiB, 33.53% compilation time)

It seems the difference in compilation times can be ignored.

I would be happy if you could fix the implementation issue and change it to a non-breaking PR

I also updated around multiplication inverse, it now returns 1.0 instead of true.

julia> π * invπ
1.0

This behavior is different from https://github.com/JuliaMath/IrrationalConstants.jl/pull/7#issuecomment-911432902.

The main issue with a package such as IrrationalConstantRules.jl is that it will contain most/only type piracy (as you already mention) which affects even packages that only use it as a indirect dependency. This will severely limit its adoption in the Julia ecosystem

I understand that. IrrationalConstantRules.jl will be my hobby use. I'm not planning to register the package, but it can visualize people who want more exact results.

Out of curiosity, what's the type piracy in ColorTypes?

julia> using ColorTypes

julia> RGB(0.0, 0.1, 0.2)
RGB{Float64}(0.0,0.1,0.2)

julia> RGB(0.0, 0.1, 0.2) + RGB(0.4, 0.1, 0.0)
ERROR: MethodError: no method matching +(::RGB{Float64}, ::RGB{Float64})
Math on colors is deliberately undefined in ColorTypes, but see the ColorVectorSpace package.
Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...) at ~/julia/julia-1.7.2/share/julia/base/operators.jl:655
Stacktrace:
 [1] top-level scope
   @ REPL[3]:1

julia> using ColorVectorSpace

julia> RGB(0.0, 0.1, 0.2) + RGB(0.4, 0.1, 0.0)
RGB{Float64}(0.4,0.2,0.2)

The ColorTypes package defines color types such as RGB, but it doesn't define mathematical operations such as +. The +(::RGB, ::RGB) is defined in ColorVectorSpace while + is a function from Base. Strictly speaking, this can be regarded as Type II piracy.

(What I was trying to do in IrrationalConstantRules.jl is a kind of Type I piracy, so that was not good example.)

[hyrodium]

Do we really need to fix tests in Julia 1.0? (x-ref: JuliaMath/Tau.jl#15) Julia v1.0.5 doesn't even support inv(π).

julia> versioninfo()
Julia Version 1.0.5
Commit 3af96bcefc (2019-09-09 19:06 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: AMD Ryzen 7 2700X Eight-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, znver1)

julia> inv(π)
ERROR: MethodError: no method matching Irrational{:π}(::Int64)
Closest candidates are:
  Irrational{:π}(::T<:Number) where T<:Number at boot.jl:725
  Irrational{:π}() where sym at irrationals.jl:18
  Irrational{:π}(::Complex) where T<:Real at complex.jl:37
  ...
Stacktrace:
 [1] convert(::Type{Irrational{:π}}, ::Int64) at ./number.jl:7
 [2] one(::Type{Irrational{:π}}) at ./number.jl:297
 [3] one(::Irrational{:π}) at ./number.jl:298
 [4] inv(::Irrational{:π}) at ./number.jl:220
 [5] top-level scope at none:0

The test failure is because of inv(sqrt3).

[JeffreySarnoff]

I see no compelling need to support any version earlier than Long-term support (LTS) release: v1.6.6 (March 28, 2021)

[devmotion]

Julia v1.0.5 doesn't even support inv(π).

That doesn't matter for us though as we don't define it? Since it's just a specific test that fails we could just not run it on Julia versions that do not support it. I don't think there's anything problematic in the package that we can control or fix but this is just a test problem.

[hyrodium]

If we skip the tests for Julia 1.0, that means the version of Julia is not well-supported in this package, I guess. Specified code for the unsupported version of Julia will be painful for me even if that was just a test problem.... sorry.

[devmotion]

I thought one could just skip the tests that fail because some method in base is not in Julia 1.0. It's still tested properly in all other Julia versions so I don't think the package is not tested properly. However, here the fix seems to be even simpler: we can just load Compat in the tests which defines one(::Irrational) and zero(::Irrational) for older Julia versions.

That is also the fix for your example above:

julia> inv(π)
ERROR: MethodError: no method matching Irrational{:π}(::Int64)
Closest candidates are:
  Irrational{:π}(::T<:Number) where T<:Number at boot.jl:725
  Irrational{:π}() where sym at irrationals.jl:18
  Irrational{:π}(::Complex) where T<:Real at complex.jl:37
  ...
Stacktrace:
 [1] convert(::Type{Irrational{:π}}, ::Int64) at ./number.jl:7
 [2] one(::Type{Irrational{:π}}) at ./number.jl:297
 [3] one(::Irrational{:π}) at ./number.jl:298
 [4] inv(::Irrational{:π}) at ./number.jl:220
 [5] top-level scope at none:0

julia> using Compat

julia> inv(π)
0.3183098861837907

[hyrodium]

I didn't know the Compat package, thanks!

[hyrodium]

Just created IrrationalConstantRules.jl here: https://github.com/hyrodium/IrrationalConstantRules.jl

[hyrodium]

One could use @time using IrrationalConstants to benchmark how long it takes to load the model

I just noticed @time_imports macro. (https://github.com/JuliaLang/julia/pull/41612)

current main branch

julia> @time_imports using IrrationalConstants
      9.6 ms  IrrationalConstants

with this PR

julia> @time_imports using IrrationalConstants
     31.5 ms  IrrationalConstants

I think we can ignore these differences because 30ms is small.