@@ -670,101 +670,6 @@ m_lgamma(double x)
670670 return r ;
671671}
672672
673- /* IEEE 754-style remainder operation: x - n*y where n*y is the nearest
674- multiple of y to x, taking n even in the case of a tie. Assuming an IEEE 754
675- binary floating-point format, the result is always exact. */
676-
677- static double
678- m_remainder (double x , double y )
679- {
680- /* Deal with most common case first. */
681- if (isfinite (x ) && isfinite (y )) {
682- double absx , absy , c , m , r ;
683-
684- if (y == 0.0 ) {
685- return Py_NAN ;
686- }
687-
688- absx = fabs (x );
689- absy = fabs (y );
690- m = fmod (absx , absy );
691-
692- /*
693- Warning: some subtlety here. What we *want* to know at this point is
694- whether the remainder m is less than, equal to, or greater than half
695- of absy. However, we can't do that comparison directly because we
696- can't be sure that 0.5*absy is representable (the multiplication
697- might incur precision loss due to underflow). So instead we compare
698- m with the complement c = absy - m: m < 0.5*absy if and only if m <
699- c, and so on. The catch is that absy - m might also not be
700- representable, but it turns out that it doesn't matter:
701-
702- - if m > 0.5*absy then absy - m is exactly representable, by
703- Sterbenz's lemma, so m > c
704- - if m == 0.5*absy then again absy - m is exactly representable
705- and m == c
706- - if m < 0.5*absy then either (i) 0.5*absy is exactly representable,
707- in which case 0.5*absy < absy - m, so 0.5*absy <= c and hence m <
708- c, or (ii) absy is tiny, either subnormal or in the lowest normal
709- binade. Then absy - m is exactly representable and again m < c.
710- */
711-
712- c = absy - m ;
713- if (m < c ) {
714- r = m ;
715- }
716- else if (m > c ) {
717- r = - c ;
718- }
719- else {
720- /*
721- Here absx is exactly halfway between two multiples of absy,
722- and we need to choose the even multiple. x now has the form
723-
724- absx = n * absy + m
725-
726- for some integer n (recalling that m = 0.5*absy at this point).
727- If n is even we want to return m; if n is odd, we need to
728- return -m.
729-
730- So
731-
732- 0.5 * (absx - m) = (n/2) * absy
733-
734- and now reducing modulo absy gives us:
735-
736- | m, if n is odd
737- fmod(0.5 * (absx - m), absy) = |
738- | 0, if n is even
739-
740- Now m - 2.0 * fmod(...) gives the desired result: m
741- if n is even, -m if m is odd.
742-
743- Note that all steps in fmod(0.5 * (absx - m), absy)
744- will be computed exactly, with no rounding error
745- introduced.
746- */
747- assert (m == c );
748- r = m - 2.0 * fmod (0.5 * (absx - m ), absy );
749- }
750- return copysign (1.0 , x ) * r ;
751- }
752-
753- /* Special values. */
754- if (isnan (x )) {
755- return x ;
756- }
757- if (isnan (y )) {
758- return y ;
759- }
760- if (isinf (x )) {
761- return Py_NAN ;
762- }
763- assert (isinf (y ));
764- return x ;
765- }
766-
767-
768673/*
769674 Various platforms (Solaris, OpenBSD) do nonstandard things for log(0),
770675 log(-ve), log(NaN). Here are wrappers for log and log10 that deal with
@@ -1298,7 +1203,7 @@ FUNC1D(log1p, m_log1p, 0,
12981203 "Return the natural logarithm of 1+x (base e).\n\n"
12991204 "The result is computed in a way which is accurate for x near zero." ,
13001205 "expected argument value > -1, got %s" )
1301- FUNC2 (remainder , m_remainder ,
1206+ FUNC2 (remainder , remainder ,
13021207 "remainder($module, x, y, /)\n--\n\n"
13031208 "Difference between x and the closest integer multiple of y.\n\n"
13041209 "Return x - n*y where n*y is the closest integer multiple of y.\n"
0 commit comments