---

def pfdintegral(f,x):
    g=(f).numerator()
    h=(f).denominator()
    
    if AA[x](h).degree()==1:
        return g*ln(abs(x+h.subs(x=0)))
    else:
        b0=g.subs(x=0)
        b1=diff(g,x)
        a0=h.subs(x=0)
        a1=diff(h,x).subs(x=0)
        s=sqrt(-a1^2 + 4*a0)
        return 1/2*b1*log(a1*x + x^2 + a0) - (a1*b1 - 2*b0)*arctan((a1 + 2*x)/s)/s
        
def ostrogradski(f,x):
    # f -- правильная дробь
    # g -- числитель f
    K=FractionField(ZZ[x])
    g=K(f).numerator()
    # h -- знаменатель f
    h=K(f).denominator()
    # Раскладываем знаменатель на множители
    L=list(ZZ[x](h).factor())
    
    # Составляем знаменатель алгебраической части
    hh=prod([p^(m-1) for (p,m) in L])
    n=hh.degree()
    # Если знаменатель алг. части -- константа, то
    # f -- уже log часть
    if n==0:
        return [0,f]
    else:
    # Составляем числитель алгебраической части 
	# с неопределенными коэффициентами
        A=var(['A'+str(i) for i in range(n)])
        gg=sum([A[i]*x^i for i in range(n)])
    
    # Составляем знаменатель log части
        h3=ZZ[x](prod([p for (p,m) in L]))
    # Составляем числитель log части 
    # с неопределенными коэффициентами
        m=h3.degree()
        B=var(['B'+str(i) for i in range(m)])
        g3=sum([B[i]*x^i for i in range(m)])

    # Составляем выражение, которое должно быть равно нулю
        F=ZZ[A+B][x]((diff(SR(gg/hh),x)+SR(g3/h3)-f).numerator())
    # Составляем СЛАУ для определения коэффициентов и решаем ее
        S=tsolve(triangulation([QQ[A+B](eq) for eq in F.coefficients()]))
    # Список: алгебраическая часть, log часть
    return [(gg).subs(S)/hh, (g3).subs(S)/h3]
def triangulation(S):
    T=[]
    while S!=[]:
        m=max([s.lm() for s in S])
        L = [s for s in S if s.lm() == m]
        S = [s for s in S if s.lm() < m]
        f=L[0]
        T.append(f)
        for g in L[1:]:
            g=f.lc()*g-g.lc()*f
            if g!=0:
                S.append(g)
        D=[s for s in S if s.degree()==0]
        if D!=[]:
            return D
    return T
    
def tsolve(T):
    T.reverse()
    D={}
    while T!=[]:
        g=T[0]
        D[g.lm()] = -(g-g.lt())/g.lc()
        T=[t.subs(D) for t in T[1:]]
    return D

def ostrogradsky_integral(f, x):
    num = QQ[x](f.numerator())
    den = QQ[x](f.denominator())
    q, r = num.quo_rem(den)  # целая часть и остаток
    
    # Интегрируем полиномиальную часть
    poly_int = sum(q[i] * x^(i+1) / (i+1) for i in range(q.degree()+1))
    
    # Если остаток нулевой, сразу возвращаем
    if r == 0:
        return poly_int
    
    # Правильная дробь
    proper = r / den
    alg_part, log_part = ostrogradski(proper, x)
    log_aa = FractionField(AA[x])(log_part)
    _, partials = log_aa.partial_fraction_decomposition()
    
    log_int = sum(pfdintegral(p, x) for p in partials)
    
    return poly_int + alg_part + log_int