Interesting point, I've made http://trac.sagemath.org/sage_trac/ticket/7253
On Oct 19, 2009, at 1:49 PM, finotti wrote:
>
> 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
>
> >
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org
-~----------~----~----~----~------~----~------~--~---
