import numpy as np
import multiprocessing as mp
import time

def sequential_determinant(matrix):
    n = matrix.shape[0]

    A = matrix.astype(np.float64).copy()

    sign = 1

    for i in range(n):
        max_row = i
        max_val = abs(A[i, i])
        
        for k in range(i + 1, n):
            if abs(A[k, i]) > max_val:
                max_row = k
                max_val = abs(A[k, i])

        if max_row != i:
            A[[i, max_row]] = A[[max_row, i]]
            sign *= -1

        if abs(A[i, i]) < 1e-15:
            return 0.0, -np.inf

        for k in range(i + 1, n):
            factor = A[k, i] / A[i, i]
            A[k, i:] -= factor * A[i, i:]

    log_abs_det = 0.0
    for i in range(n):
        if A[i, i] < 0:
            sign *= -1
        log_abs_det += np.log(abs(A[i, i]))
    
    return sign, log_abs_det

def to_float(sign, log_abs_det):
    if log_abs_det == -np.inf:
        return 0.0

    if log_abs_det > 700:
        return sign * np.inf
    
    return sign * np.exp(log_abs_det)


def worker_lu_decomposition(args):
    matrix, start, end = args
    n = matrix.shape[0]
    A = matrix.copy()
    perm_count = 0
    
    for i in range(start, min(end, n)):
        max_row = i
        for k in range(i + 1, n):
            if abs(A[k, i]) > abs(A[max_row, i]):
                max_row = k
        
        if max_row != i:
            A[[i, max_row]] = A[[max_row, i]]
            perm_count += 1
        
        if abs(A[i, i]) > 1e-15:
            for k in range(i + 1, n):
                factor = A[k, i] / A[i, i]
                A[k, i:] -= factor * A[i, i:]
    
    return A, perm_count


def parallel_determinant(matrix, num_processes = None):
    if num_processes is None:
        num_processes = mp.cpu_count()
    
    n = matrix.shape[0]
    
    if n < 100 or num_processes == 1:
        return sequential_determinant(matrix)
    
    A = matrix.astype(np.float64).copy()
    sign = 1
    
    with mp.Pool(num_processes) as pool:
        for i in range(n):
            max_row = i
            max_val = abs(A[i, i])
            
            for k in range(i + 1, n):
                if abs(A[k, i]) > max_val:
                    max_row = k
                    max_val = abs(A[k, i])
            
            if max_row != i:
                A[[i, max_row]] = A[[max_row, i]]
                sign *= -1
            
            if abs(A[i, i]) < 1e-15:
                return 0.0, -np.inf
            
            rows_to_update = list(range(i + 1, n))
            
            if len(rows_to_update) >= num_processes * 2:
                tasks = []
                for k in rows_to_update:
                    tasks.append((A[k].copy(), A[i].copy(), A[k, i], A[i, i], i))
                
                updated_rows = pool.map(_update_row_worker, tasks)

                for idx, k in enumerate(rows_to_update):
                    A[k] = updated_rows[idx]
            else:
                for k in rows_to_update:
                    factor = A[k, i] / A[i, i]
                    A[k, i:] -= factor * A[i, i:]

    log_abs_det = 0.0
    for i in range(n):
        if A[i, i] < 0:
            sign *= -1
        log_abs_det += np.log(abs(A[i, i]))
    
    return sign, log_abs_det


def _update_row_worker(args):
    row_k, row_i, a_ki, a_ii, i = args

    import os
    os.environ['OMP_NUM_THREADS'] = '1'
    os.environ['MKL_NUM_THREADS'] = '1'
    
    factor = a_ki / a_ii
    
    row_k[i:] -= factor * row_i[i:]
    
    return row_k


def generate_stable_matrix(size, scale = 1.0) -> np.ndarray:
    np.random.seed(13)

    random_matrix = np.random.randn(size, size)
    Q, R = np.linalg.qr(random_matrix)

    eigenvalues = np.random.uniform(0.5, 2.0, size) * scale
    D = np.diag(eigenvalues)

    A = Q @ D @ Q.T
    
    return A


def benchmark(sizes, thread_counts):
    results = []
    
    print("=" * 80)
    print("БЕНЧМАРК ВЫЧИСЛЕНИЯ ДЕТЕРМИНАНТА")
    print("=" * 80)
    print()
    
    for size in sizes:
        print(f"\n{'=' * 80}")
        print(f"Размер матрицы: {size}x{size}")
        print(f"{'=' * 80}")

        matrix = generate_stable_matrix(size)

        try:
            sign_numpy, logdet_numpy = np.linalg.slogdet(matrix)
            det_numpy = sign_numpy * np.exp(logdet_numpy) if logdet_numpy < 700 else sign_numpy * np.inf
            print(f"Детерминант (numpy):     {det_numpy:.6e}")
            print(f"  Sign: {sign_numpy}, Log|det|: {logdet_numpy:.6f}")
        except Exception as e:
            sign_numpy, logdet_numpy = None, None
            det_numpy = None
            print(f"Детерминант (numpy):     Ошибка: {e}")
        
        print()

        start_time = time.time()
        sign_seq, logdet_seq = sequential_determinant(matrix)
        seq_time = time.time() - start_time
        det_seq = to_float(sign_seq, logdet_seq)
        
        print(f"Последовательный алгоритм:")
        print(f"  Время:       {seq_time:.4f} сек")
        print(f"  Детерминант: {det_seq:.6e}")
        print(f"  Sign: {sign_seq}, Log|det|: {logdet_seq:.6f}")
        
        results.append({
            'size': size,
            'threads': 1,
            'time': seq_time,
            'determinant': det_seq,
            'speedup': 1.0
        })
        
        print()

        for num_threads in thread_counts:
            if num_threads == 1:
                continue
            
            start_time = time.time()
            sign_par, logdet_par = parallel_determinant(matrix, num_threads)
            par_time = time.time() - start_time
            det_par = to_float(sign_par, logdet_par)
            
            speedup = seq_time / par_time if par_time > 0 else 0
            
            print(f"Параллельный алгоритм ({num_threads} потоков):")
            print(f"  Время:       {par_time:.4f} сек")
            print(f"  Детерминант: {det_par:.6e}")
            print(f"  Ускорение:   {speedup:.2f}x")

            results.append({
                'size': size,
                'threads': num_threads,
                'time': par_time,
                'determinant': det_par,
                'speedup': speedup
            })
            
            print()
    
    return results

if __name__ == "__main__":
    sizes = [100, 300, 500]
    thread_counts = [1, 2, 4, mp.cpu_count()]

    thread_counts = sorted(list(set(thread_counts)))
    
    results = benchmark(sizes, thread_counts)

    print("\n" + "=" * 80)
    print("СВОДНАЯ ТАБЛИЦА РЕЗУЛЬТАТОВ")
    print("=" * 80)
    print()
    print(f"{'Размер':<10} {'Потоки':<10} {'Время (сек)':<15} {'Ускорение':<12}")
    print("-" * 80)
    
    for r in results:
        print(f"{r['size']:<10} {r['threads']:<10} {r['time']:<15.4f} {r['speedup']:<12.2f}")