2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
55 result
= (a
< 0.) ? ((int)(a
- half
)) : ((int)(a
+ half
));
59 real
cuberoot (real x
)
63 return (-pow(-x
,1.0/DIM
));
67 return (pow(x
,1.0/DIM
));
71 real
sign(real x
,real y
)
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
85 * ====================================================
86 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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
92 * ====================================================
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
;
105 # error ERROR: No 32 bit wide integer type found!
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 */
117 * Coefficients for approximation to erf on [0,0.84375]
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 */
132 * Coefficients for approximation to erf in [0.84375,1.25]
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 */
148 * Coefficients for approximation to erfc in [1.25,1/0.35]
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 */
167 * Coefficients for approximation to erfc in [1/.35,28]
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 */
184 double gmx_erfd(double x
)
188 double R
,S
,P
,Q
,s
,y
,z
,r
;
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)
211 i
= ((erf_u_int32_t
)hx
>>31)<<1;
212 return (double)(1-i
)+one
/x
; /* erf(+-inf)=+-1 */
222 return 0.125*(8.0*x
+efx8
*x
); /*avoid underflow */
226 r
= pp0
+z
*(pp1
+z
*(pp2
+z
*(pp3
+z
*pp4
)));
227 s
= one
+z
*(qq1
+z
*(qq2
+z
*(qq3
+z
*(qq4
+z
*qq5
))));
233 /* 0.84375 <= |x| < 1.25 */
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
;
239 if (ix
>= 0x40180000)
242 if(hx
>=0) return one
-tiny
; else return tiny
-one
;
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
)))))));
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
))))));
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)
271 r
= exp(-z
*z
-0.5625)*exp((z
-x
)*(z
+x
)+R
/S
);
279 double gmx_erfcd(double x
)
282 double R
,S
,P
,Q
,s
,y
,z
,r
;
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)
305 /* erfc(+-inf)=0,2 */
306 return (double)(((erf_u_int32_t
)hx
>>31)<<1)+one
/x
;
312 double r1
,r2
,s1
,s2
,s3
,z2
,z4
;
313 if(ix
< 0x3c700000) /* |x|<2**-56 */
316 r
= pp0
+z
*(pp1
+z
*(pp2
+z
*(pp3
+z
*pp4
)));
317 s
= one
+z
*(qq1
+z
*(qq2
+z
*(qq3
+z
*(qq4
+z
*qq5
))));
334 /* 0.84375 <= |x| < 1.25 */
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
)))));
339 z
= one
-erx
; return z
- P
/Q
;
341 z
= erx
+P
/Q
; return one
+z
;
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
)))))));
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
))))));
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)
376 r
= exp(-z
*z
-0.5625)*exp((z
-x
)*(z
+x
)+R
/S
);
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 */
401 * Coefficients for approximation to erf on [0,0.84375]
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 */
416 * Coefficients for approximation to erf in [0.84375,1.25]
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 */
432 * Coefficients for approximation to erfc in [1.25,1/0.35]
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 */
451 * Coefficients for approximation to erfc in [1/.35,28]
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 */
473 } ieee_float_shape_type
;
475 #define GET_FLOAT_WORD(i,d) \
477 ieee_float_shape_type gf_u; \
483 #define SET_FLOAT_WORD(d,i) \
485 ieee_float_shape_type sf_u; \
491 float gmx_erff(float x
)
494 float R
,S
,P
,Q
,s
,y
,z
,r
;
510 i
= ((erf_u_int32_t
)hx
>>31)<<1;
511 return (float)(1-i
)+onef
/x
; /* erf(+-inf)=+-1 */
521 return (float)0.125*((float)8.0*x
+efx8f
*x
); /*avoid underflow */
525 r
= pp0f
+z
*(pp1f
+z
*(pp2f
+z
*(pp3f
+z
*pp4f
)));
526 s
= onef
+z
*(qq1f
+z
*(qq2f
+z
*(qq3f
+z
*(qq4f
+z
*qq5f
))));
532 /* 0.84375 <= |x| < 1.25 */
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
;
538 if (ix
>= 0x40c00000)
541 if(hx
>=0) return onef
-tinyf
; else return tinyf
-onef
;
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
)))))));
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
))))));
559 conv
.i
= conv
.i
& 0xfffff000;
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
;
566 float gmx_erfcf(float x
)
569 float R
,S
,P
,Q
,s
,y
,z
,r
;
585 /* erfc(+-inf)=0,2 */
586 return (float)(((erf_u_int32_t
)hx
>>31)<<1)+onef
/x
;
593 return onef
-x
; /* |x|<2**-56 */
595 r
= pp0f
+z
*(pp1f
+z
*(pp2f
+z
*(pp3f
+z
*pp4f
)));
596 s
= onef
+z
*(qq1f
+z
*(qq2f
+z
*(qq3f
+z
*(qq4f
+z
*qq5f
))));
610 /* 0.84375 <= |x| < 1.25 */
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
)))));
615 z
= onef
-erxf
; return z
- P
/Q
;
617 z
= erxf
+P
/Q
; return onef
+z
;
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
)))))));
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
))))));
638 conv
.i
= conv
.i
& 0xfffff000;
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
;
644 if(hx
>0) return tinyf
*tinyf
; else return twof
-tinyf
;
649 gmx_bool
gmx_isfinite(real x
)
651 gmx_bool returnval
= TRUE
;
652 /* If no suitable function was found, assume the value is
656 returnval
= _finite(x
);
657 #elif defined HAVE_ISFINITE
658 returnval
= isfinite(x
);
659 #elif defined HAVE__ISFINITE
660 returnval
= _isfinite(x
);
666 check_int_multiply_for_overflow(gmx_large_int_t a
,
668 gmx_large_int_t
*result
)
670 gmx_large_int_t sign
= 1;
671 if((0 == a
) || (0 == b
))
686 if(GMX_LARGE_INT_MAX
/ b
< a
)
688 *result
= (sign
> 0) ? GMX_LARGE_INT_MAX
: GMX_LARGE_INT_MIN
;
691 *result
= sign
* a
* b
;