00001 //# Math.h: Casacore interface to <math.h> and other scalar math functions 00002 //# Copyright (C) 1993,1994,1995,1996,1997,1998,1999,2000,2001 00003 //# Associated Universities, Inc. Washington DC, USA. 00004 //# 00005 //# This library is free software; you can redistribute it and/or modify it 00006 //# under the terms of the GNU Library General Public License as published by 00007 //# the Free Software Foundation; either version 2 of the License, or (at your 00008 //# option) any later version. 00009 //# 00010 //# This library is distributed in the hope that it will be useful, but WITHOUT 00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 00012 //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public 00013 //# License for more details. 00014 //# 00015 //# You should have received a copy of the GNU Library General Public License 00016 //# along with this library; if not, write to the Free Software Foundation, 00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. 00018 //# 00019 //# Correspondence concerning AIPS++ should be addressed as follows: 00020 //# Internet email: aips2-request@nrao.edu. 00021 //# Postal address: AIPS++ Project Office 00022 //# National Radio Astronomy Observatory 00023 //# 520 Edgemont Road 00024 //# Charlottesville, VA 22903-2475 USA 00025 //# 00026 //# $Id$ 00027 00028 #ifndef CASA_MATH_H 00029 #define CASA_MATH_H 00030 00031 #include <casacore/casa/aips.h> 00032 //# The following is to get abs(int) and (is)finite. 00033 #include <casacore/casa/math.h> 00034 #include <casacore/casa/stdlib.h> 00035 00036 // On some systems the following is needed to get the finite function 00037 #if defined (AIPS_SOLARIS) || defined(AIPS_IRIX) 00038 #include <ieeefp.h> 00039 #endif 00040 00041 namespace casacore { //# NAMESPACE CASACORE - BEGIN 00042 00043 00044 // <summary> 00045 // Casacore interface to math.h and other scalar math functions 00046 // </summary> 00047 // <reviewed reviewer="UNKNOWN" date="before2004/08/25" tests="" demos=""> 00048 // </reviewed> 00049 00050 // <synopsis> 00051 00052 // Casacore interface to <src><math.h></src>. You should include this file 00053 // rather than <src><math.h></src> directly. It will be used to cover up any 00054 // deficiencies in the system <src><math.h></src>. 00055 00056 // This file does not include things like element-by-element 00057 // array operations. See the 00058 // <linkto group="ArrayMath.h#Array mathematical operations">ArrayMath</linkto> 00059 // functions for these functions. 00060 00061 // This file includes the standard math library. Hence besides the functions 00062 // defined here the following functions are also available. 00063 // <srcblock> 00064 // Double sin(Double x) Sine function 00065 // Double cos(Double x) Cosine function 00066 // Double tan(Double x) Tangent function 00067 // Double asin(Double x) Inverse sine function 00068 // Double acos(Double x) Inverse cosine function 00069 // Double atan(Double x) Inverse tangent function 00070 // Double atan2(Double y, Double x) Four quandrant inverse tangent function 00071 // Double hypot(Double y, Double x) Euclidean distance sqrt(x*x+y*y) 00072 00073 // Double sinh(Double x) Hyperbolic sine 00074 // Double cosh(Double x) Hyperbolic cosine 00075 // Double tanh(Double x) Hyperbolic tangent 00076 // Double acosh(Double x) Inverse hyperbolic sine 00077 // Double asinh(Double x) Inverse hyperbolic cosine 00078 // Double atanh(Double x) Inverse hyperbolic tangent 00079 00080 // Double sqrt(Double x) Square root 00081 // Double cbrt(Double x) Cube root 00082 00083 // Double pow(Double x, Double y) x raised to the power of y 00084 // Double exp(Double x) Exponental function 00085 // Double expm1(Double x) exp(x)-1. Use when x is small. 00086 // Double log(Double x) Natural logarithm 00087 // Double log10(Double x) Base ten logarithm 00088 // Double log1p(Double x) log(x+1). Use when x is small 00089 00090 // Double j0(Double x) Bessel function of the first kind, zeroth order 00091 // Double j1(Double x) Bessel function of the first kind, first order 00092 // Double jn(Int n, Double x) Bessel function of the first kind nth order 00093 // Double y0(Double x) Bessel function of the second kind, zeroth order 00094 // Double y1(Double x) Bessel function of the second kind, first order 00095 // Double yn(Int n, Double x) Bessel function of the second kind, nth order 00096 // 00097 // Double lgamma(Double x) Natural Log of the absolute value of the gamma 00098 // function 00099 // Double lgamma_r(Double x, Int* sign) Same as lgamma. The sign of the gamma 00100 // function is returned in the second argument. 00101 00102 // Double erf(Double x) Error function 00103 // Double erfc(Double x) Complementary error function (1 - erf(x)). 00104 // Use for large x. 00105 00106 // Double ceil(Double x) Returns the least integral value greater than or 00107 // equal to x 00108 // Double floor(Double x) Returns the least integral value than than or 00109 // equal to x 00110 // Double rint(Double x) Round to an integer using the current direction. 00111 00112 // Double fabs(Double x) Absolute value of x 00113 // Double remainder(Double x, Double y) the remainder. x - y*Int(x/y) 00114 // Double fmod(Double x, Double y) As above. May differ by +/- y 00115 // Int isNaN(Double x) Returns 1 if x is a NaN, zero otherwise 00116 // Int ilogb(Double x) Unbiased exponent of x 00117 // Double logb(Double x) As above but returns floating point result 00118 // Double scalbn(Double x, Int n) x*2**n. Uses exponent manipulation. 00119 // Double scalb(Double x, Double n) x*2**n. As above but n is a Double 00120 // Double significand(Double x) Returns the fractional part of x 00121 // (between 1 and 2) 00122 // Double copysign(Double x, Double y) returns a value with the magnitude of 00123 // x and the sign bit of y. 00124 // Double nextafter(Double x, Double y) Returns the next machine representable 00125 // number after x in the direction specified by y 00126 // </srcblock> 00127 // 00128 00129 // This file also includes the standard C library (stdlib.h). This is to obtain 00130 // a definition of the following functions. 00131 // <srcblock> 00132 // Int abs(Int x) absolute value function 00133 // </srcblock> 00134 // </synopsis> 00135 00136 // <group name="Math interface for casacore"> 00137 00138 // Returns f1**f2. The Double precision version is defined in the standard 00139 // library. But many compilers are not good enough to automatically do the type 00140 // promotion. Hence these functions are explicitly defined. 00141 // <group> 00142 inline Float pow(Float f1, Double f2) {return Float(std::pow(Double(f1), f2));} 00143 inline Float pow(Double f1, Float f2) {return Float(std::pow(f1, Double(f2)));} 00144 inline Int pow(Int f1, Int f2) {return Int(std::pow(Double(f1), Double(f2)));} 00145 // </group> 00146 00147 // Return the integer "less than" point (i.e. the one further from zero if 00148 // "point" is negative. 00149 // <group> 00150 inline Int ifloor(Float point) 00151 { if (point >= 0.0) return Int (point); else return Int(point - 1.0); } 00152 inline Int ifloor(Double point) 00153 { if (point >= 0.0) return Int(point); else return Int(point - 1.0); } 00154 // </group> 00155 00156 // Functions to get the max or min of two numbers. 00157 // <group> 00158 inline Int max(Int a, Int b) { if (a > b) return a; else return b; } 00159 inline Int min(Int a, Int b) { if (a > b) return b; else return a; } 00160 00161 inline uInt max(uInt a, uInt b){ if (a>b) return a; else return b; } 00162 inline uInt min(uInt a, uInt b){ if (a>b) return b; else return a; } 00163 00164 inline uInt64 max(uInt64 a, uInt64 b){ if (a>b) return a; else return b; } 00165 inline uInt64 min(uInt64 a, uInt64 b){ if (a>b) return b; else return a; } 00166 00167 inline Double max(Double a, Double b) { if (a > b) return a; else return b; } 00168 inline Double min(Double a, Double b) { if (a > b) return b; else return a; } 00169 inline Double max(Double a, Float b) { if (a > b) return a; else return b; } 00170 inline Double min(Double a, Float b) { if (a > b) return b; else return a; } 00171 inline Double max(Float a, Double b) { if (a > b) return a; else return b; } 00172 inline Double min(Float a, Double b) { if (a > b) return b; else return a; } 00173 00174 inline Float max(Float a, Float b) { if (a > b) return a; else return b; } 00175 inline Float min(Float a, Float b) { if (a > b) return b; else return a; } 00176 // </group> 00177 00178 // Get the absolute value of uInt. Should already be defined 00179 // for integers in <src><stdlib.h></src>. Define it for uInts so that certain 00180 // compilers can resolve the ambiguity when used in a templated class. 00181 // <group> 00182 #if defined(AIPS_BSD) 00183 inline Int64 abs(Int64 Val) {return Val;} 00184 #else 00185 inline uInt abs(uInt Val) {return Val;} 00186 #endif 00187 // </group> 00188 00189 // Return the square of a value. 00190 // <group> 00191 inline Int square(Int val) {return val*val;} 00192 inline Int64 square(Int64 val) {return val*val;} 00193 inline Float square(Float val) {return val*val;} 00194 inline Double square(Double val) {return val*val;} 00195 // </group> 00196 00197 // Return the cube of a value. 00198 // <group> 00199 inline Int cube(Int val) {return val*val*val;} 00200 inline Int64 cube(Int64 val) {return val*val*val;} 00201 inline Float cube(Float val) {return val*val*val;} 00202 inline Double cube(Double val) {return val*val*val;} 00203 // </group> 00204 00205 // Return the sign of a value. 00206 // <group> 00207 inline Int sign(Int val) {return val<0 ? -1 : (val>0 ? 1:0);} 00208 inline Int64 sign(Int64 val) {return val<0 ? -1 : (val>0 ? 1:0);} 00209 inline Float sign(Float val) {return val<0 ? -1 : (val>0 ? 1:0);} 00210 inline Double sign(Double val) {return val<0 ? -1 : (val>0 ? 1:0);} 00211 // </group> 00212 00213 // Return the floor modulo as used by Python (unlike C); divisor sign is used. 00214 // Note that function fmod can be used for C behaviour; dividend sign is used. 00215 // In Python: 5%3=2 -5%3=1 5%-3=-1 -5%-3=-2 00216 // In C: 5%3=2 -5%3=-2 5%-3=2 -5%-3=-2 00217 // <group> 00218 inline Int floormod (Int x, Int y) 00219 { 00220 Int r = x%y; 00221 if (r != 0 && (x<0) != (y<0)) r+=y; 00222 return r; 00223 } 00224 inline Int64 floormod (Int64 x, Int64 y) 00225 { 00226 Int64 r = x%y; 00227 if (r != 0 && (x<0) != (y<0)) r+=y; 00228 return r; 00229 } 00230 inline Float floormod (Float x, Float y) 00231 { 00232 Float r = fmod(x,y); 00233 if (r != 0 && (x<0) != (y<0)) r+=y; 00234 return r; 00235 } 00236 inline Double floormod (Double x, Double y) 00237 { 00238 Double r = fmod(x,y); 00239 if (r != 0 && (x<0) != (y<0)) r+=y; 00240 return r; 00241 } 00242 // </group> 00243 00244 // Functions to return whether a value is "relatively" near another. Returns 00245 // <src> tol > abs(val2 - val1)/max(abs(val1),(val2))</src>. 00246 // If tol <= 0, returns val1 == val2. If either val is 0.0, take care of area 00247 // around the minimum number that can be represented. 00248 // <group> 00249 Bool near(uInt val1, uInt val2, Double tol = 1.0e-5); 00250 Bool near(Int val1, Int val2, Double tol = 1.0e-5); 00251 Bool near(Float val1, Float val2, Double tol = 1.0e-5); 00252 Bool near(Float val1, Double val2, Double tol = 1.0e-5); 00253 Bool near(Double val1, Float val2, Double tol = 1.0e-5); 00254 Bool near(Double val1, Double val2, Double tol = 1.0e-13); 00255 // </group> 00256 00257 // The "allNear" versions are aliases for the normal "near" versions. They 00258 // exist to make template functions that work for both arrays and scalars 00259 // easier to write. These functions should be moved to ArrayMath.h 00260 // <group> 00261 inline Bool allNear(uInt val1, uInt val2, Double tol = 1.0e-5) 00262 { return near(val1, val2, tol); } 00263 inline Bool allNear(Int val1, Int val2, Double tol = 1.0e-5) 00264 { return near(val1, val2, tol); } 00265 inline Bool allNear(Float val1, Double val2, Double tol = 1.0e-5) 00266 { return near(val1, val2, tol); } 00267 inline Bool allNear(Double val1, Float val2, Double tol = 1.0e-5) 00268 { return near(val1, val2, tol); } 00269 inline Bool allNear(Float val1, Float val2, Double tol = 1.0e-5) 00270 { return near(val1, val2, tol); } 00271 inline Bool allNear(Double val1, Double val2, Double tol = 1.0e-13) 00272 { return near(val1, val2, tol); } 00273 // </group> 00274 00275 // Functions to return whether a value is "absolutely" near another. Returns 00276 // <src> tol > abs(val2 - val1)</src> 00277 // <group> 00278 Bool nearAbs(uInt val1, uInt val2, Double tol = 1.0e-5); 00279 Bool nearAbs(Int val1, Int val2, Double tol = 1.0e-5); 00280 Bool nearAbs(Float val1, Float val2, Double tol = 1.0e-5); 00281 Bool nearAbs(Float val1, Double val2, Double tol = 1.0e-5); 00282 Bool nearAbs(Double val1, Float val2, Double tol = 1.0e-5); 00283 Bool nearAbs(Double val1, Double val2, Double tol = 1.0e-13); 00284 // </group> 00285 00286 // The "allNearAbs" versions are aliases for the normal "nearAbs" 00287 // versions. They exist to make template functions that work for both arrays 00288 // and scalars easier to write. These functions should be in ArrayMath.h 00289 // <group> 00290 inline Bool allNearAbs(uInt val1, uInt val2, uInt tol = 1) 00291 { return nearAbs(val1, val2, tol); } 00292 inline Bool allNearAbs(Int val1, Int val2, Int tol = 1) 00293 { return nearAbs(val1, val2, tol); } 00294 inline Bool allNearAbs(Float val1, Float val2, Double tol = 1.0e-5) 00295 { return nearAbs(val1, val2, tol); } 00296 inline Bool allNearAbs(Float val1, Double val2, Double tol = 1.0e-5) 00297 { return nearAbs(val1, val2, tol); } 00298 inline Bool allNearAbs(Double val1, Float val2, Double tol = 1.0e-5) 00299 { return nearAbs(val1, val2, tol); } 00300 inline Bool allNearAbs(Double val1, Double val2, Double tol = 1.0e-13) 00301 { return nearAbs(val1, val2, tol); } 00302 // </group> 00303 00304 00305 // Functions to test if a floating point number is finite. 00306 // It is if it is NaN nor infinity. 00307 // <group> 00308 inline Bool isFinite (const Float& val) 00309 { 00310 #if defined(AIPS_DARWIN) 00311 return std::isfinite(val); 00312 #else 00313 return finite(val); 00314 #endif 00315 } 00316 inline Bool isFinite (const Double& val) 00317 { 00318 #if defined(AIPS_DARWIN) 00319 return std::isfinite(val); 00320 #else 00321 return finite(val); 00322 #endif 00323 } 00324 // </group> 00325 00326 // Functions to test for IEEE NaN's. The Float variant uses an in-line 00327 // Macro examining the bit pattern (for portability and efficiency). The 00328 // Double version invokes the IEEE function isnan found in ieeefp.h or math.h 00329 // <group> 00330 inline Bool isNaN (const Float& val) 00331 { 00332 return (((*(Int *)&(val) & 0x7f800000) == 0x7f800000) && 00333 ((*(Int *)&(val) & 0x007fffff) != 0x00000000)); 00334 } 00335 inline Bool isNaN(Double val) 00336 { 00337 return ( std::isnan(val) ); 00338 } 00339 // </group> 00340 00341 // Round a number to <src>ndigit</src> significant digits, usually used 00342 // for formatting for printing. 00343 // <br>A non-integer <src>ndigit=N+F<src>, with integer N and fraction F, 00344 // is interpreted as follows. 00345 // For <src>x = A*10^B</src>, where B is an integer, A is rounded to N digits 00346 // if <src>A > 10^F</src>, otherwise N+1 digits. 00347 // <br>For the default 2.5, a value of 32157 is rounded to 32000, 00348 // while 22157 is rounded to 22200. 00349 Double roundDouble(Double val, Double ndigit=2.5); 00350 00351 // Functions that return IEEE NaN's. The specific NaN returned has all bits 00352 // set. This is 'quiet' NaN, and because the sign bit is set it may be 00353 // considered a negative number (but NaN's are not numbers!). 00354 // <group> 00355 Float floatNaN(); 00356 Double doubleNaN(); 00357 void setNaN(Float& val); 00358 void setNaN(Double& val); 00359 // </group> 00360 00361 // Functions to test for IEEE Infinity's. Should work for positive or negative 00362 // infinity. 00363 // <group> 00364 Bool isInf(Float val); 00365 Bool isInf(Double val); 00366 // </group> 00367 00368 // Functions that return an IEEE Infinity, (positive infinity). 00369 // <group> 00370 Float floatInf(); 00371 Double doubleInf(); 00372 void setInf(Float& val); 00373 void setInf(Double& val); 00374 // </group> 00375 // </group> 00376 00377 00378 } //# NAMESPACE CASACORE - END 00379 00380 #endif