#This file is made by Amine EL HEND to pass the test of CMN of the third year in
ensam meknes
# I hope this code will help you as it's helping me right now
# I tried to include all the functions that you might need, and good luck <3
# Dedikas l club GADZ'IT :P
#------------------------------------------------------------------------------
#The libraries
import numpy as np
import scipy
import scipy.linalg
import struct
#------------------------------------------------------------------------------
# LU decomposition by Crout
def lu_decomposition_crout(A):
n=len(A)
L=np.zeros((n,n))
U=np.eye(n)#pour l identité
L[:,0]=A[:,0]#:touts les lignes
U[0,1:]=A[0,1:]/L[0,0]
for i in range(1,n-1):
L[i,i]=A[i,i]-np.sum(L[i,:i]*U[:i,i])
for j in range(i+1,n):
L[j,i]=A[j,i]-np.sum(L[j,:i]*U[:i,i])
U[i,j]=(A[i,j]-np.sum(L[i,:i]*U[:i,j]))/L[i,i]
L[n-1,n-1]=A[n-1,n-1]-np.sum(L[n-1,:n]*U[:n,n-1])
return L,U
#------------------------------------------------------------------------------
# LU decomposition by Doolittle
def lu_decomposition_doolittle(A):
n = len(A)
L = np.zeros((n, n))
U = np.zeros((n, n))
for i in range(n):
L[i][i] = 1 # Les éléments diagonaux de L sont 1
for j in range(i, n):
U[i][j] = A[i][j] - sum(L[i][k] * U[k][j] for k in range(i))
for j in range(i + 1, n):
L[j][i] = (A[j][i] - sum(L[j][k] * U[k][i] for k in range(i))) / U[i]
[i]
return L, U
#------------------------------------------------------------------------------
# You can as well use the scipy lib for doolittle methode
# P, L, U = scipy.linalg.lu(A)
#------------------------------------------------------------------------------
# Jacobi itterative methode ///// it takes the matrix A and b , the initial
# value called x0 ,the tolerated error called eps (epsilon) , and the maximum
# number of ittertaions possible to not get an infinite loop
def Jacobi(A,b,x0,eps,itmax):
x = x0
M = np.diag(np.diag(A))
N = M - A
err = np.linalg.norm(A.dot(x) - b)
ite = 0
while(err>eps and ite<itmax):
x=np.linalg.solve(M,N.dot(x) + b)
err=np.linalg.norm(A.dot(x)-b)
ite = ite + 1
return x,ite # it returns the approximation of the solution x ,
# and the itteration where it did stop
#------------------------------------------------------------------------------
# Gauss-Seidel methode
# same as Jacobi
def Gauss_Seidel(A,b,x0,eps,itmax):
x = x0
M = np.diag(np.diag(A))
for i in range(1,len(A)):
M = M + np.diag(np.diag(A,-i),-i)
N = M - A
err = np.linalg.norm(A.dot(x)-b)
it = 0
while err > eps and it < itmax:
x = np.linalg.solve(M,N.dot(x) + b)
err = np.linalg.norm(A.dot(x)-b)
it = it + 1
return x , it
#------------------------------------------------------------------------------
# Relaxation methode
# same as the previous two
def Relaxation(A,b,x0,omega,eps,itmax):
x = x0
M = np.diag(np.diag(A)) / omega
for i in range (1,len(A)) :
M = M + np.diag(np.diag(A,-i),-i)
N = M - A
err = np.linalg.norm(A.dot(x) - b)
it = 0
while err > eps and it < itmax :
x = np.linalg.solve(M,N.dot(x) + b)
err = np.linalg.norm(A.dot(x) - b)
it = it+1
return x , it
#------------------------------------------------------------------------------
# This function gives you the Cj and the dj matrix of jacobi
def Jacobi_Cj_dj(A,b):
n = len(b)
Cj =np.zeros((n,n), dtype=float)
dj = np.zeros((n,), dtype=float)
for i in range(n):
x = float(A[i][i])
Cj[i] = A[i] / x
dj[i] = b[i] / x
Cj = Cj - np.diag(np.diag(Cj))
Cj = -Cj
return Cj , dj
#------------------------------------------------------------------------------
# This function returns the Cg and dg matrix of Gauss-Seidel methode
def Gauss_Seidel_Cg_dg(A,b):
D,U,L = split_matrix(A)
Cg = -np.linalg.inv(D+L) @ U
dg = np.linalg.inv(D+L) @ b
return Cg , dg
#------------------------------------------------------------------------------
# This function returns the Cr and dr matrix of Relaxation methode
def Relaxation_Cr_dr(A,b,omega):
D,U,L = split_matrix(A)
Cr = np.linalg.inv( D + omega*L) @ ((1-omega)*D - omega*U )
dr = omega*np.linalg.inv( D + omega*L) @ b
return Cr , dr
# Gauss pivot methode to get the U and the c
# Note that you must not have any zeros in the diag !!!!!!
# the input b must be 2D array also
# example: b = np.array([[6],[-2],[-8],[-4]],dtype="float64")
def Gauss_pivot(A,b):
n = len(b)
i = 0
M = np.concatenate((A,b), axis=1)
while i < n:
if M[i][i] == 0.0 :
print("error 0")
return
for j in range(i+1,n):
x = M[j][i] / M[i][i]
M[j] = M[j] - ( x * M[i] )
i = i + 1
c = M[:, n-1]
U = np.delete(M, n-1, axis=1)
return U , c
#------------------------------------------------------------------------------
# Gauss pivot to get U,c and P and also the solution x
# Thanks to my collegue Ayman
def back_sub(A,b):
m = np.size(A,0)
n = np.size(A,1)
assert(n==m)
d = np.diag(A)
p = np.prod(d)
assert(p!=0)
if np.array_equal(np.triu(A),A) == False:
print("Matrice n'est pas triangulaire sup.")
return
x = np.zeros(n)
x[n-1] = b[n-1]/A[n-1,n-1]
for i in range(n-2,-1,-1):
s = 0
for k in range(i+1, n):
s = s + A[i,k]*x[k]
x[i] = (b[i]-s)/A[i,i]
return x
def permut_lign(M, L1, L2):
tmp = np.array(M[L1,:])
M[L1,:] = np.array(M[L2,:])
M[L2,:]= np.array(tmp)
return M
#Méthode d'élimination de Gauss
def Gauss(A,b):
n = len(b)
M = np.array(np.c_[A,b]) # M:matrice augmentée = concatenation de A et b
P= np.diag(np.ones(n))
for j in range(0,n-1):
if M[j,j]==0 : # pivot=0 => on fait pivotement
iligne= j+1
while M[iligne,j]==0:
iligne += 1
M= permut_lign(M,j,iligne)
P= permut_lign(P,j,iligne)
pivot = M[j,j]
for i in range(j+1,n):
M[i,:] = M[i,:]- ( M[i,j]/pivot *M[j,:] )
x= back_sub( M[:,0:n] , M[:,n])
U = M[:,0:n]
c = M[:,n]
return x,U,c,P
#------------------------------------------------------------------------------
# Cholesky methode
# We will use the easy way, and that is with the pre-existed function
# L = scipy.linalg.cholesky(A, lower=True)
# LT = scipy.linalg.cholesky(A, lower=False)
# LT is the transposit of L
#------------------------------------------------------------------------------
# A function to split the matrix to D,L and U
def split_matrix(matrix):
n = matrix.shape[0]
diagonal_matrix = np.diag(np.diag(matrix))
upper_triangular_matrix = np.triu(matrix, k=1)
lower_triangular_matrix = np.tril(matrix, k=-1)
return diagonal_matrix, upper_triangular_matrix, lower_triangular_matrix
#------------------------------------------------------------------------------
# numpy methodes
# there is a lot of numpy methodes that are usefull in our situation :
# np.linalg.solve(a,b) this methode solve a liniaire systeme
# np.linalg.inv(A) this one gives the inv of a matrix
# np.linalg.matrix_power(A,n) this one puts a matrix to the power of n
#------------------------------------------------------------------------------
# the multiplication of two matrix
# np.dot(a,b) or use @ as a multiplication operator
#------------------------------------------------------------------------------
# Convert a decimal number to ieee754 32bits forme
def float_to_bin(num):
bits, = struct.unpack('!I', struct.pack('!f', num))
return "{:032b}".format(bits)
#------------------------------------------------------------------------------
# Convert an ieee754 32 bits number to a decimal forme
def ieee754(N):
a = int(N[0])
b = int(N[1:9],2)
c = int("1"+N[9:], 2)
return (-1)**a * c /( 1<<( len(N)-9 - (b-127) ))
#------------------------------------------------------------------------------
# Convert an ieee754 64 bits number to a decimal forme
def ieee754_double_to_decimal(ieee754_binary):
if len(ieee754_binary) != 64:
raise ValueError("Input binary string must be 64 bits.")
ieee754_bytes = bytes(int(ieee754_binary[i:i+8], 2) for i in range(0, 64, 8))
decimal_number = struct.unpack('!d', ieee754_bytes)[0]
return decimal_number
#------------------------------------------------------------------------------
#convert a decimal number to ieee754 forme on 64 bits
def decimal_to_ieee754_double(decimal_number):
ieee754_bytes = struct.pack('!d', decimal_number)
ieee754_binary = ''.join(f'{byte:08b}' for byte in ieee754_bytes)
return ieee754_binary
#------------------------------------------------------------------------------
# Calculate F(b,t,e_min,e_max) and its cardinal
def floating_numbers(beta,t,L,U):
cardinal=2*(beta**t-beta**(t-1))*(U-L+1)
F=np.zeros(cardinal//2)
i=0
for e in range(L,U+1):
for m in range(beta**(t-1),beta**t):
F[i]=m*beta**(e-t)
i=i+1
F=np.concatenate((-F[-1::-1],F),axis = 0)
n = len(F)
return F , n # F is the numbers and n is the cardinal
#------------------------------------------------------------------------------
# To determine the norme of the error Ax=b
# np.linalg.norm(A.dot(x) - b))
#------------------------------------------------------------------------------
# Convert a none signed integer to a binary forme and also give the number of bits
def none_signed_integer_to_binary(n):
n = abs(n)
return bin(n)[2:] , len(bin(n)[2:])
def binary_to_decimal_str(binary_str):
decimal_value = 0
binary_str = binary_str[::-1]
for i, bit in enumerate(binary_str):
if bit == '1':
decimal_value += 2**i
return decimal_value
#------------------------------------------------------------------------------
# Convert a signed integer to a binary forme and also give the number of bits
def signed_integer_to_binary(n):
if n < 0 :
return "1" + bin(n)[3:] , len(bin(n)[3:]) + 1
elif n == 0 :
return 0 , 0
else:
return "0" + bin(n)[2:] , len(bin(n)[2:]) + 1
#------------------------------------------------------------------------------
# Convert a number to two complement forme and also the number of bits
def decimal_to_twos_complement(decimal_number):
if decimal_number >= 0:
binary_representation = bin(decimal_number)[2:]
num_bits = len(binary_representation)
else:
positive_representation = bin(abs(decimal_number) - 1)[2:]
num_bits = len(positive_representation)
inverted_representation = ''.join('0' if bit == '1' else '1' for bit in
positive_representation)
binary_representation = inverted_representation
return binary_representation, num_bits
#------------------------------------------------------------------------------
# Convert a number to one complement forme and also the number of bits
def decimal_to_ones_complement(decimal_number):
if decimal_number >= 0:
binary_representation = bin(decimal_number)[2:]
num_bits = len(binary_representation)
else:
positive_representation = bin(abs(decimal_number))[2:]
num_bits = len(positive_representation)
inverted_representation = ''.join('0' if bit == '1' else '1' for bit in
positive_representation)
binary_representation = inverted_representation
return binary_representation, num_bits
#------------------------------------------------------------------------------
# This function convert a string of a binary number to its decimal form
# Note that this function work only for after comma numbers
# We are going to use this function next
def base_of_floating_comma(n , b):
l = len(n)
x= 0
b = float(b)
power = -1
for i in range(l):
if n[i] == "0" :
power = power - 1
else:
x = x + b**power
power = power - 1
return x
#------------------------------------------------------------------------------
# Convert any ieee754 number to decimal in any number of bits
def ieee754_high_level(N,e,m,b):
signe = int(N[0])
exponent = N[1:e+1]
mantisse = N[e+1:]
decalage = 2**(e-1) -1
exponent0 = binary_to_decimal_str(exponent)
mantisse0 = base_of_floating_comma(mantisse , 2)
value = (-1)**signe * (1+mantisse0) * 2**(exponent0-decalage)
return value
#------------------------------------------------------------------------------
