AbstractAlgebra.jl icon indicating copy to clipboard operation
AbstractAlgebra.jl copied to clipboard

Optimize `AbstractAlgebra.Generic.pow_rmul`, at least use square-and-multiply

Open fingolfin opened this issue 1 year ago • 3 comments

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

fingolfin avatar Dec 26 '24 00:12 fingolfin

This is on purpose. Apparently doing square-and-multiply leads to more intermediate terms. You can ask @fieker for more details.

thofma avatar Dec 26 '24 07:12 thofma

In that case my request would be that this should be documented in the source code with a comment

fingolfin avatar Dec 26 '24 13:12 fingolfin

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_fps if we have an exact type and characteristic zero
  • the comment claims pow_fps requires an exact type or characteristic zero

So the comment or the code (or both...) are wrong...

fingolfin avatar Dec 26 '24 13:12 fingolfin