Pi raised to negative integer powers produces a strange error message, and the literal_pow obfuscation

General Usage
1

I am trying to raise pi to negative integer powers, and end up with unexpected errors:

julia> pi^-2 # works
0.10132118364233779

julia> n=2; pi^-n # doesn't work
ERROR: DomainError with -2:
Cannot raise an integer x to a negative power -2.
Convert input to float.
Stacktrace:
 [1] throw_domerr_powbysq(::Float64, ::Int64) at ./intfuncs.jl:190
 [2] power_by_squaring(::Irrational{:π}, ::Int64) at ./intfuncs.jl:214
 [3] ^(::Irrational{:π}, ::Int64) at ./intfuncs.jl:239
 [4] top-level scope at REPL[18]:1

The error message is confusing as pi is firstly not an integer, and the first line in Base.power_by_squaring is a call to Base.to_power_type that converts pi to a float.

julia> Base.to_power_type(pi)
3.141592653589793

The error message seems to originate from

@noinline throw_domerr_powbysq(::Any, p) = throw(DomainError(p,
    string("Cannot raise an integer x to a negative power ", p, '.',
           "\nConvert input to float.")))

which is a catch-all method, but the error message assumes that the base is an integer.

As I understand, Julia lowers literal integer exponentiation to literal_pow, while integer variables are treated using Base.power_by_squaring. The matter isn’t helped by the fact that @code_lowered doesn’t show you what’s happening.

julia> @code_lowered pi^-2
CodeInfo(
1 ─ %1 = Base.power_by_squaring(x, p)
└──      return %1
)

julia> n=2; @code_lowered pi^-n
CodeInfo(
1 ─ %1 = Base.power_by_squaring(x, p)
└──      return %1
)

julia> Base.power_by_squaring(pi,-n)
ERROR: DomainError with -2:
Cannot raise an integer x to a negative power -2.
Convert input to float.
Stacktrace:
 [1] throw_domerr_powbysq(::Float64, ::Int64) at ./intfuncs.jl:173
 [2] power_by_squaring(::Irrational{:π}, ::Int64) at ./intfuncs.jl:197
 [3] top-level scope at REPL[6]:1

The lowered code is identical for a literal integer and a variable. How to raise pi to a negative integer power without actually evaluating inv(pi^n)? intfuncs.jl says

# don't use the inv(x) transformation here since float^p is slightly more accurate

and I would like the accuracy to be preserved.

1 Like
2

Since pi is an abstract representation of the constant \pi, you need to choose a level of accuracy you want to work with in other representations (eg floating point). Some operations, eg 1/pi, pick that for you implicitly, but for your expression you might just use something like

julia> Float64(pi)^-2
0.10132118364233778

You may be interested in the discussion for this feature, in particular this comment

3

Thanks, but it seems to work differently with \euler

julia> ℯ
ℯ = 2.7182818284590...

julia> n=2; ℯ^-n
0.1353352832366127

The reason behind this appears to be that \euler has a special method dealing with exponentiation, and this simply calls the exponential function exp(n)

julia> @which ℯ^-n
^(::Irrational{:ℯ}, x::Integer) in Base.MathConstants at mathconstants.jl:91

julia> @which π^-n
^(x::Number, p::Integer) in Base at intfuncs.jl:222

This, however, introduces an arbitrariness in how irrational numbers behave. Is it worth defining a special exponentiation method for pi? Although looking at possible options

julia> exp(-2*log(pi))
0.10132118364233778

julia> inv(pi)^2
0.10132118364233779

It’s not clear which one is more accurate in general. One way to avoid this is to do what exp(::Real) does, by converting its argument to float. In this case we obtain

julia> pi^-2.0
0.10132118364233778

In the course of the evaluation pi is also promoted to Float64, so the exponentiation is effectively calling ^(::Float64,::Float64). This will likely be a breaking change, as

julia> float(pi)^float(5.0)
306.0196847852814

julia> pi*pi*pi*pi*pi
306.01968478528136

I can see the ambiguity here, but the present situation is pretty confusing. At the very least we need a better error message.

1 Like
4

I think this can be considered a missing method issue, which could be solved by defining a fallback method

Base.:^(x::Irrational, y) = Float64(x)^y

which would fix this for all other constants (note that I did not check method ambiguities etc). Cf

julia> n = 2; MathConstants.catalan^(-n)
ERROR: DomainError with -2:
Cannot raise an integer x to a negative power -2.
Convert input to float.

I would not worry about accuracy in the first iteration; general methods (eg power by squaring) should be fine. Having a missing method is an issue that should be resolved, please consider making a PR.

2 Likes