// write m in base n (as vector):
// [v0, v1, ..., vk] --> m = v0 + v1*n + ... + vk*n^k
base_n := function(m,n)
    if m eq 0 then
        return [];
    end if;
    tm, m0:=Quotrem(m,n);
    return [ m0 ] cat $$(tm,n);
end function;



// for pol. rings
aux_pol_p_power1 := function(f,p,k)
    res:=0;
    for mon in Monomials(f) do
        res+:=MonomialCoefficient(f,mon)^(p^k)*mon^(p^k);
    end for;
    return res;
end function;


// for fields of rational functions
aux_pol_p_power2 := function(f,p,k)
    P:=Parent(f);
    num:=Numerator(f);
    den:=Denominator(f);
    return P!(aux_pol_p_power1(num,p,k))/P!(aux_pol_p_power1(den,p,k));
end function;


// powers of powers of p
pol_p_power := function(f,p,k)
    // computes f^(p^k)
    P:=Parent(f);
    if k eq 0 then
        return f;
    end if;
    pp:=Characteristic(P);
    if pp eq p then
         // right call!
        res:=0;
        if IsFinite(P) then
            // finite field
            return f^(p^k);
        end if;
        if IsField(P) then
            // rational function field
            return aux_pol_p_power2(f,p,k);
        end if;
        // polynomial ring case
        return aux_pol_p_power1(f,p,k);
    else
        // wrong p -- regular power
        return f^(p^k);
    end if;
end function;


// compute f^n using the characteristic!
pol_power := function(f,n)
    P:=Parent(f);
    p:=Characteristic(P);
    v:=base_n(n,p); // write n in base p
    pwrs:=[ P!1, f ] cat [ P!0 : i in [2..(p-1)] ]; // save computed powers
    // pwrs[i] will have f^(i-1)
    res:=P!1;
    for i in [1..(#v)] do
        if v[i] ne 0 then
            // need this term
            if pwrs[v[i]+1] eq P!0 then
                // power not computed before -- save it!
                pwrs[v[i]+1]:=f^v[i];
            end if;
            res*:=pol_p_power(pwrs[v[i]+1],p,i-1);            
        end if;
    end for;
    return res;
end function;


// compute f^n using the characteristic!
// now precompute the powers in a different way
pol_power2 := function(f,n : pwrs:=[])
    P:=Parent(f);
    p:=Characteristic(P);
    v:=base_n(n,p); // write n in base p
    max:=Maximum(v);
    if pwrs eq [] then
        // no pwrs was given -- compute it!
        // pwrs[i] will have f^(i-1)
        pwrs:=[ P!1, f ] cat [ P!0 : i in [2..(p-1)] ]; // save computed powers
        // precompute all powers up to max
        for i in [3..(max+1)] do
            pwrs[i]:=f*pwrs[i-1];
            // now remove unecessary ones to save some memory
            if not((i-2) in v) then
                pwrs[i-1]:=P!0;
            end if;
        end for;
    end if;
    res:=P!1;
    for i in [1..(#v)] do
        if v[i] ne 0 then
            // need this term
            if pwrs[v[i]+1] eq 0 then
                // need to compute it
                pwrs[v[i]+1]:=f^(v[i]);
            end if;
            res*:=pol_p_power(pwrs[v[i]+1],p,i-1);
            // do we still need pwrs[v[i]+1]?
            // if not, remove it to save memory
            if not (v[i] in v[(i+1)..#v]) then
                pwrs[v[i]+1]:=P!0;
            end if;
        end if;
    end for;
    return res;
end function;



/*

// Example:
load 'pol_power.m';
p:=13;
F:=GF(p);
P<x,y,z>:=PolynomialRing(F,3);
nterm:=5;
pwr:=10;
f:=&+[ Random([1..(p-1)])*x^Random(pwr)*y^Random(pwr)*z^Random(pwr) : i in [1..nterm] ];
print "f = ", f;
print "";
ntry:=7;
min:=60;
max:=80;
r:=Floor((max - min)/ntry);
print "######################";
for j in [1..ntry] do
    i:=Random(min + (j-1)*r, min + j*r);
    print  "i = ", i;
    time a:=pol_power(f,i);
    time b:=pol_power2(f,i);
    time c:=f^i;
    a eq b;
    b eq c;
    print "######################";
end for;

*/
