src/gmxlib/maths.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * GROningen Mixture of Alchemy and Childrens' Stories
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39
  40 #include <math.h>
  41 #include <limits.h>
  42 #include "maths.h"
  43 #ifdef HAVE__FINITE
  44 #include "float.h"
  45 #endif
  46
  47 int gmx_nint(real a)
  48 {
  49   const real half = .5;
  50   int   result;
  51
  52   result = (a < 0.) ? ((int)(a - half)) : ((int)(a + half));
  53   return result;
  54 }
  55
  56 real cuberoot (real x)
  57 {
  58   if (x < 0)
  59     {
  60       return (-pow(-x,1.0/DIM));
  61     }
  62   else
  63     {
  64       return (pow(x,1.0/DIM));
  65     }
  66 }
  67
  68 real sign(real x,real y)
  69 {
  70   if (y < 0)
  71     return -fabs(x);
  72   else
  73     return +fabs(x);
  74 }
  75
  76 /* Double and single precision erf() and erfc() from
  77  * the Sun Freely Distributable Math Library FDLIBM.
  78  * See http://www.netlib.org/fdlibm
  79  * Specific file used: s_erf.c, version 1.3 95/01/18
  80  */
  81 /*
  82  * ====================================================
  83  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  84  *
  85  * Developed at SunSoft, a Sun Microsystems, Inc. business.
  86  * Permission to use, copy, modify, and distribute this
  87  * software is freely granted, provided that this notice
  88  * is preserved.
  89  * ====================================================
  90  */
  91
  92 #if ( (defined SIZEOF_INT && SIZEOF_INT==4) || (SIZEOF_INT_MAX == 2147483647) )
  93    typedef int erf_int32_t;
  94    typedef unsigned int erf_u_int32_t;
  95 #elif (LONG_MAX == 2147483647L)
  96    typedef long erf_int32_t;
  97    typedef unsigned long erf_u_int32_t;
  98 #elif (SHRT_MAX == 2147483647)
  99    typedef short erf_int32_t;
 100    typedef unsigned short erf_u_int32_t;
 101 #else
 102 #  error ERROR: No 32 bit wide integer type found!
 103 #endif
 104
 105
 106 #ifdef GMX_DOUBLE
 107
 108 static const double
 109 tiny        = 1e-300,
 110 half=  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
 111 one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
 112 two =  2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
 113 /* c = (float)0.84506291151 */
 114 erx =  8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
 115 /*
 116  * Coefficients for approximation to  erf on [0,0.84375]
 117  */
 118 efx =  1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
 119 efx8=  1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
 120 pp0  =  1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
 121 pp1  = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
 122 pp2  = -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
 123 pp3  = -5.77027029648944159157e-03, /* 0xBF77A291, 0x236668E4 */
 124 pp4  = -2.37630166566501626084e-05, /* 0xBEF8EAD6, 0x120016AC */
 125 qq1  =  3.97917223959155352819e-01, /* 0x3FD97779, 0xCDDADC09 */
 126 qq2  =  6.50222499887672944485e-02, /* 0x3FB0A54C, 0x5536CEBA */
 127 qq3  =  5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */
 128 qq4  =  1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */
 129 qq5  = -3.96022827877536812320e-06, /* 0xBED09C43, 0x42A26120 */
 130 /*
 131  * Coefficients for approximation to  erf  in [0.84375,1.25]
 132  */
 133 pa0  = -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */
 134 pa1  =  4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
 135 pa2  = -3.72207876035701323847e-01, /* 0xBFD7D240, 0xFBB8C3F1 */
 136 pa3  =  3.18346619901161753674e-01, /* 0x3FD45FCA, 0x805120E4 */
 137 pa4  = -1.10894694282396677476e-01, /* 0xBFBC6398, 0x3D3E28EC */
 138 pa5  =  3.54783043256182359371e-02, /* 0x3FA22A36, 0x599795EB */
 139 pa6  = -2.16637559486879084300e-03, /* 0xBF61BF38, 0x0A96073F */
 140 qa1  =  1.06420880400844228286e-01, /* 0x3FBB3E66, 0x18EEE323 */
 141 qa2  =  5.40397917702171048937e-01, /* 0x3FE14AF0, 0x92EB6F33 */
 142 qa3  =  7.18286544141962662868e-02, /* 0x3FB2635C, 0xD99FE9A7 */
 143 qa4  =  1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */
 144 qa5  =  1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */
 145 qa6  =  1.19844998467991074170e-02, /* 0x3F888B54, 0x5735151D */
 146 /*
 147  * Coefficients for approximation to  erfc in [1.25,1/0.35]
 148  */
 149 ra0  = -9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */
 150 ra1  = -6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */
 151 ra2  = -1.05586262253232909814e+01, /* 0xC0251E04, 0x41B0E726 */
 152 ra3  = -6.23753324503260060396e+01, /* 0xC04F300A, 0xE4CBA38D */
 153 ra4  = -1.62396669462573470355e+02, /* 0xC0644CB1, 0x84282266 */
 154 ra5  = -1.84605092906711035994e+02, /* 0xC067135C, 0xEBCCABB2 */
 155 ra6  = -8.12874355063065934246e+01, /* 0xC0545265, 0x57E4D2F2 */
 156 ra7  = -9.81432934416914548592e+00, /* 0xC023A0EF, 0xC69AC25C */
 157 sa1  =  1.96512716674392571292e+01, /* 0x4033A6B9, 0xBD707687 */
 158 sa2  =  1.37657754143519042600e+02, /* 0x4061350C, 0x526AE721 */
 159 sa3  =  4.34565877475229228821e+02, /* 0x407B290D, 0xD58A1A71 */
 160 sa4  =  6.45387271733267880336e+02, /* 0x40842B19, 0x21EC2868 */
 161 sa5  =  4.29008140027567833386e+02, /* 0x407AD021, 0x57700314 */
 162 sa6  =  1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */
 163 sa7  =  6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */
 164 sa8  = -6.04244152148580987438e-02, /* 0xBFAEEFF2, 0xEE749A62 */
 165 /*
 166  * Coefficients for approximation to  erfc in [1/.35,28]
 167  */
 168 rb0  = -9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */
 169 rb1  = -7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */
 170 rb2  = -1.77579549177547519889e+01, /* 0xC031C209, 0x555F995A */
 171 rb3  = -1.60636384855821916062e+02, /* 0xC064145D, 0x43C5ED98 */
 172 rb4  = -6.37566443368389627722e+02, /* 0xC083EC88, 0x1375F228 */
 173 rb5  = -1.02509513161107724954e+03, /* 0xC0900461, 0x6A2E5992 */
 174 rb6  = -4.83519191608651397019e+02, /* 0xC07E384E, 0x9BDC383F */
 175 sb1  =  3.03380607434824582924e+01, /* 0x403E568B, 0x261D5190 */
 176 sb2  =  3.25792512996573918826e+02, /* 0x40745CAE, 0x221B9F0A */
 177 sb3  =  1.53672958608443695994e+03, /* 0x409802EB, 0x189D5118 */
 178 sb4  =  3.19985821950859553908e+03, /* 0x40A8FFB7, 0x688C246A */
 179 sb5  =  2.55305040643316442583e+03, /* 0x40A3F219, 0xCEDF3BE6 */
 180 sb6  =  4.74528541206955367215e+02, /* 0x407DA874, 0xE79FE763 */
 181 sb7  = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
 182
 183 double gmx_erf(double x)
 184 {
 185
 186         erf_int32_t hx,ix,i;
 187     double R,S,P,Q,s,y,z,r;
 188
 189     union
 190     {
 191         double d;
 192         int    i[2];
 193     }
 194     conv;
 195
 196     conv.d=x;
 197
 198         /* In release-4-6 and later branches, only the test for
 199          * GMX_IEEE754_BIG_ENDIAN_WORD_ORDER will be required. */
 200 #if defined(IEEE754_BIG_ENDIAN_WORD_ORDER) || defined(GMX_IEEE754_BIG_ENDIAN_WORD_ORDER)
 201     hx=conv.i[0];
 202 #else
 203     hx=conv.i[1];
 204 #endif
 205
 206         ix = hx&0x7fffffff;
 207         if(ix>=0x7ff00000)
 208     {
 209                 /* erf(nan)=nan */
 210             i = ((erf_u_int32_t)hx>>31)<<1;
 211             return (double)(1-i)+one/x; /* erf(+-inf)=+-1 */
 212         }
 213
 214         if(ix < 0x3feb0000)
 215     {
 216                 /* |x|<0.84375 */
 217             if(ix < 0x3e300000)
 218         {
 219             /* |x|<2**-28 */
 220                 if (ix < 0x00800000)
 221                 return 0.125*(8.0*x+efx8*x);  /*avoid underflow */
 222             return x + efx*x;
 223             }
 224             z = x*x;
 225             r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
 226             s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
 227             y = r/s;
 228             return x + x*y;
 229         }
 230         if(ix < 0x3ff40000)
 231     {
 232                 /* 0.84375 <= |x| < 1.25 */
 233             s = fabs(x)-one;
 234             P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
 235             Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
 236             if(hx>=0) return erx + P/Q; else return -erx - P/Q;
 237         }
 238         if (ix >= 0x40180000)
 239     {
 240         /* inf>|x|>=6 */
 241             if(hx>=0) return one-tiny; else return tiny-one;
 242         }
 243         x = fabs(x);
 244         s = one/(x*x);
 245         if(ix< 0x4006DB6E)
 246     {
 247         /* |x| < 1/0.35 */
 248             R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))));
 249             S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+s*sa8)))))));
 250         }
 251     else
 252     {
 253         /* |x| >= 1/0.35 */
 254             R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))));
 255             S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
 256         }
 257
 258     conv.d = x;
 259
 260         /* In release-4-6 and later branches, only the test for
 261          * GMX_IEEE754_BIG_ENDIAN_WORD_ORDER will be required. */
 262 #if defined(IEEE754_BIG_ENDIAN_WORD_ORDER) || defined(GMX_IEEE754_BIG_ENDIAN_WORD_ORDER)
 263     conv.i[1] = 0;
 264 #else
 265     conv.i[0] = 0;
 266 #endif
 267
 268     z = conv.d;
 269
 270         r  =  exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
 271         if(hx>=0)
 272         return one-r/x;
 273     else
 274         return  r/x-one;
 275 }
 276
 277
 278 double gmx_erfc(double x)
 279 {
 280         erf_int32_t hx,ix;
 281         double R,S,P,Q,s,y,z,r;
 282
 283     union
 284     {
 285         double d;
 286         int    i[2];
 287     }
 288     conv;
 289
 290     conv.d = x;
 291
 292         /* In release-4-6 and later branches, only the test for
 293          * GMX_IEEE754_BIG_ENDIAN_WORD_ORDER will be required. */
 294 #if defined(IEEE754_BIG_ENDIAN_WORD_ORDER) || defined(GMX_IEEE754_BIG_ENDIAN_WORD_ORDER)
 295     hx=conv.i[0];
 296 #else
 297     hx=conv.i[1];
 298 #endif
 299
 300         ix = hx&0x7fffffff;
 301         if(ix>=0x7ff00000)
 302     {
 303         /* erfc(nan)=nan */
 304         /* erfc(+-inf)=0,2 */
 305             return (double)(((erf_u_int32_t)hx>>31)<<1)+one/x;
 306         }
 307
 308         if(ix < 0x3feb0000)
 309     {
 310                 /* |x|<0.84375 */
 311             double r1,r2,s1,s2,s3,z2,z4;
 312             if(ix < 0x3c700000)         /* |x|<2**-56 */
 313             return one-x;
 314             z = x*x;
 315             r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
 316             s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
 317             y = r/s;
 318             if(hx < 0x3fd00000)
 319         {
 320             /* x<1/4 */
 321             return one-(x+x*y);
 322             }
 323         else
 324         {
 325             r = x*y;
 326             r += (x-half);
 327                 return half - r ;
 328             }
 329         }
 330
 331     if(ix < 0x3ff40000)
 332     {
 333                 /* 0.84375 <= |x| < 1.25 */
 334         s = fabs(x)-one;
 335         P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
 336         Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
 337         if(hx>=0) {
 338             z  = one-erx; return z - P/Q;
 339         } else {
 340             z = erx+P/Q; return one+z;
 341         }
 342         }
 343         if (ix < 0x403c0000)
 344     {
 345         /* |x|<28 */
 346             x = fabs(x);
 347             s = one/(x*x);
 348             if(ix< 0x4006DB6D)
 349         {
 350             /* |x| < 1/.35 ~ 2.857143*/
 351                 R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))));
 352                 S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+s*sa8)))))));
 353             }
 354         else
 355         {
 356             /* |x| >= 1/.35 ~ 2.857143 */
 357             if(hx<0&&ix>=0x40180000)
 358                 return two-tiny; /* x < -6 */
 359                 R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))));
 360                 S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
 361             }
 362
 363         conv.d = x;
 364
 365         /* In release-4-6 and later branches, only the test for
 366          * GMX_IEEE754_BIG_ENDIAN_WORD_ORDER will be required. */
 367 #if defined(IEEE754_BIG_ENDIAN_WORD_ORDER) || defined(GMX_IEEE754_BIG_ENDIAN_WORD_ORDER)
 368         conv.i[1] = 0;
 369 #else
 370         conv.i[0] = 0;
 371 #endif
 372
 373         z = conv.d;
 374
 375         r  =  exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
 376
 377         if(hx>0)
 378             return r/x;
 379         else
 380             return two-r/x;
 381     }
 382     else
 383     {
 384             if(hx>0)
 385             return tiny*tiny;
 386         else
 387             return two-tiny;
 388         }
 389 }
 390
 391 #else /* single precision */
 392
 393
 394
 395 static const float
 396 tiny        = 1e-30,
 397 half=  5.0000000000e-01, /* 0x3F000000 */
 398 one =  1.0000000000e+00, /* 0x3F800000 */
 399 two =  2.0000000000e+00, /* 0x40000000 */
 400         /* c = (subfloat)0.84506291151 */
 401 erx =  8.4506291151e-01, /* 0x3f58560b */
 402 /*
 403  * Coefficients for approximation to  erf on [0,0.84375]
 404  */
 405 efx =  1.2837916613e-01, /* 0x3e0375d4 */
 406 efx8=  1.0270333290e+00, /* 0x3f8375d4 */
 407 pp0  =  1.2837916613e-01, /* 0x3e0375d4 */
 408 pp1  = -3.2504209876e-01, /* 0xbea66beb */
 409 pp2  = -2.8481749818e-02, /* 0xbce9528f */
 410 pp3  = -5.7702702470e-03, /* 0xbbbd1489 */
 411 pp4  = -2.3763017452e-05, /* 0xb7c756b1 */
 412 qq1  =  3.9791721106e-01, /* 0x3ecbbbce */
 413 qq2  =  6.5022252500e-02, /* 0x3d852a63 */
 414 qq3  =  5.0813062117e-03, /* 0x3ba68116 */
 415 qq4  =  1.3249473704e-04, /* 0x390aee49 */
 416 qq5  = -3.9602282413e-06, /* 0xb684e21a */
 417 /*
 418  * Coefficients for approximation to  erf  in [0.84375,1.25]
 419  */
 420 pa0  = -2.3621185683e-03, /* 0xbb1acdc6 */
 421 pa1  =  4.1485610604e-01, /* 0x3ed46805 */
 422 pa2  = -3.7220788002e-01, /* 0xbebe9208 */
 423 pa3  =  3.1834661961e-01, /* 0x3ea2fe54 */
 424 pa4  = -1.1089469492e-01, /* 0xbde31cc2 */
 425 pa5  =  3.5478305072e-02, /* 0x3d1151b3 */
 426 pa6  = -2.1663755178e-03, /* 0xbb0df9c0 */
 427 qa1  =  1.0642088205e-01, /* 0x3dd9f331 */
 428 qa2  =  5.4039794207e-01, /* 0x3f0a5785 */
 429 qa3  =  7.1828655899e-02, /* 0x3d931ae7 */
 430 qa4  =  1.2617121637e-01, /* 0x3e013307 */
 431 qa5  =  1.3637083583e-02, /* 0x3c5f6e13 */
 432 qa6  =  1.1984500103e-02, /* 0x3c445aa3 */
 433 /*
 434  * Coefficients for approximation to  erfc in [1.25,1/0.35]
 435  */
 436 ra0  = -9.8649440333e-03, /* 0xbc21a093 */
 437 ra1  = -6.9385856390e-01, /* 0xbf31a0b7 */
 438 ra2  = -1.0558626175e+01, /* 0xc128f022 */
 439 ra3  = -6.2375331879e+01, /* 0xc2798057 */
 440 ra4  = -1.6239666748e+02, /* 0xc322658c */
 441 ra5  = -1.8460508728e+02, /* 0xc3389ae7 */
 442 ra6  = -8.1287437439e+01, /* 0xc2a2932b */
 443 ra7  = -9.8143291473e+00, /* 0xc11d077e */
 444 sa1  =  1.9651271820e+01, /* 0x419d35ce */
 445 sa2  =  1.3765776062e+02, /* 0x4309a863 */
 446 sa3  =  4.3456588745e+02, /* 0x43d9486f */
 447 sa4  =  6.4538726807e+02, /* 0x442158c9 */
 448 sa5  =  4.2900814819e+02, /* 0x43d6810b */
 449 sa6  =  1.0863500214e+02, /* 0x42d9451f */
 450 sa7  =  6.5702495575e+00, /* 0x40d23f7c */
 451 sa8  = -6.0424413532e-02, /* 0xbd777f97 */
 452 /*
 453  * Coefficients for approximation to  erfc in [1/.35,28]
 454  */
 455 rb0  = -9.8649431020e-03, /* 0xbc21a092 */
 456 rb1  = -7.9928326607e-01, /* 0xbf4c9dd4 */
 457 rb2  = -1.7757955551e+01, /* 0xc18e104b */
 458 rb3  = -1.6063638306e+02, /* 0xc320a2ea */
 459 rb4  = -6.3756646729e+02, /* 0xc41f6441 */
 460 rb5  = -1.0250950928e+03, /* 0xc480230b */
 461 rb6  = -4.8351919556e+02, /* 0xc3f1c275 */
 462 sb1  =  3.0338060379e+01, /* 0x41f2b459 */
 463 sb2  =  3.2579251099e+02, /* 0x43a2e571 */
 464 sb3  =  1.5367296143e+03, /* 0x44c01759 */
 465 sb4  =  3.1998581543e+03, /* 0x4547fdbb */
 466 sb5  =  2.5530502930e+03, /* 0x451f90ce */
 467 sb6  =  4.7452853394e+02, /* 0x43ed43a7 */
 468 sb7  = -2.2440952301e+01; /* 0xc1b38712 */
 469
 470
 471 typedef union
 472 {
 473   float value;
 474   erf_u_int32_t word;
 475 } ieee_float_shape_type;
 476
 477 #define GET_FLOAT_WORD(i,d)                                     \
 478 do {                                                            \
 479   ieee_float_shape_type gf_u;                                   \
 480   gf_u.value = (d);                                             \
 481   (i) = gf_u.word;                                              \
 482 } while (0)
 483
 484
 485 #define SET_FLOAT_WORD(d,i)                                     \
 486 do {                                                            \
 487   ieee_float_shape_type sf_u;                                   \
 488   sf_u.word = (i);                                              \
 489   (d) = sf_u.value;                                             \
 490 } while (0)
 491
 492
 493 float gmx_erf(float x)
 494 {
 495         erf_int32_t hx,ix,i;
 496         float R,S,P,Q,s,y,z,r;
 497
 498     union
 499     {
 500         float  f;
 501         int    i;
 502     }
 503     conv;
 504
 505     conv.f=x;
 506     hx=conv.i;
 507
 508         ix = hx&0x7fffffff;
 509         if(ix>=0x7f800000)
 510     {
 511                 /* erf(nan)=nan */
 512             i = ((erf_u_int32_t)hx>>31)<<1;
 513             return (float)(1-i)+one/x;  /* erf(+-inf)=+-1 */
 514         }
 515
 516         if(ix < 0x3f580000)
 517     {
 518                 /* |x|<0.84375 */
 519             if(ix < 0x31800000)
 520         {
 521             /* |x|<2**-28 */
 522                 if (ix < 0x04000000)
 523                 return (float)0.125*((float)8.0*x+efx8*x);                  /*avoid underflow */
 524             return x + efx*x;
 525             }
 526         z = x*x;
 527         r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
 528             s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
 529             y = r/s;
 530             return x + x*y;
 531         }
 532         if(ix < 0x3fa00000)
 533     {
 534         /* 0.84375 <= |x| < 1.25 */
 535             s = fabs(x)-one;
 536             P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
 537             Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
 538             if(hx>=0) return erx + P/Q; else return -erx - P/Q;
 539         }
 540     if (ix >= 0x40c00000)
 541     {
 542                 /* inf>|x|>=6 */
 543         if(hx>=0) return one-tiny; else return tiny-one;
 544         }
 545         x = fabs(x);
 546         s = one/(x*x);
 547     if(ix< 0x4036DB6E)
 548     {
 549         /* |x| < 1/0.35 */
 550             R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))));
 551             S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+s*sa8)))))));
 552         }
 553     else
 554     {
 555         /* |x| >= 1/0.35 */
 556             R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))));
 557             S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
 558         }
 559
 560     conv.f = x;
 561     conv.i = conv.i & 0xfffff000;
 562     z = conv.f;
 563
 564         r  =  exp(-z*z-(float)0.5625)*exp((z-x)*(z+x)+R/S);
 565         if(hx>=0) return one-r/x; else return  r/x-one;
 566 }
 567
 568 float gmx_erfc(float x)
 569 {
 570         erf_int32_t hx,ix;
 571         float R,S,P,Q,s,y,z,r;
 572
 573     union
 574     {
 575         float  f;
 576         int    i;
 577     }
 578     conv;
 579
 580     conv.f=x;
 581     hx=conv.i;
 582
 583         ix = hx&0x7fffffff;
 584         if(ix>=0x7f800000)
 585     {
 586         /* erfc(nan)=nan */
 587         /* erfc(+-inf)=0,2 */
 588             return (float)(((erf_u_int32_t)hx>>31)<<1)+one/x;
 589         }
 590
 591         if(ix < 0x3f580000)
 592     {
 593                 /* |x|<0.84375 */
 594             if(ix < 0x23800000)
 595             return one-x;       /* |x|<2**-56 */
 596             z = x*x;
 597             r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
 598             s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
 599             y = r/s;
 600             if(hx < 0x3e800000)
 601         {
 602             /* x<1/4 */
 603             return one-(x+x*y);
 604             } else {
 605             r = x*y;
 606             r += (x-half);
 607                 return half - r ;
 608             }
 609         }
 610         if(ix < 0x3fa00000)
 611     {
 612         /* 0.84375 <= |x| < 1.25 */
 613             s = fabs(x)-one;
 614             P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
 615             Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
 616             if(hx>=0) {
 617                 z  = one-erx; return z - P/Q;
 618             } else {
 619             z = erx+P/Q; return one+z;
 620             }
 621         }
 622         if (ix < 0x41e00000)
 623     {
 624                 /* |x|<28 */
 625         x = fabs(x);
 626             s = one/(x*x);
 627             if(ix< 0x4036DB6D)
 628         {
 629             /* |x| < 1/.35 ~ 2.857143*/
 630             R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))));
 631                 S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+s*sa8)))))));
 632             } else {
 633             /* |x| >= 1/.35 ~ 2.857143 */
 634             if(hx<0&&ix>=0x40c00000) return two-tiny;/* x < -6 */
 635                 R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))));
 636                 S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
 637             }
 638
 639         conv.f = x;
 640         conv.i = conv.i & 0xfffff000;
 641         z = conv.f;
 642
 643             r  =  exp(-z*z-(float)0.5625)*exp((z-x)*(z+x)+R/S);
 644             if(hx>0) return r/x; else return two-r/x;
 645         } else {
 646             if(hx>0) return tiny*tiny; else return two-tiny;
 647         }
 648 }
 649
 650 #endif
 651
 652 gmx_bool gmx_isfinite(real x)
 653 {
 654     gmx_bool returnval = TRUE;
 655     /* If no suitable function was found, assume the value is
 656      * finite. */
 657
 658 #ifdef HAVE_ISFINITE
 659     returnval = isfinite(x);
 660 #elif defined HAVE__ISFINITE
 661     returnval = _isfinite(x);
 662 #elif defined HAVE__FINITE
 663     returnval = _finite(x);
 664 #endif
 665     return returnval;
 666 }
 667
 668 gmx_bool
 669 check_int_multiply_for_overflow(gmx_large_int_t a,
 670                                 gmx_large_int_t b,
 671                                 gmx_large_int_t *result)
 672 {
 673     gmx_large_int_t sign = 1;
 674     if((0 == a) || (0 == b))
 675     {
 676         *result = 0;
 677         return TRUE;
 678     }
 679     if(a < 0)
 680     {
 681         a = -a;
 682         sign = -sign;
 683     }
 684     if(b < 0)
 685     {
 686         b = -b;
 687         sign = -sign;
 688     }
 689     if(GMX_LARGE_INT_MAX / b < a)
 690     {
 691         *result = (sign > 0) ? GMX_LARGE_INT_MAX : GMX_LARGE_INT_MIN;
 692         return FALSE;
 693     }
 694     *result = sign * a * b;
 695     return TRUE;
 696 }