powers of polynomials in characteristic p>0

[finotti]

Dear all,

I need to compute some larger powers of polynomials in characteristic
p>0. I've noticed that Sage does not do it very efficiently, as even
f^(p^n) takes a long time. I wrote the following code then:

================================
# write m in base n (as vector):
# [v0, v1, ..., vk] --> m = v0 + v1*n + ... + vk*n^k
def base_n(m,n):
if m == 0:
return []
tm, m0 = divmod(m,n)
return [ m0 ] + base_n(tm,n)

# powers of powers of p
def pol_p_power(f,p,k):
# computes f^(p^k)
pp=parent(f).characteristic()
if pp== p:
# right call!
res=0
for mon in f.monomials():
res+=f.monomial_coefficient(mon)^(p^k)*mon^(p^k)
return res
else:
# wrong p -- regular power
f^(p^k)

# compute f^n using the characteristic!
def pol_power(f,n):
P=parent(f)
p=P.characteristic()
v=base_n(n,p) # write n in base p
pwrs={0 : P(1), 1 : f} # save computed powers
res=1
for i in range(len(v)):
if v[i]!=0:
# need this term
if v[i] in pwrs:
# power computed before
res*=pol_p_power(pwrs[v[i]],p,i)
else:
# power not computed before -- save it!
pwrs[v[i]]=f^(v[i])
res*=pol_p_power(pwrs[v[i]],p,i)
return res


# compute maximum of a vector
def my_max(v):
res=v[0]
for i in range(1,len(v)):
if v[i] > res: res=v[i]
return res


# compute f^n using the characteristic
# in this one we precompute all necessary powers
def pol_power2(f,n):
P=parent(f)
p=P.characteristic()
v=base_n(n,p)
pwrs={0 : P(1), 1 : f} # save computed powers
M=my_max(v) # maximum power
for i in range(2,M+1):
# precompute all powers up to the maximum
pwrs[i]=pwrs[i-1]*f
for i in range(M+1):
# save some memory by keeping only the ones we need
if not(i in v):
del pwrs[i]
res=1
for i in range(len(v)):
if v[i]!=0:
res*=pol_p_power(pwrs[v[i]],p,i)
return res

================================

Here is a test:

================================
boole[~/comp/sage]$ sage
----------------------------------------------------------------------
| Sage Version 4.1.2, Release Date: 2009-10-13 |
| Type notebook() for the GUI, and license() for information. |
----------------------------------------------------------------------
sage: attach 'pol_power.sage'
sage: p=13
sage: F=GF(p)
sage: P.<x,y,z>=PolynomialRing(F,3)
sage: nterm=5
sage: pwr=10
sage: f=sum( [F.random_element()*x^randint(0,pwr)*y^randint(0,pwr)
*z^randint(0,pwr) for i in\
...: range(nterm)] )
sage: print "f = ", f
f = 2*x^10*y^6*z^3 - 2*x^6*y^3*z^9 - 5*x^3*y^2*z^9 + y^7*z^6 -
x^4*y^7
sage: ntry=20
sage: min=100
sage: max=1000
sage: for i in random_sublist(range(min,max),ntry/(max-min)):
....: print "i = ", i
....: %time a=pol_power(f,i)
....: %time b=pol_power2(f,i)
....: %time c=f^i # sage's builtin
....: a==c
....: b==c
....: print "######################"
....: print ""
....:
i = 101
CPU times: user 0.09 s, sys: 0.02 s, total: 0.11 s
Wall time: 0.11 s
CPU times: user 0.10 s, sys: 0.01 s, total: 0.11 s
Wall time: 0.11 s
CPU times: user 0.89 s, sys: 0.17 s, total: 1.06 s
Wall time: 1.07 s
_18 = True
_18 = True
######################

i = 137
CPU times: user 0.12 s, sys: 0.02 s, total: 0.14 s
Wall time: 0.14 s
CPU times: user 0.13 s, sys: 0.02 s, total: 0.14 s
Wall time: 0.14 s
CPU times: user 3.36 s, sys: 0.72 s, total: 4.08 s
Wall time: 4.08 s
_21 = True
_21 = True
######################

i = 138
CPU times: user 0.17 s, sys: 0.02 s, total: 0.19 s
Wall time: 0.19 s
CPU times: user 0.17 s, sys: 0.02 s, total: 0.19 s
Wall time: 0.19 s
CPU times: user 3.49 s, sys: 0.68 s, total: 4.17 s
Wall time: 4.17 s
_24 = True
_24 = True
######################

i = 310
CPU times: user 0.91 s, sys: 0.36 s, total: 1.27 s
Wall time: 1.27 s
CPU times: user 0.91 s, sys: 0.36 s, total: 1.27 s
Wall time: 1.27 s
CPU times: user 30.18 s, sys: 6.47 s, total: 36.65 s
Wall time: 36.65 s
_27 = True
_27 = True
######################

i = 312
CPU times: user 0.23 s, sys: 0.07 s, total: 0.30 s
Wall time: 0.30 s
CPU times: user 0.21 s, sys: 0.08 s, total: 0.29 s
Wall time: 0.29 s
CPU times: user 34.07 s, sys: 7.31 s, total: 41.39 s
Wall time: 41.40 s
_30 = True
_30 = True
######################

i = 417
CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
Wall time: 0.01 s
CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
Wall time: 0.01 s
CPU times: user 74.05 s, sys: 6.54 s, total: 80.58 s
Wall time: 80.62 s
_33 = True
_33 = True
######################

i = 456
CPU times: user 0.04 s, sys: 0.00 s, total: 0.04 s
Wall time: 0.04 s
CPU times: user 0.04 s, sys: 0.00 s, total: 0.05 s
Wall time: 0.05 s
CPU times: user 91.65 s, sys: 17.02 s, total: 108.66 s
Wall time: 108.69 s
_36 = True
_36 = True
######################

i = 491
CPU times: user 4.06 s, sys: 0.94 s, total: 5.00 s
Wall time: 5.00 s
CPU times: user 4.00 s, sys: 0.95 s, total: 4.95 s
Wall time: 4.95 s
CPU times: user 144.18 s, sys: 31.28 s, total: 175.46 s
Wall time: 175.48 s
_39 = True
_39 = True
######################

i = 513
CPU times: user 0.44 s, sys: 0.19 s, total: 0.63 s
Wall time: 0.63 s
CPU times: user 0.42 s, sys: 0.21 s, total: 0.62 s
Wall time: 0.62 s
CPU times: user 213.25 s, sys: 49.80 s, total: 263.05 s
Wall time: 263.12 s
_42 = True
_42 = True
######################

i = 522
CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
Wall time: 0.00 s
CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
Wall time: 0.00 s
CPU times: user 209.37 s, sys: 48.57 s, total: 257.93 s
Wall time: 257.99 s
_46 = True
_46 = True
######################

i = 603
CPU times: user 0.24 s, sys: 0.02 s, total: 0.27 s
Wall time: 0.27 s
CPU times: user 0.22 s, sys: 0.06 s, total: 0.28 s
Wall time: 0.29 s
CPU times: user 230.83 s, sys: 57.76 s, total: 288.60 s
Wall time: 288.69 s
_50 = True
_50 = True
######################

i = 662
CPU times: user 14.26 s, sys: 3.50 s, total: 17.76 s
Wall time: 17.76 s
CPU times: user 14.60 s, sys: 3.43 s, total: 18.03 s
Wall time: 18.03 s
CPU times: user 430.12 s, sys: 106.54 s, total: 536.66 s
Wall time: 536.77 s
_54 = True
_54 = True
######################

i = 688
CPU times: user 1.95 s, sys: 0.86 s, total: 2.81 s
Wall time: 2.81 s
CPU times: user 1.98 s, sys: 0.82 s, total: 2.80 s
Wall time: 2.80 s
CPU times: user 559.77 s, sys: 128.93 s, total: 688.70 s
Wall time: 688.82 s
_58 = True
_58 = True
######################

i = 697
CPU times: user 0.06 s, sys: 0.00 s, total: 0.06 s
Wall time: 0.06 s
CPU times: user 0.04 s, sys: 0.00 s, total: 0.04 s
Wall time: 0.04 s
CPU times: user 556.29 s, sys: 123.71 s, total: 680.00 s
Wall time: 680.15 s
_62 = True
_62 = True
######################

i = 758
CPU times: user 0.17 s, sys: 0.00 s, total: 0.18 s
Wall time: 0.18 s
CPU times: user 0.20 s, sys: 0.01 s, total: 0.21 s
Wall time: 0.21 s
CPU times: user 585.18 s, sys: 144.59 s, total: 729.77 s
Wall time: 729.88 s
_65 = True
_65 = True
######################

[snip]
================================

So, it seems that these powers could be improved a lot. My second
function did not show much improvement in these cases, but I think it
would in other (maybe extreme) examples.

Unfortunately, I don't know enough to suggest a patch, but hopefully
some one can find a quick way to improve those computations.

Best to all,

Luis

[Robert Bradshaw]

Interesting point, I've made http://trac.sagemath.org/sage_trac/ticket/7253

[finotti]

Hi,

On Oct 20, 1:16 am, Robert Bradshaw <rober...@math.washington.edu>
wrote:
> Interesting point, I've madehttp://trac.sagemath.org/sage_trac/ticket/7253

I know that Sage developers have different priorities, but this is
sort of important to me... So, is there a way I can redefine how Sage
computes powers of polynomials in characteristic p>0 locally (until
this is "fixed" in Sage)? Any pointers or references?

BTW, the current way Sage performs this powers is already *much* more
efficient than MAGMA (expect when the power is p^n, in which case
MAGMA does it instantly, while sage doesn't.)

Thanks,

Luis

[Simon King]

Hi Luis!

On 20 Nov., 19:10, finotti <luis.fino...@gmail.com> wrote:
...

> I know that Sage developers have different priorities, but this is
> sort of important to me...

Definitely polynomials *are* a priority for some developers...

> So, is there a way I can redefine how Sage
> computes powers of polynomials in characteristic p>0 locally (until
> this is "fixed" in Sage)?  Any pointers or references?

First, I would produce a clone of Sage, in order to not destroy your
installation by mistake. So, in the shell, do
sage -clone work
where you can replace "work" by another word that you like (except
"main").

Then, if I am not mistaken, multivariate polynomials over finite
fields and over the rationals are implemented in
SAGE_ROOT/devel/sage-work/sage/rings/polynomial/
multi_polynomial_libsingular.pyx
if I am not mistaken. Univariate polynomials are somewhere else.

Basic arithmetic in python/cython is provided by certain special
methods -- you may consult the Python references about them. It is
__mul__ for multiplication, __add__ for sums, and __pow__ for
exponentiation. So, __pow__ is the method that you have to edit.

Once you've done, you need to compile and test your code. This is done
by starting Sage with an additional command line option,
sage -br work
"b" means "build", "r" means "run" (so, Sage will start if the
compilation worked), and "work" refers to the clone (with "main", you
would return to your original installation).

Then you can test your code.

If it works, you can produce a patch -- this is explained in the Sage
Developer's Guide. This should then be uploaded to the Sage Trac (I
think William Stein is the person to ask for an account). It will be
reviewed, and may eventually be part of the Sage distribution.

Best regards, and thanks for your contribution!
Simon

[finotti]

Dear Simon,

On Nov 20, 2:19 pm, Simon King <simon.k...@nuigalway.ie> wrote:
> Hi Luis!

>
> First, I would produce a clone of Sage, in order to not  destroy your
> installation by mistake. So, in the shell, do
>   sage -clone work
> where you can replace "work" by another word that you like (except
> "main").
>

> (...)

Thanks for the info. I will give it a try. I thought I could just
add a blob (which maybe could be loaded at start up)
redefining .__pow.__, but that I will try your suggestion.

Thanks for the quick and helpful reply.

Best,

Luis

[John Cremona]

One difficulty with implementing this is that there are lots of
different polynomial classes, and also lots of different __pow__
methods. In each case one could perhaps test to see of the
characteristic was p>0 and the exponent a power of p, and call special
code. (I think this was done a while ago for powering elements of
finite fields.)

In the case of univariate polynomials over GF(p) the type of the
polynomial is
'sage.rings.polynomial.polynomial_zmod_flint.Polynomial_zmod_flint'
which suggests to me that the powering is done within the Flint
library. In which case that is where improvements should be made (in
this case).

John Cremona