/*****************************************************************




 ****************************************************************/

Q := RationalField();
R<t> := LaurentSeriesRing(RationalField());

intrinsic formal_w(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
{}
   w := t^3 + O(t^prec);
   a1, a2, a3, a4, a6 := Explode(aInvariants(E));
   for n in [4..prec-1] do
      w := w +  t^n*(a1*Coefficient(w,n-1) + 
                a2*Coefficient(w,n-2) + 
                a3*&+[Integers()|Coefficient(w,i)*Coefficient(w,n-i) : i in [3..n-3]] + 
                a4*&+[Integers()|Coefficient(w,i)*Coefficient(w,n-1-i) : i in [3..n-4]] + 
                a6*&+[Integers()|Coefficient(w,i)*Coefficient(w,j)*Coefficient(w,n-i-j) :
                         i in [3..n-3], j in [3..n-3] | n-i-j ge 3]);
   end for;
   return w;
end intrinsic;

intrinsic formal_x(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
{}
   return t/formal_w(E,prec+5);   
end intrinsic;

intrinsic formal_y(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
{}
   return -1/formal_w(E,prec+6);   
end intrinsic;

intrinsic formal_omega(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
{Compute the formal differential in terms of t, i.e., return the function f(t) 
such that f(t)dt is the formal differerntial omega.

Use that omega = dx/(2*y + a1*x + a3) = f(t) dt.
}

   a1, a2, a3, a4, a6 := Explode(aInvariants(E));
   x := formal_x(E,prec);
   y := formal_y(E,prec);
   return Derivative(x)/(2*y + a1*x + a3);
end intrinsic;

intrinsic formal_log(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
{The formal log in terms of t.  This is sometimes written z=z(s).  It is the
integral with respect to t of the formal differential.}
   y := formal_y(E,prec);
   x := formal_x(E,prec);
   a1, a2, a3, a4, a6 := Explode(aInvariants(E));
   return Integral(1/(2*y + a1*x + a3) * Derivative(x));
end intrinsic;

intrinsic formal_inverse(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
{}
   x := formal_x(E,prec);
   y := formal_y(E,prec);
   a1, a2, a3, a4, a6 := Explode(aInvariants(E));
   return x/(y + a1*x + a3);
end intrinsic;

intrinsic formal_group(E::CrvEll, prec::RngIntElt) -> RngSerLaurElt
{}
   t1 := t;
   S<t2> := LaurentSeriesRing(R);
   w := formal_w(E,prec);
   function tsum(n)
      return &+[t2^m * t1^(n-m-1) : m in [0..n-1]];
   end function;
   lam := &+[tsum(n)*Coefficient(w,n) : n in [3..prec-1]];

//   lam  := (Evaluate(w,t2+O(t2^prec)) - Evaluate(w,t1+O(t1^prec)))/(t2-t1);

   w1 := Evaluate(w,t1+O(t^prec));
   nu := w1 - lam*(t1+O(t^prec));
   a1, a2, a3, a4, a6 := Explode(aInvariants(E));
/*   t3 := -t1 - t2 + \
         (a1*lam + a3*lam^2 - a2*nu - 2*a4*lam*nu - 3*a6*lam^2*nu)/ 
         (1 + a2*lam + a4*lam^2 + a6*lam^3) + O(t1^prec) + O(t2^prec);
*/
   t3 := -t1 - t2 + \
         (-a1*lam - a3*lam^2 - a2*nu - 2*a4*lam*nu - 3*a6*lam^2*nu)/ 
         (1 + a2*lam + a4*lam^2 + a6*lam^3) + O(t1^prec) + O(t2^prec);
   inv := formal_inverse(E, prec);
   return Evaluate(inv, t3);
end intrinsic;



intrinsic formal_mult(E::CrvEll, n::RngIntElt, prec::RngIntElt) -> RngSerLaurElt
{}
   if n eq -1 then
      return formal_inverse(E,prec);
   end if;

   // The following is less direct, but it works more quickly...
   L := FunctionField(E);
   EL := BaseExtend(E,L);
   Q := EL![L.1,L.2];
   x,y := Explode(Eltseq(n*Q));
   a := -x/y;
   xx := formal_x(E,prec);
   yy := formal_y(E,prec);
   phi := hom<Parent(x) -> Parent(xx) | xx,yy>;
   return phi(a);

// more direct method... (not used if the above isn't commented out)

   if n eq 1 then
      return t;
   end if;
   if n eq -1 then
      return formal_inverse(E,prec);
   end if;
   if n lt 0 then
      return Evaluate(formal_inverse(E,prec),formal_mult(E,-n,prec));
   end if;
   // Now n is positive and >= 2.
   F := formal_group(E,prec);
   g := t;  // 1
   for m in [2..n] do
      g := Evaluate(F,g);
   end for;
   return g;

/*
   // Use a binary powering algorithm to compute g = [n](t).
   g := t + O(t^prec);  // [1]
   ans := 0;
   while n gt 0 do
      if n mod 2 ne 0 then
         ans := Evaluate(Evaluate(F,ans), g);
      end if;
      g := Evaluate(Evaluate(F,g),g);   // [m] <-- [2m], where g(t) <--> [m].
      n := n div 2;
   end while;
   return ans;
*/
end intrinsic;

intrinsic indeterminate_sigma(prec::RngIntElt) -> RngSerElt
{}
   assert prec ge 3;
   S := PolynomialRing(Q, prec-2);
   AssignNames(~S,[Sprintf("s%o",i) : i in [1..prec-2]]);
   T<t> := LaurentSeriesRing(S);
   // the 2 ------------------------> is for special case s1=0.
   //return t + &+[ S.i*t^(i+1) : i in [2..prec-2]] + O(t^prec);
   return t + &+[ S.i*t^(i+1) : i in [1..prec-2]] + O(t^prec);
end intrinsic;

intrinsic formal_divpoly(E::CrvEll, m::RngIntElt, prec::RngIntElt) -> RngSerElt 
{}
   if m eq 1 then
      return R!1;
   end if;
   if m eq -1 then
      return R!(-1);
   end if;
   if m eq 2 then
      return Sqrt(Evaluate(DivisionPolynomial(E,m),formal_x(E,prec)));
   end if;
   assert IsOdd(m);
   return Evaluate(DivisionPolynomial(E,m),formal_x(E,prec));
end intrinsic;

intrinsic sigma_relation(E::CrvEll, m::RngIntElt, 
                  numvars::RngIntElt, prec::RngIntElt) -> SeqEnum
{}
   sig := indeterminate_sigma(numvars+2);
   T<t> := Parent(sig);
   mult := T!Eltseq(formal_mult(E,m,prec)) * t;
   fm := T!Eltseq(formal_divpoly(E,m,prec)) * t^(1-m^2);
   return (t^(m^2-1)*Evaluate(sig,mult) - t^(m^2-1)*sig^(m^2)*fm);
end intrinsic;

/* The two functions formal_sigma_in_s2 and formal_sigma_in_s2_2 compute
   the same thing in totally different ways.  The first one is more
   direct and hence vastly faster. 
*/

intrinsic WeierstrassP(E::CrvEll, prec::RngIntElt) -> FldPrElt
{}
   vprint Heegner :  "Computing WeierstrassP";
   L<w> := PolynomialRing(Rationals());
   R<t> := LaurentSeriesRing(L);
   p := t^(-2);
   c4, c6 := Explode(cInvariants(E));
   a := 4; b := -c4/12; c := -c6/216;
   for n in [1..prec] do
      ptest := p + w*t^n + O(t^(n+1));
      ptest_prime := Derivative(ptest);
      f := ptest_prime^2 - a*ptest^3 - b*ptest - c;   
      x := Coefficient(f,-4+n);
      coeff := -Coefficient(x,0)/Coefficient(x,1);
      p := p + coeff*t^n;
   end for;
   vprint Heegner :  "Done computing Weierstrass P";
   return LaurentSeriesRing(Rationals())!(p + O(t^(prec+1)));
end intrinsic;

intrinsic formal_sigma_in_s2(E::CrvEll, prec::RngIntElt) -> RngSerElt
{}
   // This is adapted from Christian's PARI code, which is also exactly
   //  what Mazur-Tate worked out on 2004-06-09.

   R<s2> := FieldOfFractions(PolynomialRing(RationalField()));
   S<z> := LaurentSeriesRing(R);
   w := S!WeierstrassP(E,prec);
   w := &+[Coefficient(w,n)*z^n : n in [0..prec-1]] + O(z^prec);
   sigma_of_z := z*Exp( s2*z^2 - Integral(Integral(w)) + O(z^prec) );
   z_of_s := S!formal_log(E, prec);
   sigma_of_s := Evaluate(sigma_of_z, z_of_s);
   R<s2> := PolynomialRing(RationalField());
   S<t> := PowerSeriesRing(R);
   return S!sigma_of_s;
end intrinsic;



/*******************************************************

Karl Rubin and I were wondering whether you could do a p-adic height
computation, the point of which is to figure out whether the p-adic
regulator of elliptic curves has a "tendency" to be not divisible by p. We
want to avoid a few degenerate cases, so the smallest p that is useful to
us is p=5. Is the strategy-- described below--, that systematically runs
through 5-adic height computations for points on quadratics twists of the
(non-3-Eisenstein) elliptic curve of conductor 37 reasonably easy to do?


%%%%%%%%%%%%%%%%%%%
Our untwisted equation is E:  y^2 = 4x^3-4x+1. For any nonzero integer D we
consider the twisted elliptic curve  E_D:  D.y^2 = 4x^3-4x+1.  Put F(x): =
4x^3-4x+1. Let x run through integers

  a=  0, \pm 1, \pm 2, \pm3, ...

and for each of these a's,

1) Factor F(a) = d.b^2,  where d is square-free;

2) throw out the ones where d is divisible by 5, or where the parity of the
rank of E_d is even;

3) for the rest, viewing P:=(a,b) as a rational point-- it will be of
infinite order--on the elliptic curve E_d, so one can compute the 5-adic
height of P on the elliptic curve E_d: this should be a 5-adic integer if
the height is normalized right, and there should be a good number of
instances where  this height is a 5-adic unit; throw them out, and

4) print the rest.
%%%%%%%%%%%%%%%%

Note that we only want to know the 5-adic height modulo 5.

Then it would be our job to figure out whether strange things occur for
this printed list of instances: e.g., the point P might be divisible by 5,
or the Mordell-Weil rank might be > 2, or other stuff..., and we want to
know what other stuff there might be.

****************************************************/

intrinsic point(a::RngIntElt) -> PtEll, RngIntElt, BoolElt
{Returns:
     * point P_d on a minimal model for E_d,
     * the integer d, 
     * a boolean -- true <==> an_rank(E_d) is odd.
}
   d, b := SquareFree(4*a^3-4*a+1);
   E := EllipticCurve([-16*d^2, 16*d^3]);
   P := E![4*d*a, 4*d^2*b];
   EM, phi := MinimalModel(E);
   P := phi(P);
   return P, d, RootNumber(E) eq 1;
end intrinsic;

function prime_divisors(n)
   return [F[1] : F in Factorization(n)];
end function;

intrinsic prep_multiple(P::PtEll, p::RngIntElt) -> PtEll, RngIntElt
{An integer n such that n*P is in the subgroup of
E that reduces to the identity mod p and lands in the connected
component of the Neron model modulo all primes.}
   require IsPrime(p) : "Argument 2 must be prime.";

   E := Curve(P);
   c := LCM([TamagawaNumber(E,q) : q in prime_divisors(Conductor(E))]);
   N := #ChangeRing(E,GF(p));
   n := LCM(c,N);
   return n;
end intrinsic;

intrinsic extended_log(x::FldPadElt) -> FldPadElt
{}
   p := Prime(Parent(x));
   u := x * p^(-Valuation(x));   
   return Log(u^(p-1))/(p-1);
end intrinsic;

intrinsic x_coord_d(Q::PtEll) -> RngIntElt
{}
   t, d := IsSquare(Denominator(Q[1]));
   assert t;
   return d;
end intrinsic;

intrinsic t_val(Q::PtEll)  -> FldRatElt
{Given input (x,y), return t = -x/y.}
   return -Q[1]/Q[2];
end intrinsic;

intrinsic formal_val_of_sigma_in_s2(E::CrvEll, t::FldRatElt, 
                     numvars::RngIntElt, prec::RngIntElt) -> .
{}
   sigma := formal_sigma_in_s2(E, numvars, prec);
   return Evaluate(sigma, t);
end intrinsic;

intrinsic mult_height_of_associated_prepared_point(P::PtEll, p::RngIntElt) -> .
{Returns multiplicative p-adic height as series in sigma.}
   n := prep_multiple(P, p);
   Q := n*P;
   d := x_coord_d(Q);  // positive sqrt of denominator
   HQ := formal_val_of_sigma_in_s2(Curve(P), t_val(Q), 5, 10)/d;
   return HQ, n;
end intrinsic;

intrinsic solve_for_s2(sigma::RngSerElt, p::RngIntElt) -> .
{}
   a := 0;   // a is what we know about s2 so far. 
   n := 0;
   S<t> := Parent(sigma);
   R<s2> := BaseRing(S);
   prec := AbsolutePrecision(sigma);
   for i in [3..prec-1] do
      c := Coefficient(sigma,i);
      d := Coefficients(c);
      m := Max([0] cat [Valuation(Denominator(x),p) : x in d]);
      if m gt 0 then
         f := p^m*c;
         // Now viewing f as a poly in s2, it must be 0 modulo p^m, so 
         // as to cancel the denominator. 
         R<x> := PolynomialRing(GF(p));     
         g := R!f;
         X := Roots(g);
         assert #X le 1;
         if #X eq 1 then
            b := Integers()!X[1][1];
            a +:= b*p^n;
            n +:= 1;
            // Now s2 = a + O(p^(n+1)).
            z := b+ p*s2;
            sigma := &+[Evaluate(Coefficient(sigma,r),z)*t^r : r in [0..prec-1]] + O(t^prec);
         end if;
      end if;
   end for;
   return pAdicField(p,n)!a;
end intrinsic;

intrinsic compute_s2(E::CrvEll, p::RngIntElt, prec::RngIntElt, s1::.) -> RngElt
{
  Compute the quantity s2, which is the coefficient of t^3 in the expansion
  of the p-adic sigma function associated to E.
}
/*
  If I computed everything correctly, we should have that

     s2 = -(a1^2+4*a2+E2(E))/24 - (a2 - a1*s1 - s1^2)/2

*/
    e2 := E2(E,p,prec);
    a1, a2, a3, a4, a6 := Explode(aInvariants(E));
    return (-(e2 + a1^2 + 4*a2)/12 -a2 +a1*s1 + s1^2)/2;
end intrinsic; 

intrinsic val_of_sigma(Q::PtEll, p::RngIntElt,
                       prec::RngIntElt) -> RngIntElt
{}
   E := Curve(Q);
   t := t_val(Q);
   time sigma := formal_sigma_in_s2(E, prec);
   s1 := pAdicField(p)!Coefficient(Coefficient(sigma,2),0);
   print "s1 = ", s1;
   time s2 := solve_for_s2(sigma, p);     
   print "s2 = ", s2;
   time s2b := compute_s2(E, p, 40, s1);
   print "s2b = ", s2b;
   print "s2 - s2b = ", Parent(s2)!(s2 - s2b);
   v := Evaluate(Evaluate(sigma, t),s2);
   return v;
end intrinsic;

intrinsic height_of_point_modp(P::PtEll, p::RngIntElt) -> .
{Returns p-adic height of point P.}
   n := prep_multiple(P, p);
   Q := n*P;
   d := x_coord_d(Q);  // positive sqrt of denominator
   HQ := extended_log(val_of_sigma(Q, p, 10) / d);
   HP := HQ/n^2;
   return GF(p)!(Integers()!HP);
end intrinsic;

intrinsic myExp(f::.) -> RngSerLaurElt
{}
   return &+[f^n/Factorial(n) : n in [0..AbsolutePrecision(f)]];
end intrinsic;

intrinsic myIsZero(a::.) -> .
{}
   return Valuation(a) ge Precision(Parent(a));
end intrinsic;

intrinsic myIntegral(f::RngSerLaurElt) -> RngSerLaurElt
{
 Integrate Laurent series, without stupid check on coefficient
 of t^(-1) being 0, which doesn't work in MAGMA, since nonzero
 for p-adic elements is messed up.
}
   require myIsZero(Coefficient(f,-1)) : "Coefficient of 1/t must be zero.";
   pr := AbsolutePrecision(f);
   t := Parent(f).1;
   return &+[(Coefficient(f,n)/(n+1))*t^(n+1) : n in [0..pr-1] | 
                not myIsZero(Coefficient(f,n)) ] + O(t^pr);   

//   return -Coefficient(f,-2)*t^(-1) + 
//        &+[(Coefficient(f,n)/(n+1))*t^(n+1) : n in [0..pr-1]] + O(t^pr);   
end intrinsic;

intrinsic sigma_alg2(E::CrvEll, p::RngIntElt, prec::RngIntElt :
                    e2prec := 50) -> RngElt
{
   Compute and return the formal series for sigma.  
}
   a1, a2, a3, a4, a6 := Explode(aInvariants(E));
   R<t> := LaurentSeriesRing(pAdicField(p));

   // Compute w such that omega = w(t) dt
   w := R!formal_omega(E, prec);
   //t := Parent(w).1;

   // Compute E2(E,omega)
   e2 := E2(E, p, e2prec);

   // Compute canonical X-function
   x := R!formal_x(E, prec);
   X := x + R!(a1^2 + 4*a2 + e2);

   /* Solve the differential equation
         -d/omega ( (d(sigma)/omega) * (1/sigma) ) = X
         subject to omega = t + higher terms
    */
   intwX := myIntegral(w*X);
print "w = ", w;
print "intwX = ", intwX;
print "w*intwX = ", w*intwX;
   f := -myIntegral(w*intwX);
print "f after int = ", f;
   // But the constant term and linear term are not well-defined just from integrating, 
   // so we set them to 0, 1, respectively.
   pr := AbsolutePrecision(f);
   f := R!([0,1] cat [Coefficient(f,n) : n in [2..pr-1]]) + O(t^(pr-5));
   print "f = ", f;
/*   S<t> := PowerSeriesRing(BaseRing(R));
   pr := AbsolutePrecision(f);
Eltseq(f);
print pr;
   f := S!([0] cat [Coefficient(f,n) : n in [1..pr-1]]) + O(t^pr);
*/
   // Finally find sigma.
print "f = ", f;
   sigma := myExp(f);

   return sigma;
end intrinsic;

intrinsic sigma_alg1(E::CrvEll, p::RngIntElt, prec::RngIntElt) -> .
{}
   time sigma_in_s2 := formal_sigma_in_s2(E, prec);
   s1 := pAdicField(p)!Coefficient(Coefficient(sigma_in_s2,2),0);
   time s2 := solve_for_s2(sigma_in_s2, p);     
   R<t> := PowerSeriesRing(pAdicField(p));
   pr := AbsolutePrecision(sigma_in_s2);
   sigma := &+[Evaluate(Coefficient(sigma_in_s2,n),s2)*t^n : n in [1..pr-1]];
   return sigma;
end intrinsic;

intrinsic sigma_using_e2(E::CrvEll, p::RngIntElt, prec::RngIntElt :
                    e2prec := 50) -> .
{
   Compute and return the formal series for sigma.  
}
   prec +:= 2;

   // 1. Compute value of p-adic E_2, using Kedlaya's MW algorithm
   e2 := E2(E,p,e2prec);

   // 2. Define some notation.
   K := pAdicField(p,Precision(Parent(e2)));
   a1, a2, a3, a4, a6 := Explode(aInvariants(E));
   S<t> := LaurentSeriesRing(K);
   
   // 3. Compute formal Weierstrass p function in terms of t.
   x  := S!formal_x(E, prec);
   wp := x + (a1^2 + 4*a2)/12;
   
   // 4. Compute wp(z), where z = FormalLog(t)
   z_of_t := S!formal_log(E, prec);
   t_of_z := Reverse(z_of_t);
   wp_of_z := Evaluate(wp, t_of_z);
   R<z> := S;
   
   // 5. Compute
   pr := AbsolutePrecision(wp_of_z);
   // Let g be    1/z^2 - \wp + E2/12.  Notice that we compute this 
   // by just ignoring the coefficients of wp_of_z for n=-2.
   g := e2/12 +  &+[z^n*Coefficient(-wp_of_z,n) : n in [-1..pr-1]] + O(z^pr);
   /* It's *really* weird that it isn't -E2/12, since the canonical
      eks function is \wp + E2/12. 

      There's a factor of -1 in my definition of E2 in the other file,
      which I think Tate and I got out of Katz's paper.   Maybe that
      sign is wrong?   There's some -1 normalization that is different
      in two papers...
   */
   
   // 6. Integrate twice, exponentiate, etc.
   sigma_of_z := z*myExp(myIntegral(myIntegral(g)));
   sigma_of_t := Evaluate(sigma_of_z, z_of_t);
   R<t> := PowerSeriesRing(K);

   return R!sigma_of_t + O(t^(prec-2));
end intrinsic;

intrinsic height_function(E::CrvEll, p::RngIntElt, prec::RngIntElt) -> .
{Returns p-adic height of point P.}
   sig := sigma_using_e2(E, p, prec : e2prec := prec) ;
   function h(P)
      assert Curve(P) eq E;
      n := prep_multiple(P, p);
      Q := n*P;
      d := x_coord_d(Q);  // positive sqrt of denominator
      t := t_val(Q);
      v := Evaluate(sig,t);
      HQ := extended_log(v/d);
      HP := HQ/n^2;
      return HP/p;
   end function;
   return h;
end intrinsic;

intrinsic regulator(E::CrvEll, p::RngIntElt) -> .
{}
    return regulator(E,p,20);
end intrinsic;

intrinsic regulator(E::CrvEll, p::RngIntElt, prec::RngIntElt) -> .
{
Computes and returns the p-adic regulator of E.
}
   h := height_function(E, p, prec);
   G, f := MordellWeilGroup(E);
   B := [f(G.i) : i in [1..Ngens(G)] | Order(G.i) eq 0];
   function pairing(P,Q)
      return (h(P+Q)-h(P)-h(Q))/2;
   end function;
   K := Parent(h(B[1]));
   M := MatrixAlgebra(K,#B)!0;
   for i in [1..#B] do
      for j in [i..#B] do
         aij := pairing(B[i],B[j]);
         M[i,j] := aij;
         M[j,i] := aij;
      end for;
   end for;
   return Determinant(M);
   //return Determinant(M), M;
end intrinsic;

intrinsic good_ordinary_primes(E::CrvEll, pmax::RngIntElt) -> SeqEnum
{}
   N := Conductor(E);
   return [p : p in [5..pmax] | IsPrime(p) and
                                N mod p ne 0 and
                                TraceOfFrobenius(ChangeRing(E,GF(p))) mod p ne 0];
end intrinsic;

intrinsic val_of_sigma_alg2(Q::PtEll, p::RngIntElt,
                       prec::RngIntElt) -> RngElt
{
  Return the value of the p-adic sigma function on Q, using the new algorithm.
}
   E := Curve(Q);
   sigma := sigma_alg2(E, p, prec);
   t := t_val(Q);
   return Evaluate(sigma, t);
end intrinsic;

intrinsic height_of_point(P::PtEll, p::RngIntElt, prec::RngIntElt) -> .
{Returns p-adic height of point P.}
   n := prep_multiple(P, p);
   Q := n*P;
   d := x_coord_d(Q);  // positive sqrt of denominator
   v := val_of_sigma(Q, p, prec);
   HQ := extended_log(v/d);
   HP := HQ/n^2;
   return HP/p;
end intrinsic;

intrinsic experiment(p::RngIntElt, a_range::SeqEnum, prec::RngIntElt) 
{}
   for a in a_range do
      P, d, t := point(a);
      if not t or d mod p eq 0 then
         continue;
      end if;
      h := height_of_point(P, p, prec) ;
      E := Curve(P);
      printf "[%o, %o], ", a, h;
   end for;
end intrinsic;

intrinsic ranks(p::RngIntElt, a_range::SeqEnum, prec::RngIntElt) 
{}
   for a in a_range do
      P, d, t := point(a);
      if not t or d mod p eq 0 then
         continue;
      end if;
      E := Curve(P);
print aInvariants(E);
      printf "a = %o\t\trank = %o\n", a, Rank(E);
   end for;
end intrinsic;

intrinsic padicprint(x::.)  -> .
{}
   z := Integers()!x;
   p := Prime(Parent(x));
   i := 0;
   s := "";
   while z ne 0 and i lt Precision(Parent(x)) do
      c := z mod p;
      if c ne 0 then
         if c ne 1 then         
            s *:= Sprintf("%o*", c);        
         end if;
         s *:= Sprintf("%o^%o + ", p, i);
      end if;
      z := z div p;
      i := i + 1;
   end while;
   s *:= Sprintf("O(%o^%o)", p, i);
   return s;
end intrinsic;


