Should power_by_squaring be defined for negative powers?

I write to tell of a disappointment I had. Perhaps people can explain to me
the rationale or the right solution.

I implemented my permutation type (I am porting GAP code and none of the
existing permutation packages quite filled my requirements). I was
pleasantly surprised that, as soon as I defined the product of 2
permutations, the power p^3 was automatically defined for a permutation p.
But I was soon disappointed: p^-2 gives an error message, even though I
defined inv( p )

I tracked this to the function power_by_squaring which raises an error
automatically on negative exponents. I think that power_by_squaring should
define p^-2 as power_by_squaring(inv( p ),2).

It would cause an incompatibility since inv is defined on integers,
yielding floats, so this would define negative powers of integers. But I
think defining inv on integers is a mistake since it makes inv
type-unstable.

What do people think?

No, inv is type-stable there. Type-stability doesn’t mean that the type returned by a function is the same as the input. It means that the type can be inferred at every step of the way. If inv of an integer is always a floating point number, then that’s type stable. The problem then is that ^ on integers wouldn’t be type-stable because it could output floats or integers, even though inv and power_by_squaring are type-stable, because of the switch where you’d check for negativity in order to apply inv which does a type change.

1 Like

Ok, I see there would be a problem for Integers. It could of course be solved by a ^ method for them
which would check that the exponent is positive before calling power_by_squaring. Still, it would be nice
to have negative powers defined automatically for types which have inverses of the same type.

Can’t you just implement it for your type if it is closed under inv? Eg

struct Foo{T <: AbstractFloat}
    x::T
end

Base.:*(a::Foo, b::Foo) = Foo(a.x * b.x)

Base.inv(a::Foo) = Foo(1/a.x)   # assume this is returns the same type

Base.one(a::Foo) = Foo(one(a.x))

function Base.:^(a::Foo, p::Integer)
    if p ≥ 0
        Base.power_by_squaring(a, p)
    else
        Base.power_by_squaring(inv(a), -p)
    end
end
1 Like

Of course I did implement it as you showed. The point is that it would be nice to have it being defined
automatically. A slightly different design for Base would achieve that.

Can you provide an example? I don’t see how you could automate it, besides relying on the return type. Eg

function Base.:^(a::T, p::Integer) where T
    if p ≥ 0
        Base.power_by_squaring(a, p)
    else
        b = inv(a)
        b isa T || error("won't calculate the inverse 'cause it is not type stable")
        Base.power_by_squaring(b, -p)
    end
end

then hope that it is optimized away.

You can add a Core.Inference call (and make place Base.@pure on it if you want to assume that the inference result won’t change) and error based on the output of that.

2 Likes

@Tamas_Papp That definition looks reasonable to me (although I might abbreviate it as b = inv(a)::T). Maybe we should add this to the definition in Base?

Core.Inference is not @pure, and it makes no sense to suggest either here, since they don’t appear to have any connection to the original question (there is no use of inference necessary here, nor are any of the results known to be pure)

3 Likes

I was thinking a call to Core.Inference would get rid of calling inv(a), but now I see that @tpapp’s way would probably just compile directly to an error in the else case so sure that works without calling Inference. And I did say that there are cases where pure is wrong there…

Note that in Julia 0.7, p^-2 lowers to inv(p)^2, so this will work and be type-stable for literal integer exponents. (make x^-n equivalent to inv(x)^n for literal n by stevengj · Pull Request #24240 · JuliaLang/julia · GitHub).

3 Likes

I don’t see any harm in it.

Thank everybody for the reactivity!
If I understand, negative powers already work in 0.7 for literal exponents. That would be an argument to
add for consistency Tamas Papp’s implementation (with Jameson’s improvement) ,so negative exponents
work more generally.