Optimize `AbstractAlgebra.Generic.pow_rmul`, at least use square-and-multiply
The function AbstractAlgebra.Generic.pow_rmul is used to implement powering for Generic.MPoly when the coefficient ring is in positive characteristic. It looks like this:
function pow_rmul(a::MPoly{T}, b::Int) where {T <: RingElement}
b < 0 && throw(DomainError(b, "exponent must be >= 0"))
if iszero(a)
return zero(a)
elseif b == 0
return one(a)
end
z = deepcopy(a)
for i = 2:b
z = mul!(z, z, a)
end
return z
end
This naive powering code clearly is suboptimal and should be improved, at the very least we can use square-and-multiply.
Example input which ends up using pow_rmul:
julia> R = residue_ring(ZZ, 4)[1]
Residue ring of integers modulo 4
julia> Rxy, (x,y) = polynomial_ring(R, [:x, :y])
(Multivariate polynomial ring in 2 variables over R, AbstractAlgebra.Generic.MPoly{EuclideanRingResidueRingElem{BigInt}}[x, y])
julia> a = 2x+2y
2*x + 2*y
julia> a^20
0
This is on purpose. Apparently doing square-and-multiply leads to more intermediate terms. You can ask @fieker for more details.
In that case my request would be that this should be documented in the source code with a comment
There is a comment before pow_fps which states this:
# Implement fps algorithm from "Sparse polynomial powering with heaps" by
# Monagan and Pearce, except that we modify the algorithm to return terms
# in ascending order and we fix some issues in the original algorithm
# http://www.cecm.sfu.ca/CAG/papers/SparsePowering.pdf
Looking into the linked PDF it explicitly mention rmul and that is much better than square-and-multiply in the sparse case, so I guess the comment I am wishing for could just reference that.
Also, at the same time this code calling pow_rmul should be corrected:
elseif !is_exact_type(T) || !iszero(characteristic(base_ring(a)))
# pow_fps requires char 0 or exact ring, so use pow_rmul if not or unsure
return pow_rmul(a, b)
else
# ... code using pow_fps ...
The comment and the condition check claim different things:
- the code uses
pow_fpsif we have an exact type and characteristic zero - the comment claims
pow_fpsrequires an exact type or characteristic zero
So the comment or the code (or both...) are wrong...