// Magma code related to "Fast Jacobian arithmetic for hyperelliptic curves of genus 3"
// by Andrew V. Sutherland, see http://arxiv.org/abs/1607.08602 for details.
// Copyright 2016, Andrew V. Sutherland, 2016.  Licensed under GPL version 3.

// This source file contains implementations of the following alogirthms given in the paper:
// * Precompute, Compose, Reduce, Adjust, Addition, Negation (these work for any genus)
// * TypicalAddition, TypicalDoubling, TypicalNegation (for typical divisor classes in genus 3)
// * TypicalAdditionFormula, TypicalDoublingFormula, TypicalNegationFormula (for typical divisor classes in genus 3)
// We also give functions for generating random hyperelliptic curves over finite fields
// and for generating random divisor classes on them, as well as functions for testing the code.

// We work thoughout with hyperelliptic curve y^2=f(x) where f is monic of degree 2g+2
// (for the Typical* functions we also assume the x^7 coeff of f is zero)
// Points at infinity are P_inf:=(1:1:0) and Pbar_inf:=(1:-1:0)
// Hyperelliptic conjugate of affine point P=(x,y) is -P=(x,-y)
// Infinity divisor is D_inf = ceil(g/2)*P_inf + floor(g/2)*Pbar_inf

// For polynomials u,v with u monic, deg(v)<deg(u) and u dividing f-v^2 we use div[u,v]
// to denote the Mumfor representations of the affine semi-reduced divisor with points (a,v(a))
// where a ranges over the roots of u.

// For u,v as above with deg(u) <= g and 0 <= n <= g-deg(u) we take
// div[u,v,n] = div[u,v] + nP_inf + (g-deg(u)-n)Pbar_inf - D_inf as the unique rep of its divisor class
// For u,v as above with deg(u) <= 2g and 0 <= n <= 2g-deg(u) we use
// div[u,v,n]* to denote div[u,v] + nP_inf + (2g-deg(u)-n)Pbar_inf

// The auxiliary polynomial V is the unique monic poly of deg g+1 for which deg(f-V^2) <= g
// It is computed by Precompute below.

Coeff := function(f,n)
    c:=Coefficients(f); return n+1 le #c select c[n+1] else 0;
end function;

// Verifies that f is a monic separable polynomial of degree 2g+2 >= 6 and returns g
ValidateHyperellipticPolynomial := function(f)
    assert Degree(f) ge 6 and IsEven(Degree(f)) and IsMonic(f) and IsSeparable(f);
    return ExactQuotient(Degree(f)-2,2);
end function;

// Given monic f of degree 2g+2 compute unique monic V of degree g+1 for which deg(f-V^2) <= g
Precompute := function(f)
    g := ValidateHyperellipticPolynomial(f);
    R<x>:=PolynomialRing(BaseRing(f));
    V := x^(g+1);
    i:=g;
    while i ge 0 do
        V +:= 1/2*Coeff(f-V^2,g+1+i)*x^i;
        i-:=1;
    end while;
    assert Degree(f-V^2) le g;
    return V;
end function;

// Verify that D=div[u,v,n] is a valid unique divisor class representative
ValidateReducedDivisor := procedure(D,f)
    g := ValidateHyperellipticPolynomial(f);
    assert IsMonic(D[1]) and Degree(D[2]) lt Degree(D[1]) and D[3] ge 0 and D[3] le g-Degree(D[1]);
    assert IsDivisibleBy(f-D[2]^2,D[1]);
end procedure;

// Verify that D=div[u,v,n]* is a valid semi-reduced divisor with deg(u)+n <= 2g
ValidateSemiReducedDivisor := procedure(D,f)
    g := ValidateHyperellipticPolynomial(f);
    assert IsMonic(D[1]) and Degree(D[2]) lt Degree(D[1]) and D[3] ge 0 and Degree(D[1])+D[3] le 2*g;
    assert IsDivisibleBy(f-D[2]^2,D[1]);
end procedure;

// Given D1=div[u1,v1,n1] and D2=div[u1,v1,n1] compute D3=div[u3,v3,n3]* with D3 ~ D1 + D2
Compose := function(D1,D2,f)
    g := ValidateHyperellipticPolynomial(f);
    ValidateReducedDivisor(D1,f);  ValidateReducedDivisor(D2,f);
    u1:=D1[1]; v1:=D1[2]; u2:=D2[1]; v2:=D2[2]; n1:=D1[3]; n2:=D2[3];
    w,c := XGCD([u1,u2,v1+v2]);
    u3 := ExactQuotient(u1*u2,w^2);
    v3 := ExactQuotient(c[1]*u1*v2+c[2]*u2*v1+c[3]*(v1*v2+f),w) mod u3;
    return <u3,v3,n1+n2+Degree(w)>;
end function;

// Given D1=div[u1,v1,n1]* with deg(u1) > g+1 compute D2=div[u2,v2,n2]* with D2 ~ D1 and deg(u2) <= deg(u1)-2
Reduce := function(D1,f)
    g := ValidateHyperellipticPolynomial(f);
    ValidateSemiReducedDivisor(D1,f);
    u1:=D1[1]; v1:=D1[2]; n1:=D1[3];
    assert Degree(u1) gt g+1;
    u2 := ExactQuotient(f-v1^2,u1);  u2 := u2 / LeadingCoefficient(u2);
    v2 := (-v1) mod u2;
    if Degree(v1) eq g+1 and LeadingCoefficient(v1) eq 1 then
        delta := Degree(u1)-(g+1);
    elif Degree(v1) eq g+1 and LeadingCoefficient(v1) eq -1 then
        delta := g+1-Degree(u2);
    else
        delta := ExactQuotient(Degree(u1)-Degree(u2),2);
    end if;
    return <u2,v2,n1+delta>;
end function;

// Given D1=div[u1,v1,n1]* with deg(u1) <= g+1 compute D2=div[u2,v1,n2] with D2 ~ D1
Adjust := function(D1,f,V)
    g := ValidateHyperellipticPolynomial(f);
    ValidateSemiReducedDivisor(D1,f);
    u1:=D1[1]; v1:=D1[2]; n1:=D1[3];
    assert Degree(u1) le g+1;
    if n1 ge Ceiling(g/2) and n1 le Ceiling(3*g/2)-Degree(u1) then
        return <u1,v1,n1-Ceiling(g/2)>;
    end if;
    if n1 lt Ceiling(g/2) then
        vhat1 := v1 - V + (V mod u1);
        u2 := ExactQuotient (f-vhat1^2,u1);  u2 := u2 / LeadingCoefficient(u2);
        v2 := -vhat1 mod u2;
        n2 := n1 + g+1 - Degree(u2);
    else
        vhat1 := v1 + V - (V mod u1);
        u2 := ExactQuotient (f-vhat1^2,u1);  u2 := u2 / LeadingCoefficient(u2);
        v2 := -vhat1 mod u2;
        n2 := n1 + Degree(u1) - (g+1);
    end if;
    assert Degree(u2) lt Degree(u1) or (n1 lt Ceiling(g/2) and n2 gt n1) or (n1 gt Ceiling(3*g/2)-Degree(u1) and n2 lt n1);
    return $$(<u2,v2,n2>,f,V);
end function;

// Given D1=div[u1,v1,n1] and D2=div[u2,v2,n2] compute D3=div[u3,v3,n3] with D3 ~ D1 + D2
Addition := function(D1,D2,f,V)
    g := ValidateHyperellipticPolynomial(f);
    D := Compose(D1,D2,f);
    while Degree(D[1]) gt g+1 do D := Reduce(D,f); end while;
    return Adjust(D,f,V);
end function;

// Given D1=div[u1,v1,n1] compute D2=div[u2,v2,n2] with D2 ~ 2*D1
Doubling := function (D1,f,V) return Addition(D1,D1,f,V); end function;

// Given D1=div[u1,v1,n1] compute D2=div[u2,v2,n2] with D2 ~ -D1
Negation := function(D1,f,V)
    g := ValidateHyperellipticPolynomial(f);
    ValidateReducedDivisor(D1,f);
    u1:=D1[1]; v1:=D1[2]; n1:=D1[3];
    if g mod 2 eq 0 then return <u1,-v1,g-Degree(u1)-n1>; end if;
    if n1 gt 0 then return <u1,-v1,g-Degree(u1)-n1+1>; end if;
    return Adjust(<u1,-v1,Ceiling(3*g/2)-Degree(u1)+1>,f,V);
end function;

// Return the representative of the trivial divisor class for the curve
Identity := function(f)
    g := ValidateHyperellipticPolynomial(f);
    return <Parent(f)!1,Parent(f)!0,Ceiling(g/2)>;
end function;

// Given an integer n and D1=div[u1,v1,n1] compute D2=div[u2,v2,n2] with D2 ~ n*D
ScalarMultiplication := function(n,D1,f,V)
    ValidateReducedDivisor(D1,f);
    if n eq 0 then return Identity(f); end if;
    if n eq 1 then return D1; end if;
    if n lt 0 then return Negation($$(-n,D1,f,V),f,V); end if;
    bits:=IntegerToSequence(n,2);
    D2 := D1;
    i := #bits-1;
    while i gt 0 do
        D2 := Doubling(D2,f,V);
        if bits[i] ne 0 then D2 := Addition(D2,D1,f,V); end if;
        i-:=1;
    end while;
    return D2;
end function;

// Given a genus g > 1 and an odd prime power q returns a random monic separable poly of degree 2g+2 in GF(q)[x] with degree 2g+1 coeff zero.
RandomCurve := function(g,q)
    assert g gt 1 and q mod 2 eq 1;
    F:=GF(q);
    R<x>:=PolynomialRing(F);
    while true do
        f:= x^(2*g+2) + R![Random(F): i in [0..2*g]];
        if IsSeparable(f) then return f; end if;
    end while;
end function;

// CRT lift via Garner's algorithm (used by ModPolySqrt below)
// Given irred relatively prime polys A=[a_1,a_2,...,a_n] over F[x] and polys B=[b1,...,bn] compute b in F[x] st b ~ b_i mod a_i
CRTLift := function(A,B)
    assert #A ge 1;
    if #A eq 1 then return B[1]; end if;
    if #A eq 2 then return (B[1]*A[2]*InverseMod(A[2],A[1]) + B[2]*A[1]*InverseMod(A[1],A[2])) mod (A[1]*A[2]); end if;
    n := Ceiling(#A/2);
    return $$([&*A[[1..n]],&*A[[n+1..#A]]],[$$(A[[1..n]],B[[1..n]]),$$(A[[n+1..#A]],B[[n+1..#B]])]);
end function;

// given separable monic poly u in F[x] and poly f in F[x], attempt to compute v such that v^2 = f mod u with deg(v) < deg(u)
// Magma can do this when u is irreducible but does not appear to know how to take sqrts in an etale algebra
// 
EtaleIsSquare := function(f,u)
    assert GCD(u,Derivative(u)) eq 1;
    R:=Parent(u);
    A := [a[1] : a in Factorization(u)];
    B := [];
    for a in A do
        S,pi:=quo<R|a>;
        b,g := IsSquare(pi(f));
        if b then b := g^2 eq pi(f); end if;    // workaround Magma bug (as of V2.22-1)
        if not b then return false,_; end if;
        Append(~B,(-1)^Random(1)*Inverse(pi)(g));  // Randomize choice of sqrt (Magam does not do this)
    end for;
    v:=CRTLift(A,B);
    assert IsDivisibleBy(v^2-f,u);
    return true,v;
end function;
    
// Returns a (reasonably) random divisor class rep D=div[u,v,n] on the curve C:y^2=f(x) over a finite field
// We only generate divisors that have order 1 at all affine points, but these clearly generate Pic^0(C)
RandomDivisor := function(f)
    g := ValidateHyperellipticPolynomial(f);
    R:=Parent(f);  F:=BaseRing(f);  assert IsFinite(F);
    while true do
        u := R![Random(F) : i in [1..g+1]];
        if GCD(u,Derivative(u)) ne 1 then continue; end if;
        u := u/LeadingCoefficient(u);
        if Degree(u) eq 0 then v:=R!0; break; end if;
        b,v := EtaleIsSquare(f,u);
        if b then break; end if;
    end while;
    return <u,v,Random(g-Degree(u))>;
end function;

// Given D1=div[u1,v1,n1] and D2=div[u2,v2,n2] on genus 3 hyperelliptic curve return D3=div[u2,v2,0] with D3 ~ D1+D2
// Implements TypicalAddition algorithm, reverts to Addition for atypical divisors
TypicalAddition := function(D1,D2,f,V)
    assert Degree(f) eq 8 and IsMonic(f) and Coeff(f,7) eq 0;
    if Degree(D1[1]) ne 3 or Degree(D2[1]) ne 3 then return Addition(D1,D2,f,V); end if;

    u1:=D1[1]; v1:=D1[2]; u2:=D2[1]; v2:=D2[2];
    if GCD(u1,u2) ne 1 then return Addition(D1,D2,f,V); end if;
    st:=((v2-v1)*InverseMod(u1,u2)) mod u2;
    if Degree(st) ne 2 then return Addition(D1,D2,f,V); end if;
    c:=1/LeadingCoefficient(st);
    s:=c*st;  z:=s*u1;
    w:=ExactQuotient(f-v1^2,u1);
    u4:=ExactQuotient(s*(z+2*c*v1)-c^2*w,u2);
    vt4 := u4+v1+1/c*(z mod u4);
    if Coeff(vt4,3) eq 0 then return Addition(D1,D2,f,V); end if;
    u5:=1/(2*Coeff(vt4,3))*ExactQuotient(vt4^2-f,u4);
    v5:=vt4 mod u5;
    return <u5,v5,3-Degree(u5)>;
end function;

// Given D1=div[u1,v1,n1] and D2=div[u2,v2,n2] return D3=div[u2,v2,0] with D3 ~ D1+D2
// Implements explicit formula for TypicalAddition, reverts to Addition for atypical inputs
TypicalAdditionFormula := function(D1,D2,f,V)
    assert Degree(f) eq 8 and IsMonic(f) and Coeff(f,7) eq 0;
    if Degree(D1[1]) ne 3 or Degree(D2[1]) ne 3 then return Addition(D1,D2,f,V); end if;

    u10:=Coeff(D1[1],0); u11:=Coeff(D1[1],1); u12:=Coeff(D1[1],2); u20:=Coeff(D2[1],0); u21:=Coeff(D2[1],1); u22:=Coeff(D2[1],2);
    v10:=Coeff(D1[2],0); v11:=Coeff(D1[2],1); v12:=Coeff(D1[2],2); v20:=Coeff(D2[2],0); v21:=Coeff(D2[2],1); v22:=Coeff(D2[2],2);
    f4:=Coeff(f,4); f5:=Coeff(f,5); f6:=Coeff(f,6);

/*
// This code can be used to verify intermediate steps below

u1:=D1[1]; v1:=D1[2]; u2:=D2[1]; v2:=D2[2]; R:=Parent(f);
if GCD(u1,u2) ne 1 then return Addition(D1,D2,f,V); end if;
st:=((v2-v1)*InverseMod(u1,u2)) mod u2;
if Degree(st) ne 2 then return Addition(D1,D2,f,V); end if;
c:=1/LeadingCoefficient(st);
s:=c*st;  z:=s*u1;
w:=ExactQuotient(f-v1^2,u1);
u4:=ExactQuotient(s*(z+2*c*v1)-c^2*w,u2);
vt4 := u4+v1+1/c*(z mod u4);
if Coeff(vt4,3) eq 0 then return Addition(D1,D2,f,V); end if;
u5:=1/(2*Coeff(vt4,3))*ExactQuotient(vt4^2-f,u4);
v5:=vt4 mod u5;
assert <u5,v5,3-Degree(u5)> eq Addition(D1,D2,f,V);
*/

    // \\textbf{\\textsc{TypicalAddition}}: $\\div[u_5,v_5,n_5]\\ \\sim\\ \\div[u_1,v_1,0]+\\div[u_2,v_2,0]$ with $\\gcd(u_1,u_2)=1$.

    // 1. Compute $r := \\Res(u_1,u_2)$ and $i(x) = i_2x^2 + i_1x + i_0 := ru_1^{-1} \\bmod u_2$ (and $w_0:=u_{11}-u_{12}$).\\opcount{15M+12A}
    t1:=u10-u20; t2:=u11-u21; w0:=u12-u22; t3:=t2-u22*w0;
    t4:=t1-u21*w0; t5:=u22*t3-t4; t6:=u20*w0+u21*t3;
    i0 := t4*t5-t3*t6; i1 := w0*t6-t2*t5; i2 := w0*t4-t2*t3;
    r := t1*i0-u20*(t3*i2+w0*i1);

//assert r eq Resultant(u1,u2) and R![i0,i1,i2] eq r*InverseMod(u1,u2);

    // 2. Compute $q(x) = q_2x^2+q_1x+q_0 := r(v_2-v_1)\\,u_1^{-1} \\bmod u_2.$\\opcount{10M+30A}
    t1:=v20-v10; t2:=v11-v21; t3:=v12-v22; t4:=t2*i1; t5:=t1*i0; t6:=t3*i2; t7:=u22*t6;
    t8:=t4+t6+t7-(t2+t3)*(i1+i2); t9:=u20+u22; t10:=(t9+u21)*(t8-t6); t11:=(t9-u21)*(t8+t6);
    q0 := t5-u20*t8;
    q1 := t4-t5+(t11-t10)/2-t7+(t1-t2)*(i0+i1);
    q2 := t6-q0-t4+(t1-t3)*(i0+i2)-(t10+t11)/2;

//assert R![q0,q1,q2] eq r*st;

    // 3. Compute $t_1 := rq_2\\tilde{v}_{43}$ via \\eqref{eq:vt43add}, and $w_1 := c^{-1} = q_2/r,\\ \\ w_2 := c = r/q_2,\\ \\ w_3:=c^2,\\ \\ w_4:=(2\\tilde{v}_{43})^{-1}$.
    // \\hspace{12pt}Then compute $s(x) = x^2+s_1x + s_0:= c(v_2-v_1)\\,u_1^{-1} \\bmod u_2$ and $\\tilde{v}_{43}$.\\opcount{I+18M+5A}
    t1:=(r+q1)^2 + q2*(r*w0+q2*u21-q1*u22-q0); t2:=2*t1; t3:=r*q2;
    // If $t_2=0$ or $t_3=0$ then abort (revert to \\textsc{Addition}).
    if t2 eq 0 or t3 eq 0 then return Addition(D1,D2,f,V); end if;
    t4 := 1/(t2*t3);  t5:=t2*t4;  t6:=r*t5;
    w1 := t5*q2^2; w2 := r*t6; w3 := w2^2; w4 := t3^2*t4;
    s0 := t6*q0; s1 := t6*q1; 
    vt43 := t1*t5;

/*
assert t1 eq r*q2*Coeff(vt4,3) and w1 eq 1/c and w2 eq c and w3 eq c^2 and w4 eq 1/(2*vt43);
assert R![s0,s1,1] eq s;
assert vt43 eq Coeff(vt4,3);
assert vt43 eq (u12-u22)+c+2*s1+(u21+s1*(s1-u22)-s0)/c; // verify identity (1)
*/

    // 4. Compute $z(x) = x^5+z_4x^4+z_3x^3+z_2x^2+z_1x+z_0:=su_1$.\\opcount{4M+15A}
    t6:=s0+s1; t1:=u10+u12; t2:=t6*(t1+u11); t3:=(t1-u11)*(s0-s1); t4:=u12*s1;
    z0:=u10*s0; z1:=(t2-t3)/2-t4; z2:=(t2+t3)/2-z0+u10; z3:=u11+s0+t4; z4:=u12+s1;

//assert R![z0,z1,z2,z3,z4,1] eq z;

    // 5. Compute $u_4(x) = x^4+u_{43}x^3+u_{42}x^2+u_{41}x+u_{40}:= (s(z+2cv_1)-c^2(f-v_1^2)/u_1)/u_2.$\\opcount{14M+31A}
    u43 := z4+s1-u22;
    t0:=s1*z4; t1:=u22*u43;
    u42 := z3+t0+s0-w3-u21-t1;
    t2:=u21*u42; t3:=(u21+u22)*(u42+u43)-t1-t2; t4:=2*w2;
    t5:=t4*v12; t6:=s0*z3; t7:=(s0+s1)*(z3+z4)-t0-t6;
    u41 := z2+t7+t5+w3*u12-u20-t3;
    u40 := z1+s1*(t5+z2)+t6+t4*v11-w3*(f6+u12^2-u11)-u20*u43-t2-u22*u41;

//assert R![u40,u41,u42,u43,1] eq u4;

    // 6. Compute $\\tilde{v}_4(x) = x^4+\\tilde{v}_{43}x^3+\\tilde{v}_{42}x^2+\\tilde{v}_{41}x+\\tilde{v}_{40}:= -\\hat{v}_4 =  v_1 + u_4 + c^{-1}(z \\bmod u_4)$.\\opcount{6M+10A}
    t1:=u43-z4+w2;
    vt40 := v10 + w1*(z0+u40*t1);
    vt41 := v11 + w1*(z1-u40+u41*t1);
    vt42 := v12 + w1*(z2-u41+u42*t1);

//assert R![vt40,vt41,vt42,vt43,1] eq vt4;

    // 7. Compute $u_5(x) = x^3+u_{52}x^2+u_{51}x+u_{50} := (2\\tilde v_{43})^{-1}(\\tilde{v}_4^2-f)/u_4$. \\opcount{9M+17A}
    u52 := vt43/2+w4*(2*vt42-f6)-u43;
    u51 := w4*(2*(vt41+vt43*vt42)-f5)-u52*u43-u42;
    u50 := w4*(vt42^2+2*(vt40+vt43*vt41)-f4)-u51*u43-u52*u42-u41;

//assert R![u50,u51,u52,1] eq u5;

    // 8. Compute $v_5(x) = v_{52}x^2+v_{51}x+v_{50}:= \\tilde{v}_4 \\bmod u_5$.\\opcount{3M+6A}
    t1 := u52 - vt43;
    v50 := vt40 + t1*u50;
    v51 := vt41 - u50 + t1*u51;
    v52 := vt42 - u51 + t1*u52;

//assert R![v50,v51,v52] eq v5;

    // 9. Output $\\div[u_5,v_5,3-\\deg u_5]$.\\opcount{Total: I+79M+126A}

    u5 := Parent(f)![u50,u51,u52,1];
    v5 := Parent(f)![v50,v51,v52];
    return <u5,v5,3-Degree(u5)>;
end function;

// Given D1=div[u1,v1,n1] on genus 3 hyperelliptic return D3=div[u2,v2,0] with D3 ~ 2*D1
// Implements TypicalDoubling algorithm, reverts to Addition for atypical input
TypicalDoubling := function(D1,f,V)
    assert Degree(f) eq 8 and IsMonic(f) and Coeff(f,7) eq 0;
    if Degree(D1[1]) ne 3 then return Addition(D1,D1,f,V); end if;
    u1:=D1[1]; v1:=D1[2]; R:=Parent(f);
    if GCD(u1,v1) ne 1 then return Addition(D1,D1,f,V); end if;
    w:=ExactQuotient(f-v1^2,u1);
    wb:=w mod u1;
    st:=(wb*InverseMod(2*v1,u1)) mod u1;
    if Degree(st) ne 2 then return Addition(D1,D1,f,V); end if;
    c:=1/LeadingCoefficient(st);
    s:=c*st;  z:=s*u1;
    w:=ExactQuotient(f-v1^2,u1);
    u4:=s^2 - ExactQuotient(c^2*w-2*c*v1*s,u1);
    vt4 := u4+v1+1/c*(z mod u4);
    if Coeff(vt4,3) eq 0 then return Addition(D1,D1,f,V); end if;
    u5:=1/(2*Coeff(vt4,3))*ExactQuotient(vt4^2-f,u4);
    v5:=vt4 mod u5;
    return <u5,v5,3-Degree(u5)>;
end function;

// Given D1=div[u1,v1,n1] on genus 3 hyperelliptic return D3=div[u2,v2,0] with D3 ~ 2*D1
// Implements explicit formula for TypicalDoubling, reverts to Addition for atypical input
TypicalDoublingFormula := function(D1,f,V)
    assert Degree(f) eq 8 and IsMonic(f) and Coeff(f,7) eq 0;
    if Degree(D1[1]) ne 3 then return Addition(D1,D1,f,V); end if;
    u10:=Coeff(D1[1],0); u11:=Coeff(D1[1],1); u12:=Coeff(D1[1],2);
    v10:=Coeff(D1[2],0); v11:=Coeff(D1[2],1); v12:=Coeff(D1[2],2);
    f3:=Coeff (f,3); f4:=Coeff(f,4); f5:=Coeff(f,5); f6:=Coeff(f,6);

/*
// This code can be used to verify intermediate steps below
u1:=D1[1]; v1:=D1[2]; R:=Parent(f);
if GCD(u1,v1) ne 1 then return Addition(D1,D1,f,V); end if;
w:=ExactQuotient(f-v1^2,u1);
wb:=w mod u1;
st:=(wb*InverseMod(2*v1,u1)) mod u1;
if Degree(st) ne 2 then return Addition(D1,D1,f,V); end if;
c:=1/LeadingCoefficient(st);
s:=c*st;  z:=s*u1;
w:=ExactQuotient(f-v1^2,u1);
u4:=s^2 - ExactQuotient(c^2*w-2*c*v1*s,u1);
assert IsMonic(u4);
vt4 := u4+v1+1/c*(z mod u4);
if Coeff(vt4,3) eq 0 then return Addition(D1,D1,f,V); end if;
u5:=1/(2*Coeff(vt4,3))*ExactQuotient(vt4^2-f,u4);
v5:=vt4 mod u5;
assert <u5,v5,3-Degree(u5)> eq Addition(D1,D1,f,V);
*/

    // \\textbf{\\textsc{TypicalDoubling}}: $\\div[u_5,v_5,n_4]\\ \\sim\\ 2\\div[u_1,v_1,0]$ with $\\gcd(u_1,v_1)=1$.

    // 1. Compute $r := \\Res(u_1,v_1)$ and $i(x)=i_2x^2+i_1x+i_0 := rv_1^{-1} \\bmod u_1$ (and $w_0:=v_{11}-u_{12}v_{12}$).\\opcount{15M+10A}
    w0:=v11-u12*v12; t2:=v10-u11*v12; t3:=t2-u12*w0; t4:=u10*v12+u11*w0;
    i0:=t2*t3+w0*t4; i1:=-(v11*t3+v12*t4); i2:=v11*w0-v12*t2;
    r:=v10*i0-u10*(w0*i2+v12*i1);

//assert R![i0,i1,i2] eq r*InverseMod(v1,u1);

    // 3. Compute $p(x)=p_2x^2+p_1x+p_0 := \\overline{w} := (f-v_1^2)/u_1 \\bmod u_1$ (and $w_1:=u_{12}^2,\\ w_2:=u_{12}^2+f_6$).\\opcount{11M+24A}
    w1:=u12^2; t2:=2*u10; t3:=3*u11; w2:=w1+f6; t5:=2*t2-f5; t6:=2*u12; t7:=t3-w2;
    p2:=f5+t6*(t7-w1)-t2;
    p1:=u12*t5-v12^2+f4-u11*(2*f6-t3)-w1*(t7+t3);
    p0:=f3-u11*(w1*t6-t5)-t2*w2-u12*p1-2*v11*v12;

//assert R![p0,p1,p2] eq wb;

    // 4. Compute $q(x)= q_2x^2+q_1x+q_0 := r\\tilde{s} = r((f-v_1^2)/u_1)v_1^{-1} \\bmod u_1$ (and $w_3:=u_{10}+u_{11}+u_{12}$). \\opcount{10M+28A}
    t1:=i1*p1; t2:=i0*p0; t3:=i2*p2; t4:=u12*t3; t5:=(i1+i2)*(p1+p2)-t1-t3-t4; t6:=u10*t5;
    t7:=u10+u12; w3:=t7+u11; t9:=t7-u11; t10:=w3*(t3+t5); t11:=t9*(t5-t3);
    q0:=t2-t6;
    q1:=t4+(i0+i1)*(p0+p1)+(t11-t10)/2-t1-t2;
    q2:=t1+t6+(i0+i2)*(p0+p2)-t2-t3-(t10+t11)/2;

//assert R![q0,q1,q2] eq r*InverseMod(v1,u1)*wb mod u1;

    // 5. Compute $t_3 := 2rq_2\\tilde{v}_{43}$ via \\eqref{eq:vt43dbl}, and $w_4 := 1/\\tilde{v}_{43},\\ \\ w_5 := 1/c,\\ \\ w_6:=c,\\ \\ w_7:=c^2$.
    // \\hspace{12pt}Then compute $s(x)=x^2+s_1x+s_0 := q/(2r)$ made monic and $\\tilde{v}_{43}$.\\opcount{I+18M+7A}
    t0:=2*r; t1:=t0^2; t2:=q2^2;  t3:=t1-q0*q2+q1*(2*t0+q1-q2*u12)+t2*u11;
    // If $r=0$ or $q_2=0$ or $t_3=0$ then abort (revert to \\textsc{Addition}).
    if t0 eq 0 or q2 eq 0 or t3 eq 0 then return Doubling(D1,f,V); end if;
    t4:=1/(t0*q2*t3);  w4:=t1*t2*t4;  t5:=t3*t4; t6:=t0*t5; 
    /// t4=1/(4r^2q2^2vhat33), w4=1/(vhat33), t5=1/(2rq2), t6=1/(q2)
    w5:=t2*t5; w6:=t1*t5; w7:=w6^2; 
    /// w5=q2/(2r), w6=2r/q2, w7=4r^2/q2^2
    s0:=t6*q0; s1:=t6*q1;
    vt43:=t3*t5;

/*
assert t3 eq 2*r*q2*vt43 and w4 eq 1/vt43 and w5 eq 1/c and w6 eq c and w7 eq c^2;
assert R![s0,s1,1] eq s;
assert vt43 eq Coeff(vt4,3);
assert vt43 eq 2*s1+c+1/c*(s1*(s1-u12)-s0+u11); // verify identity (2)
*/

    // 6. Compute $z(x)=x^5+z_4x^4+z_3x^3+z_2x^2+z_1x+z_0 := su_1$.\\opcount{4M+12A}
    t1:=w3*(s0+s1); t2:=t9*(s0-s1); t3:=u12*s1; 
    z0:=s0*u10;  z1:=(t1-t2)/2-t3;  z2:=(t1+t2)/2-z0+u10;  z3:=u11+s0+t3;  z4:=u12+s1;  

//assert R![z0,z1,z2,z3,z4,1] eq z;

    // 7. Compute $u_4(x)=x^4+u_{43}x^3+u_{42}x^2+u_{41}x+u_{40} := s^2 - (c^2(f-v_1^2)/u_1-2csv_1)/u_1$.\\opcount{8M+14A}
    t1:=v12*w6;
    u43 := 2*s1;
    u42 := 2*s0+s1^2-w7;
    u41 := 2*(s0*s1+u12*w7+t1);
    u40 := s0^2+2*(w0*w6+s1*t1)-w7*(w2+2*(w1-u11));

//assert R![u40,u41,u42,u43,1] eq u4;

    // 8. $\\tilde{v}_4(x) = \\tilde{v}_{43}x^3+\\tilde{v}_{42}x^2+\\tilde{v}_{41}x+\\tilde{v}_{30} := -\\hat{v}_4 = v_1 + u_4 + c^{-1}(z \\bmod u_4)$.  \\opcount{6M+10A}
    /// note w6=1/w5
    w7:=u43-z4+w6;
    vt40 := v10 + w5*(z0+u40*w7);
    vt41 := v11 + w5*(z1-u40+u41*w7);
    vt42 := v12 + w5*(z2-u41+u42*w7);
    /// vhat33 := w5*(z3-u32+u33*w7); computed above

//assert R![vt40,vt41,vt42,vt43,1] eq vt4;

    // 9. $u_5(x)=x^3+u_{52}x^2+u_{51}x+u_{50} := (2\\tilde{v}_{43})^{-1}(\\tilde{v}_4^2-f)/u_4$.\\opcount{7M+17A}
    u52 := vt43/2+w4*(vt42-f6/2)-u43;
    u51 := vt42+w4*(vt41-f5/2)-u52*u43-u42;
    u50 := vt41+w4*((vt42^2-f4)/2+vt40)-u51*u43-u52*u42-u41;

//assert R![u50,u51,u52,1] eq u5;

    // 10. $v_5(x)=v_{52}x^2+v_{41}x+v_{50} := \\tilde{v}_4 \\bmod u_5$. \\opcount{3M+6A}
    t1:=u52-vt43;
    v50 := vt40 + t1*u50;
    v51 := vt41 - u50 + t1*u51;
    v52 := vt42 - u51 + t1*u52;

//assert R![v50,v51,v52] eq v5;

    // 11. Output $\\div[u_4,v_4,3-\\deg u_4]$.\\opcount{Total: I+82M+128A}

    u5 := Parent(f)![u50,u51,u52,1];
    v5 := Parent(f)![v50,v51,v52];
    return <u5,v5,3-Degree(u5)>;
end function;

// Given D1=div[u1,v1,n1] return D2 ~ -D1
// Implments TypicalNegation algorithm, reverts to Negation for atypical input
TypicalNegation := function(D1,f,V)
    assert Degree(f) eq 8 and IsMonic(f) and Coeff(f,7) eq 0;
    u1:=D1[1]; v1:=D1[2]; f6:=Coeff(f,6);
    R:=Parent(f);
    vt1:=v1-V+(V mod u1);
    if f6+2*Coeff(vt1,2) eq 0 then return Negation(D1,f,V); end if;
    u2:=(f6+2*Coeff(vt1,2))^-1*ExactQuotient(f-vt1^2,u1);
    v2:=vt1 mod u2;
    return <u2,v2,0>;
end function;

// Given D1=div[u1,v1,n1] return D2 ~ -D1
// Implments explicit formula for TypicalNegation, reverts to Negation for atypical input
TypicalNegationFormula := function(D1,f,V)
    assert Degree(f) eq 8 and IsMonic(f) and Coeff(f,7) eq 0;
    u10:=Coeff(D1[1],0); u11:=Coeff(D1[1],1); u12:=Coeff(D1[1],2);
    v10:=Coeff(D1[2],0); v11:=Coeff(D1[2],1); v12:=Coeff(D1[2],2);
    f3:=Coeff(f,3); f4:=Coeff(f,4); f5:=Coeff(f,5); f6:=Coeff(f,6);

/*
// This code may be used to verify intermediate steps below
u1:=D1[1]; v1:=D1[2];
R:=Parent(f);
vt1:=v1-V+(V mod u1);
if f6+2*Coeff(vt1,2) eq 0 then return Negation(D1,f,V); end if;
u2:=(f6+2*Coeff(vt1,2))^-1*ExactQuotient(f-vt1^2,u1);
v2:=vt1 mod u2;
assert <u2,v2,0> eq Negation(D1,f,V);
*/
    // \\textbf{\\textsc{TypicalNegation}}: $\\div[u_2,v_2,0]\\ \\sim \\ -\\div[u_1,v_1,0]$.
    
    // 1. Compute $\\tilde{v}_1(x)=-x^4+\\tilde{v}_{12}x^2+\\tilde{v}_{11}x+\\tilde{v}_{10} := v_1  - V + (V \\bmod u_1)$.\\opcount{3M+5A}
    vt12 := v12 - u11 + u12^2;
    vt11 := v11 - u10 + u11*u12;
    vt10 := v10 + u10*u12;

//assert R![vt10,vt11,vt12,0,-1] eq vt1;

    // 2. Compute $u_2(x) = x^3+u_{22}x^2+u_{21}x+u_{20} := (f_6+2\\tilde{v}_{12})^{-1}(f-\\tilde{v}_1^2)/u_1$.\\hspace{16pt}\\opcount{I+8M+14A}
    t1:=2*vt12; t2:=f6+t1;
    // If $t_1 = 0$ then abort (revert to \\textsc{Negation}).
    if t2 eq 0 then return Negation(D1,f,V); end if;
    t3:=1/t2;
    u22 := t3*(f5+2*vt11) - u12;
    u21 := t3*(f4+2*vt10-vt12^2) - u11 - u12*u22;
    u20 := t3*(f3-t1*vt11) - u10 - u11*u22 - u12*u21;

//assert R![u20,u21,u22,1] eq u2;

    // 3. Compute $v_2(x) = v_{22}x^2+v_{21}x+v_{20}:= \\tilde{v}_1 \\bmod u_2$.\\opcount{3M+5A}
    v22 := vt12 - u22^2 + u21;
    v21 := vt11 - u21*u22 + u20;
    v20 := vt10 - u20*u22;

//assert R![v20,v21,v22] eq v2;

    // 4. Output $\\div[u_2,v_2,0]$.\\opcount{Total: I+14M+24A}

    return <Parent(f)![u20,u21,u22,1], Parent(f)![v20,v21,v22], 0>;
end function;

// The Test* functions below can be used to verify the excution of the algorithms above

TestRandomScalarMultiplication := function(g,q,n)
    f := RandomCurve(g,q);
    printf "Testing %o random scalar multiplications on the curve y^2 = %o\n", n,f;
    V := Precompute(f);
    i := 0;
    while true do
        N1 := Random(1000);
        N2 := Random(1000);
        D := RandomDivisor(f);
        if Addition(ScalarMultiplication(N1,D,f,V),ScalarMultiplication(N2,D,f,V),f,V) ne ScalarMultiplication(N1+N2,D,f,V) then
            print "Failed on scalar multiplication of %D on curve y^2=%o\n", D, f;
            print "%o*D + %o*D != %o*D\n", N1, N2, N1+N1;
            return false;
        end if;
        printf ".";
        i +:= 1;
        if i eq n then break; end if;
    end while;
    return true;
end function;

TestScalarMultiplication := function(g,q,n)
    if q^g gt 1000000 then print "Computing #Jac(C) could take a while...hit ctrl-c to abort."; end if;
    f := RandomCurve(g,q);
    N := Evaluate(LPolynomial(HyperellipticCurve(f)),1);
    printf "Testing scalar multiplication by #Jac(C)=%o on the curve y^2 = %o\n", N,f;
    V := Precompute(f);
    i := 0;
    I := Identity(f);
    while true do
        D:=RandomDivisor(f);
        if D eq I then continue; end if;
        if ScalarMultiplication(N,D,f,V) ne I then
            printf "Failed on D = %o\n on curve y^2=%o\n", D, f;
            printf "%o*D != 0\n", N;
            return false;
        end if;
        if ScalarMultiplication(N+1,D,f,V) eq I then
            printf "Failed on D = %o\n on curve y^2=%o\n", D, f;
            printf "(%o+1)*D = 0\n", N;
            return false;
        end if;
        printf ".";
        i +:= 1;
        if i eq n then break; end if;
    end while;
    print "";
    return true;
end function;
        
TestAdditionFormula := function(q,n)
    g:=3;
    f:=RandomCurve(g,q);
    printf "Testing %o typical additions of random generic divisors on the curve y^2 = %o\n", n,f;
    V:=Precompute(f);
    i := 0;
    while true do
        D1:=RandomDivisor(f);
        D2:=RandomDivisor(f);
        if Degree(D1[1]) lt g or Degree(D2[1]) lt g or GCD(D1[1],D2[1]) ne 1 then continue; end if;
        D3:=Addition(D1,D2,f,V);
        D4:=TypicalAdditionFormula(D1,D2,f,V);
        if D3 ne D4 then
            printf "Failed adding D1 = %o, D2 = %o\n on curve y^2=%o\n", D1, D2, f;
            printf "%o != %o\n", D4, D3;
            return false;
        end if;
        printf ".";
        i +:= 1;
        if i eq n then break; end if;
    end while;
    print "";
    return true;
end function;

TestDoublingFormula := function(q,n)
    g:=3;
    f:=RandomCurve(g,q);
    printf "Testing %o typical doublings of a random generic divisor on the curve y^2 = %o\n", n,f;
    V:=Precompute(f);
    i := 0;
    while true do
        D1:=RandomDivisor(f);
        if Degree(D1[1]) lt g or GCD(D1[1],D1[2]) ne 1 then continue; end if;
        D3:=Doubling(D1,f,V);
        D4:=TypicalDoublingFormula(D1,f,V);
        if D3 ne D4 then
            printf "Failed doubling D1 = %o on curve y^2=%o\n", D1,f;
            printf "%o != %o\n", D4, D3;
            return false;
        end if;
        printf ".";
        i +:= 1;
        if i eq n then break; end if;
    end while;
    print "";
    return true;
end function;

TestNegationFormula := function(q,n)
    g:=3;
    f:=RandomCurve(g,q);
    printf "Testing %o typical negions of a random generic divisor on the curve y^2 = %o\n", n,f;
    V:=Precompute(f);
    i := 0;
    while true do
        D1:=RandomDivisor(f);
        if Degree(D1[1]) lt g then continue; end if;
        D3:=Negation(D1,f,V);
        D4:=TypicalNegationFormula(D1,f,V);
        if D3 ne D4 then
            printf "Failed negating D1 = %o on curve y^2=%o\n", D1,f;
            printf "%o != %o\n", D4, D3;
            return false;
        end if;
        printf ".";
        i +:= 1;
        if i eq n then break; end if;
    end while;
    print "";
    return true;
end function;

// The functions below are simple schoolbook implementations from GG Ch.1
// of polynomial addition, subtraction, multiplication, division, and xgcd
// These are used solely for obtaining field operation counts for the
// unoptimized algorithm.

add := procedure(~s,a,b,~cnts);
    n:=Degree(a); m:=Degree(b); x:=Parent(a).1;
    cnts[3] +:= n lt m select n+1 else m+1;
    s := a+b;
end procedure;

sub := procedure(~d,a,b,~cnts);
    n:=Degree(a); m:=Degree(b); x:=Parent(a).1;
    cnts[3] +:= n lt m select n+1 else m+1;
    d := a-b;
end procedure;

mul := procedure(~p,a,b,~cnts)
    n:=Degree(a); m:=Degree(b); x:=Parent(a).1;
    p:=a*b;
    cnts[2] +:= (n+(IsMonic(a) select 0 else 1))*(m+(IsMonic(b) select 0 else 1));
    cnts[3] +:= (n+1)*(m+1) - (m+n+1);
end procedure;

ediv := procedure(~q,~r,a,b,~cnts)
    n:=Degree(a); m:=Degree(b); x:=Parent(a).1; c:=IsMonic(b) select 0 else 1;
    q:= Parent(a)!0; r:=a;  u:=1/LeadingCoefficient(b); cnts[1] +:= c;
    i:=n-m;
    while i ge 0 do
        if Degree(r) eq m+i then
            s := LeadingCoefficient(r)*u; r:= r-s*x^i*b; cnts[2] +:= m+1+c; cnts[3] +:=m+1;
            q +:= s*x^i;
        end if;
        i -:= 1;
    end while;
    assert q*b+r eq a;
end procedure;

xgcd := procedure(~d,~a,~b,f,g,~cnts)
    R:=Parent(f);
    r:=[f,g]; s:=[R|1,0]; t:=[R|0,1];
    i:=2; j:=1; q:=R!0; z:=R!0; w:=R!0;
    while true do
        ediv(~q,~z,r[j],r[i],~cnts);
        mul(~w,q,r[i],~cnts); sub(~r[j],r[j],w,~cnts);
        if r[j] eq 0 then break; end if;
        mul(~w,q,s[i],~cnts); sub(~s[j],s[j],w,~cnts);
        mul(~w,q,t[i],~cnts); sub(~t[j],t[j],w,~cnts);
        k:=i; i:=j; j:=k;
    end while;
    assert s[i]*f+t[i]*g eq r[i];
    d:=r[i]; a:=s[i]; b:=t[i];
    if not IsMonic(d) then
        e:=1/LeadingCoefficient(d); cnts[3]+:=1;
        d*:=e; a*:=e; b*:=e; cnts[2] +:=Degree(d)+Degree(a)+Degree(b)+2;
    end if;
    assert d eq GCD(f,g) and d eq a*f+b*g;
end procedure;

// Procedure to count field operations for addition of typical divisors
// using the standard Addition algorithm
TypicalAdditionUnoptimizedCost := function(D1,D2,f,V)
    g:=ValidateHyperellipticPolynomial(f);
    u1:=D1[1]; v1:=D1[2]; u2:=D2[1]; v2:=D2[2];
    assert Degree(u1) eq g and Degree(u2) eq g and GCD(u1,u2) eq 1;
    R:=Parent(f);
    cnts := [0,0,0];
    d:=R!0; c1:=R!0; c2:=R!0;
    xgcd(~d,~c1,~c2,u1,u2,~cnts);
    assert d eq 1;
    u3:=R!0; v3:=R!0; s1:=R!0; s2:=R!0; t1:=R!0; t2:=R!0;
    mul(~u3,u1,u2,~cnts);
    mul(~s1,u1,v2,~cnts); mul(~s2,c1,s1,~cnts);
    mul(~t1,u2,v1,~cnts); mul(~t2,c2,t1,~cnts);
    add(~t1,s2,t2,~cnts);
    v3:=R!0;
    ediv(~t2,~v3,t1,u3,~cnts);
    assert <u3,v3,0> eq Compose(D1,D2,f);
    mul(~t1,v3,v3,~cnts); sub(~t2,f,t1,~cnts);
    ediv(~t1,~s1,t2,u3,~cnts);
    assert s1 eq R!0;
    c:=R!(1/LeadingCoefficient(t1)); cnts[1]+:=1;
    u4:=R!0;
    mul(~u4,c,t1,~cnts);
    v4:=R!0;
    ediv(~t1,~v4,-v3,u4,~cnts); // don't charge for negating v3
    assert <u4,v4,1> eq Reduce(Compose(D1,D2,f),f);
    sub(~t1,v4,V,~cnts);
    ediv(~s2,~s1,V,u4,~cnts);
    vh4:=R!0;
    add(~vh4,t1,s1,~cnts);
    mul(~t1,vh4,vh4,~cnts);
    sub(~t2,f,t1,~cnts);
    ediv(~t1,~s1,t2,u4,~cnts);
    assert s1 eq R!0;
    c:=R!(1/LeadingCoefficient(t1)); cnts[1]+:=1;
    u5:=R!0;
    mul(~u5,c,t1,~cnts);
    ediv(~t1,~v5,-vh4,u5,~cnts); // don't charge for negating vh4
    D3:=<u5,v5,3-Degree(u5)>;
    assert D3 eq Addition(D1,D2,f,V);
    return cnts;
end function;

// Procedure to count field operations for doubling a typical divisor
// using the standard Addition algorithm
TypicalDoublingUnoptimizedCost := function(D1,f,V)
    g:=ValidateHyperellipticPolynomial(f);
    u1:=D1[1]; v1:=D1[2];
    assert Degree(u1) eq g and GCD(u1,v1) eq 1;
    R:=Parent(f);
    cnts := [0,0,0];
    d:=R!0; c1:=R!0; c3:=R!0;
    xgcd(~d,~c1,~c3,u1,2*v1,~cnts);
    assert d eq 1;
    u3:=R!0; v3:=R!0; s1:=R!0; s2:=R!0; t1:=R!0; t2:=R!0;
    mul(~u3,u1,u1,~cnts);
    mul(~s1,u1,v1,~cnts); mul(~s2,c1,s1,~cnts);
    mul(~t2,v1,v1,~cnts); add(~t1,t2,f,~cnts); mul(~t2,c3,t1,~cnts);
    add(~t1,s2,t2,~cnts);
    v3:=R!0;
    ediv(~t2,~v3,t1,u3,~cnts);
    assert <u3,v3,0> eq Compose(D1,D1,f);
    mul(~t1,v3,v3,~cnts); sub(~t2,f,t1,~cnts);
    ediv(~t1,~s1,t2,u3,~cnts);
    assert s1 eq R!0;
    c:=R!(1/LeadingCoefficient(t1)); cnts[1]+:=1;
    u4:=R!0;
    mul(~u4,c,t1,~cnts);
    v4:=R!0;
    ediv(~t1,~v4,-v3,u4,~cnts); // don't charge for negating v3
    assert <u4,v4,1> eq Reduce(Compose(D1,D1,f),f);
    sub(~t1,v4,V,~cnts);
    ediv(~s2,~s1,V,u4,~cnts);
    vh4:=R!0;
    add(~vh4,t1,s1,~cnts);
    mul(~t1,vh4,vh4,~cnts);
    sub(~t2,f,t1,~cnts);
    ediv(~t1,~s1,t2,u4,~cnts);
    assert s1 eq R!0;
    c:=R!(1/LeadingCoefficient(t1)); cnts[1]+:=1;
    u5:=R!0;
    mul(~u5,c,t1,~cnts);
    ediv(~t1,~v5,-vh4,u5,~cnts); // don't charge for negating vh4
    D3:=<u5,v5,3-Degree(u5)>;
    assert D3 eq Addition(D1,D1,f,V);
    return cnts;
end function;

// The code below is used to transform the code above for
// Typical*Formula into latex source

IsAlpha := function(c)
    return (c ge "a" and c le "z") or (c ge "A" and c le "Z");
end function;
        
IsDigit := function(c)
    return c ge "0" and c le "9";
end function;
        
LatexCode := procedure(str)
    lines := Split(str,"\n");
    incomment := false;
    for x in lines do
        if #x gt 0 and x[1] ne " " then continue; end if;
        line := StripWhiteSpace(x);
        if #line eq 0 then continue; end if;
        if Substring(line,1,3) eq "///" then continue; end if;
        if Substring(line,1,2) eq "//" then
            if not incomment and not "abort" in line then printf "\\hline\n"; incomment:=true; end if;
            printf "%o\\\\\n", Substring(x,8,#x-7);
            continue;
        end if;
        if Substring(line,1,2) eq "if" then continue; end if;
        if incomment then
            incomment := false;
            printf "\\hline\n";
        end if;
        printf "$";
        i:=1;
        while i le #line do
            if line[i] eq "*" then i+:=1; continue; end if;
            if not IsAlpha(line[i]) then 
                if line[i] eq ";" and i lt #line then
                    printf ";\\quad "; i +:= 1;
                else
                    printf line[i]; i+:=1;
                end if;
                continue;
            end if;
            if i lt #line-1 and line[i+1] eq "h" then
                printf "\\hat{%o}", line[i]; i +:= 2;
            elif i lt #line-1 and line[i+1] eq "t" then
                printf "\\tilde{%o}", line[i]; i +:= 2;
            else
                printf line[i]; i +:= 1;
            end if;
            if i le #line and IsDigit(line[i]) then
                printf "_{%o", line[i]; i+:=1;
                while i le #line and IsDigit(line[i]) do printf line[i]; i+:=1; end while;
                printf "}";
            end if;
        end while;
        print "$\\\\";
    end for;
    print "\\hline";
end procedure;
