#ifndef _BGC_COMPLEX_H_INCLUDED_
#define _BGC_COMPLEX_H_INCLUDED_

#include "utilities.h"
#include "angle.h"

#include <math.h>

typedef struct
{
    float real, imaginary;
} BGC_FP32_Complex;

typedef struct
{
    double real, imaginary;
} BGC_FP64_Complex;

// =================== Reset ==================== //

inline void bgc_fp32_complex_reset(BGC_FP32_Complex* complex)
{
    complex->real = 0.0f;
    complex->imaginary = 0.0f;
}

inline void bgc_fp64_complex_reset(BGC_FP64_Complex* complex)
{
    complex->real = 0.0;
    complex->imaginary = 0.0;
}

// ==================== Set ===================== //

inline void bgc_fp32_complex_make(BGC_FP32_Complex* complex, const float real, const float imaginary)
{
    complex->real = real;
    complex->imaginary = imaginary;
}

inline void bgc_fp64_complex_make(BGC_FP64_Complex* complex, const double real, const double imaginary)
{
    complex->real = real;
    complex->imaginary = imaginary;
}

// ================== Modulus =================== //

inline float bgc_fp32_complex_get_square_modulus(const BGC_FP32_Complex* number)
{
    return number->real * number->real + number->imaginary * number->imaginary;
}

inline double bgc_fp64_complex_get_square_modulus(const BGC_FP64_Complex* number)
{
    return number->real * number->real + number->imaginary * number->imaginary;
}

inline float bgc_fp32_complex_get_modulus(const BGC_FP32_Complex* number)
{
    return sqrtf(bgc_fp32_complex_get_square_modulus(number));
}

inline double bgc_fp64_complex_get_modulus(const BGC_FP64_Complex* number)
{
    return sqrt(bgc_fp64_complex_get_square_modulus(number));
}

// ================= Comparison ================= //

inline int bgc_fp32_complex_is_zero(const BGC_FP32_Complex* number)
{
    return bgc_fp32_complex_get_square_modulus(number) <= BGC_FP32_SQUARE_EPSILON;
}

inline int bgc_fp64_complex_is_zero(const BGC_FP64_Complex* number)
{
    return bgc_fp64_complex_get_square_modulus(number) <= BGC_FP64_SQUARE_EPSILON;
}

inline int bgc_fp32_complex_is_unit(const BGC_FP32_Complex* number)
{
    return bgc_fp32_is_square_unit(bgc_fp32_complex_get_square_modulus(number));
}

inline int bgc_fp64_complex_is_unit(const BGC_FP64_Complex* number)
{
    return bgc_fp64_is_square_unit(bgc_fp64_complex_get_square_modulus(number));
}

// ==================== Copy ==================== //

inline void bgc_fp32_complex_copy(BGC_FP32_Complex* destination, const BGC_FP32_Complex* source)
{
    destination->real = source->real;
    destination->imaginary = source->imaginary;
}

inline void bgc_fp64_complex_copy(BGC_FP64_Complex* destination, const BGC_FP64_Complex* source)
{
    destination->real = source->real;
    destination->imaginary = source->imaginary;
}

// ==================== Swap ==================== //

inline void bgc_fp32_complex_swap(BGC_FP32_Complex* number1, BGC_FP32_Complex* number2)
{
    const float real = number2->real;
    const float imaginary = number2->imaginary;

    number2->real = number1->real;
    number2->imaginary = number1->imaginary;

    number1->real = real;
    number1->imaginary = imaginary;
}

inline void bgc_fp64_complex_swap(BGC_FP64_Complex* number1, BGC_FP64_Complex* number2)
{
    const double real = number2->real;
    const double imaginary = number2->imaginary;

    number2->real = number1->real;
    number2->imaginary = number1->imaginary;

    number1->real = real;
    number1->imaginary = imaginary;
}

// ================== Convert =================== //

inline void bgc_fp64_complex_convert_to_fp32(BGC_FP32_Complex* destination, const BGC_FP64_Complex* source)
{
    destination->real = (float)source->real;
    destination->imaginary = (float)source->imaginary;
}

inline void bgc_fp32_complex_convert_to_fp64(BGC_FP64_Complex* destination, const BGC_FP32_Complex* source)
{
    destination->real = source->real;
    destination->imaginary = source->imaginary;
}

// ================== Negative ================== //

inline void bgc_fp32_complex_revert(BGC_FP32_Complex* number)
{
    number->real = -number->real;
    number->imaginary = -number->imaginary;
}

inline void bgc_fp64_complex_revert(BGC_FP64_Complex* number)
{
    number->real = -number->real;
    number->imaginary = -number->imaginary;
}

inline void bgc_fp32_complex_get_reverse(BGC_FP32_Complex* reverse, const BGC_FP32_Complex* number)
{
    reverse->real = -number->real;
    reverse->imaginary = -number->imaginary;
}

inline void bgc_fp64_complex_get_reverse(BGC_FP64_Complex* reverse, const BGC_FP64_Complex* number)
{
    reverse->real = -number->real;
    reverse->imaginary = -number->imaginary;
}

// ================= Normalize ================== //

inline int bgc_fp32_complex_normalize(BGC_FP32_Complex* number)
{
    const float square_modulus = bgc_fp32_complex_get_square_modulus(number);

    if (bgc_fp32_is_square_unit(square_modulus)) {
        return 1;
    }

    if (square_modulus <= BGC_FP32_SQUARE_EPSILON || isnan(square_modulus)) {
        return 0;
    }

    const float multiplicand = sqrtf(1.0f / square_modulus);

    number->real *= multiplicand;
    number->imaginary *= multiplicand;

    return 1;
}

inline int bgc_fp64_complex_normalize(BGC_FP64_Complex* number)
{
    const double square_modulus = bgc_fp64_complex_get_square_modulus(number);

    if (bgc_fp64_is_square_unit(square_modulus)) {
        return 1;
    }

    if (square_modulus <= BGC_FP64_SQUARE_EPSILON || isnan(square_modulus)) {
        return 0;
    }

    const double multiplicand = sqrt(1.0 / square_modulus);

    number->real *= multiplicand;
    number->imaginary *= multiplicand;

    return 1;
}

inline int bgc_fp32_complex_get_normalized(BGC_FP32_Complex* normalized, const BGC_FP32_Complex* number)
{
    const float square_modulus = bgc_fp32_complex_get_square_modulus(number);

    if (bgc_fp32_is_square_unit(square_modulus)) {
        normalized->real = number->real;
        normalized->imaginary = number->imaginary;
        return 1;
    }

    if (square_modulus <= BGC_FP32_SQUARE_EPSILON || isnan(square_modulus)) {
        normalized->real = 0.0f;
        normalized->imaginary = 0.0f;
        return 0;
    }

    const float multiplicand = sqrtf(1.0f / square_modulus);

    normalized->real = number->real * multiplicand;
    normalized->imaginary = number->imaginary * multiplicand;

    return 1;
}

inline int bgc_fp64_complex_get_normalized(BGC_FP64_Complex* normalized, const BGC_FP64_Complex* number)
{
    const double square_modulus = bgc_fp64_complex_get_square_modulus(number);

    if (bgc_fp64_is_square_unit(square_modulus)) {
        normalized->real = number->real;
        normalized->imaginary = number->imaginary;
        return 1;
    }

    if (square_modulus <= BGC_FP64_SQUARE_EPSILON || isnan(square_modulus)) {
        normalized->real = 0.0;
        normalized->imaginary = 0.0;
        return 0;
    }

    const double multiplicand = sqrt(1.0 / square_modulus);

    normalized->real = number->real * multiplicand;
    normalized->imaginary = number->imaginary * multiplicand;

    return 1;
}

// ================= Conjugate ================== //

inline void bgc_fp32_complex_conjugate(BGC_FP32_Complex* number)
{
    number->imaginary = -number->imaginary;
}

inline void bgc_fp64_complex_conjugate(BGC_FP64_Complex* number)
{
    number->imaginary = -number->imaginary;
}

inline void bgc_fp32_complex_get_conjugate(BGC_FP32_Complex* conjugate, const BGC_FP32_Complex* number)
{
    conjugate->real = number->real;
    conjugate->imaginary = -number->imaginary;
}

inline void bgc_fp64_complex_get_conjugate(BGC_FP64_Complex* conjugate, const BGC_FP64_Complex* number)
{
    conjugate->real = number->real;
    conjugate->imaginary = -number->imaginary;
}

// =================== Invert =================== //

inline int bgc_fp32_complex_get_inverse(BGC_FP32_Complex* inverse, const BGC_FP32_Complex* number)
{
    const float square_modulus = bgc_fp32_complex_get_square_modulus(number);

    if (square_modulus <= BGC_FP32_SQUARE_EPSILON || isnan(square_modulus)) {
        return 0;
    }

    const float multiplicand = 1.0f / square_modulus;

    inverse->real = number->real * multiplicand;
    inverse->imaginary = -number->imaginary * multiplicand;

    return 1;
}

inline int bgc_fp64_complex_get_inverse(BGC_FP64_Complex* inverse, const BGC_FP64_Complex* number)
{
    const double square_modulus = bgc_fp64_complex_get_square_modulus(number);

    if (square_modulus <= BGC_FP64_SQUARE_EPSILON || isnan(square_modulus)) {
        return 0;
    }

    const double multiplicand = 1.0 / square_modulus;

    inverse->real = number->real * multiplicand;
    inverse->imaginary = -number->imaginary * multiplicand;

    return 1;
}

inline int bgc_fp32_complex_invert(BGC_FP32_Complex* number)
{
    return bgc_fp32_complex_get_inverse(number, number);
}

inline int bgc_fp64_complex_invert(BGC_FP64_Complex* number)
{
    return bgc_fp64_complex_get_inverse(number, number);
}

// =============== Get Exponation =============== //

void bgc_fp32_complex_get_exponation(BGC_FP32_Complex* power, const BGC_FP32_Complex* base, const float real_exponent, const float imaginary_exponent);

void bgc_fp64_complex_get_exponation(BGC_FP64_Complex* power, const BGC_FP64_Complex* base, const double real_exponent, const double imaginary_exponent);

// ==================== Add ===================== //

inline void bgc_fp32_complex_add(BGC_FP32_Complex* sum, const BGC_FP32_Complex* number1, const BGC_FP32_Complex* number2)
{
    sum->real = number1->real + number2->real;
    sum->imaginary = number1->imaginary + number2->imaginary;
}

inline void bgc_fp64_complex_add(BGC_FP64_Complex* sum, const BGC_FP64_Complex* number1, const BGC_FP64_Complex* number2)
{
    sum->real = number1->real + number2->real;
    sum->imaginary = number1->imaginary + number2->imaginary;
}

// ================= Add scaled ================= //

inline void bgc_fp32_complex_add_scaled(BGC_FP32_Complex* sum, const BGC_FP32_Complex* basic_number, const BGC_FP32_Complex* scalable_number, const float scale)
{
    sum->real = basic_number->real + scalable_number->real * scale;
    sum->imaginary = basic_number->imaginary + scalable_number->imaginary * scale;
}

inline void bgc_fp64_complex_add_scaled(BGC_FP64_Complex* sum, const BGC_FP64_Complex* basic_number, const BGC_FP64_Complex* scalable_number, const double scale)
{
    sum->real = basic_number->real + scalable_number->real * scale;
    sum->imaginary = basic_number->imaginary + scalable_number->imaginary * scale;
}

// ================== Subtract ================== //

inline void bgc_fp32_complex_subtract(BGC_FP32_Complex* difference, const BGC_FP32_Complex* minuend, const BGC_FP32_Complex* subtrahend)
{
    difference->real = minuend->real - subtrahend->real;
    difference->imaginary = minuend->imaginary - subtrahend->imaginary;
}

inline void bgc_fp64_complex_subtract(BGC_FP64_Complex* difference, const BGC_FP64_Complex* minuend, const BGC_FP64_Complex* subtrahend)
{
    difference->real = minuend->real - subtrahend->real;
    difference->imaginary = minuend->imaginary - subtrahend->imaginary;
}

// ================== Multiply ================== //

inline void bgc_fp32_complex_get_product(BGC_FP32_Complex* product, const BGC_FP32_Complex* number1, const BGC_FP32_Complex* number2)
{
    const float real = number1->real * number2->real - number1->imaginary * number2->imaginary;
    const float imaginary = number1->real * number2->imaginary + number1->imaginary * number2->real;

    product->real = real;
    product->imaginary = imaginary;
}

inline void bgc_fp64_complex_get_product(BGC_FP64_Complex* product, const BGC_FP64_Complex* number1, const BGC_FP64_Complex* number2)
{
    const double real = number1->real * number2->real - number1->imaginary * number2->imaginary;
    const double imaginary = number1->real * number2->imaginary + number1->imaginary * number2->real;

    product->real = real;
    product->imaginary = imaginary;
}

// ============= Multiply By Number ============= //

inline void bgc_fp32_complex_multiply(BGC_FP32_Complex* product, const BGC_FP32_Complex* multiplicand, const float multiplier)
{
    product->real = multiplicand->real * multiplier;
    product->imaginary = multiplicand->imaginary * multiplier;
}

inline void bgc_fp64_complex_multiply(BGC_FP64_Complex* product, const BGC_FP64_Complex* multiplicand, const double multiplier)
{
    product->real = multiplicand->real * multiplier;
    product->imaginary = multiplicand->imaginary * multiplier;
}

// =================== Divide =================== //

inline int bgc_fp32_complex_get_ratio(BGC_FP32_Complex* quotient, const BGC_FP32_Complex* divident, const BGC_FP32_Complex* divisor)
{
    const float square_modulus = bgc_fp32_complex_get_square_modulus(divisor);

    if (square_modulus <= BGC_FP32_SQUARE_EPSILON) {
        return 0;
    }

    const float real = divident->real * divisor->real + divident->imaginary * divisor->imaginary;
    const float imaginary = divident->imaginary * divisor->real - divident->real * divisor->imaginary;

    const float multiplier = 1.0f / square_modulus;

    quotient->real = real * multiplier;
    quotient->imaginary = imaginary * multiplier;

    return 1;
}

inline int bgc_fp64_complex_get_ratio(BGC_FP64_Complex* quotient, const BGC_FP64_Complex* divident, const BGC_FP64_Complex* divisor)
{
    const double square_modulus = bgc_fp64_complex_get_square_modulus(divisor);

    if (square_modulus <= BGC_FP64_SQUARE_EPSILON) {
        return 0;
    }

    const double real = divident->real * divisor->real + divident->imaginary * divisor->imaginary;
    const double imaginary = divident->imaginary * divisor->real - divident->real * divisor->imaginary;

    const double multiplier = 1.0 / square_modulus;

    quotient->real = real * multiplier;
    quotient->imaginary = imaginary * multiplier;

    return 1;
}

// ============== Divide By Number ============== //

inline void bgc_fp32_complex_divide(BGC_FP32_Complex* quotient, const BGC_FP32_Complex* dividend, const float divisor)
{
    bgc_fp32_complex_multiply(quotient, dividend, 1.0f / divisor);
}

inline void bgc_fp64_complex_divide(BGC_FP64_Complex* quotient, const BGC_FP64_Complex* dividend, const double divisor)
{
    bgc_fp64_complex_multiply(quotient, dividend, 1.0 / divisor);
}

// ================== Average2 ================== //

inline void bgc_fp32_complex_get_mean2(BGC_FP32_Complex* mean, const BGC_FP32_Complex* number1, const BGC_FP32_Complex* number2)
{
    mean->real = (number1->real + number2->real) * 0.5f;
    mean->imaginary = (number1->imaginary + number2->imaginary) * 0.5f;
}

inline void bgc_fp64_complex_get_mean2(BGC_FP64_Complex* mean, const BGC_FP64_Complex* number1, const BGC_FP64_Complex* number2)
{
    mean->real = (number1->real + number2->real) * 0.5;
    mean->imaginary = (number1->imaginary + number2->imaginary) * 0.5;
}

// ================== Average3 ================== //

inline void bgc_fp32_complex_get_mean3(BGC_FP32_Complex* mean, const BGC_FP32_Complex* number1, const BGC_FP32_Complex* number2, const BGC_FP32_Complex* number3)
{
    mean->real = (number1->real + number2->real + number3->real) * BGC_FP32_ONE_THIRD;
    mean->imaginary = (number1->imaginary + number2->imaginary + number3->imaginary) * BGC_FP32_ONE_THIRD;
}

inline void bgc_fp64_complex_get_mean3(BGC_FP64_Complex* mean, const BGC_FP64_Complex* number1, const BGC_FP64_Complex* number2, const BGC_FP64_Complex* number3)
{
    mean->real = (number1->real + number2->real + number3->real) * BGC_FP64_ONE_THIRD;
    mean->imaginary = (number1->imaginary + number2->imaginary + number3->imaginary) * BGC_FP64_ONE_THIRD;
}

// =================== Linear =================== //

inline void bgc_fp32_complex_interpolate(BGC_FP32_Complex* interpolation, const BGC_FP32_Complex* number1, const BGC_FP32_Complex* number2, const float phase)
{
    const float counter_phase = 1.0f - phase;

    interpolation->real = number1->real * counter_phase + number2->real * phase;
    interpolation->imaginary = number1->imaginary * counter_phase + number2->imaginary * phase;
}

inline void bgc_fp64_complex_interpolate(BGC_FP64_Complex* interpolation, const BGC_FP64_Complex* number1, const BGC_FP64_Complex* number2, const double phase)
{
    const double counter_phase = 1.0 - phase;

    interpolation->real = number1->real * counter_phase + number2->real * phase;
    interpolation->imaginary = number1->imaginary * counter_phase + number2->imaginary * phase;
}

// ================== Are Close ================= //

inline int bgc_fp32_complex_are_close(const BGC_FP32_Complex* number1, const BGC_FP32_Complex* number2)
{
    const float square_modulus1 = bgc_fp32_complex_get_square_modulus(number1);
    const float square_modulus2 = bgc_fp32_complex_get_square_modulus(number2);

    const float d_real = number1->real - number2->real;
    const float d_imaginary = number1->imaginary - number2->imaginary;

    const float square_distance = d_real * d_real + d_imaginary * d_imaginary;

    if (square_modulus1 <= BGC_FP32_EPSILON_EFFECTIVENESS_LIMIT || square_modulus2 <= BGC_FP32_EPSILON_EFFECTIVENESS_LIMIT) {
        return square_distance <= BGC_FP32_SQUARE_EPSILON;
    }

    return square_distance <= BGC_FP32_SQUARE_EPSILON * square_modulus1 && square_distance <= BGC_FP32_SQUARE_EPSILON * square_modulus2;
}

inline int bgc_fp64_complex_are_close(const BGC_FP64_Complex* number1, const BGC_FP64_Complex* number2)
{
    const double square_modulus1 = bgc_fp64_complex_get_square_modulus(number1);
    const double square_modulus2 = bgc_fp64_complex_get_square_modulus(number2);

    const double d_real = number1->real - number2->real;
    const double d_imaginary = number1->imaginary - number2->imaginary;

    const double square_distance = d_real * d_real + d_imaginary * d_imaginary;

    if (square_modulus1 <= BGC_FP64_EPSILON_EFFECTIVENESS_LIMIT || square_modulus2 <= BGC_FP64_EPSILON_EFFECTIVENESS_LIMIT) {
        return square_distance <= BGC_FP64_SQUARE_EPSILON;
    }

    return square_distance <= BGC_FP64_SQUARE_EPSILON * square_modulus1 && square_distance <= BGC_FP64_SQUARE_EPSILON * square_modulus2;
}

#endif