/***********************************************/
/*           Copyright (c) 2025 Belmu          */
/*             All Rights Reserved             */
/***********************************************/

/*
    [References]:
        ProofWiki. (2018). Definition:Inverse Sine/Complex. https://proofwiki.org/wiki/Definition:Inverse_Sine/Complex
        Jessie. (2021). open-source-utility-code. https://github.com/Jessie-LC/open-source-utility-code/blob/main/advanced/complexNumber.glsl
        Wikipedia. (2023). Trigonométrie complexe. https://fr.m.wikipedia.org/wiki/Trigonom%C3%A9trie_complexe
*/

struct complexFloat {
    float r;
    float i;
};

const complexFloat imaginaryUnit = complexFloat(0.0,  1.0);

float complexAbs(complexFloat z) {
    return sqrt(z.r * z.r + z.i * z.i);
}

complexFloat complexArg(complexFloat z) {
    return complexFloat(atan(z.i, z.r), 0.0);
}

complexFloat complexAdd(complexFloat a, complexFloat b) {
    complexFloat c;
    c.r = a.r + b.r;
    c.i = a.i + b.i;
    return c;
}

complexFloat complexAdd(float a, complexFloat b) {
    return complexAdd(complexFloat(a, 0.0), b);
}

complexFloat complexAdd(complexFloat a, float b) {
    return complexAdd(a, complexFloat(b, 0.0));
}

complexFloat complexSub(complexFloat a, complexFloat b) {
    complexFloat c;
    c.r = a.r - b.r;
    c.i = a.i - b.i;
    return c;
}

complexFloat complexSub(float a, complexFloat b) {
    return complexSub(complexFloat(a, 0.0), b);
}

complexFloat complexSub(complexFloat a, float b) {
    return complexSub(a, complexFloat(b, 0.0));
}

complexFloat complexMul(complexFloat a, complexFloat b) {
    complexFloat c;
    c.r = a.r * b.r - a.i * b.i;
    c.i = a.r * b.i + a.i * b.r;
    return c;
}

complexFloat complexMul(float a, complexFloat b) {
    return complexMul(complexFloat(a, 0.0), b);
}

complexFloat complexMul(complexFloat a, float b) {
    return complexMul(a, complexFloat(b, 0.0));
}

complexFloat complexDiv(complexFloat a, complexFloat b) {
    complexFloat c;
    if (abs(b.r) >= abs(b.i)) { // Safe scaling
        float r = b.i / b.r;
        float d = b.r + r * b.i;
        c.r = (a.r + r * a.i) / d;
        c.i = (a.i - r * a.r) / d;
    } else {
        float r = b.r / b.i;
        float d = b.i + r * b.r;
        c.r = (a.r * r + a.i) / d;
        c.i = (a.i * r - a.r) / d;
    }
    return c;
}

complexFloat complexDiv(float a, complexFloat b) {
    return complexDiv(complexFloat(a, 0.0), b);
}

complexFloat complexDiv(complexFloat a, float b) {
    return complexDiv(a, complexFloat(b, 0.0));
}

complexFloat complexSquare(complexFloat z) {
    return complexMul(z, z);
}

complexFloat complexExp(complexFloat z) {
    return complexMul(exp(z.r), complexFloat(cos(z.i), sin(z.i)));
}

complexFloat complexLn(complexFloat z) {
    // ln(z.r² + z.i²) * 0.5 + i arg(z)
    return complexFloat(log(z.r * z.r + z.i * z.i) * 0.5, atan(z.i, z.r));
}

complexFloat complexSin(complexFloat z) {
    // (e^iz - e^(-iz)) / 2i
    return complexMul(complexSub(complexExp(complexFloat(-z.i, z.r)), complexExp(complexFloat(z.i, -z.r))), complexFloat(0.0, 0.5));
}

complexFloat complexCos(complexFloat z) {
    // (e^iz + e^(-iz)) / 2
    return complexMul(complexAdd(complexExp(complexFloat(-z.i, z.r)), complexExp(complexFloat(z.i, -z.r))), 0.5);
}

complexFloat complexArcSin(complexFloat z) {
    complexFloat zSq         = complexSquare(z);
    complexFloat oneMinusZSq = complexSub(1.0, zSq);

    // Operands inside the natural logarithm ln(a + b * c)
    complexFloat a = complexFloat(-z.i, z.r);
    complexFloat b = complexFloat(sqrt(complexAbs(oneMinusZSq)), 0.0);
    complexFloat c = complexExp(complexMul(complexFloat(0.0, 0.5), complexArg(oneMinusZSq)));

    // ln(a + b * c) / i
    return complexDiv(complexLn(complexAdd(a, complexMul(b, c))), imaginaryUnit);
}
