bgc-c/basic-geometry/complex.h

#ifndef _BGC_COMPLEX_H_
#define _BGC_COMPLEX_H_

#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(const float real, const float imaginary, BGC_FP32_Complex* complex)
{
    complex->real = real;
    complex->imaginary = imaginary;
}

inline void bgc_fp64_complex_make(const double real, const double imaginary, BGC_FP64_Complex* complex)
{
    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_EPSYLON;
}

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

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(const BGC_FP32_Complex* source, BGC_FP32_Complex* destination)
{
    destination->real = source->real;
    destination->imaginary = source->imaginary;
}

inline void bgc_fp64_complex_copy(const BGC_FP64_Complex* source, BGC_FP64_Complex* destination)
{
    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(const BGC_FP64_Complex* source, BGC_FP32_Complex* destination)
{
    destination->real = (float)source->real;
    destination->imaginary = (float)source->imaginary;
}

inline void bgc_fp32_complex_convert_to_fp64(const BGC_FP32_Complex* source, BGC_FP64_Complex* destination)
{
    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(const BGC_FP32_Complex* number, BGC_FP32_Complex* opposite)
{
    opposite->real = -number->real;
    opposite->imaginary = -number->imaginary;
}

inline void bgc_fp64_complex_get_reverse(const BGC_FP64_Complex* number, BGC_FP64_Complex* opposite)
{
    opposite->real = -number->real;
    opposite->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_EPSYLON || 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_EPSYLON || 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(const BGC_FP32_Complex* number, BGC_FP32_Complex* normalized)
{
    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_EPSYLON || 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(const BGC_FP64_Complex* number, BGC_FP64_Complex* normalized)
{
    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_EPSYLON || 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(const BGC_FP32_Complex* number, BGC_FP32_Complex* conjugate)
{
    conjugate->real = number->real;
    conjugate->imaginary = -number->imaginary;
}

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

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

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

    if (square_modulus <= BGC_FP32_SQUARE_EPSYLON || 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(const BGC_FP64_Complex* number, BGC_FP64_Complex* inverse)
{
    const double square_modulus = bgc_fp64_complex_get_square_modulus(number);

    if (square_modulus <= BGC_FP64_SQUARE_EPSYLON || 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(const BGC_FP32_Complex* base, const float real_exponent, const float imaginary_exponent, BGC_FP32_Complex* power);

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

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

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

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

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

inline void bgc_fp32_complex_add_scaled(const BGC_FP32_Complex* basic_number, const BGC_FP32_Complex* scalable_number, const float scale, BGC_FP32_Complex* sum)
{
    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(const BGC_FP64_Complex* basic_number, const BGC_FP64_Complex* scalable_number, const double scale, BGC_FP64_Complex* sum)
{
    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(const BGC_FP32_Complex* minuend, const BGC_FP32_Complex* subtrahend, BGC_FP32_Complex* difference)
{
    difference->real = minuend->real - subtrahend->real;
    difference->imaginary = minuend->imaginary - subtrahend->imaginary;
}

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

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

inline void bgc_fp32_complex_get_product(const BGC_FP32_Complex* number1, const BGC_FP32_Complex* number2, BGC_FP32_Complex* product)
{
    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(const BGC_FP64_Complex* number1, const BGC_FP64_Complex* number2, BGC_FP64_Complex* product)
{
    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(const BGC_FP32_Complex* multiplicand, const float multiplier, BGC_FP32_Complex* product)
{
    product->real = multiplicand->real * multiplier;
    product->imaginary = multiplicand->imaginary * multiplier;
}

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

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

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

    if (square_modulus <= BGC_FP32_SQUARE_EPSYLON) {
        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(const BGC_FP64_Complex* divident, const BGC_FP64_Complex* divisor, BGC_FP64_Complex* quotient)
{
    const double square_modulus = bgc_fp64_complex_get_square_modulus(divisor);

    if (square_modulus <= BGC_FP64_SQUARE_EPSYLON) {
        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(const BGC_FP32_Complex* dividend, const float divisor, BGC_FP32_Complex* quotient)
{
    bgc_fp32_complex_multiply(dividend, 1.0f / divisor, quotient);
}

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

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

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

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

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

inline void bgc_fp32_complex_get_mean3(const BGC_FP32_Complex* number1, const BGC_FP32_Complex* number2, const BGC_FP32_Complex* number3, BGC_FP32_Complex* mean)
{
    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(const BGC_FP64_Complex* number1, const BGC_FP64_Complex* number2, const BGC_FP64_Complex* number3, BGC_FP64_Complex* mean)
{
    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(const BGC_FP32_Complex* number1, const BGC_FP32_Complex* number2, const float phase, BGC_FP32_Complex* interpolation)
{
    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(const BGC_FP64_Complex* number1, const BGC_FP64_Complex* number2, const double phase, BGC_FP64_Complex* interpolation)
{
    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_EPSYLON_EFFECTIVENESS_LIMIT || square_modulus2 <= BGC_FP32_EPSYLON_EFFECTIVENESS_LIMIT) {
        return square_distance <= BGC_FP32_SQUARE_EPSYLON;
    }

    return square_distance <= BGC_FP32_SQUARE_EPSYLON * square_modulus1 && square_distance <= BGC_FP32_SQUARE_EPSYLON * 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_EPSYLON_EFFECTIVENESS_LIMIT || square_modulus2 <= BGC_FP64_EPSYLON_EFFECTIVENESS_LIMIT) {
        return square_distance <= BGC_FP64_SQUARE_EPSYLON;
    }

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

#endif