#!/usr/bin/env python
# coding: utf-8

# # Performances des nombres sous-normaux
# 
# Pour tester la représentation en base 2 : https://www.h-schmidt.net/FloatConverter/IEEE754.html

# In[ ]:


import numpy as np
import time
import math
import matplotlib.pyplot as plt


# In[ ]:


# Python utilise des nombres flottants à double précision (64 bits) :
# - signe sur 1 bit
# - mantisse de 53 bits (52 stockés)
# - exposant sur 11 bits (de -1024 à +1023)
e_max = -(1024-100)
e_min = -(1024+60) # Pour légèrement dépasser la capacité des nombres sous-normaux
N = 100000

duration = np.empty(e_max - e_min + 1)

for e in range(e_min, e_max+1):
    values = (2**e) * np.ones(N)
    tic = time.process_time()
    r = values.sum()
    duration[e - e_min] = time.process_time() - tic;

plt.figure(figsize=(12, 6))
plt.semilogy(range(e_min, e_max+1), duration)
plt.xlabel('exposant')
plt.ylabel('temps');
    


# # Approximation des réels

# In[ ]:


a = 10000.0
b = 9999.5
c = 0.1


# In[ ]:


print("{:.20f}".format(c))


# In[ ]:


a1 = a + c
print("{:.20f}".format(a1))


# In[ ]:


print("{:.20f}".format(a1 - b))


# In[ ]:


print("{:.20f}".format(0.6))


# # Unit in the Last Place

# In[ ]:


def ulp(v):
    print("{}\t{}".format(v, v - np.nextafter(v, v-1)))
    
ulp(1000)
ulp(100)
ulp(10)
ulp(1)
ulp(0.1)
ulp(0.01)
ulp(0)


# In[ ]:


(1000 + 1.5e-12) - 1000


# In[ ]:


(1000 + 1e-14) - 1000


# In[ ]:


(1 + 1e-16) - 1


# In[ ]:


(1 + 2e-16) - 1


# In[ ]:


(1 + 3e-16) - 1


# In[ ]:


(1 + 4e-16) - 1


# # Non associativité des opérations

# In[ ]:


0.1 + (0.2 - 0.3)


# In[ ]:


(0.1 + 0.2) - 0.3


# In[ ]:


1 + 1e-17 - 1


# In[ ]:


1 - 1 + 1e-17


# # Annulation catastrophique
# 
# ## Différences finies

# In[ ]:


def f(x):
    return x**3

def df(x):
    return 3 * x**2

def deriv_fd1(f, x, h):
    return (f(x+h) - f(x)) / h

def deriv_fd2(f, x, h):
    return (f(x+h) - f(x-h)) / (2*h)

x = 1
h = np.logspace(-16, -1, 101)
err1 = abs(deriv_fd1(f, x, h) - df(x))
err2 = abs(deriv_fd2(f, x, h) - df(x))

plt.figure(figsize=(12, 6))
plt.loglog(h, err1)
plt.loglog(h, err2)
plt.legend(['Ordre 1', 'Ordre 2'])
plt.xlabel('h')
plt.ylabel('erreurr');


# In[ ]:


f(x)


# In[ ]:


f(x + 1e-9) # (1 + 1e-9)^3 = 1 + 3e-9 + 3e-18 + 1e-27


# In[ ]:


f(x + 1e-9) - f(x - 1e-9) # = 6e-9 + 2e-27


# In[ ]:


f(x + 1e-9) - f(x - 1e-9) - 6e-9 # = 2e-27


# In[ ]:


(f(x + 1e-9) - f(x - 1e-9)) / 2e-9 # = 3 + 1e-18


# ## Calcul des racines d'un polynôme du second degré

# In[ ]:


a = 1.
b = 1e4
c = 1.


# In[ ]:


# Calcul des racines
delta = b*b - 4*a*c
x = (-b - math.sqrt(delta)) / (2*a)
y = (-b + math.sqrt(delta)) / (2*a)
print("x = {}".format(x))
print("y = {}".format(y))


# In[ ]:


# Vérification
print("ax² + bx + c = {}".format(a*x*x + b*x + c))
print("ay² + by + c = {}".format(a*y*y + b*y + c))


# In[ ]:


b*b


# In[ ]:


a*c


# In[ ]:


math.sqrt(b*b - 4*a*c)


# $y=(-b+\sqrt{\delta})/(2*a)$ est mal calculé (soustraction de $2$ grands nombres proches, mais $x$ est bien calculé.
# 
# On sait que $x \times y = c / a$, d'où :

# In[ ]:


y = c / (x*a)
print("y = {}".format(y))


# In[ ]:


# Vérification
print("ax² + bx + c = {}".format(a*x*x + b*x + c))
print("ay² + by + c = {}".format(a*y*y + b*y + c))


# **Exercice :**
# envisager tous les cas possibles pour le
# choix de a, b et c :
# l'écriture d'un programme numériquement robuste pour le calcul des racines d'un
# trinôme du second degré  est loin d'être simple.

# # Erreurs d'arrondi
# 
# Calcul des termes d'une suite :
# 
# 
# ## Calcul de $u_{n+1} = 4 u_n - 1$ avec $u_0 = 1/3$

# In[ ]:


x = 1/3

for i in range(50):
    x = 4*x - 1
    print(x)


# ## Calcul de $u_{n+1} = 3 u_n - 1$ avec $u_0 = 1/2$

# In[ ]:


x = 1/2

for i in range(50):
    x = 3*x - 1
    print(x)


# # Problème mal conditionné
# 
# ## Résoudre un système linéaire $Ax = b$ avec $A$ quelconque

# In[ ]:


N = 40
A = np.random.rand(N, N)
b = np.arange(1, N+1)


# In[ ]:


x = np.linalg.solve(A, b)
x


# In[ ]:


np.max(np.abs(A @ x - b))


# In[ ]:


np.linalg.cond(A)


# ## Résoudre un système linéaire $Ax = b$ avec la matrice de Hilbert
# $$A_{i,j}=\frac{1}{i+j-1}.$$

# In[ ]:


N = 40
A = np.fromfunction(lambda i, j: 1 / (i+j+1), (N, N)) # i+j+1 car i et j commencent à 0 en Python
b = np.arange(1, N+1)


# In[ ]:


x = np.linalg.solve(A, b)
x


# In[ ]:


np.max(np.abs(A @ x - b))


# In[ ]:


np.linalg.cond(A)


# ## Résoudre un système linéaire $Ax = b$ avec la matrice du Laplacien 1D
# Résoudre
# $$\Delta u = 0$$
# avec
# $$\Delta u(x) \approx \frac{u(x-h) - 2u(x) + u(x+h)}{h^2}.$$

# In[ ]:


N = 40
A = np.diag(np.ones(N-1), 1) + np.diag(np.ones(N-1), -1) - 2 * np.eye(N, N)
b = np.arange(1, N+1)


# In[ ]:


x = np.linalg.solve(A, b)
x


# In[ ]:


np.max(np.abs(A @ x - b))


# In[ ]:


np.linalg.cond(A)


# In[ ]: