#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_scalar_reset(BGC_FP32_DualScalar* number)
{
    number->real_part = 0.0f;
    number->dual_part = 0.0f;
}

inline void bgc_fp64_dual_scalar_reset(BGC_FP64_DualScalar* number)
{
    number->real_part = 0.0;
    number->dual_part = 0.0;
}

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

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

inline void bgc_fp64_dual_scalar_make(BGC_FP64_DualScalar* 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_scalar_copy(BGC_FP32_DualScalar* destination, const BGC_FP32_DualScalar* source)
{
    destination->real_part = source->real_part;
    destination->dual_part = source->dual_part;
}

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

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

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

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

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

inline void bgc_fp64_dual_scalar_convert_to_fp32(BGC_FP32_DualScalar* first, BGC_FP64_DualScalar* second)
{
    first->real_part = (float) second->real_part;
    first->dual_part = (float) second->dual_part;
}

inline void bgc_fp32_dual_scalar_convert_to_fp64(BGC_FP64_DualScalar* first, BGC_FP32_DualScalar* second)
{
    first->real_part = second->real_part;
    first->dual_part = second->dual_part;
}

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

inline void bgc_fp32_dual_scalar_revert(BGC_FP32_DualScalar* number)
{
    number->real_part = -number->real_part;
    number->dual_part = -number->dual_part;
}

inline void bgc_fp64_dual_scalar_revert(BGC_FP64_DualScalar* number)
{
    number->real_part = -number->real_part;
    number->dual_part = -number->dual_part;
}

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

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

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

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

inline void bgc_fp32_dual_scalar_conjugate(BGC_FP32_DualScalar* number)
{
    number->dual_part = -number->dual_part;
}

inline void bgc_fp64_dual_scalar_conjugate(BGC_FP64_DualScalar* number)
{
    number->dual_part = -number->dual_part;
}

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

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

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

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

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

inline void bgc_fp64_dual_scalar_add(BGC_FP64_DualScalar* sum, const BGC_FP64_DualScalar* first, const BGC_FP64_DualScalar* 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_scalar_add_scaled(BGC_FP32_DualScalar* sum, const BGC_FP32_DualScalar* base_number, const BGC_FP32_DualScalar* 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_scalar_add_scaled(BGC_FP64_DualScalar* sum, const BGC_FP64_DualScalar* base_number, const BGC_FP64_DualScalar* 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_scalar_subtract(BGC_FP32_DualScalar* difference, const BGC_FP32_DualScalar* minuend, const BGC_FP32_DualScalar* subtrahend)
{
    difference->real_part = minuend->real_part - subtrahend->real_part;
    difference->dual_part = minuend->dual_part - subtrahend->dual_part;
}

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

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

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

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

inline void bgc_fp32_dual_scalar_multiply_by_dual(BGC_FP32_DualScalar* product, const BGC_FP32_DualScalar* multiplicand, const BGC_FP32_DualScalar* 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_scalar_multiply_by_dual(BGC_FP64_DualScalar* product, const BGC_FP64_DualScalar* multiplicand, const BGC_FP64_DualScalar* 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_scalar_multiply_by_conjugate(BGC_FP32_DualScalar* product, const BGC_FP32_DualScalar* multiplicand, const BGC_FP32_DualScalar* 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_scalar_multiply_by_conjugate(BGC_FP64_DualScalar* product, const BGC_FP64_DualScalar* multiplicand, const BGC_FP64_DualScalar* 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_scalar_divide_by_real(BGC_FP32_DualScalar* quotient, const BGC_FP32_DualScalar* dividend, const float divisor)
{
    if (bgc_fp32_is_zero(divisor) || isnan(divisor)) {
        return BGC_FAILURE;
    }

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

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

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

inline int bgc_fp32_dual_scalar_divide_by_dual(BGC_FP32_DualScalar* quotient, const BGC_FP32_DualScalar* dividend, const BGC_FP32_DualScalar* 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_scalar_divide_by_dual(BGC_FP64_DualScalar* quotient, const BGC_FP64_DualScalar* dividend, const BGC_FP64_DualScalar* 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_scalar_divide_by_conjugate(BGC_FP32_DualScalar* quotient, const BGC_FP32_DualScalar* dividend, const BGC_FP32_DualScalar* 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_scalar_divide_by_conjugate(BGC_FP64_DualScalar* quotient, const BGC_FP64_DualScalar* dividend, const BGC_FP64_DualScalar* 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_scalar_get_mean2(BGC_FP32_DualScalar* mean, const BGC_FP32_DualScalar* first, const BGC_FP32_DualScalar* 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_scalar_get_mean2(BGC_FP64_DualScalar* mean, const BGC_FP64_DualScalar* first, const BGC_FP64_DualScalar* 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_scalar_get_mean3(BGC_FP32_DualScalar* mean, const BGC_FP32_DualScalar* first, const BGC_FP32_DualScalar* second, const BGC_FP32_DualScalar* 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_scalar_get_mean3(BGC_FP64_DualScalar* mean, const BGC_FP64_DualScalar* first, const BGC_FP64_DualScalar* second, const BGC_FP64_DualScalar* 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_scalar_interpolate(BGC_FP32_DualScalar* interpolation, const BGC_FP32_DualScalar* first, const BGC_FP32_DualScalar* 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_scalar_interpolate(BGC_FP64_DualScalar* interpolation, const BGC_FP64_DualScalar* first, const BGC_FP64_DualScalar* 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