src/ComplexMath.cpp

Go to the documentation of this file.

// Copyright (c) 2010 Matthew Krupcale

// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:

// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.

// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.

#ifndef COMPLEX_MATH_CPP
#define COMPLEX_MATH_CPP

//============================================================================
// Name        : ComplexMath.cpp
// Author      : Matthew Krupcale
// Version     : 0.1.2
// Copyright   : Copyright (C) 2010 Matthew Krupcale
// Description : A set of mathematical functions for complex numbers
//============================================================================

#include "ComplexNumber.h"
#include "ComplexMath.h"
#include "MathExceptions.h"
#include "Namespace.h"
#include <cmath>
#include <vector>

BEGIN_NAMESPACE(math)


BEGIN_NAMESPACE(complex)

template <class T>
T
abs(const ComplexNumber<T>& z) {
    return std::sqrt(std::pow(z.getReal(), 2) +
                     std::pow(z.getImaginary(), 2));
}

template <class T>
T
abs2(const ComplexNumber<T>& z) {
    return std::pow(z.getReal(), 2) +
           std::pow(z.getImaginary(), 2);
}

template <class T>
T
arg(const ComplexNumber<T>& z) {
    if (z.getReal() == 0 && z.getImaginary() == 0) {
        return 0;
    }
    return std::atan2(z.getImaginary(), z.getReal());
}

template <class T>
T
logabs(const ComplexNumber<T>& z) {
    if (z == 0) {
        throw LogarithmOfZeroException();
    }
    return 0.5 * std::log(abs2(z));
}

template <class T>
ComplexNumber<T>
conjugate(const ComplexNumber<T>& z) {
    return ComplexNumber<T>(z.getReal(), -z.getImaginary());
}

template <class T>
ComplexNumber<T>
exp(const ComplexNumber<T>& z) {
    return ComplexNumber<T>(std::exp(z.getReal()) * std::cos(z.getImaginary()),
                            std::exp(z.getReal()) * std::sin(z.getImaginary()));
}

template <class T>
ComplexNumber<T>
inverse(const ComplexNumber<T>& z) {
    // prevent division by zero
    if (z == 0) {
        throw DivisionByZeroException();
    }
    return (conjugate(z) / abs2(z));
}

template <class T>
ComplexNumber<T>
log(const ComplexNumber<T>& z) {
    if (z == 0) {
        throw LogarithmOfZeroException();
    }
    return ComplexNumber<T>(logabs(z), arg(z));
}

template <class T>
ComplexNumber<T>
log(const T& x) {
    if (x == 0) {
        throw LogarithmOfZeroException();
    }
    return ComplexNumber<T>(std::log(std::abs(x)), std::atan2(0, x));
}

template <class T>
ComplexNumber<T>
log10(const ComplexNumber<T>& z) {
    return (log(z) / M_LN10l);
}

template <class T>
ComplexNumber<T>
log10(const T& x) {
    if (x == 0) {
        throw LogarithmOfZeroException();
    }
    return ComplexNumber<T>(std::log(std::abs(x)) / M_LN10l,
                            std::atan2(0, x) / M_LN10l);
}

template <class T>
ComplexNumber<T>
logb(const ComplexNumber<T>& z, const ComplexNumber<T>& b) {
    return (log(z) / log(b));
}

template <class T>
ComplexNumber<T>
logb(const ComplexNumber<T>& z, const T& b) {
    return (log(z) / log(b));
}

template <class T>
ComplexNumber<T>
logb(const T& x, const ComplexNumber<T>& b) {
    return (log(x) / log(b));
}

template <class T>
ComplexNumber<T>
logb(const T& x, const T& b) {
    return (log(x) / log(b));
}

template <class T>
ComplexNumber<T>
pow(const ComplexNumber<T>& a, const ComplexNumber<T>& b) {
    if (a == 0) {
        if (b == 0) {
            throw ZeroToTheZerothPowerException();
        } else {
            throw ZeroToComplexPowerException();
        }
    }
    return exp(b * log(a));
}

template <class T>
ComplexNumber<T>
pow(const ComplexNumber<T>& a, const T& b) {
    if (a == 0 && b ==0) {
        throw ZeroToTheZerothPowerException();
    }
    return exp(b * log(a));
}

template <class T>
ComplexNumber<T>
pow(const T& a, const ComplexNumber<T>& b) {
    if (a == 0) {
        if (b == 0) {
            throw ZeroToTheZerothPowerException();
        } else {
            throw ZeroToComplexPowerException();
        }
    }
    return exp(b * log(a));
}

template <class T>
ComplexNumber<T>
sqrt(const ComplexNumber<T>& z) {
    return ComplexNumber<T>(std::pow(std::pow(z.getReal(), 2) +
                                     std::pow(z.getImaginary(), 2), 0.25) *
                            std::cos(0.5 * arg(z)),
                            std::pow(std::pow(z.getReal(), 2) +
                                     std::pow(z.getImaginary(), 2), 0.25) *
                            std::sin(0.5 * arg(z)));
}

template <class T>
ComplexNumber<T>
sqrt(const T& x) {
    if (x >= 0) {
        return ComplexNumber<T>(std::sqrt(x), 0);
    } else {
        return ComplexNumber<T>(0, std::sqrt(-x));
    }
}

template <class T>
ComplexNumber<T>
cbrt(const ComplexNumber<T>& z) {
    return ComplexNumber<T>(std::pow(std::pow(z.getReal(), 2) +
                                     std::pow(z.getImaginary(), 2), 1. / 6) *
                            std::cos(arg(z) / 3),
                            std::pow(std::pow(z.getReal(), 2) +
                                     std::pow(z.getImaginary(), 2), 1. / 6) *
                            std::sin(arg(z) / 3));
}

template <class T>
ComplexNumber<T>
nthRoot(const ComplexNumber<T>& z, const int n) {
    if (n == 0) {
        throw DivisionByZeroException();
    }
    return ComplexNumber<T>(std::pow(std::pow(z.getReal(), 2) +
                                     std::pow(z.getImaginary(), 2), 1. / (2 * n)) *
                            std::cos(arg(z) / n),
                            std::pow(std::pow(z.getReal(), 2) +
                                     std::pow(z.getImaginary(), 2), 1. / (2 * n)) *
                            std::sin(arg(z) / n));
}

template <class T>
std::vector<ComplexNumber<T> >
nthRoots(const ComplexNumber<T>& z, const int n) {
    std::vector<ComplexNumber<T> > roots(n);
    nthRoots(z, roots, n);
    return roots;
}

template <class T>
void
nthRoots(const ComplexNumber<T>& z,
         std::vector<ComplexNumber<T> >& roots, const int n) {
    if (n == 0) {
        throw DivisionByZeroException();
    }
    int absN = (n < 0) ? -n : n;
    for (int i = 0; i < absN; i++) {
        roots.push_back(ComplexNumber<T>(std::pow(std::pow(z.getReal(), 2) +
                                                  std::pow(z.getImaginary(), 2), 1. / (2 * n)) *
                                         std::cos((arg(z) + 2 * i * M_PIl) / n),
                                         std::pow(std::pow(z.getReal(), 2) +
                                                  std::pow(z.getImaginary(), 2), 1. / (2 * n)) *
                                         std::sin((arg(z) + 2 * i * M_PIl) / n)));
    }
}

template <class T>
ComplexNumber<T>
sin(const ComplexNumber<T>& z) {
    return ComplexNumber<T>(std::sin(z.getReal()) * std::cosh(z.getImaginary()),
                            std::cos(z.getReal()) * std::sinh(z.getImaginary()));
}

template <class T>
ComplexNumber<T>
cos(const ComplexNumber<T>& z) {
    return ComplexNumber<T>(std::cos(z.getReal()) * std::cosh(z.getImaginary()),
                            -std::sin(z.getReal()) * std::sinh(z.getImaginary()));
}

template <class T>
ComplexNumber<T>
tan(const ComplexNumber<T>& z) {
    return ComplexNumber<T>(std::sin(2 * z.getReal()) /
                            (std::cosh(2 * z.getImaginary()) +
                             std::cos(2 * z.getReal())),
                            std::sinh(2 * z.getImaginary()) /
                            (std::cosh(2 * z.getImaginary()) +
                             std::cos(2 * z.getReal())));
}

template <class T>
ComplexNumber<T>
csc(const ComplexNumber<T>& z) {
    return inverse(sin(z));
}

template <class T>
ComplexNumber<T>
sec(const ComplexNumber<T>& z) {
    return inverse(cos(z));
}

template <class T>
ComplexNumber<T>
cot(const ComplexNumber<T>& z) {
    return inverse(cot(z));
}

template <class T>
ComplexNumber<T>
asin(const ComplexNumber<T>& z) {
    if (z.getImaginary() == 0) {
        return asin(z.getReal());
    } else {
        return -ComplexNumber<double>::i * log(ComplexNumber<double>::i * z +
                                               sqrt(1 - pow(z, 2)));
    }
}

template <class T>
ComplexNumber<T>
asin(const T& x) {
    if (std::abs(x) <= 1) {
        return ComplexNumber<T>(::asin(x), 0);
    } else {
        if (x < 0) {
            return ComplexNumber<T>(-M_PI_2l, ::acosh(-x));
        } else {
            return ComplexNumber<T>(M_PI_2l, -::acosh(x));
        }
    }
}

template <class T>
ComplexNumber<T>
acos(const ComplexNumber<T>& z) {
    if (z.getImaginary() == 0) {
        return acos(z.getReal());
    } else {
        return -ComplexNumber<double>::i * log(z + sqrt(pow(z, 2) - 1));
    }
}

template <class T>
ComplexNumber<T>
acos(const T& x) {
    if (std::abs(x) <= 1) {
        return ComplexNumber<T>(::acos(x), 0);
    } else {
        if (x < 0) {
            return ComplexNumber<T>(M_PIl, -::acosh(-x));
        } else {
            return ComplexNumber<T>(0, ::acosh(x));
        }
    }
}

template <class T>
ComplexNumber<T>
atan(const ComplexNumber<T>& z) {
    if (z.getImaginary() == 0) {
        return atan(z.getReal());
    } else {
        return ComplexNumber<double>::i * .5 * log((ComplexNumber<double>::i + z) /
                                                   (ComplexNumber<double>::i - z));
    }
}

template <class T>
ComplexNumber<T>
atan(const T& x) {
    return ComplexNumber<T>(std::atan(x), 0);
}

template <class T>
ComplexNumber<T>
acsc(const ComplexNumber<T>& z) {
    return asin(inverse(z));
}

template <class T>
ComplexNumber<T>
acsc(const T& x) {
    return asin(1. / x);
}

template <class T>
ComplexNumber<T>
asec(const ComplexNumber<T>& z) {
    return acos(inverse(z));
}

template <class T>
ComplexNumber<T>
asec(const T& x) {
    return acos(1. / x);
}

template <class T>
ComplexNumber<T>
acot(const ComplexNumber<T>& z) {
    return atan(inverse(z));
}

template <class T>
ComplexNumber<T>
acot(const T& x) {
    return atan(1. / x);
}

template <class T>
ComplexNumber<T>
sinh(const ComplexNumber<T>& z) {
    return ComplexNumber<T>(std::sinh(z.getReal()) * std::cos(z.getImaginary()),
                            std::cosh(z.getReal()) * std::sin(z.getImaginary()));
}

template <class T>
ComplexNumber<T>
cosh(const ComplexNumber<T>& z) {
    return ComplexNumber<T>(std::cosh(z.getReal()) * std::cos(z.getImaginary()),
                            std::sinh(z.getReal()) * std::sin(z.getImaginary()));
}

template <class T>
ComplexNumber<T>
tanh(const ComplexNumber<T>& z) {
    return ComplexNumber<T>(std::sinh(2 * z.getReal()) /
                            (std::cos(2 * z.getImaginary()) +
                             std::cosh(2 * z.getReal())),
                            std::sin(2 * z.getImaginary()) /
                            (std::cos(2 * z.getImaginary()) +
                             std::cosh(2 * z.getReal())));
}

template <class T>
ComplexNumber<T>
csch(const ComplexNumber<T>& z) {
    return inverse(sinh(z));
}

template <class T>
ComplexNumber<T>
sech(const ComplexNumber<T>& z) {
    return inverse(cosh(z));
}

template <class T>
ComplexNumber<T>
coth(const ComplexNumber<T>& z) {
    return inverse(tanh(z));
}

template <class T>
ComplexNumber<T>
asinh(const ComplexNumber<T>& z) {
    if (z.getImaginary() == 0) {
        return asinh(z.getReal());
    } else {
        return log(z + sqrt(pow(z, 2) + 1));
    }
}

template <class T>
ComplexNumber<T>
asinh(const T& x) {
    return ComplexNumber<T>(::asinh(x), 0);
}

template <class T>
ComplexNumber<T>
acosh(const ComplexNumber<T>& z) {
    if (z.getImaginary() == 0) {
        return acosh(z.getReal());
    } else {
        return log(z + sqrt(z - 1) * sqrt(z + 1));
    }
}

template <class T>
ComplexNumber<T>
acosh(const T& x) {
    if (x >= 1) {
        return ComplexNumber<T>(::acosh(x), 0);
    } else {
        if (x >= -1) {
            return ComplexNumber<T>(0, ::acos(x));
        } else {
            return ComplexNumber<T>(::acosh(-x), M_PIl);
        }
    }
}

template <class T>
ComplexNumber<T>
atanh(const ComplexNumber<T>& z) {
    if (z.getImaginary() == 0) {
        return atanh(z.getReal());
    } else {
        return 0.5 * log((1 + z) / (1 - z));
    }
}

template <class T>
ComplexNumber<T>
atanh(const T& x) {
    if (std::abs(x) < 1) {
        return ComplexNumber<T>(::atanh(x), 0);
    } else {
        return ComplexNumber<T>(::atanh(1. / x), (x < 0) ? M_PI_2l : -M_PI_2l);
    }
}

template <class T>
ComplexNumber<T>
acsch(const ComplexNumber<T>& z) {
    return asinh(inverse(z));
}

template <class T>
ComplexNumber<T>
acsch(const T& x) {
    return ComplexNumber<T>(::asinh(1. / x), 0);
}

template <class T>
ComplexNumber<T>
asech(const ComplexNumber<T>& z) {
    return acosh(inverse(z));
}

template <class T>
ComplexNumber<T>
asech(const T& x) {
    if (x > 0 && x <= 1) {
        return ComplexNumber<T>(::acosh(1. / x), 0);
    } else {
        return acosh(1. / x);
    }
}

template <class T>
ComplexNumber<T>
acoth(const ComplexNumber<T>& z) {
    return atanh(inverse(z));
}

template <class T>
ComplexNumber<T>
acoth(const T& x) {
    if (std::abs(x) > 1) {
        return ComplexNumber<T>(::atanh(1. / x), 0);
    } else {
        return atanh(1. / x);
    }
}

END_NAMESPACE(complex)
END_NAMESPACE(math)

#endif /* COMPLEX_MATH_CPP */