const std = @import("std");
const vm = @import("root.zig");

/// The number of radians per one turn (τ).
pub const rad_per_turn = 1.0 * std.math.tau;

/// The number of degrees per one turn (360).
pub const deg_per_turn = 360.0;

/// The number of turns per one radian (1 / τ).
pub const turns_per_rad = 1.0 / std.math.tau;

/// The number of turns per one degree (1 / 360).
pub const turns_per_deg = 1.0 / 360.0;

/// Converts an angle in turns to radians. `@TypeOf(ang)` must be a float or
/// comptime number or a vector of floats.
pub fn turnsToRadians(ang: anytype) if (@TypeOf(ang) == comptime_int) comptime_float else @TypeOf(ang) {
    const T = @TypeOf(ang);
    switch (@typeInfo(T)) {
        .float, .comptime_float, .comptime_int => return ang * rad_per_turn,
        .vector => |V| if (@typeInfo(V.child) == .float) return ang * @as(T, @splat(rad_per_turn)),
        else => {},
    }
    @compileError("Input must be float or a comptime number, or a vector of floats.");
}

test turnsToRadians {
    try std.testing.expectEqual(0, turnsToRadians(@as(f32, 0)));
    try std.testing.expectApproxEqAbs(0.25 * std.math.tau, turnsToRadians(@as(f32, 0.25)), 0x1p-5);
    try std.testing.expectApproxEqAbs(0.5 * std.math.tau, turnsToRadians(@as(f32, 0.5)), 0x1p-5);
    try std.testing.expectApproxEqAbs(0.75 * std.math.tau, turnsToRadians(@as(f32, 0.75)), 0x1p-5);
    try std.testing.expectApproxEqAbs(std.math.tau, turnsToRadians(@as(f32, 1)), 0x1p-5);
}

/// Converts an angle in turns to degrees. `@TypeOf(ang)` must be a float or
/// comptime number or a vector of floats.
pub fn turnsToDegrees(ang: anytype) if (@TypeOf(ang) == comptime_int) comptime_float else @TypeOf(ang) {
    const T = @TypeOf(ang);
    switch (@typeInfo(T)) {
        .float, .comptime_float, .comptime_int => return ang * deg_per_turn,
        .vector => |V| if (@typeInfo(V.child) == .float) return ang * @as(T, @splat(deg_per_turn)),
        else => {},
    }
    @compileError("Input must be float or a comptime number, or a vector of floats.");
}

test turnsToDegrees {
    try std.testing.expectEqual(0, turnsToDegrees(@as(f32, 0)));
    try std.testing.expectEqual(90, turnsToDegrees(@as(f32, 0.25)));
    try std.testing.expectEqual(180, turnsToDegrees(@as(f32, 0.5)));
    try std.testing.expectEqual(270, turnsToDegrees(@as(f32, 0.75)));
    try std.testing.expectEqual(360, turnsToDegrees(@as(f32, 1)));
}

/// Converts an angle in radians to turns. `@TypeOf(ang)` must be a float or
/// comptime number or a vector of floats.
pub fn radiansToTurns(ang: anytype) if (@TypeOf(ang) == comptime_int) comptime_float else @TypeOf(ang) {
    const T = @TypeOf(ang);
    switch (@typeInfo(T)) {
        .float, .comptime_float, .comptime_int => return ang * turns_per_rad,
        .vector => |V| if (@typeInfo(V.child) == .float) return ang * @as(T, @splat(turns_per_rad)),
        else => {},
    }
    @compileError("Input must be float or a comptime number, or a vector of floats.");
}

test radiansToTurns {
    try std.testing.expectEqual(0, radiansToTurns(@as(f32, 0)));
    try std.testing.expectApproxEqAbs(0.25, radiansToTurns(@as(f32, 0.25 * std.math.tau)), 0x1p-5);
    try std.testing.expectApproxEqAbs(0.5, radiansToTurns(@as(f32, 0.5 * std.math.tau)), 0x1p-5);
    try std.testing.expectApproxEqAbs(0.75, radiansToTurns(@as(f32, 0.75 * std.math.tau)), 0x1p-5);
    try std.testing.expectApproxEqAbs(1, radiansToTurns(@as(f32, std.math.tau)), 0x1p-5);
}

/// Converts an angle in degrees to turns. `@TypeOf(ang)` must be a float or
/// comptime number or a vector of floats.
pub fn degreesToTurns(ang: anytype) if (@TypeOf(ang) == comptime_int) comptime_float else @TypeOf(ang) {
    const T = @TypeOf(ang);
    switch (@typeInfo(T)) {
        .float, .comptime_float, .comptime_int => return ang / deg_per_turn,
        .vector => |V| if (@typeInfo(V.child) == .float) return ang / @as(T, @splat(deg_per_turn)),
        else => {},
    }
    @compileError("Input must be float or a comptime number, or a vector of floats.");
}

test degreesToTurns {
    try std.testing.expectEqual(0, degreesToTurns(@as(f32, 0)));
    try std.testing.expectEqual(0.25, degreesToTurns(@as(f32, 90)));
    try std.testing.expectEqual(0.5, degreesToTurns(@as(f32, 180)));
    try std.testing.expectEqual(0.75, degreesToTurns(@as(f32, 270)));
    try std.testing.expectEqual(1, degreesToTurns(@as(f32, 360)));
}

pub fn cos(angle_turns: f32) f32 {
    return cossin(angle_turns).re;
}

test cos {
    try std.testing.expectEqual(1, cos(-1));
    try std.testing.expectEqual(0, cos(-0.75));
    try std.testing.expectEqual(-1, cos(-0.5));
    try std.testing.expectEqual(0, cos(-0.25));
    try std.testing.expectEqual(1, cos(0));
    try std.testing.expectEqual(0, cos(0.25));
    try std.testing.expectEqual(-1, cos(0.5));
    try std.testing.expectEqual(0, cos(0.75));
}

pub fn cos_x8(angle_turns: vm.f32x8) vm.f32x8 {
    return cossin_x8(angle_turns).re;
}

test cos_x8 {
    try std.testing.expectEqual(
        .{ 1, 0, -1, 0, 1, 0, -1, 0 },
        cos_x8(.{ -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75 }),
    );
}

pub fn sin(angle_turns: f32) f32 {
    return cossin(angle_turns).im;
}

test sin {
    try std.testing.expectEqual(0, sin(-1));
    try std.testing.expectEqual(1, sin(-0.75));
    try std.testing.expectEqual(0, sin(-0.5));
    try std.testing.expectEqual(-1, sin(-0.25));
    try std.testing.expectEqual(0, sin(0));
    try std.testing.expectEqual(1, sin(0.25));
    try std.testing.expectEqual(0, sin(0.5));
    try std.testing.expectEqual(-1, sin(0.75));
}

pub fn sin_x8(angle_turns: vm.f32x8) vm.f32x8 {
    return cossin_x8(angle_turns).im;
}

test sin_x8 {
    try std.testing.expectEqual(
        .{ 0, 1, 0, -1, 0, 1, 0, -1 },
        sin_x8(.{ -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75 }),
    );
}

pub fn cossin(angle_turns: f32) vm.Complex {
    @setFloatMode(.optimized);
    // Taylor series expansion for f(x)=cos(xπ/2)
    const term_cos_0: f32 = 1.0;
    const term_cos_2: f32 = -1.23370055; // -π²/8
    const term_cos_4: f32 = 0.253669508; // π⁴/384
    const term_cos_6: f32 = -0.020863481; // -π⁶/46080
    // Taylor series expansion for f(x)=sin(xπ/2)
    const term_sin_1: f32 = 1.570796327; // π/2
    const term_sin_3: f32 = -0.645964098; // -π³/48
    const term_sin_5: f32 = 0.079692626; // π⁵/3840

    const angle_01 = angle_turns - @floor(angle_turns);
    const angle_04 = 4.0 * angle_01;

    const quadrant: u32 = @intFromFloat(angle_04);
    const quadrant_odd = (quadrant & 1) * 0xFFFFFFFF;
    const sign_mask_cos = ((quadrant + 1) & 0b10) << 30;
    const sign_mask_sin = (quadrant & 0b10) << 30;

    const x = angle_04 - @floor(angle_04);
    const x2 = x * x;

    const c: u32 = @bitCast(((term_cos_6 * x2 + term_cos_4) * x2 + term_cos_2) * x2 + term_cos_0);
    const s: u32 = @bitCast(((term_sin_5 * x2 + term_sin_3) * x2 + term_sin_1) * x);

    const result_cos: f32 = @bitCast(((s & quadrant_odd) | (c & ~quadrant_odd)) ^ sign_mask_cos);
    const result_sin: f32 = @bitCast(((c & quadrant_odd) | (s & ~quadrant_odd)) ^ sign_mask_sin);

    return .init(result_cos, result_sin);
}

test cossin {
    try std.testing.expectEqual(vm.Complex.initVector2(.unit_x), cossin(-1));
    try std.testing.expectEqual(vm.Complex.initVector2(.unit_y), cossin(-0.75));
    try std.testing.expectEqual(vm.Complex.initVector2(.unit_nx), cossin(-0.5));
    try std.testing.expectEqual(vm.Complex.initVector2(.unit_ny), cossin(-0.25));
    try std.testing.expectEqual(vm.Complex.initVector2(.unit_x), cossin(0));
    try std.testing.expectEqual(vm.Complex.initVector2(.unit_y), cossin(0.25));
    try std.testing.expectEqual(vm.Complex.initVector2(.unit_nx), cossin(0.5));
    try std.testing.expectEqual(vm.Complex.initVector2(.unit_ny), cossin(0.75));
}

pub fn cossin_x8(angle_turns: vm.f32x8) vm.Complex_x8 {
    @setFloatMode(.optimized);
    // Taylor series expansion for f(x)=cos(xπ/2)
    const term_cos_0 = vm.ps(1.0);
    const term_cos_2 = vm.ps(-1.23370055); // -π²/8
    const term_cos_4 = vm.ps(0.253669508); // π⁴/384
    const term_cos_6 = vm.ps(-0.020863481); // -π⁶/46080
    // Taylor series expansion for f(x)=sin(xπ/2)
    const term_sin_1 = vm.ps(1.570796327); // π/2
    const term_sin_3 = vm.ps(-0.645964098); // -π³/48
    const term_sin_5 = vm.ps(0.079692626); // π⁵/3840

    const angle_01 = angle_turns - @floor(angle_turns);
    const angle_04 = vm.ps(4.0) * angle_01;

    const quadrant: vm.u32x8 = @intFromFloat(angle_04);
    const quadrant_odd = (quadrant & vm.epu32(1)) != vm.epu32(0);
    const sign_mask_cos = ((quadrant + vm.epu32(1)) & vm.epu32(0b10)) << @splat(30);
    const sign_mask_sin = (quadrant & vm.epu32(0b10)) << @splat(30);

    const x = angle_04 - @floor(angle_04);
    const x2 = x * x;

    const c: vm.u32x8 = @bitCast(((term_cos_6 * x2 + term_cos_4) * x2 + term_cos_2) * x2 + term_cos_0);
    const s: vm.u32x8 = @bitCast(((term_sin_5 * x2 + term_sin_3) * x2 + term_sin_1) * x);

    const result_cos: vm.f32x8 = @bitCast(@select(u32, quadrant_odd, s, c) ^ sign_mask_cos);
    const result_sin: vm.f32x8 = @bitCast(@select(u32, quadrant_odd, c, s) ^ sign_mask_sin);

    return .init(result_cos, result_sin);
}

test cossin_x8 {
    try std.testing.expectEqual(
        vm.Complex_x8.initArrayOfComplex(.{
            .initVector2(.unit_x),
            .initVector2(.unit_y),
            .initVector2(.unit_nx),
            .initVector2(.unit_ny),
            .initVector2(.unit_x),
            .initVector2(.unit_y),
            .initVector2(.unit_nx),
            .initVector2(.unit_ny),
        }),
        cossin_x8(.{ -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75 }),
    );
}