# Berechnung der Chebyshev-Punkte mittels Tridiagonalmatrix
# Aufgabe 6.3

import numpy as np

def chebyshev_tri_matrix(n):
    J = np.zeros((n,n))
    for i in range(n-1):
        J[i,i+1] = 1/2
        J[i+1,i] = 1/2
    J[0,1]=1/np.sqrt(2)
    J[1,0]=1/np.sqrt(2)
    return J

def mysign(x): #Konvention sign(0)=1
    y=np.sign(x)
    if abs(x)<1e-15:
        y=1
    return y

def qr(A):
    n_row, n_col = A.shape
    P = np.eye(n_row,n_row)
    for k in range(n_col):
        T1= A[k:,k]
        alpha = -mysign(T1[0])*np.linalg.norm(T1)
        v=np.array(T1.copy())
        v[0] = v[0]-alpha
        if np.linalg.norm(v)<1e-15:
            c=0
        else:     
            c = 2/np.linalg.norm(v)**2
        P[k:,:] = P[k:,:]-np.outer(c*v,v@P[k:,:])
        A[k:,k:] = A[k:,k:]-np.outer(c*v,v@A[k:,k:])
    return P.transpose(), A


n=8
for i in range(2,n):
    J=chebyshev_tri_matrix(i+1)
    # QR-Iteration
    y_alt=np.zeros(J.shape[0])
    maxiter=1000
    for j in range(maxiter):
        if j>maxiter:
            print("n=",n," QR konvergiert nicht")
            break
        Q,R=qr(J)
        J=R@Q
        y=np.sort(np.diag(R))
        if np.linalg.norm(y-y_alt)<1e-10:
           print("n=",n," konvergiert nach ",j," QR-Schritten")
           print("Chebyshev-Punkte:",y)
           break
        y_alt=y.copy()