import numpy as np
import numpy.linalg as alg
import scipy.special as sc
import scipy.integrate as scint
import numpy.random as rd
import matplotlib.pyplot as plt
from math import *

## PLANCHE 1
    
# Q1

def Res(P,Q) :
    M = np.zeros((4,4))
    for i in range(2) :
        M[i,i] = P[0]
        M[i+1,i] = P[1]
        M[i+2,i] = P[2]
    for j in range(2) :
        M[j,2+j] = Q[0]
        M[j+1,2+j] = Q[1]
        M[j+2,2+j] = Q[2]
    print(M)
    return alg.det(M)

print(Res([2,-3,1],[3,-4,1]))
print(Res([2,-3,1],[12,-7,1]))
print(Res([1,-2,1],[-2,0,1]))
# Le résultant de P et Q semble valoir 0 quand P et Q ont une racine commune.


## PLANCHE 2

# Q1

def mat(n) :
    M = np.zeros((n+1,n+1))
    for i in range(n+1) :
        for j in range(n+1) :
            M[i,j] = 1/(n-i+j+1)*sc.comb(n,i)
    return M

def diag() :
    for i in range(2,7) :
        print(alg.eig(mat(i))[0])
        
diag()
# Valeurs propres distinctes dans chaque cas : An est diagonalisable
    
# Q2a

def prod_scal(x,y,n) :
    S = 0
    for i in range(n+1) :
        for j in range(n+1) :
            S += x[i,0]*y[j,0]/(i+j+1)
    return S

# Q2b

def ortho(n) :
    A = alg.eig(mat(n))[1]
    for i in range(n+1) :
        for j in range(i) :
            print(prod_scal(A[:,i:i+1],A[:,j:j+1],n))

# Les sous-espaces propres sont de dimension 1, on calcule les produits scalaires
# des vecteurs propres, deux à deux.

ortho(8)

## PLANCHE 3

#Q1

def Test1(beta,n) :
    return -n*(n**(beta)/(n+1)**(beta)-1)

print(Test1(4,500))

# Il semble que lambda = beta
# On connaît le critère de convergence d'une série de Riemann...

# Q2

def Test2(n) :
    return -n*((n*log(n)**2/((n+1)*log(n+1)**2))-1)

print(Test2(100000))
# Il semble que lambda = 1

def ConvSerie(N):
    S = 0
    for n in range(2,N) :
        S += 1/(n*log(n)**2)
    return S

print(ConvSerie(100))

## PLANCHE 5

# Q1a

def f(x) :
    return sin(x)**2/x**2

print(scint.quad(f,0,10000))
# L'intégrale semble converger

# Q1b

def integrale(m,seuil) :
    def f(x) :
        return sin(x)**(2*m+1)/x
    return scint.quad(f,0,seuil)

print(integrale(5000,500))
# Les intégrales semblent converger et tendre vers 0 quand m tend vers l'infini
    
## PLANCHE 6

# Q1a

def M(n,p) :
    mat = np.zeros((n,n))
    for i in range(n) :
        for j in range(n) :
            if i == j :
                mat[i,j] = 1
            else :
                mat[i,j] = rd.binomial(1,p)
    return mat

print(M(4,0.5))

# Q1b

def NbIsoles(n,p) :
    mat = M(n,p)
    N = 0
    for i in range(n) :
        L = mat[i:i+1,:]
        S = 0
        for j in range(n) :
            S += L[0,j]
        N += int(S == 1)
    return N

print(NbIsoles(4,0.5))

# Q1c

def Appro(n,p) :
    S = 0
    for i in range(3000) :
        S += NbIsoles(n,p)
    return S/3000

print(Appro(4,0.5))

# Q1d

def Graphe(p) :
    N = range(2,11)
    O = []
    for i in range(2,11) :
        print(i)
        O += [ log(Appro(i,p)/i) ]
        print(i,"concat")
    plt.plot(N,O,'-*r')
    plt.show()

Graphe(0.7)
# On obtient des tracés proches d'une droite de coefficient directeur négatif.
# Ainsi E(S_n) est de la forme K*n*exp(-a*n) avec a > 0.

## PLANCHE 7

# Q1

def Conjecture(p) :
    S = 0
    for j in range(p,2*p+1) :
        S += sc.comb(j,p)/2**j
    return S

print(Conjecture(120))
# Cette quantité semble valoir 1 pour tout entier p.

# Q2a

def Simu(N) :
    P1 = N
    P2 = N
    while P1 != 0 and P2 != 0 :
        paquet = rd.randint(1,3)
        if paquet == 1 :
            P1 += -1
        else :
            P2 += -1
    return P1 + P2

print(Simu(5))

# Q2b

def Esp(N,seuil) :
    S = 0
    for i in range(seuil) :
        S += Simu(N)
    return S/seuil

print(Esp(5,2000))

# Q2c

def Graphique() :
    X = range(1,101)
    Y = []
    for N in X :
        Y += [ Esp(N,400) ]
    plt.plot(X,Y)
    plt.show()
    
Graphique()