newlib-cygwin/newlib/libm/common/s_scalbn.c

110 lines
2.8 KiB
C
Raw Permalink Normal View History

2000-02-18 03:39:52 +08:00
/* @(#)s_scalbn.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
FUNCTION
<<scalbn>>, <<scalbnf>>, <<scalbln>>, <<scalblnf>>---scale by power of FLT_RADIX (=2)
2000-02-18 03:39:52 +08:00
INDEX
scalbn
INDEX
scalbnf
2009-03-25 Craig Howland <howland@LGSInnovations.com> * libc/include/math.h: (llround, llroundf): Declare. * libm/common/s_llround.c: New file, implementing llround(). * libm/common/sf_llround.c: New file, implementing llroundf(). * libm/common/sf_lround.c: Remove spurious cast in _DOUBLE_IS_32BITS version of function. * libm/common/sf_lrint.c: Ditto. * libm/common/sf_logb.c: Corrected return for subnormal argument by replacing existing function with a version created from sf_ilogb.c. * libm/common/s_logb.c: Ditto, except starting point s_ilogb.c. Also added documentation for logb() and logbf(). * libm/common/s_signbit.c: Add signbit() documentation. * libm/common/s_log2.c: Update return values to match what w_log2.c has, since log2 uses log(); add note about being derived instead of direct. * libm/common/sf_fma.c: Add casts to attempt to get correct results, as well as comments pointing out problems with the implementation. * libm/common/s_fma.c: Add fma() and fmaf() documentation. * libm/common/sf_remquo.c: Incorrect quotient returns for large values corrected by discarding existing function and replacing with Sun verion, with some enhancements. * libm/common/s_remquo.c: Ditto. Add remquo() and remquof() documentation. * libm/common/s_fmax.c: Add fmax() and fmaxf() documentation. * libm/common/s_fmin.c: Add fmin() and fminf() documentation. * libm/common/s_fdim.c: Return NAN for NAN arg, add fdim() and fdimf() documentation. * libm/common/sf_fdim.c: Return NAN for NAN arg, HUGE_VALF for inf arg. * libm/common/s_trunc.c: Add trunc() and truncf() documentation. * libm/common/s_rint.c: Add rint() and rintf() documentation. * libm/common/s_round.c: Add round() and roundf() documentation. * libm/common/s_scalbn.c: Add scalbln() and scalblnf() documentation. * libm/common/s_infinity.c: Add infinity() and infinityf() documentation. * libm/common/s_lround.c: Add lround(), lroundf(), llround(), and llroundf() documentation. * libm/common/s_lrint.c: Add lrint(), lrintf(), llrint(), and llrintf() documentation. * libm/common/isgreater.c: New file for documenting math.h function-like macros isgreater(), isgreaterequal(), isless(), islessequal(), islessgreater(), and isunordered(). * libm/common/s_isnan.c: Add documentation for function-like macros fpclassify(), isfinite(), isinf(), isnan(), and isnormal(). * libm/common/s_nearbyint.c: Add nearbyint() and nearbyintf() documentation. * libm/common/Makefile.am: Add s_llround.c (src); sf_llround.c (fsrc); s_fdim.def, s_fma.def, s_fmax.def, s_fmin.def, s_logb.def, s_lrint.def, s_lround.def, s_nearbyint.def, s_remquo.def, s_rint.def, s_round.def, s_signbit.def, s_trunc.def, and isgreater.def (chobj); re-name all existing chew files (chobj) to match source file base names (put in underscores), delete all special targets for chew files (leaving all to be generated by rule). * libm/common/Makefile.in: regenerate. * libm/math/w_exp2.c: Add "base 2" to documentation description (and delete TRAD_SYNOPSIS). * libm/math/w_gamma.c: Add tgamma() and tgammaf() documentation, along with some history behind the function names. * libm/math/math.tex: Add includes for newly-added documentation (see .def additions to common/Makefile.am and math/Makefile.am in this ChangeLog list), adjusted existing .def file names to match source file base names (added underscores); add mention of HUGE_VALF; rename "Version of library" section to "Error Handling" and add some text about floating-point exception; added section "Standards Compliance And Portability". * libm/math/Makefile.am: Add w_exp2.def (chobj); re-name all existing chew files (chobj) to match source file base names, delete all special targets for chew files (leaving all to be generated by rule). * libm/math/Makefile.in: regenerated * doc/makedoc.c: Change silent ignoring of commands < 5 characters to a failure when reading macro file for commands < 4 characters; add -v (verbose) option for printing some debugging information; get rid of spurious translation of "@*" to "*" (no source files used @*, so no existing doc pages were affected); clean up some compiler warnings. * doc/doc.str: add BUGS and SEEALSO sections (to match texi2pod.pl which has them); Remove ITEM command (redundant with makedoc built-in "o", not used in any present source file so nothing is lost, anyway). * HOWTO: New file to hold information for maintainers regarding how to do things. Initial sections on documentation and ELIX levels.
2009-03-26 03:13:24 +08:00
INDEX
scalbln
INDEX
scalblnf
2000-02-18 03:39:52 +08:00
SYNOPSIS
2000-02-18 03:39:52 +08:00
#include <math.h>
2009-03-25 Craig Howland <howland@LGSInnovations.com> * libc/include/math.h: (llround, llroundf): Declare. * libm/common/s_llround.c: New file, implementing llround(). * libm/common/sf_llround.c: New file, implementing llroundf(). * libm/common/sf_lround.c: Remove spurious cast in _DOUBLE_IS_32BITS version of function. * libm/common/sf_lrint.c: Ditto. * libm/common/sf_logb.c: Corrected return for subnormal argument by replacing existing function with a version created from sf_ilogb.c. * libm/common/s_logb.c: Ditto, except starting point s_ilogb.c. Also added documentation for logb() and logbf(). * libm/common/s_signbit.c: Add signbit() documentation. * libm/common/s_log2.c: Update return values to match what w_log2.c has, since log2 uses log(); add note about being derived instead of direct. * libm/common/sf_fma.c: Add casts to attempt to get correct results, as well as comments pointing out problems with the implementation. * libm/common/s_fma.c: Add fma() and fmaf() documentation. * libm/common/sf_remquo.c: Incorrect quotient returns for large values corrected by discarding existing function and replacing with Sun verion, with some enhancements. * libm/common/s_remquo.c: Ditto. Add remquo() and remquof() documentation. * libm/common/s_fmax.c: Add fmax() and fmaxf() documentation. * libm/common/s_fmin.c: Add fmin() and fminf() documentation. * libm/common/s_fdim.c: Return NAN for NAN arg, add fdim() and fdimf() documentation. * libm/common/sf_fdim.c: Return NAN for NAN arg, HUGE_VALF for inf arg. * libm/common/s_trunc.c: Add trunc() and truncf() documentation. * libm/common/s_rint.c: Add rint() and rintf() documentation. * libm/common/s_round.c: Add round() and roundf() documentation. * libm/common/s_scalbn.c: Add scalbln() and scalblnf() documentation. * libm/common/s_infinity.c: Add infinity() and infinityf() documentation. * libm/common/s_lround.c: Add lround(), lroundf(), llround(), and llroundf() documentation. * libm/common/s_lrint.c: Add lrint(), lrintf(), llrint(), and llrintf() documentation. * libm/common/isgreater.c: New file for documenting math.h function-like macros isgreater(), isgreaterequal(), isless(), islessequal(), islessgreater(), and isunordered(). * libm/common/s_isnan.c: Add documentation for function-like macros fpclassify(), isfinite(), isinf(), isnan(), and isnormal(). * libm/common/s_nearbyint.c: Add nearbyint() and nearbyintf() documentation. * libm/common/Makefile.am: Add s_llround.c (src); sf_llround.c (fsrc); s_fdim.def, s_fma.def, s_fmax.def, s_fmin.def, s_logb.def, s_lrint.def, s_lround.def, s_nearbyint.def, s_remquo.def, s_rint.def, s_round.def, s_signbit.def, s_trunc.def, and isgreater.def (chobj); re-name all existing chew files (chobj) to match source file base names (put in underscores), delete all special targets for chew files (leaving all to be generated by rule). * libm/common/Makefile.in: regenerate. * libm/math/w_exp2.c: Add "base 2" to documentation description (and delete TRAD_SYNOPSIS). * libm/math/w_gamma.c: Add tgamma() and tgammaf() documentation, along with some history behind the function names. * libm/math/math.tex: Add includes for newly-added documentation (see .def additions to common/Makefile.am and math/Makefile.am in this ChangeLog list), adjusted existing .def file names to match source file base names (added underscores); add mention of HUGE_VALF; rename "Version of library" section to "Error Handling" and add some text about floating-point exception; added section "Standards Compliance And Portability". * libm/math/Makefile.am: Add w_exp2.def (chobj); re-name all existing chew files (chobj) to match source file base names, delete all special targets for chew files (leaving all to be generated by rule). * libm/math/Makefile.in: regenerated * doc/makedoc.c: Change silent ignoring of commands < 5 characters to a failure when reading macro file for commands < 4 characters; add -v (verbose) option for printing some debugging information; get rid of spurious translation of "@*" to "*" (no source files used @*, so no existing doc pages were affected); clean up some compiler warnings. * doc/doc.str: add BUGS and SEEALSO sections (to match texi2pod.pl which has them); Remove ITEM command (redundant with makedoc built-in "o", not used in any present source file so nothing is lost, anyway). * HOWTO: New file to hold information for maintainers regarding how to do things. Initial sections on documentation and ELIX levels.
2009-03-26 03:13:24 +08:00
double scalbn(double <[x]>, int <[n]>);
float scalbnf(float <[x]>, int <[n]>);
double scalbln(double <[x]>, long int <[n]>);
float scalblnf(float <[x]>, long int <[n]>);
2000-02-18 03:39:52 +08:00
DESCRIPTION
2009-03-25 Craig Howland <howland@LGSInnovations.com> * libc/include/math.h: (llround, llroundf): Declare. * libm/common/s_llround.c: New file, implementing llround(). * libm/common/sf_llround.c: New file, implementing llroundf(). * libm/common/sf_lround.c: Remove spurious cast in _DOUBLE_IS_32BITS version of function. * libm/common/sf_lrint.c: Ditto. * libm/common/sf_logb.c: Corrected return for subnormal argument by replacing existing function with a version created from sf_ilogb.c. * libm/common/s_logb.c: Ditto, except starting point s_ilogb.c. Also added documentation for logb() and logbf(). * libm/common/s_signbit.c: Add signbit() documentation. * libm/common/s_log2.c: Update return values to match what w_log2.c has, since log2 uses log(); add note about being derived instead of direct. * libm/common/sf_fma.c: Add casts to attempt to get correct results, as well as comments pointing out problems with the implementation. * libm/common/s_fma.c: Add fma() and fmaf() documentation. * libm/common/sf_remquo.c: Incorrect quotient returns for large values corrected by discarding existing function and replacing with Sun verion, with some enhancements. * libm/common/s_remquo.c: Ditto. Add remquo() and remquof() documentation. * libm/common/s_fmax.c: Add fmax() and fmaxf() documentation. * libm/common/s_fmin.c: Add fmin() and fminf() documentation. * libm/common/s_fdim.c: Return NAN for NAN arg, add fdim() and fdimf() documentation. * libm/common/sf_fdim.c: Return NAN for NAN arg, HUGE_VALF for inf arg. * libm/common/s_trunc.c: Add trunc() and truncf() documentation. * libm/common/s_rint.c: Add rint() and rintf() documentation. * libm/common/s_round.c: Add round() and roundf() documentation. * libm/common/s_scalbn.c: Add scalbln() and scalblnf() documentation. * libm/common/s_infinity.c: Add infinity() and infinityf() documentation. * libm/common/s_lround.c: Add lround(), lroundf(), llround(), and llroundf() documentation. * libm/common/s_lrint.c: Add lrint(), lrintf(), llrint(), and llrintf() documentation. * libm/common/isgreater.c: New file for documenting math.h function-like macros isgreater(), isgreaterequal(), isless(), islessequal(), islessgreater(), and isunordered(). * libm/common/s_isnan.c: Add documentation for function-like macros fpclassify(), isfinite(), isinf(), isnan(), and isnormal(). * libm/common/s_nearbyint.c: Add nearbyint() and nearbyintf() documentation. * libm/common/Makefile.am: Add s_llround.c (src); sf_llround.c (fsrc); s_fdim.def, s_fma.def, s_fmax.def, s_fmin.def, s_logb.def, s_lrint.def, s_lround.def, s_nearbyint.def, s_remquo.def, s_rint.def, s_round.def, s_signbit.def, s_trunc.def, and isgreater.def (chobj); re-name all existing chew files (chobj) to match source file base names (put in underscores), delete all special targets for chew files (leaving all to be generated by rule). * libm/common/Makefile.in: regenerate. * libm/math/w_exp2.c: Add "base 2" to documentation description (and delete TRAD_SYNOPSIS). * libm/math/w_gamma.c: Add tgamma() and tgammaf() documentation, along with some history behind the function names. * libm/math/math.tex: Add includes for newly-added documentation (see .def additions to common/Makefile.am and math/Makefile.am in this ChangeLog list), adjusted existing .def file names to match source file base names (added underscores); add mention of HUGE_VALF; rename "Version of library" section to "Error Handling" and add some text about floating-point exception; added section "Standards Compliance And Portability". * libm/math/Makefile.am: Add w_exp2.def (chobj); re-name all existing chew files (chobj) to match source file base names, delete all special targets for chew files (leaving all to be generated by rule). * libm/math/Makefile.in: regenerated * doc/makedoc.c: Change silent ignoring of commands < 5 characters to a failure when reading macro file for commands < 4 characters; add -v (verbose) option for printing some debugging information; get rid of spurious translation of "@*" to "*" (no source files used @*, so no existing doc pages were affected); clean up some compiler warnings. * doc/doc.str: add BUGS and SEEALSO sections (to match texi2pod.pl which has them); Remove ITEM command (redundant with makedoc built-in "o", not used in any present source file so nothing is lost, anyway). * HOWTO: New file to hold information for maintainers regarding how to do things. Initial sections on documentation and ELIX levels.
2009-03-26 03:13:24 +08:00
The <<scalbn>> and <<scalbln>> functions compute
@ifnottex
<[x]> times FLT_RADIX to the power <[n]>.
@end ifnottex
@tex
$x \cdot FLT\_RADIX^n$.
@end tex
efficiently. The result is computed by manipulating the exponent, rather than
by actually performing an exponentiation or multiplication. In this
floating-point implementation FLT_RADIX=2, which makes the <<scalbn>>
functions equivalent to the <<ldexp>> functions.
2000-02-18 03:39:52 +08:00
RETURNS
2009-03-25 Craig Howland <howland@LGSInnovations.com> * libc/include/math.h: (llround, llroundf): Declare. * libm/common/s_llround.c: New file, implementing llround(). * libm/common/sf_llround.c: New file, implementing llroundf(). * libm/common/sf_lround.c: Remove spurious cast in _DOUBLE_IS_32BITS version of function. * libm/common/sf_lrint.c: Ditto. * libm/common/sf_logb.c: Corrected return for subnormal argument by replacing existing function with a version created from sf_ilogb.c. * libm/common/s_logb.c: Ditto, except starting point s_ilogb.c. Also added documentation for logb() and logbf(). * libm/common/s_signbit.c: Add signbit() documentation. * libm/common/s_log2.c: Update return values to match what w_log2.c has, since log2 uses log(); add note about being derived instead of direct. * libm/common/sf_fma.c: Add casts to attempt to get correct results, as well as comments pointing out problems with the implementation. * libm/common/s_fma.c: Add fma() and fmaf() documentation. * libm/common/sf_remquo.c: Incorrect quotient returns for large values corrected by discarding existing function and replacing with Sun verion, with some enhancements. * libm/common/s_remquo.c: Ditto. Add remquo() and remquof() documentation. * libm/common/s_fmax.c: Add fmax() and fmaxf() documentation. * libm/common/s_fmin.c: Add fmin() and fminf() documentation. * libm/common/s_fdim.c: Return NAN for NAN arg, add fdim() and fdimf() documentation. * libm/common/sf_fdim.c: Return NAN for NAN arg, HUGE_VALF for inf arg. * libm/common/s_trunc.c: Add trunc() and truncf() documentation. * libm/common/s_rint.c: Add rint() and rintf() documentation. * libm/common/s_round.c: Add round() and roundf() documentation. * libm/common/s_scalbn.c: Add scalbln() and scalblnf() documentation. * libm/common/s_infinity.c: Add infinity() and infinityf() documentation. * libm/common/s_lround.c: Add lround(), lroundf(), llround(), and llroundf() documentation. * libm/common/s_lrint.c: Add lrint(), lrintf(), llrint(), and llrintf() documentation. * libm/common/isgreater.c: New file for documenting math.h function-like macros isgreater(), isgreaterequal(), isless(), islessequal(), islessgreater(), and isunordered(). * libm/common/s_isnan.c: Add documentation for function-like macros fpclassify(), isfinite(), isinf(), isnan(), and isnormal(). * libm/common/s_nearbyint.c: Add nearbyint() and nearbyintf() documentation. * libm/common/Makefile.am: Add s_llround.c (src); sf_llround.c (fsrc); s_fdim.def, s_fma.def, s_fmax.def, s_fmin.def, s_logb.def, s_lrint.def, s_lround.def, s_nearbyint.def, s_remquo.def, s_rint.def, s_round.def, s_signbit.def, s_trunc.def, and isgreater.def (chobj); re-name all existing chew files (chobj) to match source file base names (put in underscores), delete all special targets for chew files (leaving all to be generated by rule). * libm/common/Makefile.in: regenerate. * libm/math/w_exp2.c: Add "base 2" to documentation description (and delete TRAD_SYNOPSIS). * libm/math/w_gamma.c: Add tgamma() and tgammaf() documentation, along with some history behind the function names. * libm/math/math.tex: Add includes for newly-added documentation (see .def additions to common/Makefile.am and math/Makefile.am in this ChangeLog list), adjusted existing .def file names to match source file base names (added underscores); add mention of HUGE_VALF; rename "Version of library" section to "Error Handling" and add some text about floating-point exception; added section "Standards Compliance And Portability". * libm/math/Makefile.am: Add w_exp2.def (chobj); re-name all existing chew files (chobj) to match source file base names, delete all special targets for chew files (leaving all to be generated by rule). * libm/math/Makefile.in: regenerated * doc/makedoc.c: Change silent ignoring of commands < 5 characters to a failure when reading macro file for commands < 4 characters; add -v (verbose) option for printing some debugging information; get rid of spurious translation of "@*" to "*" (no source files used @*, so no existing doc pages were affected); clean up some compiler warnings. * doc/doc.str: add BUGS and SEEALSO sections (to match texi2pod.pl which has them); Remove ITEM command (redundant with makedoc built-in "o", not used in any present source file so nothing is lost, anyway). * HOWTO: New file to hold information for maintainers regarding how to do things. Initial sections on documentation and ELIX levels.
2009-03-26 03:13:24 +08:00
<[x]> times 2 to the power <[n]>. A range error may occur.
2000-02-18 03:39:52 +08:00
PORTABILITY
2009-03-25 Craig Howland <howland@LGSInnovations.com> * libc/include/math.h: (llround, llroundf): Declare. * libm/common/s_llround.c: New file, implementing llround(). * libm/common/sf_llround.c: New file, implementing llroundf(). * libm/common/sf_lround.c: Remove spurious cast in _DOUBLE_IS_32BITS version of function. * libm/common/sf_lrint.c: Ditto. * libm/common/sf_logb.c: Corrected return for subnormal argument by replacing existing function with a version created from sf_ilogb.c. * libm/common/s_logb.c: Ditto, except starting point s_ilogb.c. Also added documentation for logb() and logbf(). * libm/common/s_signbit.c: Add signbit() documentation. * libm/common/s_log2.c: Update return values to match what w_log2.c has, since log2 uses log(); add note about being derived instead of direct. * libm/common/sf_fma.c: Add casts to attempt to get correct results, as well as comments pointing out problems with the implementation. * libm/common/s_fma.c: Add fma() and fmaf() documentation. * libm/common/sf_remquo.c: Incorrect quotient returns for large values corrected by discarding existing function and replacing with Sun verion, with some enhancements. * libm/common/s_remquo.c: Ditto. Add remquo() and remquof() documentation. * libm/common/s_fmax.c: Add fmax() and fmaxf() documentation. * libm/common/s_fmin.c: Add fmin() and fminf() documentation. * libm/common/s_fdim.c: Return NAN for NAN arg, add fdim() and fdimf() documentation. * libm/common/sf_fdim.c: Return NAN for NAN arg, HUGE_VALF for inf arg. * libm/common/s_trunc.c: Add trunc() and truncf() documentation. * libm/common/s_rint.c: Add rint() and rintf() documentation. * libm/common/s_round.c: Add round() and roundf() documentation. * libm/common/s_scalbn.c: Add scalbln() and scalblnf() documentation. * libm/common/s_infinity.c: Add infinity() and infinityf() documentation. * libm/common/s_lround.c: Add lround(), lroundf(), llround(), and llroundf() documentation. * libm/common/s_lrint.c: Add lrint(), lrintf(), llrint(), and llrintf() documentation. * libm/common/isgreater.c: New file for documenting math.h function-like macros isgreater(), isgreaterequal(), isless(), islessequal(), islessgreater(), and isunordered(). * libm/common/s_isnan.c: Add documentation for function-like macros fpclassify(), isfinite(), isinf(), isnan(), and isnormal(). * libm/common/s_nearbyint.c: Add nearbyint() and nearbyintf() documentation. * libm/common/Makefile.am: Add s_llround.c (src); sf_llround.c (fsrc); s_fdim.def, s_fma.def, s_fmax.def, s_fmin.def, s_logb.def, s_lrint.def, s_lround.def, s_nearbyint.def, s_remquo.def, s_rint.def, s_round.def, s_signbit.def, s_trunc.def, and isgreater.def (chobj); re-name all existing chew files (chobj) to match source file base names (put in underscores), delete all special targets for chew files (leaving all to be generated by rule). * libm/common/Makefile.in: regenerate. * libm/math/w_exp2.c: Add "base 2" to documentation description (and delete TRAD_SYNOPSIS). * libm/math/w_gamma.c: Add tgamma() and tgammaf() documentation, along with some history behind the function names. * libm/math/math.tex: Add includes for newly-added documentation (see .def additions to common/Makefile.am and math/Makefile.am in this ChangeLog list), adjusted existing .def file names to match source file base names (added underscores); add mention of HUGE_VALF; rename "Version of library" section to "Error Handling" and add some text about floating-point exception; added section "Standards Compliance And Portability". * libm/math/Makefile.am: Add w_exp2.def (chobj); re-name all existing chew files (chobj) to match source file base names, delete all special targets for chew files (leaving all to be generated by rule). * libm/math/Makefile.in: regenerated * doc/makedoc.c: Change silent ignoring of commands < 5 characters to a failure when reading macro file for commands < 4 characters; add -v (verbose) option for printing some debugging information; get rid of spurious translation of "@*" to "*" (no source files used @*, so no existing doc pages were affected); clean up some compiler warnings. * doc/doc.str: add BUGS and SEEALSO sections (to match texi2pod.pl which has them); Remove ITEM command (redundant with makedoc built-in "o", not used in any present source file so nothing is lost, anyway). * HOWTO: New file to hold information for maintainers regarding how to do things. Initial sections on documentation and ELIX levels.
2009-03-26 03:13:24 +08:00
ANSI C, POSIX
SEEALSO
<<ldexp>>
2000-02-18 03:39:52 +08:00
*/
/*
* scalbn (double x, int n)
* scalbn(x,n) returns x* 2**n computed by exponent
* manipulation rather than by actually performing an
* exponentiation or a multiplication.
*/
#include "fdlibm.h"
#ifndef _DOUBLE_IS_32BITS
#ifdef __STDC__
static const double
#else
static double
#endif
two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
huge = 1.0e+300,
tiny = 1.0e-300;
#ifdef __STDC__
double scalbn (double x, int n)
#else
double scalbn (x,n)
double x; int n;
#endif
{
__int32_t k,hx,lx;
EXTRACT_WORDS(hx,lx,x);
k = (hx&0x7ff00000)>>20; /* extract exponent */
if (k==0) { /* 0 or subnormal x */
if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
x *= two54;
GET_HIGH_WORD(hx,x);
k = ((hx&0x7ff00000)>>20) - 54;
if (n< -50000) return tiny*x; /*underflow*/
}
if (k==0x7ff) return x+x; /* NaN or Inf */
libm: Fixing overflow handling issue for scalbnf and scalbn cc Aldy Hernandez <aldyh@redhat.com> and Andrew MacLeod <amacleod@redhat.com>, they are author of new VRP analysis for GCC, just to make sure I didn't mis-understanding or mis-interpreting anything on GCC site. GCC 11 have better value range analysis, that give GCC more confidence to perform more aggressive optimization, but it cause scalbn/scalbnf get wrong result. Using scalbn to demostrate what happened on GCC 11, see comments with VRP prefix: ```c double scalbn (double x, int n) { /* VRP RESULT: n = [-INF, +INF] */ __int32_t k,hx,lx; ... k = (hx&0x7ff00000)>>20; /* VRP RESULT: k = [0, 2047] */ if (k==0) { /* VRP RESULT: k = 0 */ ... k = ((hx&0x7ff00000)>>20) - 54; if (n< -50000) return tiny*x; /*underflow*/ /* VRP RESULT: k = -54 */ } /* VRP RESULT: k = [-54, 2047] */ if (k==0x7ff) return x+x; /* NaN or Inf */ /* VRP RESULT: k = [-54, 2046] */ k = k+n; if (k > 0x7fe) return huge*copysign(huge,x); /* overflow */ /* VRP RESULT: k = [-INF, 2046] */ /* VRP RESULT: n = [-INF, 2100], because k + n <= 0x7fe is false, so: 1. -INF < [-54, 2046] + n <= 0x7fe(2046) < INF 2. -INF < [-54, 2046] + n <= 2046 < INF 3. -INF < n <= 2046 - [-54, 2046] < INF 4. -INF < n <= [0, 2100] < INF 5. n = [-INF, 2100] */ if (k > 0) /* normal result */ {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;} if (k <= -54) { /* VRP OPT: Evaluate n > 50000 as true...*/ if (n > 50000) /* in case integer overflow in n+k */ return huge*copysign(huge,x); /*overflow*/ else return tiny*copysign(tiny,x); /*underflow*/ } k += 54; /* subnormal result */ SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x*twom54; } ``` However give the input n = INT32_MAX, k = k+n will overflow, and then we expect got `huge*copysign(huge,x)`, but new VRP optimization think `n > 50000` is never be true, so optimize that into `tiny*copysign(tiny,x)`. so the solution here is to moving the overflow handle logic before `k = k + n`.
2021-07-19 16:50:22 +08:00
if (n > 50000) /* in case integer overflow in n+k */
return huge*copysign(huge,x); /*overflow*/
2000-02-18 03:39:52 +08:00
k = k+n;
if (k > 0x7fe) return huge*copysign(huge,x); /* overflow */
if (k > 0) /* normal result */
{SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
libm: Fixing overflow handling issue for scalbnf and scalbn cc Aldy Hernandez <aldyh@redhat.com> and Andrew MacLeod <amacleod@redhat.com>, they are author of new VRP analysis for GCC, just to make sure I didn't mis-understanding or mis-interpreting anything on GCC site. GCC 11 have better value range analysis, that give GCC more confidence to perform more aggressive optimization, but it cause scalbn/scalbnf get wrong result. Using scalbn to demostrate what happened on GCC 11, see comments with VRP prefix: ```c double scalbn (double x, int n) { /* VRP RESULT: n = [-INF, +INF] */ __int32_t k,hx,lx; ... k = (hx&0x7ff00000)>>20; /* VRP RESULT: k = [0, 2047] */ if (k==0) { /* VRP RESULT: k = 0 */ ... k = ((hx&0x7ff00000)>>20) - 54; if (n< -50000) return tiny*x; /*underflow*/ /* VRP RESULT: k = -54 */ } /* VRP RESULT: k = [-54, 2047] */ if (k==0x7ff) return x+x; /* NaN or Inf */ /* VRP RESULT: k = [-54, 2046] */ k = k+n; if (k > 0x7fe) return huge*copysign(huge,x); /* overflow */ /* VRP RESULT: k = [-INF, 2046] */ /* VRP RESULT: n = [-INF, 2100], because k + n <= 0x7fe is false, so: 1. -INF < [-54, 2046] + n <= 0x7fe(2046) < INF 2. -INF < [-54, 2046] + n <= 2046 < INF 3. -INF < n <= 2046 - [-54, 2046] < INF 4. -INF < n <= [0, 2100] < INF 5. n = [-INF, 2100] */ if (k > 0) /* normal result */ {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;} if (k <= -54) { /* VRP OPT: Evaluate n > 50000 as true...*/ if (n > 50000) /* in case integer overflow in n+k */ return huge*copysign(huge,x); /*overflow*/ else return tiny*copysign(tiny,x); /*underflow*/ } k += 54; /* subnormal result */ SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x*twom54; } ``` However give the input n = INT32_MAX, k = k+n will overflow, and then we expect got `huge*copysign(huge,x)`, but new VRP optimization think `n > 50000` is never be true, so optimize that into `tiny*copysign(tiny,x)`. so the solution here is to moving the overflow handle logic before `k = k + n`.
2021-07-19 16:50:22 +08:00
if (k <= -54)
return tiny*copysign(tiny,x); /*underflow*/
2000-02-18 03:39:52 +08:00
k += 54; /* subnormal result */
SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
return x*twom54;
}
#endif /* _DOUBLE_IS_32BITS */