def decomposition_g0_pth_power(fct):
    C = fct.curve
    Fxy, Rxy, xy, y = C.fct_field
    if fct.function in Rxy:
        return (fct, 0*C.x)
    '''Decompose fct as g0 + A^p, if possible. Output: (g0, A).'''
    omega = fct.diffn().regular_form()
    g0 = omega.integral()
    A = (fct - g0).pth_root()
    return (g0, A)

def decomposition_g0_p2th_power(fct):
    '''Decompose fct as g0 + A^(p^2), if possible. Output: (g0, A).'''
    C = fct.curve
    p = C.characteristic
    g0, A = decomposition_g0_pth_power(fct)
    A0, A1 = decomposition_g0_pth_power(A)
    return (g0 + A0^p, A1)

def decomposition_omega0_hpdh(omega):
    '''Decompose omega = (regular on U0) + h^(p-1) dh, so that Cartier(omega) = (regular on U0) + dh.
    Result: (regular on U0, h)'''
    C = omega.curve
    if omega.is_regular_on_U0():
        return (omega, 0*C.x)
    omega1 = omega.cartier().cartier()
    omega1 = omega1.inv_cartier().inv_cartier()
    fct = (omega.cartier() - omega1.cartier()).integral()
    return (omega1, fct)

def decomposition_omega8_hpdh(omega, prec = 50):
    '''Decompose omega = (regular on U8) + h^(p-1) dh, so that Cartier(omega) = (regular on U8) + dh.
    Result: (regular on U8, h)'''
    C = omega.curve
    g = C.genus()
    Fxy, Rxy, x, y = C.fct_field
    F = C.base_ring
    p = C.characteristic
    C = omega.curve
    if omega.is_regular_on_Uinfty():
        return (omega, 0*C.x)
    Rt.<t> = LaurentSeriesRing(F)
    omega_analytic = Rt(laurent_analytic_part(omega.expansion_at_infty(prec = prec)))
    Cv = C.uniformizer()
    v = Fxy(Cv.function)
    omega_analytic = Fxy(omega_analytic(t = v))
    omega_analytic = superelliptic_function(C, omega_analytic)*Cv.diffn()
    omega8 = omega - omega_analytic
    dh = omega.cartier() - omega8.cartier()
    h = dh.integral()
    return (omega8, h)

def decomposition_g8_pth_power(fct, prec = 50):
    '''Decompose fct as g8 + A^p, if possible. Output: (g8, A).'''
    C = fct.curve
    F = C.base_ring
    Rt.<t> = LaurentSeriesRing(F)
    Fxy, Rxy, x, y = C.fct_field
    if fct.expansion_at_infty().valuation() >= 0:
        return (fct, 0*C.x)
    A = laurent_analytic_part(fct.expansion_at_infty(prec = prec))
    Cv = C.uniformizer()
    v = Cv.function
    A = A(t = v)
    A = superelliptic_function(C, A)
    A = A.pth_root()
    g8 = fct - A^p
    return (g8, A)

def decomposition_g8_p2th_power(fct):
    '''Decompose fct as g8 + A^(p^2), if possible. Output: (g8, A).'''
    g0, A = decomposition_g8_pth_power(fct)
    A0, A1 = decomposition_g8_pth_power(A)
    return (g0 + A0^p, A1)