dead code in lgamma_r[f]?
Steve Kargl
sgk at troutmask.apl.washington.edu
Thu Dec 5 18:23:25 UTC 2013
On Thu, Dec 05, 2013 at 09:37:00AM -0800, Steve Kargl wrote:
> In reviewing the implementation details for lgamma_r[f],
> I noticed that in src/e_lgamma.c (and similar code in
> e_lgammaf.c):
These should be e_lgamma_r.c and e_lgammaf_r.c.
I've only done some limited testing, but the
patch following my sig seems to work.
Note the first part of the diff allows us to remove the
cast to float in e_lgamma_r.c, which minimizes the diff
between e_lgamma_r.c and e_lgammaf_r.c (and e_lgammal_r.c :).
--
Steve
Index: src/e_lgamma_r.c
===================================================================
--- src/e_lgamma_r.c (revision 1427)
+++ src/e_lgamma_r.c (working copy)
@@ -173,19 +173,15 @@
*/
z = floor(y);
if(z!=y) { /* inexact anyway */
- y *= 0.5;
- y = 2.0*(y - floor(y)); /* y = |x| mod 2.0 */
- n = (int) (y*4.0);
+ y /= 2;
+ y = 2*(y - floor(y)); /* y = |x| mod 2.0 */
+ n = (int)(y*4);
} else {
- if(ix>=0x43400000) {
- y = zero; n = 0; /* y must be even */
- } else {
- if(ix<0x43300000) z = y+two52; /* exact */
- GET_LOW_WORD(n,z);
- n &= 1;
- y = n;
- n<<= 2;
- }
+ z = y+two52; /* exact */
+ GET_LOW_WORD(n,z);
+ n &= 1;
+ y = n;
+ n<<= 2;
}
switch (n) {
case 0: y = __kernel_sin(pi*y,zero,0); break;
Index: src/e_lgammaf_r.c
===================================================================
--- src/e_lgammaf_r.c (revision 1427)
+++ src/e_lgammaf_r.c (working copy)
@@ -106,19 +106,15 @@
*/
z = floorf(y);
if(z!=y) { /* inexact anyway */
- y *= (float)0.5;
- y = (float)2.0*(y - floorf(y)); /* y = |x| mod 2.0 */
- n = (int) (y*(float)4.0);
+ y /= 2;
+ y = 2*(y - floorf(y)); /* y = |x| mod 2.0 */
+ n = (int)(y*4);
} else {
- if(ix>=0x4b800000) {
- y = zero; n = 0; /* y must be even */
- } else {
- if(ix<0x4b000000) z = y+two23; /* exact */
- GET_FLOAT_WORD(n,z);
- n &= 1;
- y = n;
- n<<= 2;
- }
+ z = y+two23; /* exact */
+ GET_FLOAT_WORD(n,z);
+ n &= 1;
+ y = n;
+ n<<= 2;
}
switch (n) {
case 0: y = __kernel_sindf(pi*y); break;
More information about the freebsd-numerics
mailing list