Ray: /ray/src/lib/numerics/num

00001 /*
00002  * lib/numerics/num_math.h
00003  * 
00004  * Numerics library general maths header. 
00005  * 
00006  * Copyright (c) 2003--2004 by Wolfgang Wieser ] wwieser (a) gmx <*> de [ 
00007  * 
00008  * This file may be distributed and/or modified under the terms of the 
00009  * GNU General Public License version 2 as published by the Free Software 
00010  * Foundation. (See COPYING.GPL for details.)
00011  * 
00012  * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
00013  * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
00014  * 
00015  */
00016 
00017 #ifndef _LIB_NUMERICS_MATH_H_
00018 #define _LIB_NUMERICS_MATH_H_
00019 
00020 #include <lib/sconfig.h>    /* MUST be first */
00021 
00022 
00047 #include <math.h>
00048 
00071 //     HUGE_VAL, INFINITY, ...
00072 
00073 
00074 
00076 namespace NUM  // numerics
00077 {
00078 
00079 // Make sure to cancel the most common one :)
00080 #undef M_PI
00081 
00082 // This is somewhat ugly to support compilers which do not support 
00083 // static const floating point members directly assigned. 
00084 // In case we do not provide the FP numbers here, they are provided 
00085 // externally in fpconst.cc. 
00086 #if USE_STATIC_CONST_FP_MEMB || defined(DOXYGEN_IS_SCANNING)
00087 #  define SCFMinit(x)   =x
00088 #else
00089 #  define SCFMinit(x)
00090 #endif
00091 
00098 template<typename T>struct FPConst
00099 {
00101     T macheps1;
00103     T sqrt_macheps1;
00104     
00107     T macheps1_prec;
00108     
00109     // Some more constants which need not be computed like macheps1. 
00110     // NOTE: If you add somehing here, DO NOT FORGET to update the list 
00111     //       in fpconst.cc. 
00112     static const T pi       SCFMinit(3.1415926535897932384626433832795029);   
00113     static const T div_pi_2 SCFMinit(1.5707963267948966192313216916397514);   
00114     static const T div_pi_3 SCFMinit(1.047197551196597746154214461093168);    
00115     static const T div_pi_4 SCFMinit(0.7853981633974483096156608458198757);   
00116     static const T mul_pi_2 SCFMinit(6.283185307179586476925286766559006);    
00117     static const T mul_pi_3 SCFMinit(9.424777960769379715387930149838509);    
00118     static const T mul_pi_4 SCFMinit(12.56637061435917295385057353311801);    
00119     static const T mul_pi_2_3 SCFMinit(2.094395102393195492308428922186335);  
00120     static const T div_1_pi SCFMinit(0.3183098861837906715377675267450287);   
00121     static const T div_2_pi SCFMinit(0.6366197723675813430755350534900574);   
00122     static const T div_2_sqrt_pi SCFMinit(1.1283791670955125738961589031215452); 
00123     
00124     static const T sqrt_2  SCFMinit(1.4142135623730950488016887242096981);       
00125     static const T div_1_sqrt_2 SCFMinit(0.7071067811865475244008443621048490);  
00126     
00127     static const T e       SCFMinit(2.7182818284590452353602874713526625);    
00128     static const T log2_e  SCFMinit(1.4426950408889634073599246810018921);    
00129     static const T log10_e SCFMinit(0.4342944819032518276511289189166051);    
00130     static const T ln_2    SCFMinit(0.6931471805599453094172321214581766);    
00131     static const T ln_10   SCFMinit(2.3025850929940456840179914546843642);    
00132     
00133     static const T div_1_pow_2_31   SCFMinit(1.0/2147483648.0);        
00134     static const T div_1_pow_2_31m1 SCFMinit(1.0/2147483647.0);        
00135     static const T div_1_pow_2_32   SCFMinit(1.0/4294967296.0);        
00136     static const T div_1_pow_2_32m1 SCFMinit(1.0/4294967295.0);        
00137     static const T div_1_pow_2_53   SCFMinit(1.0/9007199254740992.0);  
00138     static const T div_1_pow_2_53m1 SCFMinit(1.0/9007199254740991.0);  
00139     
00140     //static const T div_1_31 SCFMinit(1.0/31.0);      ///< 1/31
00141     //static const T div_1_63 SCFMinit(1.0/63.0);      ///< 1/63
00142     //static const T div_1_255 SCFMinit(1.0/255.0);     ///< 1/255
00143     //static const T div_1_65535 SCFMinit(1.0/65535.0);   ///< 1/65535
00144     
00150     static const T huge_val;
00151     
00153     static const T max_val;
00155     static const T min_val;
00156 };
00157 
00158 // Get rid of that again...
00159 #undef SCFMinit
00160 
00162 extern const FPConst<dbl> fpconst_dbl;
00164 extern const FPConst<flt> fpconst_flt;
00165 
00166 
00167 //==============================================================================
00168 // BASIC FUNCTIONS START HERE.
00169 //==============================================================================
00170 
00171 #ifdef sgn
00172 #error sgn defined
00173 #endif
00174 
00180 template<typename T>inline _constfn_ int sgn(T x)
00181     {  return(x>0 ? 1 : (x<0 ? -1 : 0));  }
00182 
00191 template<typename T>inline _constfn_ int sgn11(T x)
00192     {  return(x<0 ? -1 : 1);  }
00193 
00194 
00195 #ifndef DOXYGEN_IS_SCANNING
00196 
00197 // Test sign of a floating point number, and also of integers. 
00198 // (Fast) specialisation for FP numbers immediately below. 
00199 template<typename T>inline _constfn_ bool is_negative(T x)
00200     {  return(x<0);  }
00201 // FIXME: Tests indicate that "x<0" is faster than "signbit". WTH?
00202 // Ultra-fast (factor 5 over both!) sign checking can be done by using 
00203 // a more clever algo than glibc signbit() implementation like 
00204 //   return(((int32*)&x)[1]<0);  <-- for 8byte double
00205 // Similar holds for finite(), isnan(), isinf(). 
00206 inline _constfn_ bool is_negative(double x)
00207     {  return(signbit(x));  }  // <-- NOTE: Not ::signbit() as it is a macro. 
00208 inline _constfn_ bool is_negative(float x)
00209     {  return(signbit(x));  }  // <-- NOTE: Not ::signbit() as it is a macro. 
00210 inline _constfn_ bool is_negative(long double x)
00211     {  return(signbit(x));  }  // <-- NOTE: Not ::signbit() as it is a macro. 
00212 
00213 
00214 // Test for infinity of floating point numbers. 
00215 #if HAVE_ISINF
00216     template<typename T>inline _constfn_ bool is_inf(T x)
00217         {  return(isinf(x));  }  // <-- NOTE: actually a macro. 
00218 #elif HAVE_FPCLASSIFY
00219     template<typename T>inline _constfn_ bool is_inf(T x)
00220         {  return(fpclassify(x)==FP_INFINITE);  }  // <-NOTE: actually a macro. 
00221 #else
00222 #  error "You lack isinf()/fpclassify()."
00223 #endif
00224 
00225 // Test for NaN of floating point numbers. 
00226 #if HAVE_ISNAN
00227     template<typename T>inline _constfn_ bool is_nan(T x)
00228         {  return(isnan(x));  }  // <-- NOTE: actually a macro. 
00229 #elif HAVE_FPCLASSIFY
00230     template<typename T>inline _constfn_ bool is_nan(T x)
00231         {  return(fpclassify(x)==FP_NAN);  }  // <-NOTE: actually a macro. 
00232 #else
00233 #  error "You lack isnan()/fpclassify()."
00234 #endif
00235 
00236 // Test for finite floating point numbers. 
00237 #if HAVE_FINITE
00238     template<typename T>inline _constfn_ bool is_finite(T x)
00239         {  return(finite(x));  }  // <-- NOTE: actually a macro. 
00240 #elif HAVE_ISFINITE
00241     template<typename T>inline _constfn_ bool is_finite(T x)
00242         {  return(isfinite(x));  }  // <-- NOTE: actually a macro. 
00243 #elif HAVE_FPCLASSIFY
00244     template<typename T>inline _constfn_ bool is_finite(T x)
00245         {  return(fpclassify(x)!=FP_NAN && fpclassify(x)!=FP_INFINITE);  }
00246 #else
00247 #  error "You lack finite()/isfinite()/fpclassify()."
00248 #endif
00249 
00250 
00251 #if defined(sqr) || defined(SQR)
00252 #error sqr defined
00253 #endif
00254 
00255 template<typename T>inline _constfn_ T sqr(T x)  {  return(x*x);  }
00256 
00257 #ifdef fabs
00258 #error fabs defined
00259 #endif
00260 // Return absolute value of x; this is the general template; specialisation 
00261 // using NUM_OVERLOADED(fabs) below. 
00262 template<typename T>inline T fabs(T x)  {  return(x<0 ? -x : x);  }
00263 
00264 #if defined(max) || defined(min)
00265 #error max or min defined
00266 #endif
00267 
00268 template<typename T>inline T max(T a,T b)
00269     {  return(a<b ? b : a);  }
00271 template<typename T>inline T min(T a,T b)
00272     {  return(a<b ? a : b);  }
00273 
00274 
00275 #define NUM_OVERLOADED_R(func,cfunc) \
00276     inline double func(double x)  {  return(::cfunc(x));  } \
00277     inline float func(float x)    {  return(::cfunc ## f(x));  } \
00278     inline long double func(long double x)  {  return(::cfunc ## l(x));  }
00279 #define NUM_OVERLOADED(func) NUM_OVERLOADED_R(func,func)
00280 
00281 #define NUM_OVERLOADED2_R(func,cfunc) \
00282     inline double func(double x,double y)  {  return(::cfunc(x,y));  } \
00283     inline float func(float x,float y)     {  return(::cfunc ## f(x,y));  } \
00284     inline long double func(long double x,long double y) \
00285         {  return(::cfunc ## l(x,y));  }
00286 #define NUM_OVERLOADED2(func) NUM_OVERLOADED2_R(func,func)
00287 
00288 #if defined(sqrt) || defined(sin) || defined(cos) || defined(tan) || defined(exp) || defined(log)
00289 #error defined math functions
00290 #endif
00291 
00292 NUM_OVERLOADED(fabs)
00293 NUM_OVERLOADED(sqrt)  NUM_OVERLOADED2(pow)
00294 NUM_OVERLOADED(sin)   NUM_OVERLOADED(asin)
00295 NUM_OVERLOADED(cos)   NUM_OVERLOADED(acos)
00296 NUM_OVERLOADED(sinh)  NUM_OVERLOADED(asinh)
00297 NUM_OVERLOADED(cosh)  NUM_OVERLOADED(acosh)
00298 NUM_OVERLOADED(tan)   NUM_OVERLOADED(atan)  NUM_OVERLOADED2(atan2)
00299 NUM_OVERLOADED(tanh)  NUM_OVERLOADED(atanh)
00300 NUM_OVERLOADED(exp)   NUM_OVERLOADED(log)
00301 
00302 NUM_OVERLOADED(ceil)  NUM_OVERLOADED(floor)
00303 NUM_OVERLOADED(trunc)
00304 
00305 NUM_OVERLOADED2(remainder)
00306 
00307 // FP specialisation of the min/max functions if available. Otherwise, use 
00308 // general version above. 
00309 #if HAVE_FMAX
00310     NUM_OVERLOADED2_R(max,fmax)
00311 #endif
00312 #if HAVE_FMIN
00313     NUM_OVERLOADED2_R(min,fmin)
00314 #endif
00315 // Do not use these; use plain max(), min() instead. 
00316 extern _deprecated_ dbl fmax(dbl x,dbl y);
00317 extern _deprecated_ flt fmaxf(flt x,flt y);
00318 extern _deprecated_ dbl fmin(dbl x,dbl y);
00319 extern _deprecated_ flt fminf(flt x,flt y);
00320 
00321 // These are currently not used. 
00322 // There are HAVE_xxx #defines for all 3 of them because some systems do not 
00323 // support them. 
00324 //  NUM_OVERLOADED(rint)
00325 //  NUM_OVERLOADED(nearbyint)
00326 //  NUM_OVERLOADED(round)
00327 extern _deprecated_ dbl rint(dbl x);  // do not use (currently)
00328 extern _deprecated_ dbl nearbyint(dbl x);  // do not use (currently)
00329 extern _deprecated_ dbl round(dbl x);  // do not use (currently)
00330 
00331 // Round to nearest integer value using CURRENT rounding direction. 
00332 // Specialisations below. 
00333 template<typename I,typename FP>inline I iroundc(FP x);
00334 // FIXME: We should change that for 64 bit systems. 
00335 template<> inline int32 iroundc<int32,double>     (double x)      { return(::lrint(x));   }
00336 template<> inline int32 iroundc<int32,float>      (float x)       { return(::lrintf(x));  }
00337 template<> inline int32 iroundc<int32,long double>(long double x) { return(::lrintl(x));  }
00338 template<> inline int64 iroundc<int64,double>     (double x)      { return(::llrint(x));  }
00339 template<> inline int64 iroundc<int64,float>      (float x)       { return(::llrintf(x)); }
00340 template<> inline int64 iroundc<int64,long double>(long double x) { return(::llrintl(x)); }
00341 // "Overloaded" for 8 and 16 bit integers internally using the 32bit version. 
00342 #define IROUND_OVERLOADED(func,itype) \
00343     template<> inline itype func<itype,double>(double x) \
00344         { return(NUM::func<int32>(x)); } \
00345     template<> inline itype func<itype,float>(float x) \
00346         { return(NUM::func<int32>(x)); } \
00347     template<> inline itype func<itype,long double>(long double x) \
00348         { return(NUM::func<int32>(x)); }
00349 IROUND_OVERLOADED(iroundc,int8);
00350 IROUND_OVERLOADED(iroundc,uint8);
00351 IROUND_OVERLOADED(iroundc,int16);
00352 IROUND_OVERLOADED(iroundc,uint16);
00353 // To allow for compiler selection: 
00354 template<typename I,typename FP>inline I iroundc(I *d,FP x)
00355     {  return(*d=NUM::iroundc<I,FP>(x));  }
00356 
00357 // Round to nearest integer value rounding away from zero. 
00358 // Specialisations below. 
00359 template<typename I,typename FP>inline I iround(FP x);
00360 // FIXME: We should change that for 64 bit systems. 
00361 template<> inline int32 iround<int32,double>     (double x)      { return(::lround(x));   }
00362 template<> inline int32 iround<int32,float>      (float x)       { return(::lroundf(x));  }
00363 template<> inline int32 iround<int32,long double>(long double x) { return(::lroundl(x));  }
00364 template<> inline int64 iround<int64,double>     (double x)      { return(::llround(x));  }
00365 template<> inline int64 iround<int64,float>      (float x)       { return(::llroundf(x)); }
00366 template<> inline int64 iround<int64,long double>(long double x) { return(::llroundl(x)); }
00367 // "Overloaded" for 8 and 16 bit integersinternally using the 32bit version. 
00368 IROUND_OVERLOADED(iround,int8);
00369 IROUND_OVERLOADED(iround,uint8);
00370 IROUND_OVERLOADED(iround,int16);
00371 IROUND_OVERLOADED(iround,uint16);
00372 #undef IROUND_OVERLOADED
00373 // To allow for compiler selection: 
00374 template<typename I,typename FP>inline I iround(I *d,FP x)
00375     {  return(*d=NUM::iround<I,FP>(x));  }
00376 
00377 
00378 #if HAVE_CBRT
00379     NUM_OVERLOADED(cbrt)
00380 #else
00381     template<typename T>inline T cbrt(T x)
00382         {  return( NUM::is_negative(x) ? 
00383             -NUM::pow(-x,0.3333333333333333333333333333333333) : 
00384              NUM::pow( x,0.3333333333333333333333333333333333) );  }
00385 #endif
00386 
00387 #if HAVE_EXP10
00388     NUM_OVERLOADED(exp10)
00389 #else
00390     template<typename T>inline T exp10(T x)
00391         {  return(NUM::pow(T(10),x));  }
00392 #endif
00393 
00394 #if HAVE_LOG10
00395     NUM_OVERLOADED(log10)
00396 #else
00397     template<typename T>inline T log10(T x)
00398         {  return(NUM::log(x)/FPConst<T>::ln_10);  }
00399 #endif
00400 
00401 #if HAVE_EXP2
00402     NUM_OVERLOADED(exp2)
00403 #else
00404     template<typename T>inline T exp2(T x)
00405         {  return(NUM::pow(T(2),x));  }
00406 #endif
00407 
00408 #if HAVE_LOG2
00409     NUM_OVERLOADED(log2)
00410 #elif HAVE_LOGB
00411     NUM_OVERLOADED_R(log2,logb)
00412 #else
00413     template<typename T>inline T log2(T x)
00414         {  return(NUM::log(x)/FPConst<T>::ln_2);  }
00415 #endif
00416 // Do not use logb(), use log2() instead. 
00417 extern _deprecated_ dbl logb(dbl x);
00418 extern _deprecated_ flt logbf(flt x);
00419 
00420 #if HAVE_ILOGB
00421     inline int ilog2(double x)       {  return(::ilogb(x));  }
00422     inline int ilog2(float x)        {  return(::ilogbf(x));  }
00423     inline int ilog2(long double x)  {  return(::ilogbl(x));  }
00424 #else
00425     // Verified to work correctly (WW, 12/2004). 
00426     template<typename T>inline int ilog2(T x)
00427         {  return((int)NUM::floor(NUM::log2(x)));  }
00428 #endif
00429 // Do not use ilogb(), use ilog2() instead. 
00430 extern _deprecated_ int ilogb(dbl x);
00431 extern _deprecated_ int ilogbf(flt x);
00432 
00433 #if HAVE_EXPM1
00434     NUM_OVERLOADED(expm1)
00435 #else
00436     template<typename T>inline T expm1(T x)
00437         {  return(NUM::exp(x)-1);  }
00438 #endif
00439 
00440 #if HAVE_LOG1P
00441     NUM_OVERLOADED(log1p)
00442 #else
00443     template<typename T>inline T log1p(T x)
00444         {  return(NUM::log(1+x));  }
00445 #endif
00446 
00447 #if HAVE_COPYSIGN
00448     NUM_OVERLOADED2(copysign)
00449 #else
00450     template<typename T>inline T copysign(T x,T y)
00451         {  return(NUM::is_negative(x)==NUM::is_negative(y) ? x : -x);  }
00452 #endif
00453 
00454 #if HAVE_FDIM
00455     NUM_OVERLOADED2(fdim)
00456 #else
00457     template<typename T>inline T fdim(T x,T y)
00458         {  T tmp=x-y;  return(is_negative(tmp) ? T(0) : tmp);  }
00459 #endif
00460 
00461 #if HAVE_FMA
00462     inline double fma(double x,double y,double z)  {  return(::fma(x,y,z));  }
00463     inline float  fma(float x,float y,float z)     {  return(::fmaf(x,y,z));  }
00464     inline long double fma(long double x,long double y,long double z)
00465         {  return(::fmal(x,y,z));  }
00466 #else
00467     template<typename T>inline T fma(T x,T y,T z)
00468         {  return(x*y+z));  }
00469 #endif
00470 
00471 #if defined(sincos)
00472 #error defined math function sincos
00473 #endif
00474 #if HAVE_SINCOS
00475     inline void sincos(double x,double *sinx,double *cosx)
00476         {  ::sincos(x,sinx,cosx);  }
00477     inline void sincos(float x,float *sinx,float *cosx)
00478         {  ::sincosf(x,sinx,cosx);  }
00479     inline void sincos(long double x,long double *sinx,long double *cosx)
00480         {  ::sincosl(x,sinx,cosx);  }
00481 #else
00482     template<typename T>inline void sincos(T x,T *sinx,T *cosx)
00483         {  *sinx=NUM::sin(x);  *cosx=NUM::cos(x);  }
00484 #endif
00485 
00486 // The ldexp() function returns the result of multiplying the 
00487 // floating-point number x by 2 raised to the power exp.
00488 inline double      ldexp(double x,int exp)       {  return(::ldexp(x,exp));  }
00489 inline float       ldexp(float x,int exp)        {  return(::ldexpf(x,exp));  }
00490 inline long double ldexp(long double x,int exp)  {  return(::ldexpl(x,exp));  }
00491 // Do not use scalbn() or scalb(), use ldexp() instead. 
00492 extern _deprecated_ dbl scalbn(dbl x,int exp);
00493 extern _deprecated_ dbl scalb(dbl x,int exp);
00494 
00495 // The modf() function breaks the argument x into an integral part and 
00496 // a fractional part, each of which has the same sign as x. 
00497 // The integral part is stored in iptr.
00498 inline double modf(double x,double *iptr)
00499     {  return(::modf(x,iptr));  }
00500 inline float modf(float x,float *iptr)
00501     {  return(::modff(x,iptr));  }
00502 inline long double modf(long double x,long double *iptr)
00503     {  return(::modfl(x,iptr));  }
00504 
00505 // The frexp() function is used to split the number x into a normalized 
00506 // fraction and an exponent which is stored in exp.
00507 inline double      frexp(double x,int *exp)       {  return(::frexp(x,exp));  }
00508 inline float       frexp(float x,int *exp)        {  return(::frexpf(x,exp));  }
00509 inline long double frexp(long double x,int *exp)  {  return(::frexpl(x,exp));  }
00510 
00511 // This is not to be used. Use remainder() instead. 
00512 extern _deprecated_ dbl drem(dbl x,dbl y);   // do not use
00513 extern _deprecated_ flt drem(flt x,flt y);   // do not use
00514 extern _deprecated_ flt dremf(flt x,flt y);  // do not use
00515 
00516 #undef NUM_OVERLOADED
00517 #undef NUM_OVERLOADED2
00518 #undef NUM_OVERLOADED_R
00519 #undef NUM_OVERLOADED2_R
00520 
00521 #else  /* defined(DOXYGEN_IS_SCANNING) */
00522 
00523 // THIS IS ONLY DOCU FOR DOXYGEN. 
00524 #error "This should not be compiled; only for doxygen."
00525 
00535 template<typename T>inline bool is_inf(T x);  // only for doxygen
00536 
00546 template<typename T>inline bool is_nan(T x);  // only for doxygen
00547 
00557 template<typename T>inline bool is_finite(T x);  // only for doxygen
00558 
00559 
00569 template<typename T>inline bool is_negative(T x);  // only for doxygen
00570 
00581 template<typename T>inline T fabs(T x);   // only for doxygen
00582 
00593 template<typename T>inline T copysign(T x,T y);   // only for doxygen
00594 
00602 template<typename T>inline T max(T a,T b);  // only for doxygen
00603 
00611 template<typename T>inline T min(T a,T b);  // only for doxygen
00612 
00614 template<typename T>inline T asin(T x);  // only for doxygen
00616 template<typename T>inline T acos(T x);  // only for doxygen
00618 template<typename T>inline T tan(T x);  // only for doxygen
00619 
00621 template<typename T>inline T sinh(T x);  // only for doxygen
00623 template<typename T>inline T cosh(T x);  // only for doxygen
00625 template<typename T>inline T tanh(T x);  // only for doxygen
00626 
00628 template<typename T>inline T asinh(T x);  // only for doxygen
00630 template<typename T>inline T acosh(T x);  // only for doxygen
00632 template<typename T>inline T atanh(T x);  // only for doxygen
00633 
00642 template<typename T>inline T sin(T x);  // only for doxygen
00643 
00652 template<typename T>inline T cos(T x);  // only for doxygen
00653 
00660 template<typename T>inline T atan(T x);  // only for doxygen
00661 
00669 template<typename T>inline T atan2(T x);  // only for doxygen
00670 
00672 template<typename T>inline T sqr(T x);   // only for doxygen
00674 template<typename T>inline T sqrt(T x);  // only for doxygen
00676 template<typename T>inline T cbrt(T x);  // only for doxygen
00677 
00683 template<typename T>inline T pow(T x,T y);  // only for doxygen
00684 
00691 template<typename T>inline T exp(T x);  // only for doxygen
00692 
00698 template<typename T>inline T log(T x);  // only for doxygen
00699 
00706 template<typename T>inline T expm1(T x);  // only for doxygen
00707 
00713 template<typename T>inline T log1p(T x);  // only for doxygen
00714 
00716 template<typename T>inline T exp10(T x);  // only for doxygen
00718 template<typename T>inline T log10(T x);  // only for doxygen
00719 
00727 template<typename T>inline T fdim(T x,T y);  // only for doxygen
00728 
00734 template<typename T>inline T exp2(T x);  // only for doxygen
00735 
00741 template<typename T>inline T log2(T x);  // only for doxygen
00742 
00750 template<typename T>inline int ilog2(T x);  // only for doxygen
00751 
00760 template<typename T>inline T ldexp(T x,int exp);  // only for doxygen
00761 
00771 template<typename T>inline T modf(T x,T *iptr);  // only for doxygen
00772 
00784 template<typename T>inline T frexp(T x,int *exp);  // only for doxygen
00785 
00787 template<typename T>inline void sincos(T x,T *sinx,T *cosx); // only for doxygen
00788 
00798 template<typename T>inline T remainder(T x,T y);  // only for doxygen
00799 
00805 template<typename T>inline T fma(T x,T y,T z);  // only for doxygen
00806 
00808 template<typename T>inline T ceil(T x);  // only for doxygen
00810 template<typename T>inline T floor(T x);  // only for doxygen
00812 template<typename T>inline T trunc(T x);  // only for doxygen
00813 
00828 template<typename I,typename FP>inline I iroundc(FP x);  // only for doxygen
00829 
00841 template<typename I,typename FP>inline I iroundc(I *d,FP x); // only for doxygen
00842 
00857 template<typename I,typename FP>inline I iround(FP x);  // only for doxygen
00858 
00870 template<typename I,typename FP>inline I iround(I *d,FP x);  // only for doxygen
00871 
00872 
00873 #endif  /* DOXYGEN_IS_SCANNING */
00874 
00882 template<typename T>inline T dabs(T x,T y)
00883     {  return(NUM::fabs(x-y));  }
00884 
00885 //==============================================================================
00886 // ALL BASIC FUNCTIONS ARE ABOVE THIS LINE.
00887 //==============================================================================
00888 
00889 #ifndef DOXYGEN_IS_SCANNING
00890 
00891 // Now, as all basic functions are provided, we can undef some of the ugly 
00892 // defines. 
00893 #undef signbit  /* Use is_negative() instead. */
00894 extern _deprecated_ int signbit(dbl x);  // Do not use; use is_negative(). 
00895 
00896 #undef isinf  /* Use is_inf() instead. */
00897 extern _deprecated_ int isinf(dbl x);  // Do not use; use is_inf(). 
00898 
00899 #undef isnan  /* Use is_nan() instead. */
00900 extern _deprecated_ int isnan(dbl x);  // Do not use; use is_nan(). 
00901 
00902 #undef finite    /* Use is_finite() instead. */
00903 #undef isfinite  /* Use is_finite() instead. */
00904 extern _deprecated_ int finite(dbl x);    // Do not use; use is_finite(). 
00905 extern _deprecated_ int isfinite(dbl x);  // Do not use; use is_finite(). 
00906 
00907 // Do not use fpclassify directly. Use is_inf(), is_nan(), is_finite(). 
00908 #undef fpclassify
00909 extern _deprecated_ int fpclassify(dbl x);
00910 
00911 // Do not use significand() unless we really need it; use frexp() instead 
00912 // which also eliminates the corresponding call to ilog2(). 
00913 #undef significand
00914 extern _deprecated_ dbl significand(dbl x);
00915 
00916 #endif  /* !defined(DOXYGEN_IS_SCANNING) */
00917 
00918 //==============================================================================
00919 // ALL HIGHER-LEVEL FUNCTIONS ARE BELOW THIS LINE.
00920 //==============================================================================
00921 
00923 template<typename T>inline T max(T a,T b,T c)
00924     {  return(b<a ? NUM::max(a,c) : NUM::max(b,c));  }
00926 template<typename T>inline T min(T a,T b,T c)
00927     {  return(a<b ? NUM::min(a,c) : NUM::min(b,c));  }
00928 
00930 inline _constfn_ int powii(int x,int y)
00931 {
00932     if(y<0) return(-1);
00933     int res=1;
00934     for(;y;x*=x)  {  if(y&1) res*=x;  y/=2;  }
00935     return(res);
00936 }
00938 template<typename T>inline _constfn_ T powi(T x,int y)
00939 {
00940     #if 0
00941         bool neg=(y<0 ? (y=-y,1) : 0);
00942         T res=1;
00943         for(;y;x*=x)  {  if(y&1) res*=x;  y/=2;  }
00944         return(neg ? 1/res : res);
00945     #else
00946         // This is repeated squaring method, based on GSL implementation. 
00947         // Returns 0.0^0 = 1.0, so continuous in x
00948         if(y<0)  {  x=1/x;  y=-y;  }
00949         T res=1;
00950         do { if(y&1) res*=x;  x*=x; } while((y/=2));
00951         return(res);
00952     #endif
00953 }
00954 
00955 
00962 template<typename T>inline _constfn_ T clamp2(T x,T low,T high)
00963     {  return(x<low ? low : (x>high ? high : x));  }
00964 
00970 template<typename T>inline T clamp1(T x,T b)
00971     {  return(NUM::fabs(x)>b ? (x>b ? b : -b) : x);  }
00972 
00980 template<typename T>inline T hypot2(T a,T b)
00981 {
00982     a=NUM::fabs(a); b=NUM::fabs(b);
00983     return(a>b ? (a * sqrt(1 + sqr(b/a))) : 
00984         (b==0 ? 0 : b * sqrt(1 + sqr(a/b))) );
00985 }
00986 
00987 
01001 template<typename T,typename P>inline T lerp2(T l,T r,P p,P p1)
01002     {  return( p*r + p1*l );  }
01003 
01013 template<typename T,typename P>inline T lerp(T l,T r,P p)
01014     {  return(NUM::lerp2<T,P>(l,r,p,1-p));  }
01015 
01016 
01036 template<typename T,typename P>inline T bilerp2(
01037     T lt,T rt,T lb,T rb,
01038     P px,P py,P px1,P py1)
01039     {  return( py * (px*rb + px1*lb) + py1 * (px*rt + px1*lt) );  }
01040 
01050 template<typename T,typename P>inline T bilerp(
01051     T lt,T rt,T lb,T rb,
01052     P px,P py)
01053     {  return(NUM::bilerp2<T,P>(lt,rt,lb,rb,px,py,1-px,1-py));  }
01054 
01055 
01064 template<typename T>inline T FermiDirac(T x)
01065     {  return(1/(1+NUM::exp(-x)));  }
01066 
01085 template<typename T>inline T FermiDirac(T x,T t)
01086     {  return(NUM::FermiDirac(x*t));  }
01087 
01088 
01099 template<typename T>inline T StepLinear(T x,T a,T b)
01100     {  return( x<=a ? 0 : (x>=b ? 1 : (x-a)/(b-a)) );  }
01101 
01112 template<typename T>inline T StepCubic(T x,T a,T b)
01113     {  T v=NUM::StepLinear(x,a,b);  return((3-2*v)*v*v);  }
01114 
01115 
01125 template<typename T>inline T NormEuclidic2(const T *x,int dim)
01126 {
01127     T sum=0;
01128     for(int i=0; i<dim; i++)  sum+=NUM::sqr(x[i]);
01129     return(sum);
01130 }
01131 
01141 template<typename T>inline T NormEuclidic(const T *x,int dim)
01142 {
01143     if(dim==1)  return(NUM::fabs(*x));
01144     T sum=0;
01145     for(int i=0; i<dim; i++)  sum+=NUM::sqr(x[i]);
01146     return(sqrt(sum));
01147 }
01148 
01156 template<typename T>inline T DistEuclidic2(const T *x,const T *y,int dim)
01157 {
01158     T sum=0;
01159     for(int i=0; i<dim; i++)  sum+=NUM::sqr(y[i]-x[i]);
01160     return(sum);
01161 }
01162 
01170 template<typename T>inline T DistEuclidic(const T *x,const T *y,int dim)
01171 {
01172     if(dim==1)  return(NUM::fabs(*y-*x));
01173     T sum=0;
01174     for(int i=0; i<dim; i++)  sum+=sqr(y[i]-x[i]);
01175     return(sqrt(sum));
01176 }
01177 
01178 
01186 template<typename T>inline T FindMaximum(const T *a,int dim)
01187 {
01188     if(!dim) return(0);
01189     T max=a[0];
01190     for(int i=1; i<dim; i++)
01191         if(max<a[i]) max=a[i];
01192     return(max);
01193 }
01194 
01200 template<class T>inline T FindMinimum(const T *a,int dim)
01201 {
01202     if(!dim) return(0);
01203     T min=a[0];
01204     for(int i=1; i<dim; i++)
01205         if(min>a[i]) min=a[i];
01206     return(min);
01207 }
01208 
01216 template<class T>inline int FindMinimumIdx(const T *a,int dim)
01217 {
01218     if(!dim) return(-1);
01219     int midx=0;
01220     T min=a[0];
01221     for(int i=1; i<dim; i++)
01222         if(min>a[i])
01223         {  min=a[i];  midx=i;  }
01224     return(midx);
01225 }
01226 
01227 
01239 template<typename T>inline T NormalizeVector(T *x,int dim)
01240 {
01241     T len=NormEuclidic(x,dim),fact=T(1)/len;
01242     for(int i=0; i<dim; i++)
01243     {  x[i]*=fact;  }
01244     return(len);
01245 }
01246 
01256 template<typename T>inline T NormalizeVector(T *dest,const T *src,int dim)
01257 {
01258     T len=NormEuclidic(src,dim),fact=T(1)/len;
01259     for(int i=0; i<dim; i++)
01260     {  dest[i]=src[i]*fact;  }
01261     return(len);
01262 }
01263 
01265 template<typename T>inline T NormalizeVector2D(T *x)
01266 { T len=hypot(x[0],x[1]),fact=T(1)/len;  x[0]*=fact; x[1]*=fact; return(len); }
01267 
01269 template<typename T>inline void CrossProduct(T *dest,const T *a,const T *b)
01270 {
01271     dest[0] = a[1]*b[2] - a[2]*b[1];
01272     dest[1] = a[2]*b[0] - a[0]*b[2];
01273     dest[2] = a[0]*b[1] - a[1]*b[0];
01274 }
01275 
01276 
01283 template<typename T>inline T DotProduct(const T *a,const T*b,int dim)
01284 {
01285     T res=0;
01286     for(int i=0; i<dim; i++)  res+=a[i]*b[i];
01287     return(res);
01288 }
01290 template<typename T>inline T DotProduct3D(const T *a,const T*b)
01291     {  return(a[0]*b[0]+a[1]*b[1]+a[2]*b[2]);  }
01292 
01293 
01304 template<typename T>inline T PolvEval(T x,const T *c,int n)
01305 {
01306     T r=c[n];
01307     while(n--)  r=r*x+c[n];
01308     return(r);
01309 }
01310 
01311 
01326 template<typename T>inline int ArrayBinsearchInterval(T *array,int size,T value)
01327 {
01328     if(!size)  return(0);
01329     // Prevent us from crossing the array borders below: 
01330     if(value<array[0])  return(-1);
01331     if(value>=array[size-1])  return(size);
01332 Assert(size>=2);
01333     int a=0,b=size-1;
01334     while(b-a>1)
01335     {
01336         int m=(a+b)/2;
01337         if(value>=array[m])  a=m;
01338         else  b=m;
01339     }
01340 Assert(array[a]<=value);
01341 Assert(value<array[b]);
01342     return(a);
01343 }
01344 
01345 
01357 template<typename T>inline int TestOrderAscenting(const T *array,int n,
01358     T val_eps=1e-7)
01359 {
01360     if(!array) return(0);
01361     for(int i=1; i<n; i++)
01362         if(array[i]-array[i-1]<=val_eps)
01363             return(i);
01364     return(0);
01365 }
01366 
01367 }  // end of namespace NUM
01368 
01369 #endif  /* _LIB_NUMERICS_MATH_H_ */
/ray/src/lib/numerics/num_math.h
