#!/usr/bin/env python3

import scipy.linalg as la
import numpy as np

def Main():
    print('N\tS_m\t\tS_M\t\tS')

    N = 1000

    # Assymptotic solution (N -> ∞)
    S = np.log(N+2) + np.log(2) - 1

    E, S_k = Solution(N)

    i_m = np.argmin(S_k)
    i_M = np.argmax(S_k)

    print(f'{N}\t{S_k[i_m]:.8f}\t{S_k[i_M]:.8f}\t{S:.8f}')

def Solution(N):
    n = np.arange(1, N+2)

    # Quasi-momentum
    K = n*np.pi/(N+2)

    # Energy
    E = 2*np.cos(K)/N

    S_k = []

    for k in K:
        # Wave-function
        Psi = np.sin(k*n)/np.sqrt((N+2)/2)
        s = [p**2*np.log(p**2) for p in Psi if abs(p) > 1e-12]
        S_k.append( -np.sum(s) )

    return E, S_k

Main()