/* -*- mode: C -*-  */
/*
   IGraph library.
   Copyright (C) 2010-2012  Gabor Csardi <csardi.gabor@gmail.com>
   334 Harvard street, Cambridge, MA 02139 USA

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
   02110-1301 USA

*/

#include "igraph_complex.h"
#include "igraph_math.h"
#include <math.h>

/**
 * \example igraph_complex.c
 */

igraph_complex_t igraph_complex(igraph_real_t x, igraph_real_t y) {
    igraph_complex_t res;
    IGRAPH_REAL(res) = x;
    IGRAPH_IMAG(res) = y;
    return res;
}

igraph_complex_t igraph_complex_polar(igraph_real_t r, igraph_real_t theta) {
    igraph_complex_t res;
    IGRAPH_REAL(res) = r * cos(theta);
    IGRAPH_IMAG(res) = r * sin(theta);
    return res;
}

igraph_bool_t igraph_complex_eq_tol(igraph_complex_t z1,
                                    igraph_complex_t z2,
                                    igraph_real_t tol) {
    if (fabs(IGRAPH_REAL(z1) - IGRAPH_REAL(z2)) > tol ||
        fabs(IGRAPH_IMAG(z1) - IGRAPH_IMAG(z2)) > tol) {
        return 0;
    }
    return 1;
}

igraph_real_t igraph_complex_mod(igraph_complex_t z) {
    igraph_real_t x = IGRAPH_REAL(z);
    igraph_real_t y = IGRAPH_IMAG(z);
    return hypot(x, y);
}

igraph_real_t igraph_complex_arg(igraph_complex_t z) {
    igraph_real_t x = IGRAPH_REAL(z);
    igraph_real_t y = IGRAPH_IMAG(z);
    if (x == 0.0 && y == 0.0) {
        return 0.0;
    }
    return atan2(y, x);
}

igraph_complex_t igraph_complex_add(igraph_complex_t z1,
                                    igraph_complex_t z2) {
    igraph_complex_t res;
    IGRAPH_REAL(res) = IGRAPH_REAL(z1) + IGRAPH_REAL(z2);
    IGRAPH_IMAG(res) = IGRAPH_IMAG(z1) + IGRAPH_IMAG(z2);
    return res;
}

igraph_complex_t igraph_complex_sub(igraph_complex_t z1,
                                    igraph_complex_t z2) {
    igraph_complex_t res;
    IGRAPH_REAL(res) = IGRAPH_REAL(z1) - IGRAPH_REAL(z2);
    IGRAPH_IMAG(res) = IGRAPH_IMAG(z1) - IGRAPH_IMAG(z2);
    return res;
}

igraph_complex_t igraph_complex_mul(igraph_complex_t z1,
                                    igraph_complex_t z2) {
    igraph_complex_t res;
    IGRAPH_REAL(res) = IGRAPH_REAL(z1) * IGRAPH_REAL(z2) -
                       IGRAPH_IMAG(z1) * IGRAPH_IMAG(z2);
    IGRAPH_IMAG(res) = IGRAPH_REAL(z1) * IGRAPH_IMAG(z2) +
                       IGRAPH_IMAG(z1) * IGRAPH_REAL(z2);
    return res;
}

igraph_complex_t igraph_complex_div(igraph_complex_t z1,
                                    igraph_complex_t z2) {
    igraph_complex_t res;
    igraph_real_t z1r = IGRAPH_REAL(z1), z1i = IGRAPH_IMAG(z1);
    igraph_real_t z2r = IGRAPH_REAL(z2), z2i = IGRAPH_IMAG(z2);
    igraph_real_t s = 1.0 / igraph_complex_abs(z2);
    igraph_real_t sz2r = s * z2r;
    igraph_real_t sz2i = s * z2i;
    IGRAPH_REAL(res) = (z1r * sz2r + z1i * sz2i) * s;
    IGRAPH_IMAG(res) = (z1i * sz2r - z1r * sz2i) * s;
    return res;
}

igraph_complex_t igraph_complex_add_real(igraph_complex_t z,
        igraph_real_t x) {
    igraph_complex_t res;
    IGRAPH_REAL(res) = IGRAPH_REAL(z) + x;
    IGRAPH_IMAG(res) = IGRAPH_IMAG(z);
    return res;
}

igraph_complex_t igraph_complex_add_imag(igraph_complex_t z,
        igraph_real_t y) {
    igraph_complex_t res;
    IGRAPH_REAL(res) = IGRAPH_REAL(z);
    IGRAPH_IMAG(res) = IGRAPH_IMAG(z) + y;
    return res;
}

igraph_complex_t igraph_complex_sub_real(igraph_complex_t z,
        igraph_real_t x) {
    igraph_complex_t res;
    IGRAPH_REAL(res) = IGRAPH_REAL(z) - x;
    IGRAPH_IMAG(res) = IGRAPH_IMAG(z);
    return res;
}

igraph_complex_t igraph_complex_sub_imag(igraph_complex_t z,
        igraph_real_t y) {
    igraph_complex_t res;
    IGRAPH_REAL(res) = IGRAPH_REAL(z);
    IGRAPH_IMAG(res) = IGRAPH_IMAG(z) - y;
    return res;
}

igraph_complex_t igraph_complex_mul_real(igraph_complex_t z,
        igraph_real_t x) {
    igraph_complex_t res;
    IGRAPH_REAL(res) = IGRAPH_REAL(z) * x;
    IGRAPH_IMAG(res) = IGRAPH_IMAG(z) * x;
    return res;
}

igraph_complex_t igraph_complex_mul_imag(igraph_complex_t z,
        igraph_real_t y) {
    igraph_complex_t res;
    IGRAPH_REAL(res) = - IGRAPH_IMAG(z) * y;
    IGRAPH_IMAG(res) =   IGRAPH_REAL(z) * y;
    return res;
}

igraph_complex_t igraph_complex_div_real(igraph_complex_t z,
        igraph_real_t x) {
    igraph_complex_t res;
    IGRAPH_REAL(res) = IGRAPH_REAL(z) / x;
    IGRAPH_IMAG(res) = IGRAPH_IMAG(z) / x;
    return res;
}

igraph_complex_t igraph_complex_div_imag(igraph_complex_t z,
        igraph_real_t y) {
    igraph_complex_t res;
    IGRAPH_REAL(res) =   IGRAPH_IMAG(z) / y;
    IGRAPH_IMAG(res) = - IGRAPH_REAL(z) / y;
    return res;
}

igraph_complex_t igraph_complex_conj(igraph_complex_t z) {
    igraph_complex_t res;
    IGRAPH_REAL(res) =   IGRAPH_REAL(z);
    IGRAPH_IMAG(res) = - IGRAPH_IMAG(z);
    return res;
}

igraph_complex_t igraph_complex_neg(igraph_complex_t z) {
    igraph_complex_t res;
    IGRAPH_REAL(res) = - IGRAPH_REAL(z);
    IGRAPH_IMAG(res) = - IGRAPH_IMAG(z);
    return res;
}

igraph_complex_t igraph_complex_inv(igraph_complex_t z) {
    igraph_complex_t res;
    igraph_real_t s = 1.0 / igraph_complex_abs(z);
    IGRAPH_REAL(res) =   (IGRAPH_REAL(z) * s) * s;
    IGRAPH_IMAG(res) = - (IGRAPH_IMAG(z) * s) * s;
    return res;
}

igraph_real_t igraph_complex_abs(igraph_complex_t z) {
    return hypot(IGRAPH_REAL(z), IGRAPH_IMAG(z));
}

igraph_real_t igraph_complex_logabs(igraph_complex_t z) {
    igraph_real_t xabs = fabs(IGRAPH_REAL(z));
    igraph_real_t yabs = fabs(IGRAPH_IMAG(z));
    igraph_real_t max, u;
    if (xabs >= yabs) {
        max = xabs;
        u = yabs / xabs;
    } else {
        max = yabs;
        u = xabs / yabs;
    }
    return log (max) + 0.5 * log1p (u * u);
}

igraph_complex_t igraph_complex_sqrt(igraph_complex_t z) {
    igraph_complex_t res;

    if (IGRAPH_REAL(z) == 0.0 && IGRAPH_IMAG(z) == 0.0) {
        IGRAPH_REAL(res) = IGRAPH_IMAG(res) = 0.0;
    } else {
        igraph_real_t x = fabs (IGRAPH_REAL(z));
        igraph_real_t y = fabs (IGRAPH_IMAG(z));
        igraph_real_t w;
        if (x >= y)  {
            igraph_real_t t = y / x;
            w = sqrt (x) * sqrt (0.5 * (1.0 + sqrt (1.0 + t * t)));
        } else {
            igraph_real_t t = x / y;
            w = sqrt (y) * sqrt (0.5 * (t + sqrt (1.0 + t * t)));
        }

        if (IGRAPH_REAL(z) >= 0.0) {
            igraph_real_t ai = IGRAPH_IMAG(z);
            IGRAPH_REAL(res) = w;
            IGRAPH_IMAG(res) = ai / (2.0 * w);
        } else {
            igraph_real_t ai = IGRAPH_IMAG(z);
            igraph_real_t vi = (ai >= 0) ? w : -w;
            IGRAPH_REAL(res) = ai / (2.0 * vi);
            IGRAPH_IMAG(res) = vi;
        }
    }

    return res;
}

igraph_complex_t igraph_complex_sqrt_real(igraph_real_t x) {
    igraph_complex_t res;
    if (x >= 0) {
        IGRAPH_REAL(res) = sqrt(x);
        IGRAPH_IMAG(res) = 0.0;
    } else {
        IGRAPH_REAL(res) = 0.0;
        IGRAPH_IMAG(res) = sqrt(-x);
    }
    return res;
}

igraph_complex_t igraph_complex_exp(igraph_complex_t z) {
    igraph_real_t rho   = exp(IGRAPH_REAL(z));
    igraph_real_t theta = IGRAPH_IMAG(z);
    igraph_complex_t res;
    IGRAPH_REAL(res) = rho * cos(theta);
    IGRAPH_IMAG(res) = rho * sin(theta);
    return res;
}

igraph_complex_t igraph_complex_pow(igraph_complex_t z1,
                                    igraph_complex_t z2) {
    igraph_complex_t res;

    if (IGRAPH_REAL(z1) == 0 && IGRAPH_IMAG(z1) == 0.0) {
        if (IGRAPH_REAL(z2) == 0 && IGRAPH_IMAG(z2) == 0.0) {
            IGRAPH_REAL(res) = 1.0;
            IGRAPH_IMAG(res) = 0.0;
        } else {
            IGRAPH_REAL(res) = IGRAPH_IMAG(res) = 0.0;
        }
    } else if (IGRAPH_REAL(z2) == 1.0 && IGRAPH_IMAG(z2) == 0.0) {
        IGRAPH_REAL(res) = IGRAPH_REAL(z1);
        IGRAPH_IMAG(res) = IGRAPH_IMAG(z1);
    } else if (IGRAPH_REAL(z2) == -1.0 && IGRAPH_IMAG(z2) == 0.0) {
        res = igraph_complex_inv(z1);
    } else {
        igraph_real_t logr = igraph_complex_logabs (z1);
        igraph_real_t theta = igraph_complex_arg (z1);
        igraph_real_t z2r = IGRAPH_REAL(z2), z2i = IGRAPH_IMAG(z2);
        igraph_real_t rho = exp (logr * z2r - z2i * theta);
        igraph_real_t beta = theta * z2r + z2i * logr;
        IGRAPH_REAL(res) = rho * cos(beta);
        IGRAPH_IMAG(res) = rho * sin(beta);
    }

    return res;
}

igraph_complex_t igraph_complex_pow_real(igraph_complex_t z,
        igraph_real_t x) {
    igraph_complex_t res;
    if (IGRAPH_REAL(z) == 0.0 && IGRAPH_IMAG(z) == 0.0) {
        if (x == 0) {
            IGRAPH_REAL(res) = 1.0;
            IGRAPH_IMAG(res) = 0.0;
        } else {
            IGRAPH_REAL(res) = IGRAPH_IMAG(res) = 0.0;
        }
    } else {
        igraph_real_t logr = igraph_complex_logabs(z);
        igraph_real_t theta = igraph_complex_arg(z);
        igraph_real_t rho = exp (logr * x);
        igraph_real_t beta = theta * x;
        IGRAPH_REAL(res) = rho * cos(beta);
        IGRAPH_IMAG(res) = rho * sin(beta);
    }
    return res;
}

igraph_complex_t igraph_complex_log(igraph_complex_t z) {
    igraph_complex_t res;
    IGRAPH_REAL(res) = igraph_complex_logabs(z);
    IGRAPH_IMAG(res) = igraph_complex_arg(z);
    return res;
}

igraph_complex_t igraph_complex_log10(igraph_complex_t z) {
    return igraph_complex_mul_real(igraph_complex_log(z), 1 / log(10.0));
}

igraph_complex_t igraph_complex_log_b(igraph_complex_t z,
                                      igraph_complex_t b) {
    return igraph_complex_div (igraph_complex_log(z), igraph_complex_log(b));
}

igraph_complex_t igraph_complex_sin(igraph_complex_t z) {
    igraph_real_t zr = IGRAPH_REAL(z);
    igraph_real_t zi = IGRAPH_IMAG(z);
    igraph_complex_t res;
    if (zi == 0.0) {
        IGRAPH_REAL(res) = sin(zr);
        IGRAPH_IMAG(res) = 0.0;
    } else {
        IGRAPH_REAL(res) = sin(zr) * cosh(zi);
        IGRAPH_IMAG(res) = cos(zr) * sinh(zi);
    }
    return res;
}

igraph_complex_t igraph_complex_cos(igraph_complex_t z) {
    igraph_real_t zr = IGRAPH_REAL(z);
    igraph_real_t zi = IGRAPH_IMAG(z);
    igraph_complex_t res;
    if (zi == 0.0) {
        IGRAPH_REAL(res) = cos(zr);
        IGRAPH_IMAG(res) = 0.0;
    } else {
        IGRAPH_REAL(res) = cos(zr) * cosh(zi);
        IGRAPH_IMAG(res) = sin(zr) * sinh(-zi);
    }
    return res;
}

igraph_complex_t igraph_complex_tan(igraph_complex_t z) {
    igraph_real_t zr = IGRAPH_REAL(z);
    igraph_real_t zi = IGRAPH_IMAG(z);
    igraph_complex_t res;
    if (fabs (zi) < 1) {
        igraph_real_t D = pow (cos (zr), 2.0) + pow (sinh (zi), 2.0);
        IGRAPH_REAL(res) = 0.5 * sin (2 * zr) / D;
        IGRAPH_IMAG(res) = 0.5 * sinh (2 * zi) / D;
    } else {
        igraph_real_t u = exp (-zi);
        igraph_real_t C = 2 * u / (1 - pow (u, 2.0));
        igraph_real_t D = 1 + pow (cos (zr), 2.0) * pow (C, 2.0);
        igraph_real_t S = pow (C, 2.0);
        igraph_real_t T = 1.0 / tanh (zi);
        IGRAPH_REAL(res) = 0.5 * sin (2 * zr) * S / D;
        IGRAPH_IMAG(res) = T / D;
    }
    return res;
}

igraph_complex_t igraph_complex_sec(igraph_complex_t z) {
    return igraph_complex_inv(igraph_complex_cos(z));
}

igraph_complex_t igraph_complex_csc(igraph_complex_t z) {
    return igraph_complex_inv(igraph_complex_sin(z));
}

igraph_complex_t igraph_complex_cot(igraph_complex_t z) {
    return igraph_complex_inv(igraph_complex_tan(z));
}