#################################
########### AGRO 2023 ###########
#################################

## Importations

import math as m

## PARTIE 2

## 1. Coefficients binomiaux

def C(n,k) :
    Num = 1
    for i in range(2,n+1) :
        Num = Num*i
    Den = 1
    for i in range(2,k+1) :
        Den = Den*i
    for i in range(2,n-k+1) :
        Den = Den*i
    return Num/Den

# la variable locale 'Num' est initialisée à 1
# et la boucle des lignes 3 et 4 la multiplie
# successivement par 2, 3 ... n
# donc à l'issue de cette boucle, Num vaut n!
# De même, à l'issue de la boucle lignes 6 et 7,
# 'Den' vaut k! puis à l'issue de la boucle
# lignes 8 et 9, Den vaut k!*(n-k)!
# On renvoie le quotient Num/Den, soit n!/(k!(n-k)!)
# ce qui est le coefficient binomial k parmi n.

# Il y a 3 boucles, de taille n-1, k-1 et n-k-1
# donc au total 2n-3 opérations, soit environ 2n.
# De plus, les flottants 'Num' et 'Den' (surtout 'Num')
# risquent d'être très grands (risque de dépassement
# de capacité).

## 2. Optimisation des coefficients binomiaux

def Copt(n,k) :
    Num = 1
    if k > n-k :
        k = n-k
    for i in range(n-k+1,n+1) :
        Num *= i
    Den = 1
    for i in range(2,k+1) :
        Den *= i
    return Num/Den

# on a ici simplifié n!/(n-k)!, ce qui donne le produit
# de tous les entiers de n-k+1 à n.
# on a également utilisé la relation de symétrie :
# k parmi n = (n-k) parmi n, donc on choisit le plus
# petit entre k et n-k.

## 3. Fonction 'Feuilles'

def Feuilles(n,u,s0,setoile) :
    F = []
    for k in range(n+1) :
        F.append(max(u**(2*k-n)*s0-setoile , 0))
    return F

# Le fonction 'Feuilles' renvoie une liste de flottants
# de taille (n+1). Les éléments de la liste renvoyées
# sont les termes (u**(2k-n)s0 - s*) lorsque ceux-ci
# sont positifs, et zéro sinon.

## 4. Fonction 'call'

def call(n,u,s0,setoile,p) :
    cn = 0
    L = Feuilles(n,u,s0,setoile)
    for k in range(n+1) :
        cn += L[k]*Copt(n,k)*p**k*(1-p)**(n-k)
    return cn


## PARTIE 3

## 1. Méthode de Newton

nmax = 100

def Newton(u0, epsilon, F, Fprime) :
    u = u0
    v = u - F(u)/Fprime(u)
    compteur = 1
    while abs(v-u) > epsilon and compteur < nmax :
        u = v
        v = u - F(u)/Fprime(u)
        compteur += 1
    return u

# pour tester
def f(t) :
    return (t-3)**2

def fprime(t) :
    return 2*(t-3)

## 2. Choix de epsilon

# on peut choisir epsilon = 10**(-10) pour avoir
# une grande précision. En effet, la méthode de Newton
# quand elle s'applique, converge très rapidement.

## 3. Cas d'erreur

# f'(un) = 0 conduit à une division par zéro, donc
# à un message d'erreur. On peut modifier ainsi :

def Newton_2(u0, epsilon, F, Fprime) :
    u = u0
    if Fprime(u) == 0 :
        return 'méthode non applicable'
    v = u - F(u)/Fprime(u)
    compteur = 1
    while abs(v-u) > epsilon and compteur < nmax :
        u = v
        if Fprime(u) == 0 :
            return 'méthode non applicable'
        v = u - F(u)/Fprime(u)
        compteur += 1
    return u

# rappel : un test d'égalité entre flottants n'est
# pas toujours fiable.


## 4. Fonction mystère

def mystere(F,Fprime,s0,setoile,epsilon) :
    a = n
    u = (setoile/s0)**(1/(2*a-n))
    while F(u) < 0 and a >= m.floor(1+n/2) :
        a -= 1
        u = (setoile/s0)**(1/(2*a-n))
    a += 1
    u = (setoile/s0)**(1/(2*a-n))
    u = Newton(u,epsilon,F,Fprime)
    return u

# La fonction 'mystere' renvoie une valeur approchée
# à epsilon près d'une solution de l'équation F(u) = 0
# donc une valeur de u telle que : call(u) = c*
# On peut aussi répondre que la fonction 'mystere'
# renvoie un message d'erreur : name 'n' is not defined...

## 5. Volatilité

def Volatilite(F,Fprime,s0,setoile,epsilon) :
    u = mystere(F,Fprime,s0,setoile,epsilon)
    return m.log(u)