GRALE
realtype.h
1 /*
2 
3  This file is a part of GRALE, a library to facilitate the simulation
4  and inversion of gravitational lenses.
5 
6  Copyright (C) 2008-2012 Jori Liesenborgs
7 
8  Contact: jori.liesenborgs@gmail.com
9 
10  This program is free software; you can redistribute it and/or modify
11  it under the terms of the GNU General Public License as published by
12  the Free Software Foundation; either version 2 of the License, or
13  (at your option) any later version.
14 
15  This program is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  GNU General Public License for more details.
19 
20  You should have received a copy of the GNU General Public License
21  along with this program; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
23 
24 */
25 
26 #ifndef REALTYPE_H
27 
28 #define REALTYPE_H
29 
30 #include <string.h>
31 #include <cmath>
32 #ifdef GRALECONFIG_SUPPORT_INTELIPP
33  #include <ipps.h>
34 #endif // GRALECONFIG_SUPPORT_INTELIPP
35 
36 namespace grale
37 {
38 #ifndef WIN32
39  inline float asinh(float x)
40  {
41  return ::asinhf(x);
42  }
43 
44  inline double asinh(double x)
45  {
46  return ::asinh(x);
47  }
48 
49  inline long double asinh(long double x)
50  {
51  return ::asinhl(x);
52  }
53 
54  inline float acosh(float x)
55  {
56  return ::acoshf(x);
57  }
58 
59  inline double acosh(double x)
60  {
61  return ::acosh(x);
62  }
63 
64  inline long double acosh(long double x)
65  {
66  return ::acoshl(x);
67  }
68 
69  inline float atanh(float x)
70  {
71  return ::atanhf(x);
72  }
73 
74  inline double atanh(double x)
75  {
76  return ::atanh(x);
77  }
78 
79  inline long double atanh(long double x)
80  {
81  return ::atanhl(x);
82  }
83 
84  template<class T>
85  inline T min(T x, T y)
86  {
87  return (x < y)?x:y;
88  }
89 
90  template<class T>
91  inline T max(T x, T y)
92  {
93  return (x > y)?x:y;
94  }
95 #else
96  template<class T>
97  T asinh(T x)
98  {
99  return (std::exp(x)-std::exp(-x))/(T)2;
100  }
101 
102  template<class T>
103  T acosh(T x)
104  {
105  return (std::exp(x)+std::exp(-x))/(T)2;
106  }
107 
108  template<class T>
109  T atanh(T x)
110  {
111  return (std::exp(x)-std::exp(-x))/(std::exp(x)+std::exp(-x));
112  }
113 #endif
114 
115 #ifndef GRALECONFIG_SUPPORT_INTELIPP
116 
117  template<class T>
118  inline T CalculateDotProduct(const T *a, const T *b, int num)
119  {
120  T v = 0;
121  for (int i = 0 ; i < num ; i++)
122  v += a[i]*b[i];
123  return v;
124  }
125 
126  template<class T>
127  inline void CalculateAddProductC(T *pDest, const T *pSrc, T C, int num, T *tmpbuf)
128  {
129  for (int i = 0 ; i < num ; i++)
130  pDest[i] = pDest[i] + pSrc[i]*C;
131  }
132 
133  template<class T>
134  inline void CopyVector(T *pDest, const T *pSrc, int num)
135  {
136  memcpy(pDest, pSrc, num*sizeof(T));
137  }
138 
139  template<class T>
140  inline void SquareVector(T *pVect, int num)
141  {
142  for (int i = 0 ; i < num ; i++)
143  pVect[i] *= pVect[i];
144  }
145 
146  template<class T>
147  inline void SquareVector(T *pDest, const T *pSrc, int num)
148  {
149  for (int i = 0 ; i < num ; i++)
150  pDest[i] = pSrc[i]*pSrc[i];
151  }
152  template<class T>
153  inline void MultiplyVector(T *pSrcDst, const T *pSrc, int num)
154  {
155  for (int i = 0 ; i < num ; i++)
156  pSrcDst[i] *= pSrc[i];
157  }
158 
159  template<class T>
160  inline void MultiplyVector(T *pDst, const T *pSrc, const T factor, int num)
161  {
162  for (int i = 0 ; i < num ; i++)
163  pDst[i] = pSrc[i]*factor;
164  }
165 
166  template<class T>
167  inline void MultiplyVector(T *pSrcDst, const T factor, int num)
168  {
169  for (int i = 0 ; i < num ; i++)
170  pSrcDst[i] *= factor;
171  }
172 
173  template<class T>
174  inline void MultiplyVector(T *pDst, const T *pSrc1, const T *pSrc2, int num)
175  {
176  for (int i = 0 ; i < num ; i++)
177  pDst[i] = pSrc1[i]*pSrc2[i];
178  }
179 
180  template<class T>
181  inline void AbsVector(T *pVect, int num)
182  {
183  for (int i = 0 ; i < num ; i++)
184  pVect[i] = std::abs(pVect[i]);
185  }
186 
187  template<class T>
188  inline void AddVector(T *pSrcDst, const T *pSrc, int num)
189  {
190  for (int i = 0 ; i < num ; i++)
191  pSrcDst[i] += pSrc[i];
192  }
193 
194  template<class T>
195  inline void AddVector(T *pSrcDst, T value, int num)
196  {
197  for (int i = 0 ; i < num ; i++)
198  pSrcDst[i] += value;
199  }
200 
201  template<class T>
202  inline void AddVector(T *pDst, const T *pSrc1, const T *pSrc2, int num)
203  {
204  for (int i = 0 ; i < num ; i++)
205  pDst[i] = pSrc1[i] + pSrc2[i];
206  }
207 
208  template<class T>
209  inline void SubVector(T *pSrcDst, const T *pSrc, int num)
210  {
211  for (int i = 0 ; i < num ; i++)
212  pSrcDst[i] -= pSrc[i];
213  }
214 
215  template<class T>
216  inline void SubVector(T *pDst, const T *pSrc1, const T *pSrc2, int num)
217  {
218  for (int i = 0 ; i < num ; i++)
219  pDst[i] = pSrc1[i] - pSrc2[i];
220  }
221 
222 
223 #else
224 
225 //#error "TODO: add float versions"
226 
227  inline double CalculateDotProduct(const double *a, const double *b, int num)
228  {
229  double val;
230  ippsDotProd_64f(a, b, num, &val);
231  return val;
232  }
233 
234  inline void CalculateAddProductC(double *pDest, const double *pSrc, double C, int num, double *tmpbuf)
235  {
236  ippsMulC_64f(pSrc, C, tmpbuf, num);
237  ippsAdd_64f_I(tmpbuf, pDest, num);
238  }
239 
240  inline void CopyVector(double *pDest, const double *pSrc, int num)
241  {
242  ippsCopy_64f(pSrc, pDest, num);
243  }
244 
245  inline void SquareVector(double *pVect, int num)
246  {
247  ippsSqr_64f_I(pVect, num);
248  }
249 
250  inline void MultiplyVector(double *pSrcDst, const double *pSrc, int num)
251  {
252  ippsMul_64f_I(pSrc, pSrcDst, num);
253  }
254 
255  inline void MultiplyVector(double *pSrcDst, double factor, int num)
256  {
257  ippsMulC_64f_I(factor, pSrcDst, num);
258  }
259 
260  inline void AbsVector(double *pVect, int num)
261  {
262  ippsAbs_64f_I(pVect, num);
263  }
264 
265  inline void AddVector(double *pSrcDst, const double *pSrc, int num)
266  {
267  ippsAdd_64f_I(pSrc, pSrcDst, num);
268  }
269 
270  inline void SubVector(double *pSrcDst, const double *pSrc, int num)
271  {
272  ippsSub_64f_I(pSrc, pSrcDst, num);
273  }
274 
275  inline void SubVector(double *pDst, const double *pSrc1, const double *pSrc2, int num)
276  {
277  ippsSub_64f(pSrc2, pSrc1, pDst, num);
278  }
279 
280  // float versions
281 
282  inline float CalculateDotProduct(const float *a, const float *b, int num)
283  {
284  float val;
285  ippsDotProd_32f(a, b, num, &val);
286  return val;
287  }
288 
289  inline void CalculateAddProductC(float *pDest, const float *pSrc, float C, int num, float *tmpbuf)
290  {
291  ippsMulC_32f(pSrc, C, tmpbuf, num);
292  ippsAdd_32f_I(tmpbuf, pDest, num);
293  }
294 
295  inline void CopyVector(float *pDest, const float *pSrc, int num)
296  {
297  ippsCopy_32f(pSrc, pDest, num);
298  }
299 
300  inline void SquareVector(float *pVect, int num)
301  {
302  ippsSqr_32f_I(pVect, num);
303  }
304 
305  inline void SquareVector(float *pDst, const float *pSrc, int num)
306  {
307  ippsSqr_32f(pSrc, pDst, num);
308  }
309 
310  inline void MultiplyVector(float *pSrcDst, const float *pSrc, int num)
311  {
312  ippsMul_32f_I(pSrc, pSrcDst, num);
313  }
314 
315  inline void MultiplyVector(float *pDst, const float *pSrc, float factor, int num)
316  {
317  ippsMulC_32f(pSrc, factor, pDst, num);
318  }
319 
320  inline void MultiplyVector(float *pDst, const float *pSrc1, const float *pSrc2, int num)
321  {
322  ippsMul_32f(pSrc1, pSrc2, pDst, num);
323  }
324 
325  inline void AbsVector(float *pVect, int num)
326  {
327  ippsAbs_32f_I(pVect, num);
328  }
329 
330  inline void AddVector(float *pSrcDst, const float *pSrc, int num)
331  {
332  ippsAdd_32f_I(pSrc, pSrcDst, num);
333  }
334 
335  inline void AddVector(float *pDst, const float *pSrc1, const float *pSrc2, int num)
336  {
337  ippsAdd_32f(pSrc1, pSrc2, pDst, num);
338  }
339 
340  inline void AddVector(float *pSrcDst, float value, int num)
341  {
342  ippsAddC_32f_I(value, pSrcDst, num);
343  }
344 
345  inline void SubVector(float *pSrcDst, const float *pSrc, int num)
346  {
347  ippsSub_32f_I(pSrc, pSrcDst, num);
348  }
349 
350  inline void SubVector(float *pDst, const float *pSrc1, const float *pSrc2, int num)
351  {
352  ippsSub_32f(pSrc2, pSrc1, pDst, num);
353  }
354 
355  inline void MultiplyVector(float *pSrcDst, float factor, int num)
356  {
357  ippsMulC_32f_I(factor, pSrcDst, num);
358  }
359 
360 #endif // !GRALECONFIG_SUPPORT_INTELIPP
361 
362 
363 }
364 
365 #define ABS(x) std::abs(x)
366 #define SQRT(x) std::sqrt(x)
367 #define COS(x) std::cos(x)
368 #define SIN(x) std::sin(x)
369 #define TAN(x) std::tan(x)
370 #define ACOS(x) std::acos(x)
371 #define ASIN(x) std::asin(x)
372 #define ATAN(x) std::atan(x)
373 #define LN(x) std::log(x)
374 #define ASINH(x) grale::asinh(x)
375 #define ACOSH(x) grale::acosh(x)
376 #define ATANH(x) grale::atanh(x)
377 #define POW(x,y) std::pow(x,y)
378 #define EXP(x) std::exp(x)
379 #define SINH(x) std::sinh(x)
380 #define COSH(x) std::cosh(x)
381 #define TANH(x) std::tanh(x)
382 
383 #ifndef WIN32
384 #define MIN(x,y) grale::min(x,y)
385 #define MAX(x,y) grale::max(x,y)
386 #else
387 #define MIN(x,y) (((x)<(y))?(x):(y))
388 #define MAX(x,y) (((x)>(y))?(x):(y))
389 #endif
390 
391 
392 #endif // REALTYPE_H
393