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 * Copyright (c) 2013,2014,2015, 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.
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.
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.
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.
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.
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.
39 #include "utilities.h"
50 #ifndef GMX_NATIVE_WINDOWS
59 result = (a < 0.) ? ((int)(a - half)) : ((int)(a + half));
67 return (-pow(-x, 1.0/3.0));
71 return (pow(x, 1.0/3.0));
75 real sign(real x, real y)
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
93 * ====================================================
94 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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
100 * ====================================================
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 */
111 * Coefficients for approximation to erf on [0,0.84375]
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 */
126 * Coefficients for approximation to erf in [0.84375,1.25]
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 */
142 * Coefficients for approximation to erfc in [1.25,1/0.35]
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 */
161 * Coefficients for approximation to erfc in [1/.35,28]
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 */
178 double gmx_erfd(double x)
180 #if GMX_FLOAT_FORMAT_IEEE754
181 gmx_int32_t hx, ix, i;
182 double R, S, P, Q, s, y, z, r;
193 #if GMX_IEEE754_BIG_ENDIAN_WORD_ORDER
200 if (ix >= 0x7ff00000)
203 i = ((gmx_uint32_t)hx>>31)<<1;
204 return (double)(1-i)+one/x; /* erf(+-inf)=+-1 */
215 return 0.125*(8.0*x+efx8*x); /*avoid underflow */
220 r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
221 s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
227 /* 0.84375 <= |x| < 1.25 */
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)))));
240 if (ix >= 0x40180000)
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)))))));
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))))));
269 #if GMX_IEEE754_BIG_ENDIAN_WORD_ORDER
277 r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
287 /* No IEEE754 information. We need to trust that the OS provides erf(). */
293 double gmx_erfcd(double x)
295 #if GMX_FLOAT_FORMAT_IEEE754
297 double R, S, P, Q, s, y, z, r;
308 #if GMX_IEEE754_BIG_ENDIAN_WORD_ORDER
315 if (ix >= 0x7ff00000)
318 /* erfc(+-inf)=0,2 */
319 return (double)(((gmx_uint32_t)hx>>31)<<1)+one/x;
325 double r1, r2, s1, s2, s3, z2, z4;
326 if (ix < 0x3c700000) /* |x|<2**-56 */
331 r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
332 s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
349 /* 0.84375 <= |x| < 1.25 */
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)))));
355 z = one-erx; return z - P/Q;
359 z = erx+P/Q; return one+z;
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)))))));
375 /* |x| >= 1/.35 ~ 2.857143 */
376 if (hx < 0 && ix >= 0x40180000)
378 return two-tiny; /* x < -6 */
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))))));
386 #if GMX_IEEE754_BIG_ENDIAN_WORD_ORDER
394 r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
417 /* No IEEE754 information. We need to trust that the OS provides erfc(). */
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 */
431 * Coefficients for approximation to erf on [0,0.84375]
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 */
446 * Coefficients for approximation to erf in [0.84375,1.25]
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 */
462 * Coefficients for approximation to erfc in [1.25,1/0.35]
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 */
481 * Coefficients for approximation to erfc in [1/.35,28]
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 */
503 } ieee_float_shape_type;
505 #define GET_FLOAT_WORD(i, d) \
507 ieee_float_shape_type gf_u; \
513 #define SET_FLOAT_WORD(d, i) \
515 ieee_float_shape_type sf_u; \
521 float gmx_erff(float x)
523 gmx_int32_t hx, ix, i;
524 float R, S, P, Q, s, y, z, r;
537 if (ix >= 0x7f800000)
540 i = ((gmx_uint32_t)hx>>31)<<1;
541 return (float)(1-i)+onef/x; /* erf(+-inf)=+-1 */
552 return (float)0.125*((float)8.0*x+efx8f*x); /*avoid underflow */
557 r = pp0f+z*(pp1f+z*(pp2f+z*(pp3f+z*pp4f)));
558 s = onef+z*(qq1f+z*(qq2f+z*(qq3f+z*(qq4f+z*qq5f))));
564 /* 0.84375 <= |x| < 1.25 */
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)))));
577 if (ix >= 0x40c00000)
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)))))));
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))))));
605 conv.i = conv.i & 0xfffff000;
608 r = exp(-z*z-(float)0.5625)*exp((z-x)*(z+x)+R/S);
619 float gmx_erfcf(float x)
622 float R, S, P, Q, s, y, z, r;
635 if (ix >= 0x7f800000)
638 /* erfc(+-inf)=0,2 */
639 return (float)(((gmx_uint32_t)hx>>31)<<1)+onef/x;
647 return onef-x; /* |x|<2**-56 */
650 r = pp0f+z*(pp1f+z*(pp2f+z*(pp3f+z*pp4f)));
651 s = onef+z*(qq1f+z*(qq2f+z*(qq3f+z*(qq4f+z*qq5f))));
667 /* 0.84375 <= |x| < 1.25 */
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)))));
673 z = onef-erxf; return z - P/Q;
677 z = erxf+P/Q; return onef+z;
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)))))));
693 /* |x| >= 1/.35 ~ 2.857143 */
694 if (hx < 0 && ix >= 0x40c00000)
696 return twof-tinyf; /* x < -6 */
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))))));
703 conv.i = conv.i & 0xfffff000;
706 r = exp(-z*z-(float)0.5625)*exp((z-x)*(z+x)+R/S);
730 gmx_bool gmx_isfinite(real gmx_unused x)
735 returnval = _finite(x);
736 #elif defined HAVE_ISFINITE
737 returnval = isfinite(x);
738 #elif defined HAVE__ISFINITE
739 returnval = _isfinite(x);
741 /* If no suitable function was found, assume the value is
748 gmx_bool gmx_isnan(real x)
754 gmx_within_tol(double f1,
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)) )
770 gmx_numzero(double a)
772 return gmx_within_tol(a, 0.0, GMX_REAL_MIN/GMX_REAL_EPS);
776 gmx_log2i(unsigned int n)
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
786 _BitScanReverse(&i, n);
789 #elif defined(__xlC__)
790 return 31 - __cntlz4(n);
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)
801 unsigned int r; /* r will be lg(n) */
802 unsigned int t, tt; /* temporaries */
804 if ((tt = n >> 16) != 0)
806 r = ((t = tt >> 8) != 0) ? 24 + LogTable256[t] : 16 + LogTable256[tt];
810 r = ((t = n >> 8) != 0) ? 8 + LogTable256[t] : LogTable256[n];
817 check_int_multiply_for_overflow(gmx_int64_t a,
821 gmx_int64_t sign = 1;
822 if ((0 == a) || (0 == b))
837 if (GMX_INT64_MAX / b < a)
839 *result = (sign > 0) ? GMX_INT64_MAX : GMX_INT64_MIN;
842 *result = sign * a * b;
846 int gmx_greatest_common_divisor(int p, int q)
858 int gmx_feenableexcept()
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
866 * Might also work on non-Apple Unix. But should be tested
869 unsigned int excepts = FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW;
871 unsigned int new_excepts = excepts & FE_ALL_EXCEPT,
872 old_excepts; // previous masks
874 if (fegetenv (&fenv) )
878 old_excepts = fenv.__control & FE_ALL_EXCEPT;
881 fenv.__control &= ~new_excepts;
882 fenv.__mxcsr &= ~(new_excepts << 7);
884 return ( fesetenv (&fenv) ? -1 : old_excepts );