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.