1 // -*- C++ -*-
2 /***************************************************************************
3 * blitz/array/funcs.h Math functions on arrays
4 *
5 * $Id$
6 *
7 * Copyright (C) 1997-2011 Todd Veldhuizen <tveldhui@acm.org>
8 *
9 * This file is a part of Blitz.
10 *
11 * Blitz is free software: you can redistribute it and/or modify
12 * it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation, either version 3
14 * of the License, or (at your option) any later version.
15 *
16 * Blitz is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU Lesser General Public License for more details.
20 *
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with Blitz. If not, see <http://www.gnu.org/licenses/>.
23 *
24 * Suggestions: blitz-devel@lists.sourceforge.net
25 * Bugs: blitz-support@lists.sourceforge.net
26 *
27 * For more information, please see the Blitz++ Home Page:
28 * https://sourceforge.net/projects/blitz/
29 *
30 ****************************************************************************/
31 #ifndef BZ_ARRAY_FUNCS_H
32 #define BZ_ARRAY_FUNCS_H
33
34 #include <blitz/funcs.h>
35 #include <blitz/array/newet-macros.h>
36 #include <blitz/array/reduce.h>
37 #include <blitz/levicivita.h>
38
39 namespace blitz {
40
41 /* This file contains definitions of the math functions for ETBase
42 types. Note that these definitions are in the blitz namespace and
43 will hide the definitions of the builtin math functions in the
44 global namespace. For this reason, the calls to the builtin math
45 functions in the functors' apply() methods must explicitly qualify
46 the namespace.
47 */
48
49
50 // unary functions
51
BZ_DECLARE_ARRAY_ET_UNARY(abs,Fn_abs)52 BZ_DECLARE_ARRAY_ET_UNARY(abs, Fn_abs)
53 BZ_DECLARE_ARRAY_ET_UNARY(acos, Fn_acos)
54 BZ_DECLARE_ARRAY_ET_UNARY(asin, Fn_asin)
55 BZ_DECLARE_ARRAY_ET_UNARY(atan, Fn_atan)
56 BZ_DECLARE_ARRAY_ET_UNARY(ceil, Fn_ceil)
57 BZ_DECLARE_ARRAY_ET_UNARY(cexp, Fn_exp)
58 BZ_DECLARE_ARRAY_ET_UNARY(cos, Fn_cos)
59 BZ_DECLARE_ARRAY_ET_UNARY(cosh, Fn_cosh)
60 BZ_DECLARE_ARRAY_ET_UNARY(csqrt, Fn_sqrt)
61 BZ_DECLARE_ARRAY_ET_UNARY(cube, Fn_cube)
62 BZ_DECLARE_ARRAY_ET_UNARY(exp, Fn_exp)
63 BZ_DECLARE_ARRAY_ET_UNARY(fabs, Fn_fabs)
64 BZ_DECLARE_ARRAY_ET_UNARY(floor, Fn_floor)
65 BZ_DECLARE_ARRAY_ET_UNARY(log, Fn_log)
66 BZ_DECLARE_ARRAY_ET_UNARY(log10, Fn_log10)
67 BZ_DECLARE_ARRAY_ET_UNARY(pow2, Fn_sqr)
68 BZ_DECLARE_ARRAY_ET_UNARY(pow3, Fn_cube)
69 BZ_DECLARE_ARRAY_ET_UNARY(pow4, Fn_pow4)
70 BZ_DECLARE_ARRAY_ET_UNARY(pow5, Fn_pow5)
71 BZ_DECLARE_ARRAY_ET_UNARY(pow6, Fn_pow6)
72 BZ_DECLARE_ARRAY_ET_UNARY(pow7, Fn_pow7)
73 BZ_DECLARE_ARRAY_ET_UNARY(pow8, Fn_pow8)
74 BZ_DECLARE_ARRAY_ET_UNARY(sin, Fn_sin)
75 BZ_DECLARE_ARRAY_ET_UNARY(sinh, Fn_sinh)
76 BZ_DECLARE_ARRAY_ET_UNARY(sqr, Fn_sqr)
77 BZ_DECLARE_ARRAY_ET_UNARY(sqrt, Fn_sqrt)
78 BZ_DECLARE_ARRAY_ET_UNARY(tan, Fn_tan)
79 BZ_DECLARE_ARRAY_ET_UNARY(tanh, Fn_tanh)
80
81 #ifdef BZ_HAVE_COMPLEX_FCNS
82 BZ_DECLARE_ARRAY_ET_UNARY(arg, Fn_arg)
83 BZ_DECLARE_ARRAY_ET_UNARY(conj, Fn_conj)
84 BZ_DECLARE_ARRAY_ET_UNARY(imag, Fn_imag)
85 BZ_DECLARE_ARRAY_ET_UNARY(norm, Fn_norm)
86 BZ_DECLARE_ARRAY_ET_UNARY(real, Fn_real)
87 #endif
88
89 #ifdef BZ_HAVE_IEEE_MATH
90 // finite and trunc omitted: blitz-bugs/archive/0189.html
91 BZ_DECLARE_ARRAY_ET_UNARY(acosh, Fn_acosh)
92 BZ_DECLARE_ARRAY_ET_UNARY(asinh, Fn_asinh)
93 BZ_DECLARE_ARRAY_ET_UNARY(atanh, Fn_atanh)
94 BZ_DECLARE_ARRAY_ET_UNARY(cbrt, Fn_cbrt)
95 BZ_DECLARE_ARRAY_ET_UNARY(erf, Fn_erf)
96 BZ_DECLARE_ARRAY_ET_UNARY(erfc, Fn_erfc)
97 BZ_DECLARE_ARRAY_ET_UNARY(expm1, Fn_expm1)
98 // BZ_DECLARE_ARRAY_ET_UNARY(finite, Fn_finite)
99 BZ_DECLARE_ARRAY_ET_UNARY(ilogb, Fn_ilogb)
100 BZ_DECLARE_ARRAY_ET_UNARY(blitz_isnan, Fn_isnan)
101 BZ_DECLARE_ARRAY_ET_UNARY(j0, Fn_j0)
102 BZ_DECLARE_ARRAY_ET_UNARY(j1, Fn_j1)
103 BZ_DECLARE_ARRAY_ET_UNARY(lgamma, Fn_lgamma)
104 BZ_DECLARE_ARRAY_ET_UNARY(logb, Fn_logb)
105 BZ_DECLARE_ARRAY_ET_UNARY(log1p, Fn_log1p)
106 BZ_DECLARE_ARRAY_ET_UNARY(rint, Fn_rint)
107 // BZ_DECLARE_ARRAY_ET_UNARY(trunc, Fn_trunc)
108 BZ_DECLARE_ARRAY_ET_UNARY(y0, Fn_y0)
109 BZ_DECLARE_ARRAY_ET_UNARY(y1, Fn_y1)
110 #endif
111
112 #ifdef BZ_HAVE_SYSTEM_V_MATH
113 BZ_DECLARE_ARRAY_ET_UNARY(_class, Fn__class)
114 BZ_DECLARE_ARRAY_ET_UNARY(itrunc, Fn_itrunc)
115 BZ_DECLARE_ARRAY_ET_UNARY(nearest, Fn_nearest)
116 BZ_DECLARE_ARRAY_ET_UNARY(rsqrt, Fn_rsqrt)
117 BZ_DECLARE_ARRAY_ET_UNARY(uitrunc, Fn_uitrunc)
118 #endif
119
120 // cast() function
121
122 template<typename T_cast, typename T1>
123 _bz_inline_et
124 _bz_ArrayExpr<_bz_ArrayExprUnaryOp<_bz_typename asExpr<T1>::T_expr,
125 Cast<_bz_typename asExpr<T1>::T_expr::T_numtype, T_cast> > >
126 cast(const ETBase<T1>& expr)
127 {
128 return _bz_ArrayExpr<_bz_ArrayExprUnaryOp<
129 _bz_typename asExpr<T1>::T_expr,
130 Cast<_bz_typename asExpr<T1>::T_expr::T_numtype,T_cast> > >
131 (expr.unwrap());
132 }
133
134 // binary functions
135
BZ_DECLARE_ARRAY_ET_BINARY(atan2,Fn_atan2)136 BZ_DECLARE_ARRAY_ET_BINARY(atan2, Fn_atan2)
137 BZ_DECLARE_ARRAY_ET_BINARY(fmod, Fn_fmod)
138 BZ_DECLARE_ARRAY_ET_BINARY(pow, Fn_pow)
139
140 #ifdef BZ_HAVE_COMPLEX_FCNS
141 BZ_DECLARE_ARRAY_ET_BINARY(polar, Fn_polar)
142 #endif
143
144 #ifdef BZ_HAVE_SYSTEM_V_MATH
145 BZ_DECLARE_ARRAY_ET_BINARY(copysign, Fn_copysign)
146 BZ_DECLARE_ARRAY_ET_BINARY(drem, Fn_drem)
147 BZ_DECLARE_ARRAY_ET_BINARY(hypot, Fn_hypot)
148 BZ_DECLARE_ARRAY_ET_BINARY(nextafter, Fn_nextafter)
149 BZ_DECLARE_ARRAY_ET_BINARY(remainder, Fn_remainder)
150 BZ_DECLARE_ARRAY_ET_BINARY(scalb, Fn_scalb)
151 BZ_DECLARE_ARRAY_ET_BINARY(unordered, Fn_unordered)
152 #endif
153
154 BZ_DECLARE_ARRAY_ET_BINARY((min), Min)
155 BZ_DECLARE_ARRAY_ET_BINARY((max), Max)
156
157 #ifdef BZ_HAVE_SYSTEM_V_MATH
158
159 #define BZ_DECLARE_ARRAY_ET_SCALAR_FUNCS(sca) \
160 \
161 BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(atan2, Fn_atan2, sca) \
162 BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(fmod, Fn_fmod, sca) \
163 BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(pow, Fn_pow, sca) \
164 BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(copysign, Fn_copysign, sca) \
165 BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(drem, Fn_drem, sca) \
166 BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(hypot, Fn_hypot, sca) \
167 BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(nextafter, Fn_nextafter, sca) \
168 BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(remainder, Fn_remainder, sca) \
169 BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(scalb, Fn_scalb, sca) \
170 BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(unordered, Fn_unordered, sca) \
171
172 #else
173
174 #define BZ_DECLARE_ARRAY_ET_SCALAR_FUNCS(sca) \
175 \
176 BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(atan2, Fn_atan2, sca) \
177 BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(fmod, Fn_fmod, sca) \
178 BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(pow, Fn_pow, sca) \
179
180 #endif
181
182 BZ_DECLARE_ARRAY_ET_SCALAR_FUNCS(int)
183 BZ_DECLARE_ARRAY_ET_SCALAR_FUNCS(float)
184 BZ_DECLARE_ARRAY_ET_SCALAR_FUNCS(double)
185 BZ_DECLARE_ARRAY_ET_SCALAR_FUNCS(long double)
186
187 #ifdef BZ_HAVE_COMPLEX_FCNS
188 BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(polar, Fn_polar, int)
189 BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(polar, Fn_polar, float)
190 BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(polar, Fn_polar, double)
191 BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(polar, Fn_polar, long double)
192
193 template<typename T1, typename T2>
194 inline _bz_ArrayExprBinaryOp<
195 typename asExpr<complex<T1> >::T_expr,
196 typename asExpr<T2>::T_expr,
197 Fn_pow<complex<T1>,typename asExpr<T2>::T_expr::T_numtype> >
198 pow(const complex<T1> d1, const ETBase<T2>& d2)
199 {
200 return _bz_ArrayExprBinaryOp<
201 typename asExpr<complex<T1> >::T_expr,
202 typename asExpr<T2>::T_expr,
203 Fn_pow<complex<T1>,typename asExpr<T2>::T_expr::T_numtype> >
204 (asExpr<complex<T1> >::getExpr(d1),
205 asExpr<T2>::getExpr(d2.unwrap()));
206 }
207
208
209 // global functions that don't fit anywhere else
210
211 // we define a generalized dot product for all classes as sum(a*b)
212 template<typename T1, typename T2>
213 inline
214 _bz_typename ReduceSum<_bz_typename blitz::BzBinaryExprResult<Multiply,T1,T2>::T_result::T_numtype
215 >::T_resulttype
dot(const ETBase<T1> & d1,const ETBase<T2> & d2)216 dot(const ETBase<T1>& d1, const ETBase<T2>& d2)
217 {
218 return sum(d1 * d2);
219 }
220
221 // we define a generalized cross product for all classes using the
222 // Levi-Civita symbol. Return type is nice (ever heard of "write-once
223 // code")... it took 10 times longer to figure out how to write the
224 // return type than to do everything else...
225
226 template<typename T1, typename T2>
227 inline
228
229 _bz_ArrayExpr<
230 _bz_ArrayExprReduce<
231 _bz_ArrayExpr<
232 _bz_ArrayExprReduce<
233 _bz_typename BzBinaryExprResult<
234 Multiply,
235 _bz_typename BzBinaryExprResult<
236 Multiply,
237 _bz_ArrayExpr<LeviCivita>,
238 _bz_ArrayExpr<
239 ArrayIndexMapping<
240 _bz_typename asExpr<T1>::T_expr, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0>
241 >
242 >::T_result,
243 _bz_ArrayExpr<
244 ArrayIndexMapping<
245 _bz_typename asExpr<T2>::T_expr, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0>
246 >
247 >::T_result,
248 2,
249 ReduceSum<
250 _bz_typename BzBinaryExprResult<Multiply,_bz_typename BzBinaryExprResult<Multiply,_bz_ArrayExpr<LeviCivita>,_bz_ArrayExpr<ArrayIndexMapping<_bz_typename asExpr<T1>::T_expr, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0> > >::T_result,_bz_ArrayExpr<ArrayIndexMapping<_bz_typename asExpr<T2>::T_expr, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0> > > ::T_result::T_numtype,
251 BZ_SUMTYPE(bzCC(_bz_typename BzBinaryExprResult<Multiply,_bz_typename BzBinaryExprResult<Multiply,_bz_ArrayExpr<LeviCivita>,_bz_ArrayExpr<ArrayIndexMapping<_bz_typename asExpr<T1>::T_expr, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0> > >::T_result,_bz_ArrayExpr<ArrayIndexMapping<_bz_typename asExpr<T2>::T_expr, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0> > > ::T_result::T_numtype))>
252 >
253 >,
254 1,
255 ReduceSum<BZ_SUMTYPE(bzCC(_bz_typename BzBinaryExprResult<Multiply,_bz_typename BzBinaryExprResult<Multiply,_bz_ArrayExpr<LeviCivita>,_bz_ArrayExpr<ArrayIndexMapping<_bz_typename asExpr<T1>::T_expr, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0> > >::T_result,_bz_ArrayExpr<ArrayIndexMapping<_bz_typename asExpr<T2>::T_expr, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0> > > ::T_result::T_numtype))>
256 >
257 >
258
259 //int
cross(const ETBase<T1> & d1,const ETBase<T2> & d2)260 cross(const ETBase<T1>& d1, const ETBase<T2>& d2)
261 {
262 _bz_ArrayExpr<
263 _bz_ArrayExprReduce<
264 _bz_ArrayExpr<
265 _bz_ArrayExprReduce<
266 _bz_typename BzBinaryExprResult<
267 Multiply,
268 _bz_typename BzBinaryExprResult<
269 Multiply,
270 _bz_ArrayExpr<LeviCivita>,
271 _bz_ArrayExpr<
272 ArrayIndexMapping<
273 _bz_typename asExpr<T1>::T_expr, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0>
274 >
275 >::T_result,
276 _bz_ArrayExpr<
277 ArrayIndexMapping<
278 _bz_typename asExpr<T2>::T_expr, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0>
279 >
280 >::T_result,
281 2,
282 ReduceSum<
283 _bz_typename BzBinaryExprResult<Multiply,_bz_typename BzBinaryExprResult<Multiply,_bz_ArrayExpr<LeviCivita>,_bz_ArrayExpr<ArrayIndexMapping<_bz_typename asExpr<T1>::T_expr, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0> > >::T_result,_bz_ArrayExpr<ArrayIndexMapping<_bz_typename asExpr<T2>::T_expr, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0> > > ::T_result::T_numtype,
284 BZ_SUMTYPE(bzCC(_bz_typename BzBinaryExprResult<Multiply,_bz_typename BzBinaryExprResult<Multiply,_bz_ArrayExpr<LeviCivita>,_bz_ArrayExpr<ArrayIndexMapping<_bz_typename asExpr<T1>::T_expr, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0> > >::T_result,_bz_ArrayExpr<ArrayIndexMapping<_bz_typename asExpr<T2>::T_expr, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0> > > ::T_result::T_numtype))>
285 >
286 >,
287 1,
288 ReduceSum<BZ_SUMTYPE(bzCC(_bz_typename BzBinaryExprResult<Multiply,_bz_typename BzBinaryExprResult<Multiply,_bz_ArrayExpr<LeviCivita>,_bz_ArrayExpr<ArrayIndexMapping<_bz_typename asExpr<T1>::T_expr, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0> > >::T_result,_bz_ArrayExpr<ArrayIndexMapping<_bz_typename asExpr<T2>::T_expr, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0> > > ::T_result::T_numtype))>
289 >
290 >* x;
291 //int a=*x;
292 return sum(sum(LeviCivita()*d1.unwrap()(tensor::j)*d2.unwrap()(tensor::k),
293 tensor::k),
294 tensor::j);
295 }
296
297
298 // "scalar" function converts something into a scalar expression. this
299 // is used to prevent parsing of multicomponent expressions as bitwise
300 // expressions between dissimilar containers.
301
302 template <typename T>
303 inline
304 _bz_ArrayExpr<_bz_ArrayExprConstant<T> >
scalar(const T & x)305 scalar(const T& x) {
306 return _bz_ArrayExpr<_bz_ArrayExprConstant<T> >(x);
307 }
308
309 #endif
310
311 }
312
313 #endif // BZ_ARRAY_FUNCS_H
314