src/gromacs/math/utilities.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "utilities.h"
  40
  41 #include "config.h"
  42
  43 #include <assert.h>
  44 #include <limits.h>
  45 #include <math.h>
  46
  47 #ifdef HAVE__FINITE
  48 #include <float.h>
  49 #endif
  50 #ifndef GMX_NATIVE_WINDOWS
  51 #include <fenv.h>
  52 #endif
  53
  54 int gmx_nint(real a)
  55 {
  56     const real half = .5;
  57     int        result;
  58
  59     result = (a < 0.) ? ((int)(a - half)) : ((int)(a + half));
  60     return result;
  61 }
  62
  63 real cuberoot(real x)
  64 {
  65     if (x < 0)
  66     {
  67         return (-pow(-x, 1.0/3.0));
  68     }
  69     else
  70     {
  71         return (pow(x, 1.0/3.0));
  72     }
  73 }
  74
  75 real sign(real x, real y)
  76 {
  77     if (y < 0)
  78     {
  79         return -fabs(x);
  80     }
  81     else
  82     {
  83         return +fabs(x);
  84     }
  85 }
  86
  87 /* Double and single precision erf() and erfc() from
  88  * the Sun Freely Distributable Math Library FDLIBM.
  89  * See http://www.netlib.org/fdlibm
  90  * Specific file used: s_erf.c, version 1.3 95/01/18
  91  */
  92 /*
  93  * ====================================================
  94  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  95  *
  96  * Developed at SunSoft, a Sun Microsystems, Inc. business.
  97  * Permission to use, copy, modify, and distribute this
  98  * software is freely granted, provided that this notice
  99  * is preserved.
 100  * ====================================================
 101  */
 102
 103 static const double
 104     tiny        = 1e-300,
 105     half        =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
 106     one         =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
 107     two         =  2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
 108 /* c = (float)0.84506291151 */
 109     erx =  8.45062911510467529297e-01,         /* 0x3FEB0AC1, 0x60000000 */
 110 /*
 111  * Coefficients for approximation to  erf on [0,0.84375]
 112  */
 113     efx  =  1.28379167095512586316e-01,        /* 0x3FC06EBA, 0x8214DB69 */
 114     efx8 =  1.02703333676410069053e+00,        /* 0x3FF06EBA, 0x8214DB69 */
 115     pp0  =  1.28379167095512558561e-01,        /* 0x3FC06EBA, 0x8214DB68 */
 116     pp1  = -3.25042107247001499370e-01,        /* 0xBFD4CD7D, 0x691CB913 */
 117     pp2  = -2.84817495755985104766e-02,        /* 0xBF9D2A51, 0xDBD7194F */
 118     pp3  = -5.77027029648944159157e-03,        /* 0xBF77A291, 0x236668E4 */
 119     pp4  = -2.37630166566501626084e-05,        /* 0xBEF8EAD6, 0x120016AC */
 120     qq1  =  3.97917223959155352819e-01,        /* 0x3FD97779, 0xCDDADC09 */
 121     qq2  =  6.50222499887672944485e-02,        /* 0x3FB0A54C, 0x5536CEBA */
 122     qq3  =  5.08130628187576562776e-03,        /* 0x3F74D022, 0xC4D36B0F */
 123     qq4  =  1.32494738004321644526e-04,        /* 0x3F215DC9, 0x221C1A10 */
 124     qq5  = -3.96022827877536812320e-06,        /* 0xBED09C43, 0x42A26120 */
 125 /*
 126  * Coefficients for approximation to  erf  in [0.84375,1.25]
 127  */
 128     pa0  = -2.36211856075265944077e-03,        /* 0xBF6359B8, 0xBEF77538 */
 129     pa1  =  4.14856118683748331666e-01,        /* 0x3FDA8D00, 0xAD92B34D */
 130     pa2  = -3.72207876035701323847e-01,        /* 0xBFD7D240, 0xFBB8C3F1 */
 131     pa3  =  3.18346619901161753674e-01,        /* 0x3FD45FCA, 0x805120E4 */
 132     pa4  = -1.10894694282396677476e-01,        /* 0xBFBC6398, 0x3D3E28EC */
 133     pa5  =  3.54783043256182359371e-02,        /* 0x3FA22A36, 0x599795EB */
 134     pa6  = -2.16637559486879084300e-03,        /* 0xBF61BF38, 0x0A96073F */
 135     qa1  =  1.06420880400844228286e-01,        /* 0x3FBB3E66, 0x18EEE323 */
 136     qa2  =  5.40397917702171048937e-01,        /* 0x3FE14AF0, 0x92EB6F33 */
 137     qa3  =  7.18286544141962662868e-02,        /* 0x3FB2635C, 0xD99FE9A7 */
 138     qa4  =  1.26171219808761642112e-01,        /* 0x3FC02660, 0xE763351F */
 139     qa5  =  1.36370839120290507362e-02,        /* 0x3F8BEDC2, 0x6B51DD1C */
 140     qa6  =  1.19844998467991074170e-02,        /* 0x3F888B54, 0x5735151D */
 141 /*
 142  * Coefficients for approximation to  erfc in [1.25,1/0.35]
 143  */
 144     ra0  = -9.86494403484714822705e-03,        /* 0xBF843412, 0x600D6435 */
 145     ra1  = -6.93858572707181764372e-01,        /* 0xBFE63416, 0xE4BA7360 */
 146     ra2  = -1.05586262253232909814e+01,        /* 0xC0251E04, 0x41B0E726 */
 147     ra3  = -6.23753324503260060396e+01,        /* 0xC04F300A, 0xE4CBA38D */
 148     ra4  = -1.62396669462573470355e+02,        /* 0xC0644CB1, 0x84282266 */
 149     ra5  = -1.84605092906711035994e+02,        /* 0xC067135C, 0xEBCCABB2 */
 150     ra6  = -8.12874355063065934246e+01,        /* 0xC0545265, 0x57E4D2F2 */
 151     ra7  = -9.81432934416914548592e+00,        /* 0xC023A0EF, 0xC69AC25C */
 152     sa1  =  1.96512716674392571292e+01,        /* 0x4033A6B9, 0xBD707687 */
 153     sa2  =  1.37657754143519042600e+02,        /* 0x4061350C, 0x526AE721 */
 154     sa3  =  4.34565877475229228821e+02,        /* 0x407B290D, 0xD58A1A71 */
 155     sa4  =  6.45387271733267880336e+02,        /* 0x40842B19, 0x21EC2868 */
 156     sa5  =  4.29008140027567833386e+02,        /* 0x407AD021, 0x57700314 */
 157     sa6  =  1.08635005541779435134e+02,        /* 0x405B28A3, 0xEE48AE2C */
 158     sa7  =  6.57024977031928170135e+00,        /* 0x401A47EF, 0x8E484A93 */
 159     sa8  = -6.04244152148580987438e-02,        /* 0xBFAEEFF2, 0xEE749A62 */
 160 /*
 161  * Coefficients for approximation to  erfc in [1/.35,28]
 162  */
 163     rb0  = -9.86494292470009928597e-03,        /* 0xBF843412, 0x39E86F4A */
 164     rb1  = -7.99283237680523006574e-01,        /* 0xBFE993BA, 0x70C285DE */
 165     rb2  = -1.77579549177547519889e+01,        /* 0xC031C209, 0x555F995A */
 166     rb3  = -1.60636384855821916062e+02,        /* 0xC064145D, 0x43C5ED98 */
 167     rb4  = -6.37566443368389627722e+02,        /* 0xC083EC88, 0x1375F228 */
 168     rb5  = -1.02509513161107724954e+03,        /* 0xC0900461, 0x6A2E5992 */
 169     rb6  = -4.83519191608651397019e+02,        /* 0xC07E384E, 0x9BDC383F */
 170     sb1  =  3.03380607434824582924e+01,        /* 0x403E568B, 0x261D5190 */
 171     sb2  =  3.25792512996573918826e+02,        /* 0x40745CAE, 0x221B9F0A */
 172     sb3  =  1.53672958608443695994e+03,        /* 0x409802EB, 0x189D5118 */
 173     sb4  =  3.19985821950859553908e+03,        /* 0x40A8FFB7, 0x688C246A */
 174     sb5  =  2.55305040643316442583e+03,        /* 0x40A3F219, 0xCEDF3BE6 */
 175     sb6  =  4.74528541206955367215e+02,        /* 0x407DA874, 0xE79FE763 */
 176     sb7  = -2.24409524465858183362e+01;        /* 0xC03670E2, 0x42712D62 */
 177
 178 double gmx_erfd(double x)
 179 {
 180 #ifdef GMX_FLOAT_FORMAT_IEEE754
 181     gmx_int32_t hx, ix, i;
 182     double      R, S, P, Q, s, y, z, r;
 183
 184     union
 185     {
 186         double d;
 187         int    i[2];
 188     }
 189     conv;
 190
 191     conv.d = x;
 192
 193 #ifdef GMX_IEEE754_BIG_ENDIAN_WORD_ORDER
 194     hx = conv.i[0];
 195 #else
 196     hx = conv.i[1];
 197 #endif
 198
 199     ix = hx&0x7fffffff;
 200     if (ix >= 0x7ff00000)
 201     {
 202         /* erf(nan)=nan */
 203         i = ((gmx_uint32_t)hx>>31)<<1;
 204         return (double)(1-i)+one/x; /* erf(+-inf)=+-1 */
 205     }
 206
 207     if (ix < 0x3feb0000)
 208     {
 209         /* |x|<0.84375 */
 210         if (ix < 0x3e300000)
 211         {
 212             /* |x|<2**-28 */
 213             if (ix < 0x00800000)
 214             {
 215                 return 0.125*(8.0*x+efx8*x);  /*avoid underflow */
 216             }
 217             return x + efx*x;
 218         }
 219         z = x*x;
 220         r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
 221         s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
 222         y = r/s;
 223         return x + x*y;
 224     }
 225     if (ix < 0x3ff40000)
 226     {
 227         /* 0.84375 <= |x| < 1.25 */
 228         s = fabs(x)-one;
 229         P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
 230         Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
 231         if (hx >= 0)
 232         {
 233             return erx + P/Q;
 234         }
 235         else
 236         {
 237             return -erx - P/Q;
 238         }
 239     }
 240     if (ix >= 0x40180000)
 241     {
 242         /* inf>|x|>=6 */
 243         if (hx >= 0)
 244         {
 245             return one-tiny;
 246         }
 247         else
 248         {
 249             return tiny-one;
 250         }
 251     }
 252     x = fabs(x);
 253     s = one/(x*x);
 254     if (ix < 0x4006DB6E)
 255     {
 256         /* |x| < 1/0.35 */
 257         R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))));
 258         S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+s*sa8)))))));
 259     }
 260     else
 261     {
 262         /* |x| >= 1/0.35 */
 263         R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))));
 264         S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
 265     }
 266
 267     conv.d = x;
 268
 269 #ifdef GMX_IEEE754_BIG_ENDIAN_WORD_ORDER
 270     conv.i[1] = 0;
 271 #else
 272     conv.i[0] = 0;
 273 #endif
 274
 275     z = conv.d;
 276
 277     r  =  exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
 278     if (hx >= 0)
 279     {
 280         return one-r/x;
 281     }
 282     else
 283     {
 284         return r/x-one;
 285     }
 286 #else
 287     /* No IEEE754 information. We need to trust that the OS provides erf(). */
 288     return erf(x);
 289 #endif
 290 }
 291
 292
 293 double gmx_erfcd(double x)
 294 {
 295 #ifdef GMX_FLOAT_FORMAT_IEEE754
 296     gmx_int32_t hx, ix;
 297     double      R, S, P, Q, s, y, z, r;
 298
 299     union
 300     {
 301         double d;
 302         int    i[2];
 303     }
 304     conv;
 305
 306     conv.d = x;
 307
 308 #ifdef GMX_IEEE754_BIG_ENDIAN_WORD_ORDER
 309     hx = conv.i[0];
 310 #else
 311     hx = conv.i[1];
 312 #endif
 313
 314     ix = hx&0x7fffffff;
 315     if (ix >= 0x7ff00000)
 316     {
 317         /* erfc(nan)=nan */
 318         /* erfc(+-inf)=0,2 */
 319         return (double)(((gmx_uint32_t)hx>>31)<<1)+one/x;
 320     }
 321
 322     if (ix < 0x3feb0000)
 323     {
 324         /* |x|<0.84375 */
 325         double r1, r2, s1, s2, s3, z2, z4;
 326         if (ix < 0x3c700000)     /* |x|<2**-56 */
 327         {
 328             return one-x;
 329         }
 330         z = x*x;
 331         r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
 332         s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
 333         y = r/s;
 334         if (hx < 0x3fd00000)
 335         {
 336             /* x<1/4 */
 337             return one-(x+x*y);
 338         }
 339         else
 340         {
 341             r  = x*y;
 342             r += (x-half);
 343             return half - r;
 344         }
 345     }
 346
 347     if (ix < 0x3ff40000)
 348     {
 349         /* 0.84375 <= |x| < 1.25 */
 350         s = fabs(x)-one;
 351         P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
 352         Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
 353         if (hx >= 0)
 354         {
 355             z  = one-erx; return z - P/Q;
 356         }
 357         else
 358         {
 359             z = erx+P/Q; return one+z;
 360         }
 361     }
 362     if (ix < 0x403c0000)
 363     {
 364         /* |x|<28 */
 365         x = fabs(x);
 366         s = one/(x*x);
 367         if (ix < 0x4006DB6D)
 368         {
 369             /* |x| < 1/.35 ~ 2.857143*/
 370             R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))));
 371             S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+s*sa8)))))));
 372         }
 373         else
 374         {
 375             /* |x| >= 1/.35 ~ 2.857143 */
 376             if (hx < 0 && ix >= 0x40180000)
 377             {
 378                 return two-tiny; /* x < -6 */
 379             }
 380             R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))));
 381             S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
 382         }
 383
 384         conv.d = x;
 385
 386 #ifdef GMX_IEEE754_BIG_ENDIAN_WORD_ORDER
 387         conv.i[1] = 0;
 388 #else
 389         conv.i[0] = 0;
 390 #endif
 391
 392         z = conv.d;
 393
 394         r  =  exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
 395
 396         if (hx > 0)
 397         {
 398             return r/x;
 399         }
 400         else
 401         {
 402             return two-r/x;
 403         }
 404     }
 405     else
 406     {
 407         if (hx > 0)
 408         {
 409             return tiny*tiny;
 410         }
 411         else
 412         {
 413             return two-tiny;
 414         }
 415     }
 416 #else
 417     /* No IEEE754 information. We need to trust that the OS provides erfc(). */
 418     return erfc(x);
 419 #endif
 420 }
 421
 422
 423 static const float
 424     tinyf =  1e-30,
 425     halff =  5.0000000000e-01, /* 0x3F000000 */
 426     onef  =  1.0000000000e+00, /* 0x3F800000 */
 427     twof  =  2.0000000000e+00, /* 0x40000000 */
 428 /* c = (subfloat)0.84506291151 */
 429     erxf =  8.4506291151e-01,  /* 0x3f58560b */
 430 /*
 431  * Coefficients for approximation to  erf on [0,0.84375]
 432  */
 433     efxf  =  1.2837916613e-01, /* 0x3e0375d4 */
 434     efx8f =  1.0270333290e+00, /* 0x3f8375d4 */
 435     pp0f  =  1.2837916613e-01, /* 0x3e0375d4 */
 436     pp1f  = -3.2504209876e-01, /* 0xbea66beb */
 437     pp2f  = -2.8481749818e-02, /* 0xbce9528f */
 438     pp3f  = -5.7702702470e-03, /* 0xbbbd1489 */
 439     pp4f  = -2.3763017452e-05, /* 0xb7c756b1 */
 440     qq1f  =  3.9791721106e-01, /* 0x3ecbbbce */
 441     qq2f  =  6.5022252500e-02, /* 0x3d852a63 */
 442     qq3f  =  5.0813062117e-03, /* 0x3ba68116 */
 443     qq4f  =  1.3249473704e-04, /* 0x390aee49 */
 444     qq5f  = -3.9602282413e-06, /* 0xb684e21a */
 445 /*
 446  * Coefficients for approximation to  erf  in [0.84375,1.25]
 447  */
 448     pa0f = -2.3621185683e-03,  /* 0xbb1acdc6 */
 449     pa1f =  4.1485610604e-01,  /* 0x3ed46805 */
 450     pa2f = -3.7220788002e-01,  /* 0xbebe9208 */
 451     pa3f =  3.1834661961e-01,  /* 0x3ea2fe54 */
 452     pa4f = -1.1089469492e-01,  /* 0xbde31cc2 */
 453     pa5f =  3.5478305072e-02,  /* 0x3d1151b3 */
 454     pa6f = -2.1663755178e-03,  /* 0xbb0df9c0 */
 455     qa1f =  1.0642088205e-01,  /* 0x3dd9f331 */
 456     qa2f =  5.4039794207e-01,  /* 0x3f0a5785 */
 457     qa3f =  7.1828655899e-02,  /* 0x3d931ae7 */
 458     qa4f =  1.2617121637e-01,  /* 0x3e013307 */
 459     qa5f =  1.3637083583e-02,  /* 0x3c5f6e13 */
 460     qa6f =  1.1984500103e-02,  /* 0x3c445aa3 */
 461 /*
 462  * Coefficients for approximation to  erfc in [1.25,1/0.35]
 463  */
 464     ra0f = -9.8649440333e-03,  /* 0xbc21a093 */
 465     ra1f = -6.9385856390e-01,  /* 0xbf31a0b7 */
 466     ra2f = -1.0558626175e+01,  /* 0xc128f022 */
 467     ra3f = -6.2375331879e+01,  /* 0xc2798057 */
 468     ra4f = -1.6239666748e+02,  /* 0xc322658c */
 469     ra5f = -1.8460508728e+02,  /* 0xc3389ae7 */
 470     ra6f = -8.1287437439e+01,  /* 0xc2a2932b */
 471     ra7f = -9.8143291473e+00,  /* 0xc11d077e */
 472     sa1f =  1.9651271820e+01,  /* 0x419d35ce */
 473     sa2f =  1.3765776062e+02,  /* 0x4309a863 */
 474     sa3f =  4.3456588745e+02,  /* 0x43d9486f */
 475     sa4f =  6.4538726807e+02,  /* 0x442158c9 */
 476     sa5f =  4.2900814819e+02,  /* 0x43d6810b */
 477     sa6f =  1.0863500214e+02,  /* 0x42d9451f */
 478     sa7f =  6.5702495575e+00,  /* 0x40d23f7c */
 479     sa8f = -6.0424413532e-02,  /* 0xbd777f97 */
 480 /*
 481  * Coefficients for approximation to  erfc in [1/.35,28]
 482  */
 483     rb0f = -9.8649431020e-03,  /* 0xbc21a092 */
 484     rb1f = -7.9928326607e-01,  /* 0xbf4c9dd4 */
 485     rb2f = -1.7757955551e+01,  /* 0xc18e104b */
 486     rb3f = -1.6063638306e+02,  /* 0xc320a2ea */
 487     rb4f = -6.3756646729e+02,  /* 0xc41f6441 */
 488     rb5f = -1.0250950928e+03,  /* 0xc480230b */
 489     rb6f = -4.8351919556e+02,  /* 0xc3f1c275 */
 490     sb1f =  3.0338060379e+01,  /* 0x41f2b459 */
 491     sb2f =  3.2579251099e+02,  /* 0x43a2e571 */
 492     sb3f =  1.5367296143e+03,  /* 0x44c01759 */
 493     sb4f =  3.1998581543e+03,  /* 0x4547fdbb */
 494     sb5f =  2.5530502930e+03,  /* 0x451f90ce */
 495     sb6f =  4.7452853394e+02,  /* 0x43ed43a7 */
 496     sb7f = -2.2440952301e+01;  /* 0xc1b38712 */
 497
 498
 499 typedef union
 500 {
 501     float         value;
 502     gmx_uint32_t  word;
 503 } ieee_float_shape_type;
 504
 505 #define GET_FLOAT_WORD(i, d)                 \
 506     do {                                \
 507         ieee_float_shape_type gf_u;                   \
 508         gf_u.value = (d);                     \
 509         (i)        = gf_u.word;                      \
 510     } while (0)
 511
 512
 513 #define SET_FLOAT_WORD(d, i)                 \
 514     do {                                \
 515         ieee_float_shape_type sf_u;                   \
 516         sf_u.word = (i);                      \
 517         (d)       = sf_u.value;                     \
 518     } while (0)
 519
 520
 521 float gmx_erff(float x)
 522 {
 523     gmx_int32_t hx, ix, i;
 524     float       R, S, P, Q, s, y, z, r;
 525
 526     union
 527     {
 528         float  f;
 529         int    i;
 530     }
 531     conv;
 532
 533     conv.f = x;
 534     hx     = conv.i;
 535
 536     ix = hx&0x7fffffff;
 537     if (ix >= 0x7f800000)
 538     {
 539         /* erf(nan)=nan */
 540         i = ((gmx_uint32_t)hx>>31)<<1;
 541         return (float)(1-i)+onef/x; /* erf(+-inf)=+-1 */
 542     }
 543
 544     if (ix < 0x3f580000)
 545     {
 546         /* |x|<0.84375 */
 547         if (ix < 0x31800000)
 548         {
 549             /* |x|<2**-28 */
 550             if (ix < 0x04000000)
 551             {
 552                 return (float)0.125*((float)8.0*x+efx8f*x);             /*avoid underflow */
 553             }
 554             return x + efxf*x;
 555         }
 556         z = x*x;
 557         r = pp0f+z*(pp1f+z*(pp2f+z*(pp3f+z*pp4f)));
 558         s = onef+z*(qq1f+z*(qq2f+z*(qq3f+z*(qq4f+z*qq5f))));
 559         y = r/s;
 560         return x + x*y;
 561     }
 562     if (ix < 0x3fa00000)
 563     {
 564         /* 0.84375 <= |x| < 1.25 */
 565         s = fabs(x)-onef;
 566         P = pa0f+s*(pa1f+s*(pa2f+s*(pa3f+s*(pa4f+s*(pa5f+s*pa6f)))));
 567         Q = onef+s*(qa1f+s*(qa2f+s*(qa3f+s*(qa4f+s*(qa5f+s*qa6f)))));
 568         if (hx >= 0)
 569         {
 570             return erxf + P/Q;
 571         }
 572         else
 573         {
 574             return -erxf - P/Q;
 575         }
 576     }
 577     if (ix >= 0x40c00000)
 578     {
 579         /* inf>|x|>=6 */
 580         if (hx >= 0)
 581         {
 582             return onef-tinyf;
 583         }
 584         else
 585         {
 586             return tinyf-onef;
 587         }
 588     }
 589     x = fabs(x);
 590     s = onef/(x*x);
 591     if (ix < 0x4036DB6E)
 592     {
 593         /* |x| < 1/0.35 */
 594         R = ra0f+s*(ra1f+s*(ra2f+s*(ra3f+s*(ra4f+s*(ra5f+s*(ra6f+s*ra7f))))));
 595         S = onef+s*(sa1f+s*(sa2f+s*(sa3f+s*(sa4f+s*(sa5f+s*(sa6f+s*(sa7f+s*sa8f)))))));
 596     }
 597     else
 598     {
 599         /* |x| >= 1/0.35 */
 600         R = rb0f+s*(rb1f+s*(rb2f+s*(rb3f+s*(rb4f+s*(rb5f+s*rb6f)))));
 601         S = onef+s*(sb1f+s*(sb2f+s*(sb3f+s*(sb4f+s*(sb5f+s*(sb6f+s*sb7f))))));
 602     }
 603
 604     conv.f = x;
 605     conv.i = conv.i & 0xfffff000;
 606     z      = conv.f;
 607
 608     r  =  exp(-z*z-(float)0.5625)*exp((z-x)*(z+x)+R/S);
 609     if (hx >= 0)
 610     {
 611         return onef-r/x;
 612     }
 613     else
 614     {
 615         return r/x-onef;
 616     }
 617 }
 618
 619 float gmx_erfcf(float x)
 620 {
 621     gmx_int32_t hx, ix;
 622     float       R, S, P, Q, s, y, z, r;
 623
 624     union
 625     {
 626         float  f;
 627         int    i;
 628     }
 629     conv;
 630
 631     conv.f = x;
 632     hx     = conv.i;
 633
 634     ix = hx&0x7fffffff;
 635     if (ix >= 0x7f800000)
 636     {
 637         /* erfc(nan)=nan */
 638         /* erfc(+-inf)=0,2 */
 639         return (float)(((gmx_uint32_t)hx>>31)<<1)+onef/x;
 640     }
 641
 642     if (ix < 0x3f580000)
 643     {
 644         /* |x|<0.84375 */
 645         if (ix < 0x23800000)
 646         {
 647             return onef-x;  /* |x|<2**-56 */
 648         }
 649         z = x*x;
 650         r = pp0f+z*(pp1f+z*(pp2f+z*(pp3f+z*pp4f)));
 651         s = onef+z*(qq1f+z*(qq2f+z*(qq3f+z*(qq4f+z*qq5f))));
 652         y = r/s;
 653         if (hx < 0x3e800000)
 654         {
 655             /* x<1/4 */
 656             return onef-(x+x*y);
 657         }
 658         else
 659         {
 660             r  = x*y;
 661             r += (x-halff);
 662             return halff - r;
 663         }
 664     }
 665     if (ix < 0x3fa00000)
 666     {
 667         /* 0.84375 <= |x| < 1.25 */
 668         s = fabs(x)-onef;
 669         P = pa0f+s*(pa1f+s*(pa2f+s*(pa3f+s*(pa4f+s*(pa5f+s*pa6f)))));
 670         Q = onef+s*(qa1f+s*(qa2f+s*(qa3f+s*(qa4f+s*(qa5f+s*qa6f)))));
 671         if (hx >= 0)
 672         {
 673             z  = onef-erxf; return z - P/Q;
 674         }
 675         else
 676         {
 677             z = erxf+P/Q; return onef+z;
 678         }
 679     }
 680     if (ix < 0x41e00000)
 681     {
 682         /* |x|<28 */
 683         x = fabs(x);
 684         s = onef/(x*x);
 685         if (ix < 0x4036DB6D)
 686         {
 687             /* |x| < 1/.35 ~ 2.857143*/
 688             R = ra0f+s*(ra1f+s*(ra2f+s*(ra3f+s*(ra4f+s*(ra5f+s*(ra6f+s*ra7f))))));
 689             S = onef+s*(sa1f+s*(sa2f+s*(sa3f+s*(sa4f+s*(sa5f+s*(sa6f+s*(sa7f+s*sa8f)))))));
 690         }
 691         else
 692         {
 693             /* |x| >= 1/.35 ~ 2.857143 */
 694             if (hx < 0 && ix >= 0x40c00000)
 695             {
 696                 return twof-tinyf;                     /* x < -6 */
 697             }
 698             R = rb0f+s*(rb1f+s*(rb2f+s*(rb3f+s*(rb4f+s*(rb5f+s*rb6f)))));
 699             S = onef+s*(sb1f+s*(sb2f+s*(sb3f+s*(sb4f+s*(sb5f+s*(sb6f+s*sb7f))))));
 700         }
 701
 702         conv.f = x;
 703         conv.i = conv.i & 0xfffff000;
 704         z      = conv.f;
 705
 706         r  =  exp(-z*z-(float)0.5625)*exp((z-x)*(z+x)+R/S);
 707         if (hx > 0)
 708         {
 709             return r/x;
 710         }
 711         else
 712         {
 713             return twof-r/x;
 714         }
 715     }
 716     else
 717     {
 718         if (hx > 0)
 719         {
 720             return tinyf*tinyf;
 721         }
 722         else
 723         {
 724             return twof-tinyf;
 725         }
 726     }
 727 }
 728
 729
 730 gmx_bool gmx_isfinite(real gmx_unused x)
 731 {
 732     gmx_bool returnval;
 733
 734 #ifdef HAVE__FINITE
 735     returnval = _finite(x);
 736 #elif defined HAVE_ISFINITE
 737     returnval = isfinite(x);
 738 #elif defined HAVE__ISFINITE
 739     returnval = _isfinite(x);
 740 #else
 741     /* If no suitable function was found, assume the value is
 742      * finite. */
 743     returnval = TRUE;
 744 #endif
 745     return returnval;
 746 }
 747
 748 gmx_bool gmx_isnan(real x)
 749 {
 750     return x != x;
 751 }
 752
 753 int
 754 gmx_within_tol(double   f1,
 755                double   f2,
 756                double   tol)
 757 {
 758     /* The or-equal is important - otherwise we return false if f1==f2==0 */
 759     if (fabs(f1-f2) <= tol*0.5*(fabs(f1)+fabs(f2)) )
 760     {
 761         return 1;
 762     }
 763     else
 764     {
 765         return 0;
 766     }
 767 }
 768
 769 int
 770 gmx_numzero(double a)
 771 {
 772     return gmx_within_tol(a, 0.0, GMX_REAL_MIN/GMX_REAL_EPS);
 773 }
 774
 775 unsigned int
 776 gmx_log2i(unsigned int n)
 777 {
 778     assert(n != 0); /* behavior differs for 0 */
 779 #if defined(__INTEL_COMPILER)
 780     return _bit_scan_reverse(n);
 781 #elif defined(__GNUC__) && UINT_MAX == 4294967295U /*also for clang*/
 782     return __builtin_clz(n) ^ 31U;                 /* xor gets optimized out */
 783 #elif defined(_MSC_VER) && _MSC_VER >= 1400
 784     {
 785         unsigned long i;
 786         _BitScanReverse(&i, n);
 787         return i;
 788     }
 789 #elif defined(__xlC__)
 790     return 31 - __cntlz4(n);
 791 #else
 792     /* http://graphics.stanford.edu/~seander/bithacks.html#IntegerLogLookup */
 793 #define LT(n) n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n
 794     static const char     LogTable256[256] = {
 795         -1, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
 796         LT(4), LT(5), LT(5), LT(6), LT(6), LT(6), LT(6),
 797         LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7)
 798     };
 799 #undef LT
 800
 801     unsigned int r;     /* r will be lg(n) */
 802     unsigned int t, tt; /* temporaries */
 803
 804     if ((tt = n >> 16) != 0)
 805     {
 806         r = ((t = tt >> 8) != 0) ? 24 + LogTable256[t] : 16 + LogTable256[tt];
 807     }
 808     else
 809     {
 810         r = ((t = n >> 8) != 0) ? 8 + LogTable256[t] : LogTable256[n];
 811     }
 812     return r;
 813 #endif
 814 }
 815
 816 gmx_bool
 817 check_int_multiply_for_overflow(gmx_int64_t  a,
 818                                 gmx_int64_t  b,
 819                                 gmx_int64_t *result)
 820 {
 821     gmx_int64_t sign = 1;
 822     if ((0 == a) || (0 == b))
 823     {
 824         *result = 0;
 825         return TRUE;
 826     }
 827     if (a < 0)
 828     {
 829         a    = -a;
 830         sign = -sign;
 831     }
 832     if (b < 0)
 833     {
 834         b    = -b;
 835         sign = -sign;
 836     }
 837     if (GMX_INT64_MAX / b < a)
 838     {
 839         *result = (sign > 0) ? GMX_INT64_MAX : GMX_INT64_MIN;
 840         return FALSE;
 841     }
 842     *result = sign * a * b;
 843     return TRUE;
 844 }
 845
 846 int gmx_greatest_common_divisor(int p, int q)
 847 {
 848     int tmp;
 849     while (q != 0)
 850     {
 851         tmp = q;
 852         q   = p % q;
 853         p   = tmp;
 854     }
 855     return p;
 856 }
 857
 858 int gmx_feenableexcept()
 859 {
 860 #ifdef HAVE_FEENABLEEXCEPT
 861     return feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
 862 #elif (defined(__i386__) || defined(__x86_64__)) && defined(__APPLE__)
 863     /* Author:  David N. Williams
 864      * License:  Public Domain
 865      *
 866      * Might also work on non-Apple Unix. But should be tested
 867      * before enabling.
 868      */
 869     unsigned int  excepts = FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW;
 870     static fenv_t fenv;
 871     unsigned int  new_excepts = excepts & FE_ALL_EXCEPT,
 872                   old_excepts; // previous masks
 873
 874     if (fegetenv (&fenv) )
 875     {
 876         return -1;
 877     }
 878     old_excepts = fenv.__control & FE_ALL_EXCEPT;
 879
 880     // unmask
 881     fenv.__control &= ~new_excepts;
 882     fenv.__mxcsr   &= ~(new_excepts << 7);
 883
 884     return ( fesetenv (&fenv) ? -1 : old_excepts );
 885 #else
 886     return -1;
 887 #endif
 888 }