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