This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
[PATCH] Demystify the magic number 134217729.0
- From: Siddhesh Poyarekar <siddhesh at redhat dot com>
- To: libc-alpha at sourceware dot org
- Date: Fri, 28 Dec 2012 20:25:14 +0530
- Subject: [PATCH] Demystify the magic number 134217729.0
Hi,
The number 134217729.0 gets used in various places in e_pow.c but
there is no explanation of what that number is. Attached patch adds
the explanation in dla.h and uses the macro in the code instead of the
bare value. Verified that there are no regressions resulting from
this. OK to commit?
Siddhesh
* sysdeps/ieee754/dbl-64/branred.h: Include dla.h.
(split): Use macro CN instead of the bare value.
* sysdeps/ieee754/dbl-64/dla.h: Add comment to explain why CN
could be used.
* sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Use CN
instead of the bare value.
(power1): Likewise.
diff --git a/sysdeps/ieee754/dbl-64/branred.h b/sysdeps/ieee754/dbl-64/branred.h
index c1e7354..2e4c69d 100644
--- a/sysdeps/ieee754/dbl-64/branred.h
+++ b/sysdeps/ieee754/dbl-64/branred.h
@@ -26,6 +26,7 @@
#ifndef BRANRED_H
#define BRANRED_H
+#include <dla.h>
#ifdef BIG_ENDI
static const mynumber
@@ -74,6 +75,6 @@ static const double toverp[75] = { /* 2/ PI base 24*/
12618859.0, 4703257.0, 12806093.0, 14477321.0, 2786137.0,
12875403.0, 9837734.0, 14528324.0, 13719321.0, 343717.0 };
-static const double split = 134217729.0;
+static const double split = CN; /* 2^27 + 1 */
#endif
diff --git a/sysdeps/ieee754/dbl-64/dla.h b/sysdeps/ieee754/dbl-64/dla.h
index c92fcfe..d3d8892 100644
--- a/sysdeps/ieee754/dbl-64/dla.h
+++ b/sysdeps/ieee754/dbl-64/dla.h
@@ -34,7 +34,8 @@
/* IEEE double. */
/***********************************************************************/
-/* CN = 1+2**27 = '41a0000002000000' IEEE double format */
+/* CN = 1+2**27 = '41a0000002000000' IEEE double format. Use it to split a
+ double for better accuracy. */
#define CN 134217729.0
diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c
index 5131718..1d535cb 100644
--- a/sysdeps/ieee754/dbl-64/e_pow.c
+++ b/sysdeps/ieee754/dbl-64/e_pow.c
@@ -95,10 +95,10 @@ __ieee754_pow(double x, double y) {
if (ABS (y) < 0x1p-64)
y = y < 0 ? -0x1p-64 : 0x1p-64;
z = log1(x,&aa,&error); /* x^y =e^(y log (X)) */
- t = y*134217729.0;
+ t = y*CN;
y1 = t - (t-y);
y2 = y - y1;
- t = z*134217729.0;
+ t = z*CN;
a1 = t - (t-z);
a2 = (z - a1)+aa;
a = y1*a1;
@@ -182,10 +182,10 @@ SECTION
power1(double x, double y) {
double z,a,aa,error, t,a1,a2,y1,y2;
z = my_log2(x,&aa,&error);
- t = y*134217729.0;
+ t = y*CN;
y1 = t - (t-y);
y2 = y - y1;
- t = z*134217729.0;
+ t = z*CN;
a1 = t - (t-z);
a2 = z - a1;
a = y*z;