|
Scythe-1.0.3
|
00001 /* 00002 * Scythe Statistical Library Copyright (C) 2000-2002 Andrew D. Martin 00003 * and Kevin M. Quinn; 2002-present Andrew D. Martin, Kevin M. Quinn, 00004 * and Daniel Pemstein. All Rights Reserved. 00005 * 00006 * This program is free software; you can redistribute it and/or 00007 * modify under the terms of the GNU General Public License as 00008 * published by Free Software Foundation; either version 2 of the 00009 * License, or (at your option) any later version. See the text files 00010 * COPYING and LICENSE, distributed with this source code, for further 00011 * information. 00012 * -------------------------------------------------------------------- 00013 * scythestat/smath.h 00014 * 00015 */ 00016 00032 #ifndef SCYTHE_MATH_H 00033 #define SCYTHE_MATH_H 00034 00035 #ifdef SCYTHE_COMPILE_DIRECT 00036 #include "matrix.h" 00037 #include "algorithm.h" 00038 #include "error.h" 00039 #else 00040 #include "scythestat/matrix.h" 00041 #include "scythestat/algorithm.h" 00042 #include "scythestat/error.h" 00043 #endif 00044 00045 #include <cmath> 00046 #include <numeric> 00047 #include <set> 00048 00049 namespace scythe { 00050 00051 namespace { 00052 typedef unsigned int uint; 00053 } 00054 00055 /* Almost every function in this file follows one of the two patterns 00056 * described by these macros. The first macro handles single-argument 00057 * functions. The second handles two-matrix-argument functions (or 00058 * scalar-matrix, matrix-scalar. The second macro also permits 00059 * cross-type operations (these are limited only by the capabilities 00060 * of the underlying functions). 00061 */ 00062 #define SCYTHE_MATH_OP(NAME, OP) \ 00063 template <matrix_order RO, matrix_style RS, typename T, \ 00064 matrix_order PO, matrix_style PS> \ 00065 Matrix<T,RO,RS> \ 00066 NAME (const Matrix<T,PO,PS>& A) \ 00067 { \ 00068 Matrix<T,RO,RS> res(A.rows(), A.cols(), false); \ 00069 std::transform(A.begin_f(), A.end_f(), res.begin_f(), OP); \ 00070 return res; \ 00071 } \ 00072 \ 00073 template <typename T, matrix_order O, matrix_style S> \ 00074 Matrix<T,O,Concrete> \ 00075 NAME (const Matrix<T,O,S>& A) \ 00076 { \ 00077 return NAME<O,Concrete>(A); \ 00078 } 00079 00080 #define SCYTHE_MATH_OP_2ARG(NAME, OP) \ 00081 template <matrix_order RO, matrix_style RS, typename T, \ 00082 matrix_order PO1, matrix_style PS1, \ 00083 matrix_order PO2, matrix_style PS2, typename S> \ 00084 Matrix<T,RO,RS> \ 00085 NAME (const Matrix<T,PO1,PS1>& A, const Matrix<S,PO2,PS2>& B) \ 00086 { \ 00087 SCYTHE_CHECK_10 (A.size() != 1 && B.size() != 1 && \ 00088 A.size() != B.size(), scythe_conformation_error, \ 00089 "Matrices with dimensions (" << A.rows() \ 00090 << ", " << A.cols() \ 00091 << ") and (" << B.rows() << ", " << B.cols() \ 00092 << ") are not conformable"); \ 00093 \ 00094 Matrix<T,RO,RS> res; \ 00095 \ 00096 if (A.size() == 1) { \ 00097 res.resize2Match(B); \ 00098 std::transform(B.template begin_f<RO>(), B.template end_f<RO>(),\ 00099 res.begin_f(), std::bind1st(OP, A(0))); \ 00100 } else if (B.size() == 1) { \ 00101 res.resize2Match(A); \ 00102 std::transform(A.template begin_f<RO>(), A.template end_f<RO>(),\ 00103 res.begin_f(), std::bind2nd(OP, B(0))); \ 00104 } else { \ 00105 res.resize2Match(A); \ 00106 std::transform(A.template begin_f<RO>(), A.template end_f<RO>(),\ 00107 B.template begin_f<RO>(), res.begin_f(), OP); \ 00108 } \ 00109 \ 00110 return res; \ 00111 } \ 00112 \ 00113 template <typename T, matrix_order PO1, matrix_style PS1, \ 00114 matrix_order PO2, matrix_style PS2, \ 00115 typename S> \ 00116 Matrix<T,PO1,Concrete> \ 00117 NAME (const Matrix<T,PO1,PS1>& A, const Matrix<S,PO2,PS2>& B) \ 00118 { \ 00119 return NAME<PO1,Concrete>(A, B); \ 00120 } \ 00121 \ 00122 template<matrix_order RO, matrix_style RS, typename T, \ 00123 matrix_order PO, matrix_style PS, typename S> \ 00124 Matrix<T,RO,RS> \ 00125 NAME (const Matrix<T,PO,PS>& A, S b) \ 00126 { \ 00127 return NAME<RO,RS>(A, Matrix<S,RO,Concrete>(b)); \ 00128 } \ 00129 \ 00130 template <typename T, typename S, matrix_order PO, matrix_style PS> \ 00131 Matrix<T,PO,Concrete> \ 00132 NAME (const Matrix<T,PO,PS>& A, S b) \ 00133 { \ 00134 return NAME<PO,Concrete>(A, Matrix<S,PO,Concrete>(b)); \ 00135 } \ 00136 \ 00137 template<matrix_order RO, matrix_style RS, typename T, \ 00138 matrix_order PO, matrix_style PS, typename S> \ 00139 Matrix<T,RO,RS> \ 00140 NAME (T a, const Matrix<S,PO,PS>& B) \ 00141 { \ 00142 return NAME<RO,RS>(Matrix<S, RO,Concrete>(a), B); \ 00143 } \ 00144 \ 00145 template <typename T, typename S, matrix_order PO, matrix_style PS> \ 00146 Matrix<T,PO,Concrete> \ 00147 NAME (T a, const Matrix<S,PO,PS>& B) \ 00148 { \ 00149 return NAME<PO,Concrete>(Matrix<S,PO,Concrete>(a), B); \ 00150 } 00151 00152 00153 /* calc the inverse cosine of each element of a Matrix */ 00154 00176 SCYTHE_MATH_OP(acos, ::acos) 00177 00178 /* calc the inverse hyperbolic cosine of each element of a Matrix */ 00201 SCYTHE_MATH_OP(acosh, ::acosh) 00202 00203 /* calc the inverse sine of each element of a Matrix */ 00204 00227 SCYTHE_MATH_OP(asin, ::asin) 00228 00229 /* calc the inverse hyperbolic sine of each element of a Matrix */ 00230 00253 SCYTHE_MATH_OP(asinh, ::asinh) 00254 00255 /* calc the inverse tangent of each element of a Matrix */ 00256 00279 SCYTHE_MATH_OP(atan, ::atan) 00280 00281 /* calc the inverse hyperbolic tangent of each element of a Matrix */ 00304 SCYTHE_MATH_OP(atanh, ::atanh) 00305 00306 /* calc the angle whose tangent is y/x */ 00307 00332 SCYTHE_MATH_OP_2ARG(atan2, std::ptr_fun(::atan2)) 00333 00334 /* calc the cube root of each element of a Matrix */ 00346 SCYTHE_MATH_OP(cbrt, ::cbrt) 00347 00348 /* calc the ceil of each element of a Matrix */ 00360 SCYTHE_MATH_OP(ceil, ::ceil) 00361 00362 /* create a matrix containing the absval of the first input and the 00363 * sign of the second 00364 */ 00376 SCYTHE_MATH_OP_2ARG(copysign, std::ptr_fun(::copysign)) 00377 00378 /* calc the cosine of each element of a Matrix */ 00379 00401 SCYTHE_MATH_OP(cos, ::cos) 00402 00403 /* calc the hyperbolic cosine of each element of a Matrix */ 00425 SCYTHE_MATH_OP(cosh, ::cosh) 00426 00427 /* calc the error function of each element of a Matrix */ 00438 SCYTHE_MATH_OP(erf, ::erf) 00439 00440 /* calc the complementary error function of each element of a Matrix */ 00452 SCYTHE_MATH_OP(erfc, ::erfc) 00453 00454 /* calc the vaue e^x of each element of a Matrix */ 00466 SCYTHE_MATH_OP(exp, ::exp) 00467 00468 /* calc the exponent - 1 of each element of a Matrix */ 00480 SCYTHE_MATH_OP(expm1, ::expm1) 00481 00482 /* calc the absval of each element of a Matrix */ 00491 SCYTHE_MATH_OP(fabs, ::fabs) 00492 00493 /* calc the floor of each element of a Matrix */ 00505 SCYTHE_MATH_OP(floor, ::floor) 00506 00507 /* calc the remainder of the division of each matrix element */ 00518 SCYTHE_MATH_OP_2ARG(fmod, std::ptr_fun(::fmod)) 00519 00520 /* calc the fractional val of input and return exponents in int 00521 * matrix reference 00522 */ 00523 00526 template <matrix_order RO, matrix_style RS, typename T, 00527 matrix_order PO1, matrix_style PS1, 00528 matrix_order PO2, matrix_style PS2> 00529 Matrix<T,RO,RS> 00530 frexp (const Matrix<T,PO1,PS1>& A, Matrix<int,PO2,PS2>& ex) 00531 { 00532 SCYTHE_CHECK_10(A.size() != ex.size(), scythe_conformation_error, 00533 "The input matrix sizes do not match"); 00534 Matrix<T,PO1,Concrete> res(A.rows(), A.cols()); 00535 00536 typename Matrix<T,PO1,PS1>::const_forward_iterator it; 00537 typename Matrix<T,PO1,Concrete>::forward_iterator rit 00538 = res.begin_f(); 00539 typename Matrix<int,PO2,PS2>::const_forward_iterator it2 00540 = ex.begin_f(); 00541 for (it = A.begin_f(); it != A.end_f(); ++it) { 00542 *rit = ::frexp(*it, &(*it2)); 00543 ++it2; ++rit; 00544 } 00545 00546 return res; 00547 } 00548 00549 template <typename T, matrix_order PO1, matrix_style PS1, 00550 matrix_order PO2, matrix_style PS2> 00551 Matrix<T,PO1,Concrete> 00552 frexp (Matrix<T,PO1,PS1>& A, Matrix<int,PO2,PS2>& ex) 00553 { 00554 return frexp<PO1,Concrete>(A,ex); 00555 } 00556 00557 /* calc the euclidean distance between the two inputs */ 00568 SCYTHE_MATH_OP_2ARG(hypot, std::ptr_fun(::hypot)) 00569 00570 /* return (int) logb */ 00571 SCYTHE_MATH_OP(ilogb, ::ilogb) 00572 00573 /* compute the bessel func of the first kind of the order 0 */ 00589 SCYTHE_MATH_OP(j0, ::j0) 00590 00591 /* compute the bessel func of the first kind of the order 1 */ 00607 SCYTHE_MATH_OP(j1, ::j1) 00608 00609 /* compute the bessel func of the first kind of the order n 00610 * TODO: This definition causes the compiler to issue some warnings. 00611 * Fix 00612 */ 00629 SCYTHE_MATH_OP_2ARG(jn, std::ptr_fun(::jn)) 00630 00631 /* calc x * 2 ^ex */ 00641 SCYTHE_MATH_OP_2ARG(ldexp, std::ptr_fun(::ldexp)) 00642 00643 /* compute the natural log of the absval of gamma function */ 00644 00656 SCYTHE_MATH_OP(lgamma, ::lgamma) 00657 00658 /* calc the natural log of each element of a Matrix */ 00671 SCYTHE_MATH_OP(log, ::log) 00672 00673 /* calc the base-10 log of each element of a Matrix */ 00686 SCYTHE_MATH_OP(log10, ::log10) 00687 00688 /* calc the natural log of 1 + each element of a Matrix */ 00701 SCYTHE_MATH_OP(log1p, ::log1p) 00702 00703 /* calc the logb of each element of a Matrix */ 00716 SCYTHE_MATH_OP(logb, ::logb) 00717 00718 /* x = frac + i, return matrix of frac and place i in 2nd matrix 00719 */ 00720 template <matrix_order RO, matrix_style RS, typename T, 00721 matrix_order PO1, matrix_style PS1, 00722 matrix_order PO2, matrix_style PS2> 00723 Matrix<T,RO,RS> 00724 modf (const Matrix<T,PO1,PS1>& A, Matrix<double,PO2,PS2>& ipart) 00725 { 00726 SCYTHE_CHECK_10(A.size() != ipart.size(), scythe_conformation_error, 00727 "The input matrix sizes do not match"); 00728 Matrix<T,PO1,Concrete> res(A.rows(), A.cols()); 00729 00730 typename Matrix<T,PO1,PS1>::const_forward_iterator it; 00731 typename Matrix<T,PO1,Concrete>::forward_iterator rit 00732 = res.begin_f(); 00733 typename Matrix<double,PO2,PS2>::const_forward_iterator it2 00734 = ipart.begin_f(); 00735 for (it = A.begin_f(); it != A.end_f(); ++it) { 00736 *rit = ::modf(*it, &(*it2)); 00737 ++it2; ++rit; 00738 } 00739 00740 return res; 00741 } 00742 00743 template <typename T, matrix_order PO1, matrix_style PS1, 00744 matrix_order PO2, matrix_style PS2> 00745 Matrix<T,PO1,Concrete> 00746 modf (Matrix<T,PO1,PS1>& A, Matrix<double,PO2,PS2>& ipart) 00747 { 00748 return modf<PO1,Concrete>(A,ipart); 00749 } 00750 00751 /* calc x^ex of each element of a Matrix */ 00752 00762 SCYTHE_MATH_OP_2ARG(pow, std::ptr_fun(::pow)) 00763 00764 /* calc rem == x - n * y */ 00765 SCYTHE_MATH_OP_2ARG(remainder, std::ptr_fun(::remainder)) 00766 00767 /* return x rounded to nearest int */ 00768 00778 SCYTHE_MATH_OP(rint, ::rint) 00779 00780 /* returns x * FLT_RADIX^ex */ 00781 SCYTHE_MATH_OP_2ARG(scalbn, std::ptr_fun(::scalbn)) 00782 00783 /* calc the sine of x */ 00784 00806 SCYTHE_MATH_OP(sin, ::sin) 00807 00808 /* calc the hyperbolic sine of x */ 00830 SCYTHE_MATH_OP(sinh, ::sinh) 00831 00832 /* calc the sqrt of x */ 00844 SCYTHE_MATH_OP(sqrt, ::sqrt) 00845 00846 /* calc the tangent of x */ 00847 00869 SCYTHE_MATH_OP(tan, ::tan) 00870 00871 /* calc the hyperbolic tangent of x */ 00893 SCYTHE_MATH_OP(tanh, ::tanh) 00894 00895 /* bessel function of the second kind of order 0*/ 00911 SCYTHE_MATH_OP(y0, ::y0) 00912 00913 /* bessel function of the second kind of order 1*/ 00929 SCYTHE_MATH_OP(y1, ::y1) 00930 00931 /* bessel function of the second kind of order n 00932 * TODO: This definition causes the compiler to issue some warnings. 00933 * Fix 00934 */ 00951 SCYTHE_MATH_OP_2ARG(yn, std::ptr_fun(::yn)) 00952 00953 } // end namespace scythe 00954 00955 #endif /* SCYTHE_MATH_H */
1.7.6.1