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