Math.h

Go to the documentation of this file.

// # Math.h: Casacore interface to <math.h> and other scalar math functions
// # Copyright (C) 1993,1994,1995,1996,1997,1998,1999,2000,2001
// # Associated Universities, Inc. Washington DC, USA.
// #
// # This library is free software; you can redistribute it and/or modify it
// # under the terms of the GNU Library General Public License as published by
// # the Free Software Foundation; either version 2 of the License, or (at your
// # option) any later version.
// #
// # This library 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 Library General Public
// # License for more details.
// #
// # You should have received a copy of the GNU Library General Public License
// # along with this library; if not, write to the Free Software Foundation,
// # Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
// #
// # Correspondence concerning AIPS++ should be addressed as follows:
// #        Internet email: casa-feedback@nrao.edu.
// #        Postal address: AIPS++ Project Office
// #                        National Radio Astronomy Observatory
// #                        520 Edgemont Road
// #                        Charlottesville, VA 22903-2475 USA
 
#ifndef CASA_MATH_H
#define CASA_MATH_H
 
#include <casacore/casa/aips.h>
// # The following is to get abs(int) and (is)finite.
#include <casacore/casa/math.h>
#include <casacore/casa/stdlib.h>
 
// On some systems the following is needed to get the finite function
#if defined(AIPS_SOLARIS) || defined(AIPS_IRIX)
#include <ieeefp.h>
#endif
 
namespace casacore {  // # NAMESPACE CASACORE - BEGIN
 
// <summary>
// Casacore interface to math.h and other scalar math functions
// </summary>
// <reviewed reviewer="UNKNOWN" date="before2004/08/25" tests="" demos="">
// </reviewed>
 
// <synopsis>
 
// Casacore interface to <src><math.h></src>. You should include this file
// rather than <src><math.h></src> directly.  It will be used to cover up any
// deficiencies in the system <src><math.h></src>.
 
// This file does not include things like element-by-element
// array operations. See the
// <linkto group="ArrayMath.h#Array mathematical operations">ArrayMath</linkto>
// functions for these functions.
 
// This file includes the standard math library. Hence besides the functions
// defined here the following functions are also available.
// <srcblock>
// Double sin(Double x)       Sine function
// Double cos(Double x)       Cosine function
// Double tan(Double x)       Tangent function
// Double asin(Double x)      Inverse sine function
// Double acos(Double x)      Inverse cosine function
// Double atan(Double x)      Inverse tangent function
// Double atan2(Double y, Double x) Four quandrant inverse tangent function
// Double hypot(Double y, Double x) Euclidean distance sqrt(x*x+y*y)
 
// Double sinh(Double x)      Hyperbolic sine
// Double cosh(Double x)      Hyperbolic cosine
// Double tanh(Double x)      Hyperbolic tangent
// Double acosh(Double x)     Inverse hyperbolic sine
// Double asinh(Double x)     Inverse hyperbolic cosine
// Double atanh(Double x)     Inverse hyperbolic tangent
 
// Double sqrt(Double x)      Square root
// Double cbrt(Double x)      Cube root
 
// Double pow(Double x, Double y) x raised to the power of y
// Double exp(Double x)       Exponental function
// Double expm1(Double x)     exp(x)-1. Use when x is small.
// Double log(Double x)       Natural logarithm
// Double log10(Double x)     Base ten logarithm
// Double log1p(Double x)     log(x+1). Use when x is small
 
// Double j0(Double x)        Bessel function of the first kind, zeroth order
// Double j1(Double x)        Bessel function of the first kind, first order
// Double jn(Int n, Double x) Bessel function of the first kind nth order
// Double y0(Double x)        Bessel function of the second kind, zeroth order
// Double y1(Double x)        Bessel function of the second kind, first order
// Double yn(Int n, Double x) Bessel function of the second kind, nth order
//
// Double lgamma(Double x)    Natural Log of the absolute value of the gamma
//                            function
// Double lgamma_r(Double x, Int* sign) Same as lgamma. The sign of the gamma
//                            function is returned in the second argument.
 
// Double erf(Double x)       Error function
// Double erfc(Double x)      Complementary error function (1 - erf(x)).
//                            Use for large x.
 
// Double ceil(Double x)      Returns the least integral value greater than or
//                            equal to x
// Double floor(Double x)     Returns the least integral value than than or
//                            equal to x
// Double rint(Double x)      Round to an integer using the current direction.
 
// Double fabs(Double x)      Absolute value of x
// Double remainder(Double x, Double y) the remainder. x - y*Int(x/y)
// Double fmod(Double x, Double y) As above. May differ by +/- y
// Int isNaN(Double x) Returns 1 if x is a NaN, zero otherwise
// Int ilogb(Double x) Unbiased exponent of x
// Double logb(Double x)      As above but returns floating point result
// Double scalbn(Double x, Int n) x*2**n. Uses exponent manipulation.
// Double scalb(Double x, Double n) x*2**n. As above but n is a Double
// Double significand(Double x) Returns the fractional part of x
//                            (between 1 and 2)
// Double copysign(Double x, Double y)  returns a value with the magnitude of
//                                      x and the sign bit of y.
// Double nextafter(Double x, Double y) Returns the next machine representable
//                            number after x in the direction specified by y
// </srcblock>
//
 
// This file also includes the standard C library (stdlib.h). This is to obtain
// a definition of the following functions.
// <srcblock>
// Int abs(Int x)             absolute value function
// </srcblock>
// </synopsis>
 
// <group name="Math interface for casacore">
 
// Returns f1**f2. The Double precision version is defined in the standard
// library. But many compilers are not good enough to automatically do the type
// promotion. Hence these functions are explicitly defined.
// <group>
inline Float pow(Float f1, Double f2) { return Float(std::pow(Double(f1), f2)); }
inline Float pow(Double f1, Float f2) { return Float(std::pow(f1, Double(f2))); }
inline Int pow(Int f1, Int f2) { return Int(std::pow(Double(f1), Double(f2))); }
// </group>
 
// Return the integer "less than" point (i.e. the one further from zero if
// "point" is negative.
// <group>
inline Int ifloor(Float point) {
  if (point >= 0.0)
    return Int(point);
  else
    return Int(point - 1.0);
}
inline Int ifloor(Double point) {
  if (point >= 0.0)
    return Int(point);
  else
    return Int(point - 1.0);
}
// </group>
 
//  Functions to get the max or min of two numbers.
// <group>
inline Int max(Int a, Int b) {
  if (a > b)
    return a;
  else
    return b;
}
inline Int min(Int a, Int b) {
  if (a > b)
    return b;
  else
    return a;
}
 
inline uInt max(uInt a, uInt b) {
  if (a > b)
    return a;
  else
    return b;
}
inline uInt min(uInt a, uInt b) {
  if (a > b)
    return b;
  else
    return a;
}
 
inline uInt64 max(uInt64 a, uInt64 b) {
  if (a > b)
    return a;
  else
    return b;
}
inline uInt64 min(uInt64 a, uInt64 b) {
  if (a > b)
    return b;
  else
    return a;
}
 
inline Double max(Double a, Double b) {
  if (a > b)
    return a;
  else
    return b;
}
inline Double min(Double a, Double b) {
  if (a > b)
    return b;
  else
    return a;
}
inline Double max(Double a, Float b) {
  if (a > b)
    return a;
  else
    return b;
}
inline Double min(Double a, Float b) {
  if (a > b)
    return b;
  else
    return a;
}
inline Double max(Float a, Double b) {
  if (a > b)
    return a;
  else
    return b;
}
inline Double min(Float a, Double b) {
  if (a > b)
    return b;
  else
    return a;
}
 
inline Float max(Float a, Float b) {
  if (a > b)
    return a;
  else
    return b;
}
inline Float min(Float a, Float b) {
  if (a > b)
    return b;
  else
    return a;
}
// </group>
 
// Return the square of a value.
// <group>
inline Int square(Int val) { return val * val; }
inline Int64 square(Int64 val) { return val * val; }
inline Float square(Float val) { return val * val; }
inline Double square(Double val) { return val * val; }
// </group>
 
// Return the cube of a value.
// <group>
inline Int cube(Int val) { return val * val * val; }
inline Int64 cube(Int64 val) { return val * val * val; }
inline Float cube(Float val) { return val * val * val; }
inline Double cube(Double val) { return val * val * val; }
// </group>
 
// Return the sign of a value.
// <group>
inline Int sign(Int val) { return val < 0 ? -1 : (val > 0 ? 1 : 0); }
inline Int64 sign(Int64 val) { return val < 0 ? -1 : (val > 0 ? 1 : 0); }
inline Float sign(Float val) { return val < 0 ? -1 : (val > 0 ? 1 : 0); }
inline Double sign(Double val) { return val < 0 ? -1 : (val > 0 ? 1 : 0); }
// </group>
 
// Return the floor modulo as used by Python (unlike C); divisor sign is used.
// Note that function fmod can be used for C behaviour; dividend sign is used.
// In Python:   5%3=2  -5%3=1   5%-3=-1  -5%-3=-2
// In C:        5%3=2  -5%3=-2  5%-3=2   -5%-3=-2
// <group>
inline Int floormod(Int x, Int y) {
  Int r = x % y;
  if (r != 0 && (x < 0) != (y < 0)) r += y;
  return r;
}
inline Int64 floormod(Int64 x, Int64 y) {
  Int64 r = x % y;
  if (r != 0 && (x < 0) != (y < 0)) r += y;
  return r;
}
inline Float floormod(Float x, Float y) {
  Float r = fmod(x, y);
  if (r != 0 && (x < 0) != (y < 0)) r += y;
  return r;
}
inline Double floormod(Double x, Double y) {
  Double r = fmod(x, y);
  if (r != 0 && (x < 0) != (y < 0)) r += y;
  return r;
}
// </group>
 
// Functions to return whether a value is "relatively" near another. Returns
// <src> tol > abs(val2 - val1)/max(abs(val1),(val2))</src>.
// If tol <= 0, returns val1 == val2. If either val is 0.0, take care of area
// around the minimum number that can be represented.
// <group>
Bool near(uInt val1, uInt val2, Double tol = 1.0e-5);
Bool near(Int val1, Int val2, Double tol = 1.0e-5);
Bool near(Float val1, Float val2, Double tol = 1.0e-5);
Bool near(Float val1, Double val2, Double tol = 1.0e-5);
Bool near(Double val1, Float val2, Double tol = 1.0e-5);
Bool near(Double val1, Double val2, Double tol = 1.0e-13);
// </group>
 
// The "allNear" versions are aliases for the normal "near" versions. They
// exist to make template functions that work for both arrays and scalars
// easier to write. These functions should be moved to ArrayMath.h
// <group>
inline Bool allNear(uInt val1, uInt val2, Double tol = 1.0e-5) { return near(val1, val2, tol); }
inline Bool allNear(Int val1, Int val2, Double tol = 1.0e-5) { return near(val1, val2, tol); }
inline Bool allNear(Float val1, Double val2, Double tol = 1.0e-5) { return near(val1, val2, tol); }
inline Bool allNear(Double val1, Float val2, Double tol = 1.0e-5) { return near(val1, val2, tol); }
inline Bool allNear(Float val1, Float val2, Double tol = 1.0e-5) { return near(val1, val2, tol); }
inline Bool allNear(Double val1, Double val2, Double tol = 1.0e-13) {
  return near(val1, val2, tol);
}
// </group>
 
// Functions to return whether a value is "absolutely" near another. Returns
// <src> tol > abs(val2 - val1)</src>
// <group>
Bool nearAbs(uInt val1, uInt val2, Double tol = 1.0e-5);
Bool nearAbs(Int val1, Int val2, Double tol = 1.0e-5);
Bool nearAbs(Float val1, Float val2, Double tol = 1.0e-5);
Bool nearAbs(Float val1, Double val2, Double tol = 1.0e-5);
Bool nearAbs(Double val1, Float val2, Double tol = 1.0e-5);
Bool nearAbs(Double val1, Double val2, Double tol = 1.0e-13);
// </group>
 
// The "allNearAbs" versions are aliases for the normal "nearAbs"
// versions. They exist to make template functions that work for both arrays
// and scalars easier to write. These functions should be in ArrayMath.h
// <group>
inline Bool allNearAbs(uInt val1, uInt val2, uInt tol = 1) { return nearAbs(val1, val2, tol); }
inline Bool allNearAbs(Int val1, Int val2, Int tol = 1) { return nearAbs(val1, val2, tol); }
inline Bool allNearAbs(Float val1, Float val2, Double tol = 1.0e-5) {
  return nearAbs(val1, val2, tol);
}
inline Bool allNearAbs(Float val1, Double val2, Double tol = 1.0e-5) {
  return nearAbs(val1, val2, tol);
}
inline Bool allNearAbs(Double val1, Float val2, Double tol = 1.0e-5) {
  return nearAbs(val1, val2, tol);
}
inline Bool allNearAbs(Double val1, Double val2, Double tol = 1.0e-13) {
  return nearAbs(val1, val2, tol);
}
// </group>
 
// Functions to test if a floating point number is finite.
// It is if it is NaN nor infinity.
// <group>
inline Bool isFinite(const Float& val) {
#if defined(AIPS_DARWIN)
  return std::isfinite(val);
#else
  return finite(val);
#endif
}
inline Bool isFinite(const Double& val) {
#if defined(AIPS_DARWIN)
  return std::isfinite(val);
#else
  return finite(val);
#endif
}
// </group>
 
// Functions to test for IEEE NaN's.  The Float variant uses an in-line
// Macro examining the bit pattern (for portability and efficiency). The
// Double version invokes the IEEE function isnan found in ieeefp.h or math.h
// <group>
inline Bool isNaN(const Float& val) {
  return (((*(Int*)&(val) & 0x7f800000) == 0x7f800000) &&
          ((*(Int*)&(val) & 0x007fffff) != 0x00000000));
}
inline Bool isNaN(Double val) { return (std::isnan(val)); }
// </group>
 
// Round a number to <src>ndigit</src> significant digits, usually used
// for formatting for printing.
// <br>A non-integer <src>ndigit=N+F<src>, with integer N and fraction F,
// is interpreted as follows.
// For <src>x = A*10^B</src>, where B is an integer, A is rounded to N digits
// if <src>A > 10^F</src>, otherwise N+1 digits.
// <br>For the default 2.5, a value of 32157 is rounded to 32000,
// while 22157 is rounded to 22200.
Double roundDouble(Double val, Double ndigit = 2.5);
 
// Functions that return IEEE NaN's. The specific NaN returned has all bits
// set. This is 'quiet' NaN, and because the sign bit is set it may be
// considered a negative number (but NaN's are not numbers!).
// <group>
Float floatNaN();
Double doubleNaN();
void setNaN(Float& val);
void setNaN(Double& val);
// </group>
 
// Functions to test for IEEE Infinity's. Should work for positive or negative
// infinity.
// <group>
Bool isInf(Float val);
Bool isInf(Double val);
// </group>
 
// Functions that return an IEEE Infinity,  (positive infinity).
// <group>
Float floatInf();
Double doubleInf();
void setInf(Float& val);
void setInf(Double& val);
// </group>
// </group>
 
}  // namespace casacore
 
#endif