#ifndef _BGC_DUAL_NUMBER_H_INCLUDED_
#define _BGC_DUAL_NUMBER_H_INCLUDED_

#include <math.h>

#include "./utilities.h"
#include "./types.h"

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

inline void bgc_fp32_dual_number_reset(BGC_FP32_DualNumber* number)
{
    number->real_part = 0.0f;
    number->dual_part = 0.0f;
}

inline void bgc_fp64_dual_number_reset(BGC_FP64_DualNumber* number)
{
    number->real_part = 0.0;
    number->dual_part = 0.0;
}

// ==================== Make ==================== //

inline void bgc_fp32_dual_number_make(BGC_FP32_DualNumber* number, const float real_part, const float dual_part)
{
    number->real_part = real_part;
    number->dual_part = dual_part;
}

inline void bgc_fp64_dual_number_make(BGC_FP64_DualNumber* number, const double real_part, const double dual_part)
{
    number->real_part = real_part;
    number->dual_part = dual_part;
}

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

inline void bgc_fp32_dual_number_copy(BGC_FP32_DualNumber* destination, const BGC_FP32_DualNumber* source)
{
    destination->real_part = source->real_part;
    destination->dual_part = source->dual_part;
}

inline void bgc_fp64_dual_number_copy(BGC_FP64_DualNumber* destination, const BGC_FP64_DualNumber* source)
{
    destination->real_part = source->real_part;
    destination->dual_part = source->dual_part;
}

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

inline void bgc_fp32_dual_number_swap(BGC_FP32_DualNumber* first, BGC_FP32_DualNumber* second)
{
    first->real_part = second->real_part;
    first->dual_part = second->dual_part;
}

inline void bgc_fp64_dual_number_swap(BGC_FP64_DualNumber* first, BGC_FP64_DualNumber* second)
{
    first->real_part = second->real_part;
    first->dual_part = second->dual_part;
}

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

inline double bgc_fp32_dual_number_get_modulus(const BGC_FP32_DualNumber* number)
{
    return fabsf(number->real_part);
}

inline double bgc_fp64_dual_number_get_modulus(const BGC_FP64_DualNumber* number)
{
    return fabs(number->real_part);
}

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

inline void bgc_fp64_dual_number_convert_to_fp32(BGC_FP32_DualNumber* first, BGC_FP64_DualNumber* second)
{
    first->real_part = (float) second->real_part;
    first->dual_part = (float) second->dual_part;
}

inline void bgc_fp32_dual_number_convert_to_fp64(BGC_FP64_DualNumber* first, BGC_FP32_DualNumber* second)
{
    first->real_part = second->real_part;
    first->dual_part = second->dual_part;
}

// =================== Revert =================== //

inline void bgc_fp32_dual_number_revert(BGC_FP32_DualNumber* number)
{
    number->real_part = -number->real_part;
    number->dual_part = -number->dual_part;
}

inline void bgc_fp64_dual_number_revert(BGC_FP64_DualNumber* number)
{
    number->real_part = -number->real_part;
    number->dual_part = -number->dual_part;
}

// ================ Get Reverse ================= //

inline void bgc_fp32_dual_number_get_reverse(BGC_FP32_DualNumber* reverse, const BGC_FP32_DualNumber* number)
{
    reverse->real_part = -number->real_part;
    reverse->dual_part = -number->dual_part;
}

inline void bgc_fp64_dual_number_get_reverse(BGC_FP64_DualNumber* reverse, const BGC_FP64_DualNumber* number)
{
    reverse->real_part = -number->real_part;
    reverse->dual_part = -number->dual_part;
}

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

inline int bgc_fp32_dual_number_invert(BGC_FP32_DualNumber* number)
{
    if (bgc_fp32_is_zero(number->real_part) || isnan(number->real_part)) {
        return BGC_FAILURE;
    }

    const float multiplicator = 1.0f / (number->real_part * number->real_part);

    number->real_part =  number->real_part * multiplicator;
    number->dual_part = -number->dual_part * multiplicator;

    return BGC_SUCCESS;
}

inline int bgc_fp64_dual_number_invert(BGC_FP64_DualNumber* number)
{
    if (bgc_fp64_is_zero(number->real_part) || isnan(number->real_part)) {
        return BGC_FAILURE;
    }

    const double multiplicator = 1.0 / (number->real_part * number->real_part);

    number->real_part =  number->real_part * multiplicator;
    number->dual_part = -number->dual_part * multiplicator;

    return BGC_SUCCESS;
}

// ================ Get Inverse ================= //

inline int bgc_fp32_dual_number_get_inverse(BGC_FP32_DualNumber* inverse, const BGC_FP32_DualNumber* number)
{
    if (bgc_fp32_is_zero(number->real_part) || isnan(number->real_part)) {
        inverse->real_part = 0.0f;
        inverse->dual_part = 0.0f;
        return BGC_FAILURE;
    }

    const float multiplicator = 1.0f / (number->real_part * number->real_part);

    inverse->real_part =  number->real_part * multiplicator;
    inverse->dual_part = -number->dual_part * multiplicator;

    return BGC_SUCCESS;
}

inline int bgc_fp64_dual_number_get_inverse(BGC_FP64_DualNumber* inverse, const BGC_FP64_DualNumber* number)
{
    if (bgc_fp64_is_zero(number->real_part) || isnan(number->real_part)) {
        inverse->real_part = 0.0;
        inverse->dual_part = 0.0;
        return BGC_FAILURE;
    }

    const double multiplicator = 1.0 / (number->real_part * number->real_part);

    inverse->real_part =  number->real_part * multiplicator;
    inverse->dual_part = -number->dual_part * multiplicator;

    return BGC_SUCCESS;
}

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

inline void bgc_fp32_dual_number_conjugate(BGC_FP32_DualNumber* number)
{
    number->dual_part = -number->dual_part;
}

inline void bgc_fp64_dual_number_conjugate(BGC_FP64_DualNumber* number)
{
    number->dual_part = -number->dual_part;
}

// =============== Get Conjugate ================ //

inline void bgc_fp32_dual_number_get_conjugate(BGC_FP32_DualNumber* conjugate, const BGC_FP32_DualNumber* number)
{
    conjugate->real_part =  number->real_part;
    conjugate->dual_part = -number->dual_part;
}

inline void bgc_fp64_dual_number_get_conjugate(BGC_FP64_DualNumber* conjugate, const BGC_FP64_DualNumber* number)
{
    conjugate->real_part =  number->real_part;
    conjugate->dual_part = -number->dual_part;
}

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

inline void bgc_fp32_dual_number_add(BGC_FP32_DualNumber* sum, const BGC_FP32_DualNumber* first, const BGC_FP32_DualNumber* second)
{
    sum->real_part = first->real_part + second->real_part;
    sum->dual_part = first->dual_part + second->dual_part;
}

inline void bgc_fp64_dual_number_add(BGC_FP64_DualNumber* sum, const BGC_FP64_DualNumber* first, const BGC_FP64_DualNumber* second)
{
    sum->real_part = first->real_part + second->real_part;
    sum->dual_part = first->dual_part + second->dual_part;
}

// ================= Add Scaled ================= //

inline void bgc_fp32_dual_number_add_scaled(BGC_FP32_DualNumber* sum, const BGC_FP32_DualNumber* base_number, const BGC_FP32_DualNumber* scalable_number, const float scale)
{
    sum->real_part = base_number->real_part + scalable_number->real_part * scale;
    sum->dual_part = base_number->dual_part + scalable_number->dual_part * scale;
}

inline void bgc_fp64_dual_number_add_scaled(BGC_FP64_DualNumber* sum, const BGC_FP64_DualNumber* base_number, const BGC_FP64_DualNumber* scalable_number, const double scale)
{
    sum->real_part = base_number->real_part + scalable_number->real_part * scale;
    sum->dual_part = base_number->dual_part + scalable_number->dual_part * scale;
}

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

inline void bgc_fp32_dual_number_subtract(BGC_FP32_DualNumber* difference, const BGC_FP32_DualNumber* minuend, const BGC_FP32_DualNumber* subtrahend)
{
    difference->real_part = minuend->real_part - subtrahend->real_part;
    difference->dual_part = minuend->dual_part - subtrahend->dual_part;
}

inline void bgc_fp64_dual_number_subtract(BGC_FP64_DualNumber* difference, const BGC_FP64_DualNumber* minuend, const BGC_FP64_DualNumber* subtrahend)
{
    difference->real_part = minuend->real_part - subtrahend->real_part;
    difference->dual_part = minuend->dual_part - subtrahend->dual_part;
}

// ============== Subtract Scaled =============== //

inline void bgc_fp32_dual_number_subtract_scaled(BGC_FP32_DualNumber* difference, const BGC_FP32_DualNumber* base_number, const BGC_FP32_DualNumber* scalable_number, const float scale)
{
    difference->real_part = base_number->real_part - scalable_number->real_part * scale;
    difference->dual_part = base_number->dual_part - scalable_number->dual_part * scale;
}

inline void bgc_fp64_dual_number_subtract_scaled(BGC_FP64_DualNumber* difference, const BGC_FP64_DualNumber* base_number, const BGC_FP64_DualNumber* scalable_number, const double scale)
{
    difference->real_part = base_number->real_part - scalable_number->real_part * scale;
    difference->dual_part = base_number->dual_part - scalable_number->dual_part * scale;
}

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

inline void bgc_fp32_dual_number_multiply_by_real(BGC_FP32_DualNumber* product, const BGC_FP32_DualNumber* multiplicand, const float multiplier)
{
    product->real_part = multiplicand->real_part * multiplier;
    product->dual_part = multiplicand->dual_part * multiplier;
}

inline void bgc_fp64_dual_number_multiply_by_real(BGC_FP64_DualNumber* product, const BGC_FP64_DualNumber* multiplicand, const double multiplier)
{
    product->real_part = multiplicand->real_part * multiplier;
    product->dual_part = multiplicand->dual_part * multiplier;
}

inline void bgc_fp32_dual_number_multiply_by_dual(BGC_FP32_DualNumber* product, const BGC_FP32_DualNumber* multiplicand, const BGC_FP32_DualNumber* multiplier)
{
    const float real_part = multiplicand->real_part * multiplier->real_part;
    const float dual_part = multiplicand->dual_part * multiplier->real_part + multiplicand->real_part * multiplier->dual_part;

    product->real_part = real_part;
    product->dual_part = dual_part;
}

inline void bgc_fp64_dual_number_multiply_by_dual(BGC_FP64_DualNumber* product, const BGC_FP64_DualNumber* multiplicand, const BGC_FP64_DualNumber* multiplier)
{
    const double real_part = multiplicand->real_part * multiplier->real_part;
    const double dual_part = multiplicand->dual_part * multiplier->real_part + multiplicand->real_part * multiplier->dual_part;

    product->real_part = real_part;
    product->dual_part = dual_part;
}

inline void bgc_fp32_dual_number_multiply_by_conjugate(BGC_FP32_DualNumber* product, const BGC_FP32_DualNumber* multiplicand, const BGC_FP32_DualNumber* multiplier_to_conjugate)
{
    const float real_part = multiplicand->real_part * multiplier_to_conjugate->real_part;
    const float dual_part = multiplicand->dual_part * multiplier_to_conjugate->real_part - multiplicand->real_part * multiplier_to_conjugate->dual_part;

    product->real_part = real_part;
    product->dual_part = dual_part;
}

inline void bgc_fp64_dual_number_multiply_by_conjugate(BGC_FP64_DualNumber* product, const BGC_FP64_DualNumber* multiplicand, const BGC_FP64_DualNumber* multiplier_to_conjugate)
{
    const double real_part = multiplicand->real_part * multiplier_to_conjugate->real_part;
    const double dual_part = multiplicand->dual_part * multiplier_to_conjugate->real_part - multiplicand->real_part * multiplier_to_conjugate->dual_part;

    product->real_part = real_part;
    product->dual_part = dual_part;
}
// =================== Divide =================== //

inline int bgc_fp32_dual_number_divide_by_real(BGC_FP32_DualNumber* quotient, const BGC_FP32_DualNumber* dividend, const float divisor)
{
    if (bgc_fp32_is_zero(divisor) || isnan(divisor)) {
        return BGC_FAILURE;
    }

    bgc_fp32_dual_number_multiply_by_real(quotient, dividend, 1.0f / divisor);
    return BGC_SUCCESS;
}

inline int bgc_fp64_dual_number_divide_by_real(BGC_FP64_DualNumber* quotient, const BGC_FP64_DualNumber* dividend, const double divisor)
{
    if (bgc_fp64_is_zero(divisor) || isnan(divisor)) {
        return BGC_FAILURE;
    }

    bgc_fp64_dual_number_multiply_by_real(quotient, dividend, 1.0 / divisor);
    return BGC_SUCCESS;
}

inline int bgc_fp32_dual_number_divide_by_dual(BGC_FP32_DualNumber* quotient, const BGC_FP32_DualNumber* dividend, const BGC_FP32_DualNumber* divisor)
{
    if (bgc_fp32_is_zero(divisor->real_part)) {
        return BGC_FAILURE;
    }

    const float multiplier = 1.0f / divisor->real_part;

    const float real_part = dividend->real_part * multiplier;
    const float dual_part = dividend->dual_part * multiplier - (dividend->real_part * multiplier) * (divisor->dual_part * multiplier);

    quotient->real_part = real_part;
    quotient->dual_part = dual_part;

    return BGC_SUCCESS;
}

inline int bgc_fp64_dual_number_divide_by_dual(BGC_FP64_DualNumber* quotient, const BGC_FP64_DualNumber* dividend, const BGC_FP64_DualNumber* divisor)
{
    if (bgc_fp64_is_zero(divisor->real_part)) {
        return BGC_FAILURE;
    }

    const double multiplier = 1.0 / divisor->real_part;

    const double real_part = dividend->real_part * multiplier;
    const double dual_part = dividend->dual_part * multiplier - (dividend->real_part * multiplier) * (divisor->dual_part * multiplier);

    quotient->real_part = real_part;
    quotient->dual_part = dual_part;

    return BGC_SUCCESS;
}

inline int bgc_fp32_dual_number_divide_by_conjugate(BGC_FP32_DualNumber* quotient, const BGC_FP32_DualNumber* dividend, const BGC_FP32_DualNumber* divisor_to_conjugate)
{
    if (bgc_fp32_is_zero(divisor_to_conjugate->real_part)) {
        return BGC_FAILURE;
    }

    const float multiplier = 1.0f / divisor_to_conjugate->real_part;

    const float real_part = dividend->real_part * multiplier;
    const float dual_part = dividend->dual_part * multiplier + (dividend->real_part * multiplier) * (divisor_to_conjugate->dual_part * multiplier);

    quotient->real_part = real_part;
    quotient->dual_part = dual_part;

    return BGC_SUCCESS;
}

inline int bgc_fp64_dual_number_divide_by_conjugate(BGC_FP64_DualNumber* quotient, const BGC_FP64_DualNumber* dividend, const BGC_FP64_DualNumber* divisor_to_conjugate)
{
    if (bgc_fp64_is_zero(divisor_to_conjugate->real_part)) {
        return BGC_FAILURE;
    }

    const double multiplier = 1.0 / divisor_to_conjugate->real_part;

    const double real_part = dividend->real_part * multiplier;
    const double dual_part = dividend->dual_part * multiplier + (dividend->real_part * multiplier) * (divisor_to_conjugate->dual_part * multiplier);

    quotient->real_part = real_part;
    quotient->dual_part = dual_part;

    return BGC_SUCCESS;
}

// ================ Mean of Two ================= //

inline void bgc_fp32_dual_number_get_mean2(BGC_FP32_DualNumber* mean, const BGC_FP32_DualNumber* first, const BGC_FP32_DualNumber* second)
{
    mean->real_part = (first->real_part + second->real_part) * 0.5f;
    mean->dual_part = (first->dual_part + second->dual_part) * 0.5f;
}

inline void bgc_fp64_dual_number_get_mean2(BGC_FP64_DualNumber* mean, const BGC_FP64_DualNumber* first, const BGC_FP64_DualNumber* second)
{
    mean->real_part = (first->real_part + second->real_part) * 0.5;
    mean->dual_part = (first->dual_part + second->dual_part) * 0.5;
}

// =============== Mean of Three ================ //

inline void bgc_fp32_dual_number_get_mean3(BGC_FP32_DualNumber* mean, const BGC_FP32_DualNumber* first, const BGC_FP32_DualNumber* second, const BGC_FP32_DualNumber* third)
{
    mean->real_part = (first->real_part + second->real_part + third->real_part) * BGC_FP32_ONE_THIRD;
    mean->dual_part = (first->dual_part + second->dual_part + third->dual_part) * BGC_FP32_ONE_THIRD;
}

inline void bgc_fp64_dual_number_get_mean3(BGC_FP64_DualNumber* mean, const BGC_FP64_DualNumber* first, const BGC_FP64_DualNumber* second, const BGC_FP64_DualNumber* third)
{
    mean->real_part = (first->real_part + second->real_part + third->real_part) * BGC_FP64_ONE_THIRD;
    mean->dual_part = (first->dual_part + second->dual_part + third->dual_part) * BGC_FP64_ONE_THIRD;
}

// ============ Linear Interpolation ============ //

inline void bgc_fp32_dual_number_interpolate(BGC_FP32_DualNumber* interpolation, const BGC_FP32_DualNumber* first, const BGC_FP32_DualNumber* second, const float phase)
{
    const float counter_phase = 1.0f - phase;

    interpolation->real_part = first->real_part * counter_phase + second->real_part * phase;
    interpolation->dual_part = first->dual_part * counter_phase + second->dual_part * phase;
}

inline void bgc_fp64_dual_number_interpolate(BGC_FP64_DualNumber* interpolation, const BGC_FP64_DualNumber* first, const BGC_FP64_DualNumber* second, const double phase)
{
    const double counter_phase = 1.0 - phase;

    interpolation->real_part = first->real_part * counter_phase + second->real_part * phase;
    interpolation->dual_part = first->dual_part * counter_phase + second->dual_part * phase;
}

#endif