.. _program_listing_file_src_utilities.hpp:

Program Listing for File utilities.hpp
======================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_utilities.hpp>` (``src/utilities.hpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   /* This file is part of brille.
   
   Copyright © 2019,2020 Greg Tucker <greg.tucker@stfc.ac.uk>
   
   brille is free software: you can redistribute it and/or modify it under the
   terms of the GNU Affero General Public License as published by the Free
   Software Foundation, either version 3 of the License, or (at your option)
   any later version.
   
   brille 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 Affero General Public License for more details.
   
   You should have received a copy of the GNU Affero General Public License
   along with brille. If not, see <https://www.gnu.org/licenses/>.            */
   
   #ifndef BRILLE_UTILITIES_HPP_
   #define BRILLE_UTILITIES_HPP_
   #include <cmath>
   #include <type_traits>
   #include <limits>
   #include <stdexcept>  // for std::overflow_error
   #include <string>
   #include <math.h>
   #include <numeric>
   #include <iostream>
   #include "debug.hpp"
   #include "approx.hpp"
   
   
   namespace brille{
     const double             pi = 3.14159265358979323846;
     const double         halfpi = 1.57079632679489661923;
     // const double      quarterpi = 0.785398163397448309616;
     // const double      inversepi = 0.318309886183790671538;
     // const double      twooverpi = 0.636619772367581343076;
     // const double  twooversqrtpi = 1.12837916709551257390;
     // const double        sqrttwo = 1.41421356237309504880;
     // const double sqrttwoovertwo = 0.707106781186547524401;
     // const double              e = 2.71828182845904523536;
     // const double          log2e = 1.44269504088896340736;
     // const double         log10e = 0.434294481903251827651;
     // const double            ln2 = 0.693147180559945309417;
     // const double           ln10 = 2.30258509299404568402;
     namespace utils{
   
   template<typename T, int N=3> T trace(const T *M);
   
   template<typename T, int N, int M> void copy_array(T *D, const T *S);
   template<typename T, int N=3> void copy_matrix(T *M, const T *A);
   template<typename T, int N=3> void copy_vector(T *V, const T *A);
   
   template<typename T, int N, int M> bool equal_array(const T *A, const T *B, const T tol=0);
   template<typename T, int N=3> bool equal_matrix(const T *A, const T *B, const T tol=0);
   template<typename T, int N=3> bool equal_vector(const T *A, const T *B, const T tol=0);
   
   template<typename T, typename R, typename S, size_t N, size_t I, size_t M> void multiply_arrays(T *C, const R *A, const S *B);
   template<typename T, typename R, typename S, size_t N=3> void multiply_matrix_matrix(T *C, const R *A, const S *B);
   template<typename T, typename R, typename S, size_t N=3> void multiply_matrix_vector(T *C, const R *A, const S *b);
   template<typename T, typename R, typename S, size_t N=3> void multiply_vector_matrix(T *C, const R *a, const S *B);
   
   // array multiplication with specializations for complex array type(s)
   template<class T, class R, class S, class I> void mul_arrays(T* C, const I n, const I l, const I m, const R* A, const S* B);
   template<class T, class R, class S, class I> void mul_arrays(std::complex<T>* C, const I n, const I l, const I m, const R* A, const std::complex<S>* B);
   template<class T, class R, class S, class I> void mul_arrays(std::complex<T>* C, const I n, const I l, const I m, const std::complex<R>* A, const S* B);
   template<class T, class R, class S, class I> void mul_arrays(std::complex<T>* C, const I n, const I l, const I m, const std::complex<R>* A, const std::complex<S>* B);
   // specializations for matrix*matrix, matrix*vector, and vector*matrix
   template<class T, class R, class S, class I> void mul_mat_mat(T* C, const I n, const R* A, const S* B){ mul_arrays(C, n, n, n, A, B);}
   template<class T, class R, class S, class I> void mul_mat_vec(T* C, const I n, const R* A, const S* B){ mul_arrays(C, n, n, I(1), A, B);}
   template<class T, class R, class S, class I> void mul_vec_mat(T* C, const I n, const R* A, const S* B){ mul_arrays(C, I(1), n, n, A, B);}
   
   // template<class T,class R,class I,class S=std::common_type_t<T,R>>
   // std::enable_if_t<std::is_unsigned_v<I>&&std::is_same_v<S,R>>
   // mul_mat_vec_inplace(const I n, const T* A, R* B){
   //   S* C = new S[n]();
   //   mul_arrays<S,T,R,I>(C, n, n, I(1), A, B);
   //   for (I i=0; i<n; ++i) B[n] = C[n];
   //   delete[] C;
   // }
   template<class T,class R,class I,class S=std::common_type_t<T,R>>
   std::enable_if_t<std::is_unsigned_v<I>>
   mul_mat_vec_inplace(const I n, const T* A, R* B){
     S* C = new S[n]();
     mul_arrays(C, n, n, I(1), A, B);
     //for (I i=0; i<n; ++i) B[n] = C[n];
     std::copy(C,C+n,B);
     delete[] C;
   }
   
   template<typename T, typename R, typename S, int N, int M> void add_arrays(T *C, const R *A, const S *B);
   template<typename T, typename R, typename S, int N=3> void add_matrix(T *C, const R *A, const S *B);
   template<typename T, typename R, typename S, int N, int M> void subtract_arrays(T *C, const R *A, const S *B);
   template<typename T, typename R, typename S, int N> void subtract_matrix(T *C, const R *A, const S *B);
   
   template<typename T,typename R> T my_cast(const R a);
   template<typename T, typename R, int N, int M> void cast_array(T *A, const R *B);
   template<typename T, typename R, int N> void cast_matrix(T *A, const R *B);
   template<typename T, typename R, int N> void cast_vector(T *A, const R *B);
   
   template<typename R> void array_cofactor(R *C, const R *A, const int i, const int j, const int N=3, const int M=3);
   template<typename R> void matrix_cofactor(R *C, const R *A, const int i, const int j, const int N=3);
   
   template<typename R> R matrix_determinant(const R *M, const int N=3);
   template<typename R> void matrix_adjoint(R *A, const R *M, const int N=3);
   template<typename R> R matrix_determinant_and_inverse(R *invM, const R *M, const R tol=0, const int N=3);
   template<typename R> bool matrix_inverse(R *invM, const R *M, R tol=0, const int N=3);
   
   template<typename T, int N=3> bool similar_matrix(T *M, const T *A, const T *B, const T tol=0);
   
   template<typename R, int N, int M> void array_transpose(R *D, const R *S);
   template<typename R, int N=3> void matrix_transpose(R *D, const R *S);
   template<typename R, int N=3> void matrix_transpose(R *B);
   
   template<typename R, int N=3> void matrix_metric(R *M, const R *L);
   
   template<typename R, int N=3> R vector_norm_squared(const R * v);
   
   template<typename T, typename R, typename S, int N> class vector_cross_impl {
   public:
       static void vector_cross(T*, const R*, const S*) {
           throw std::runtime_error("The cross product is only defined for 3-vectors");
       }
   };
   template<typename T, typename R, typename S> class vector_cross_impl<T,R,S,3>{
   public:
       static void vector_cross(T * c, const R * a, const S * b) {
           c[0] = static_cast<T>(a[1])* static_cast<T>(b[2]) - static_cast<T>(a[2])* static_cast<T>(b[1]);
           c[1] = static_cast<T>(a[2])* static_cast<T>(b[0]) - static_cast<T>(a[0])* static_cast<T>(b[2]);
           c[2] = static_cast<T>(a[0])* static_cast<T>(b[1]) - static_cast<T>(a[1])* static_cast<T>(b[0]);
       }
   };
   template<typename T, typename R, typename S, int N = 3>
   void vector_cross(T* c, const R* a, const S* b) {
       vector_cross_impl<T, R, S, N>::vector_cross(c, a, b);
   }
   template<typename S,typename T, typename R> S vector_dot(const size_t, const T* a, const R* b);
   template<typename R, int N=3> R vector_dot(const R *a, const R *b);
   
   template<typename T> T mod1(const T a);
   
   template<typename T, typename R, int N=3> bool is_int_matrix(const T * A, const R tol);
   template<typename R> bool is_int_matrix(const int *, const R);
   
   template<class I, class T> T frobenius_distance(const I n, const T* A, const T* B);
   template<class I, class T> T frobenius_distance(const I n, const std::complex<T>* A, const std::complex<T>* B);
   
   template<class I, class T> T vector_angle(const I n, const T* A, const T* B);
   template<class I, class T> T vector_angle(const I n, const std::complex<T>* A, const std::complex<T>* B);
   template<class I, class T> T euclidean_angle(const I n, const std::complex<T>* A, const std::complex<T>* B);
   template<class I, class T> T hermitian_angle(const I n, const std::complex<T>* A, const std::complex<T>* B);
   template<class I, class T> T hermitian_angle(const I n, const T* A, const T* B);
   template<class I, class T> T hermitian_product(const I n, const T* a, const T* b);
   template<class I, class T> std::complex<T> hermitian_product(const I n, const T* a, const std::complex<T>* b);
   template<class I, class T> std::complex<T> hermitian_product(const I n, const std::complex<T>* a, const T* b);
   template<class I, class T> std::complex<T> hermitian_product(const I n, const std::complex<T>* a, const std::complex<T>* b);
   
   template<class I, class T> T vector_distance(const I n, const T* a, const T* b);
   template<class I, class T> T vector_distance(const I n, const std::complex<T>* a, const std::complex<T>* b);
   template<class I, class T> T vector_product(const I n, const T* a, const T* b);
   template<class I, class T> T vector_product(const I n, const T* a, const std::complex<T>* b);
   template<class I, class T> T vector_product(const I n, const std::complex<T>* a, const T* b);
   template<class I, class T> T vector_product(const I n, const std::complex<T>* a, const std::complex<T>* b);
   
   template<class I, class T> T inner_product(const I n, const T* a, const T* b);
   template<class I, class T> T inner_product(const I n, const std::complex<T>* a, const std::complex<T>* b);
   
   template<class T> T squared_distance(const T&A, const T& B);
   template<class T> T squared_distance(const std::complex<T>& A, const std::complex<T>& B);
   template<class T> T magnitude(const T a);
   template<class T> T magnitude(const std::complex<T> a);
   
   template<class T> size_t encode_array_signs(const size_t n, const T* a);
   template<class T> size_t encode_array_signs(const size_t n, const std::complex<T>* a);
   
   template<class T> int make_eigenvectors_equivalent(const size_t n, const T* v0, T* v1);
   template<class T> int make_eigenvectors_equivalent(const size_t n, const std::complex<T>* v0, std::complex<T> v1);
   
   template<class T, class R>
   enable_if_t<std::is_unsigned<T>::value&&std::is_unsigned<R>::value, unsigned long long>
   binomial_coefficient(const T n, const R k);
   
   template<typename S,typename U> S u2s(const U u);
   template<typename U,typename S> U s2u(const S s);
   
   #include "utilities.tpp"
   
     } // namespace utils
   } // namespace brille
   
   
   #endif