4234.sagews

F(x) = integrate(2*sqrt(1-x^2)/pi,x)

Fpaul = pi*F/2

def skip_factorial(n):
    assert n in Integers()
    if n <= 0: 
        return 1
    return prod(i for i in [1..n] if i % 2 == n % 2)
    
def Fpaul_coeff(n):
    if is_even(n):
        return 0
    if n == 1:
        return skip_factorial(n-4)*skip_factorial(n-2)/factorial(n)
    else:
        return -skip_factorial(n-4)*skip_factorial(n-2)/factorial(n)

assert Fpaul.taylor(x,0,100).coefficients() == [[Fpaul_coeff(n),n] for n in range(100) if n % 2 == 1]

T.<t> = PowerSeriesRing(Rationals())

def Fpaul_series(MaxExp):
    return sum(Fpaul_coeff(n)*t^n for n in range(MaxExp))+O(t^MaxExp)

# Compare to series at http://www.wolframalpha.com/input/?i=1%2F2+%28x+sqrt%281-x^2%29%2Bsin^%28-1%29%28x%29%29
Fpaul_series(10)

t - 1/6*t^3 - 1/40*t^5 - 1/112*t^7 - 5/1152*t^9 + O(t^10)

Fpaul(1)

1/4*pi

plot(Fpaul_series(10),(t,-1,1))

Fpaul_reversion = Fpaul_series(100).reversion()

reversion_coefficients = list(enumerate(Fpaul_reversion.list()))

def partial_reversion(MaxPrec):
    return lambda x: sum(coeff*x^exponent for (exponent,coeff) in reversion_coefficients[:MaxPrec])

print Fpaul_series(25)
print partial_reversion(25)(t) + O(t^25)
print partial_reversion(25)(Fpaul_series(25))

t - 1/6*t^3 - 1/40*t^5 - 1/112*t^7 - 5/1152*t^9 - 7/2816*t^11 - 21/13312*t^13 - 11/10240*t^15 - 429/557056*t^17 - 715/1245184*t^19 - 2431/5505024*t^21 - 4199/12058624*t^23 + O(t^25) t + 1/6*t^3 + 13/120*t^5 + 493/5040*t^7 + 37369/362880*t^9 + 4732249/39916800*t^11 + 901188997/6227020800*t^13 + 240798388357/1307674368000*t^15 + 6611433892597/27360571392000*t^17 + 39504564917358001/121645100408832000*t^19 + 22726779729476308093/51090942171709440000*t^21 + 941059359881411415629/1520706866993233920000*t^23 + O(t^25) t + O(t^25)

# A comparer avec ta courbe rouge, tu verras que c est pas les memes
plot(partial_reversion(45),(-Fpaul(1),Fpaul(1)))

X = Fpaul(1)
# Ceci est une expansion autour de 0, donc pres de Fpaul(1), on met plus de terms pour etre proche. 
plot(partial_reversion(45),(-X,X)) + plot(partial_reversion(5),(-X,X), color = "red") + plot(partial_reversion(99),(-X,X), color = "green")