CORDIC

import Pkg; Pkg.add("Printf")

using Printf

function cordic(alpha, iteration = 24)
    # Если угол больше 2π, может быть сбой знака при определении нового знака
    newsgn = isodd(Int(floor(alpha / 2π))) ? 1 : -1
    
    # Приведение угла к диапазону [-π/2, π/2]
    alpha < -π/2 && return newsgn * cordic(alpha + π, iteration)
    alpha > π/2 && return newsgn * cordic(alpha - π, iteration)

    # Таблица заранее вычисленных арктангенсов 2^(-j), где j = 0, 1, 2, ...
    angles = [
        0.78539816339745,   0.46364760900081,   0.24497866312686,   0.12435499454676,
        0.06241880999596,   0.03123983343027,   0.01562372862048,   0.00781234106010,
        0.00390623013197,   0.00195312251648,   0.00097656218956,   0.00048828121119,
        0.00024414062015,   0.00012207031189,   0.00006103515617,   0.00003051757812,
        0.00001525878906,   0.00000762939453,   0.00000381469727,   0.00000190734863,
        0.00000095367432,   0.00000047683716,   0.00000023841858,   0.00000011920929,
        0.00000005960464,   0.00000002980232,   0.00000001490116,   0.00000000745058]

    # Таблица коэффициентов масштабирования Kn, учитываемая в конце
    Kvalues = [
        0.70710678118655,   0.63245553203368,   0.61357199107790,   0.60883391251775,
        0.60764825625617,   0.60735177014130,   0.60727764409353,   0.60725911229889,
        0.60725447933256,   0.60725332108988,   0.60725303152913,   0.60725295913894,
        0.60725294104140,   0.60725293651701,   0.60725293538591,   0.60725293510314,
        0.60725293503245,   0.60725293501477,   0.60725293501035,   0.60725293500925,
        0.60725293500897,   0.60725293500890,   0.60725293500889,   0.60725293500888]

    # Выбираем значение коэффициента масштабирования в зависимости от числа итераций
    Kn = Kvalues[min(iteration, length(Kvalues))]  

    # Начальное значение вектора [cos(0), sin(0)] = [1, 0]
    v = [1, 0]  
    poweroftwo = 1  # степень двойки, 2^0 = 1, изменяется на каждой итерации
    angle = angles[1]  # текущий угол поворота
    
    # Основной цикл итераций CORDIC
    for j = 0:iteration-1
        if alpha < 0
            sigma = -1  # поворачиваем против часовой стрелки
        else
            sigma = 1   # поворачиваем по часовой стрелке
        end
        
        factor = sigma * poweroftwo  # коэффициент для матрицы поворота

        # Матрица поворота на угол atan(2^(-j)), записанная в CORDIC-виде
        R = [1 -factor
             factor  1]

        # Умножаем вектор v на матрицу R (поворачиваем вектор)
        v = R * v  

        # Обновляем значение оставшегося угла для следующей итерации
        alpha -= sigma * angle 

        # Переход к следующей степени двойки
        poweroftwo /= 2  

        # Обновляем угол для следующей итерации
        if j + 2 > length(angles)
            angle /= 2  # если таблица закончилась, уменьшаем угол вдвое
        else
            angle = angles[j + 2]  # иначе берем следующий угол из таблицы
        end
    end

    # Применяем масштабирующий коэффициент, чтобы получить истинные значения
    v .*= Kn
    return v  # возвращаем [cos(alpha), sin(alpha)]
end

function test_cordic()
    println("   x      sin(x)        Δsin         cos(x)        Δcos ")

    # Проверка для углов от -90° до +90° с шагом 15°
    for θ in -90:15:90 
        # Преобразуем градусы в радианы и запускаем функцию CORDIC
        cosθ, sinθ = cordic(deg2rad(θ))

        # Печатаем значения в формате: угол, sin, ошибка синуса, cos, ошибка косинуса
        @printf("%+05.1f°  %+.8f (%+.8f) %+.8f (%+.8f)\n",
           θ, sinθ, sinθ - sind(θ), cosθ, cosθ - cosd(θ))
    end

    println("\nx(рад)    sin(x)        Δsin         cos(x)        Δcos ")

    # Проверка для некоторых углов в радианах
    for θr in [-9, 0, 1.5, 6]
        cosθ, sinθ = cordic(θr)
        @printf("%+3.1f    %+.8f (%+.8f) %+.8f (%+.8f)\n",
           θr, sinθ, sinθ - sin(θr), cosθ, cosθ - cos(θr))
    end
end

test_cordic()