[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gawk-diffs] [SCM] gawk branch, long-double, updated. 3e6f4db8c8724164a5
|
From: |
John Haque |
|
Subject: |
[gawk-diffs] [SCM] gawk branch, long-double, updated. 3e6f4db8c8724164a52b12216746c6bb3e050b69 |
|
Date: |
Mon, 04 Feb 2013 12:18:44 +0000 |
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "gawk".
The branch, long-double has been updated
via 3e6f4db8c8724164a52b12216746c6bb3e050b69 (commit)
from fddb32fc890cabd6175ab022302ca3d852a4731d (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
http://git.sv.gnu.org/cgit/gawk.git/commit/?id=3e6f4db8c8724164a52b12216746c6bb3e050b69
commit 3e6f4db8c8724164a52b12216746c6bb3e050b69
Author: John Haque <address@hidden>
Date: Mon Feb 4 06:17:09 2013 -0600
Added replacement trig functions for long double.
diff --git a/misc/ChangeLog b/misc/ChangeLog
index 4dd9c3f..d668726 100644
--- a/misc/ChangeLog
+++ b/misc/ChangeLog
@@ -1,3 +1,10 @@
+2013-02-04 John Haque <address@hidden>
+
+ * gawk_math.c (gawk_sinl, gawk_cosl, gawk_atan2l):
+ Replacement long double math functions.
+ * ldbl-tests/*.awk: Tests for the math functions.
+ * floatcmp.awk: Added diff mode.
+
2013-02-02 John Haque <address@hidden>
* float80.c: New file. Provides 80-bit long-double support
diff --git a/misc/Makefile b/misc/Makefile
index bafec60..2a07457 100644
--- a/misc/Makefile
+++ b/misc/Makefile
@@ -42,9 +42,14 @@ LDBL_TESTS = \
exp64 \
exp113 \
pow64 \
- pow113
+ pow113 \
+ sin64 \
+ sin113 \
+ cos64 \
+ cos113
INFILE = ldbl-tests/data.in
+INFILE32 = ldbl-tests/sincos.in
# An attempt to print something that can be grepped for in build logs
pass-fail:
@@ -61,6 +66,10 @@ ldbl-tests: $(LDBL_TESTS)
$(INFILE):
@$(AWK) -M -vPREC=quad -f ldblin.awk > $(INFILE) 2>&1
+$(INFILE32): $(INFILE)
+ @$(AWK) -M -vPREC=quad -f ldblin.awk | \
+$(AWK) -M -vPREC=quad '($$1 > -2^32 && $$1 < 2^32){ print $$0 }' > $(INFILE32)
2>&1
+
ldblint64:
@echo $@
@$(AWK) -B0 -f ldbl-tests/address@hidden > ldbl-tests/_$@ 2>&1
@@ -73,48 +82,72 @@ ldblint128:
sqrt64: $(INFILE)
@echo $@
- @$(AWK) -B0 -vDIG=17 -f ldbl-tests/sqrt.awk > ldbl-tests/_$@ $(INFILE)
2>&1
- @$(AWK) -M -vDIG=17 -vPREC=64 -f ldbl-tests/sqrt.awk >
ldbl-tests/address@hidden $(INFILE) 2>&1
+ @$(AWK) -B0 -vDIG=17 -f ldbl-tests/sqrt.awk $(INFILE) > ldbl-tests/_$@
2>&1
+ @$(AWK) -M -vDIG=17 -vPREC=64 -f ldbl-tests/sqrt.awk $(INFILE) >
ldbl-tests/address@hidden 2>&1
@-$(AWK) -M -f floatcmp.awk ldbl-tests/_$@ ldbl-tests/address@hidden &&
rm ldbl-tests/_$@
sqrt113: $(INFILE)
@echo $@
- @$(AWK) -B1 -vDIG=32 -f ldbl-tests/sqrt.awk > ldbl-tests/_$@ $(INFILE)
2>&1
- @$(AWK) -M -vDIG=32 -vPREC=113 -f ldbl-tests/sqrt.awk >
ldbl-tests/address@hidden $(INFILE) 2>&1
+ @$(AWK) -B1 -vDIG=32 -f ldbl-tests/sqrt.awk $(INFILE) > ldbl-tests/_$@
2>&1
+ @$(AWK) -M -vDIG=32 -vPREC=113 -f ldbl-tests/sqrt.awk $(INFILE) >
ldbl-tests/address@hidden 2>&1
@-$(AWK) -M -f floatcmp.awk ldbl-tests/_$@ ldbl-tests/address@hidden &&
rm ldbl-tests/_$@
log64: $(INFILE)
@echo $@
- @$(AWK) -B0 -vDIG=17 -f ldbl-tests/log.awk > ldbl-tests/_$@ $(INFILE)
2>&1
- @$(AWK) -M -vDIG=17 -vPREC=64 -f ldbl-tests/log.awk >
ldbl-tests/address@hidden $(INFILE) 2>&1
+ @$(AWK) -B0 -vDIG=17 -f ldbl-tests/log.awk $(INFILE) > ldbl-tests/_$@
2>&1
+ @$(AWK) -M -vDIG=17 -vPREC=64 -f ldbl-tests/log.awk $(INFILE) >
ldbl-tests/address@hidden 2>&1
@-$(AWK) -M -f floatcmp.awk ldbl-tests/_$@ ldbl-tests/address@hidden &&
rm ldbl-tests/_$@
log113: $(INFILE)
@echo $@
- @$(AWK) -B1 -vDIG=32 -f ldbl-tests/log.awk > ldbl-tests/_$@ $(INFILE)
2>&1
- @$(AWK) -M -vDIG=32 -vPREC=113 -f ldbl-tests/log.awk >
ldbl-tests/address@hidden $(INFILE) 2>&1
+ @$(AWK) -B1 -vDIG=32 -f ldbl-tests/log.awk $(INFILE) > ldbl-tests/_$@
2>&1
+ @$(AWK) -M -vDIG=32 -vPREC=113 -f ldbl-tests/log.awk $(INFILE) >
ldbl-tests/address@hidden 2>&1
@-$(AWK) -M -f floatcmp.awk ldbl-tests/_$@ ldbl-tests/address@hidden &&
rm ldbl-tests/_$@
exp64: $(INFILE)
@echo $@
- @$(AWK) -B0 -vDIG=17 -f ldbl-tests/exp.awk > ldbl-tests/_$@ $(INFILE)
2>&1
- @$(AWK) -M -vDIG=17 -vPREC=64 -f ldbl-tests/exp.awk >
ldbl-tests/address@hidden $(INFILE) 2>&1
+ @$(AWK) -B0 -vDIG=17 -f ldbl-tests/exp.awk $(INFILE) > ldbl-tests/_$@
2>&1
+ @$(AWK) -M -vDIG=17 -vPREC=64 -f ldbl-tests/exp.awk $(INFILE) >
ldbl-tests/address@hidden 2>&1
@-$(AWK) -M -f floatcmp.awk ldbl-tests/_$@ ldbl-tests/address@hidden &&
rm ldbl-tests/_$@
exp113: $(INFILE)
@echo $@
- @$(AWK) -B1 -vDIG=32 -f ldbl-tests/exp.awk > ldbl-tests/_$@ $(INFILE)
2>&1
- @$(AWK) -M -vDIG=32 -vPREC=113 -f ldbl-tests/exp.awk >
ldbl-tests/address@hidden $(INFILE) 2>&1
+ @$(AWK) -B1 -vDIG=32 -f ldbl-tests/exp.awk $(INFILE) > ldbl-tests/_$@
2>&1
+ @$(AWK) -M -vDIG=32 -vPREC=113 -f ldbl-tests/exp.awk $(INFILE) >
ldbl-tests/address@hidden 2>&1
@-$(AWK) -M -f floatcmp.awk ldbl-tests/_$@ ldbl-tests/address@hidden &&
rm ldbl-tests/_$@
pow64: $(INFILE)
@echo $@
- @$(AWK) -B0 -vDIG=17 -f ldbl-tests/pow.awk > ldbl-tests/_$@ $(INFILE)
2>&1
- @$(AWK) -M -vDIG=17 -vPREC=64 -f ldbl-tests/pow.awk >
ldbl-tests/address@hidden $(INFILE) 2>&1
+ @$(AWK) -B0 -vDIG=17 -f ldbl-tests/pow.awk $(INFILE) > ldbl-tests/_$@
2>&1
+ @$(AWK) -M -vDIG=17 -vPREC=64 -f ldbl-tests/pow.awk $(INFILE) >
ldbl-tests/address@hidden 2>&1
@-$(AWK) -M -vTOL=10 -f floatcmp.awk ldbl-tests/_$@
ldbl-tests/address@hidden && rm ldbl-tests/_$@
pow113: $(INFILE)
@echo $@
- @$(AWK) -B1 -vDIG=32 -f ldbl-tests/pow.awk > ldbl-tests/_$@ $(INFILE)
2>&1
- @$(AWK) -M -vDIG=32 -vPREC=113 -f ldbl-tests/pow.awk >
ldbl-tests/address@hidden $(INFILE) 2>&1
- @-$(AWK) -M -vTOL=10 -f floatcmp.awk ldbl-tests/_$@
ldbl-tests/address@hidden && rm ldbl-tests/_$@
+ @$(AWK) -B1 -vDIG=32 -f ldbl-tests/pow.awk $(INFILE) > ldbl-tests/_$@
2>&1
+ @$(AWK) -M -vDIG=32 -vPREC=113 -f ldbl-tests/pow.awk $(INFILE) >
ldbl-tests/address@hidden 2>&1
+ @-$(AWK) -M -vTOL=20 -f floatcmp.awk ldbl-tests/_$@
ldbl-tests/address@hidden && rm ldbl-tests/_$@
+
+sin64:
+ @echo $@
+ @$(AWK) -B0 -vDIG=17 -f ldbl-tests/sin.awk $(INFILE32) > ldbl-tests/_$@
2>&1
+ @$(AWK) -M -vDIG=17 -vPREC=64 -f ldbl-tests/sin.awk $(INFILE32) >
ldbl-tests/address@hidden 2>&1
+ @-$(AWK) -M -f floatcmp.awk ldbl-tests/_$@ ldbl-tests/address@hidden &&
rm ldbl-tests/_$@
+
+sin113: $(INFILE32)
+ @echo $@
+ @$(AWK) -B1 -vDIG=32 -f ldbl-tests/sin.awk $(INFILE32) > ldbl-tests/_$@
2>&1
+ @$(AWK) -M -vDIG=32 -vPREC=113 -f ldbl-tests/sin.awk $(INFILE32) >
ldbl-tests/address@hidden 2>&1
+ @-$(AWK) -M -f floatcmp.awk ldbl-tests/_$@ ldbl-tests/address@hidden &&
rm ldbl-tests/_$@
+
+cos64:
+ @echo $@
+ @$(AWK) -B0 -vDIG=17 -f ldbl-tests/cos.awk $(INFILE32) > ldbl-tests/_$@
2>&1
+ @$(AWK) -M -vDIG=17 -vPREC=64 -f ldbl-tests/cos.awk $(INFILE32) >
ldbl-tests/address@hidden 2>&1
+ @-$(AWK) -M -f floatcmp.awk ldbl-tests/_$@ ldbl-tests/address@hidden &&
rm ldbl-tests/_$@
+
+cos113: $(INFILE32)
+ @echo $@
+ @$(AWK) -B1 -vDIG=32 -f ldbl-tests/cos.awk $(INFILE32) > ldbl-tests/_$@
2>&1
+ @$(AWK) -M -vDIG=32 -vPREC=113 -f ldbl-tests/cos.awk $(INFILE32) >
ldbl-tests/address@hidden 2>&1
+ @-$(AWK) -M -f floatcmp.awk ldbl-tests/_$@ ldbl-tests/address@hidden &&
rm ldbl-tests/_$@
diff --git a/misc/floatcmp.awk b/misc/floatcmp.awk
index 87c8ea2..ac2a08e 100644
--- a/misc/floatcmp.awk
+++ b/misc/floatcmp.awk
@@ -1,9 +1,14 @@
-# floatcmp.awk --- check floating-point numbers for differences in
-# sig digs.
+# floatcmp.awk --- check floating-point numbers for differences in sig digs.
+# -vMODE=[diff|plot] -- diff mode or output data suitable fro plotting
+# -vTOL=tol -- ignore sig dig difference <= tol
BEGIN {
+ if (! ("mpfr_version" in PROCINFO)) {
+ print "floatcmp.awk: requires -M option" > "/dev/stderr"
+ exit(1)
+ }
if (ARGC < 3) {
- printf("Usage: gawk -M [-vTOL=1] -f floatcmp.awk file1 file2\n")
+ print "Usage: gawk -M [-vTOL=tol] [-vMODE=diff|plot] -f
floatcmp.awk file1 file2" > "/dev/stderr"
exit(1)
}
@@ -20,10 +25,11 @@ BEGIN {
ret1 = (getline v1 < file1)
ret2 = (getline v2 < file2)
+ f1 = v1; f2 = v2; # save originals for diff mode
if (ret1 > 0 && ret2 > 0) {
line++
if ((v1 "") == (v2 ""))
- continue;
+ continue
e1 = index(v1, "e")
e2 = index(v2, "e")
if (e1 > 0 && e2 > 0 && # exclude nans and
infinities
@@ -41,8 +47,17 @@ BEGIN {
continue
}
- printf("%s %s differ: byte ?, line %d\n", file1, file2,
line)
- exit(1)
+ if (! MODE) {
+ printf("%s %s differ: byte ?, line %d\n",
file1, file2, line)
+ exit(1)
+ }
+ if (MODE == "plot")
+ printf("%d\t%d\n", line, diff)
+ else {
+ dig = length(v1)
+ printf("%-*.*s %-*.*s %+*.*d\n", dig+7, dig+7,
f1, dig+7, dig+7, f2, dig, dig, diff)
+ }
+ continue
}
if (ret1 == 0 && ret2 == 0)
diff --git a/misc/fp_math.awk b/misc/fp_math.awk
index e2847e6..60048b1 100644
--- a/misc/fp_math.awk
+++ b/misc/fp_math.awk
@@ -39,7 +39,7 @@ function euler_atan_one_over(x, \
do {
term *= (i + 2) / (i + 3)
- err = term /= xpow2_plus_one
+ term /= xpow2_plus_one
i += 2
sum += term
err = term / sum
diff --git a/misc/gawk_math.c b/misc/gawk_math.c
index 6981725..27ca215 100644
--- a/misc/gawk_math.c
+++ b/misc/gawk_math.c
@@ -27,40 +27,250 @@
#define _0L LDC(0.0)
#define _1L LDC(1.0)
#define _2L LDC(2.0)
+#define _3L LDC(3.0)
+#define _4L LDC(4.0)
/*
* Constants for computation using long doubles with enough digits for the
128-bit quad.
*/
-#define GAWK_LOG2 LDC(0.693147180559945309417232121458176568) /* log 2
(base e) */
+#define GAWK_LOG2 LDC(0.693147180559945309417232121458176568) /*
log 2 (base e) */
#define GAWK_LOG2_HIGH LDC(0.6931471801362931728363037109375) /*
high 32 bits (exact representation) */
-#define GAWK_LOG2_LOW LDC(4.236521365809284105206765680755001344e-10) /*
variable precision low bits */
+#define GAWK_LOG2_LOW LDC(4.236521365809284105206765680755001344e-10)
/* variable precision low bits */
-#define GAWK_SQRT2 LDC(1.414213562373095048801688724209698079) /*
sqrt(2) */
-#define GAWK_SQRT1_2 LDC(0.707106781186547524400844362104849039) /*
1/sqrt(2) */
+#define GAWK_SQRT2 LDC(1.414213562373095048801688724209698079)
/* sqrt(2) */
+#define GAWK_SQRT1_2 LDC(0.707106781186547524400844362104849039)
/* 1/sqrt(2) */
+#define GAWK_PI LDC(3.141592653589793238462643383279502884)
/* pi */
+#define GAWK_PI_2 LDC(1.570796326794896619231321691639751442)
/* pi/2 */
+#define GAWK_PI_4 LDC(0.785398163397448309615660845819875721)
/* pi/4 */
+#define GAWK_PI_4_HIGH LDC(0.78539816313423216342926025390625)
/* 32-bit 1st part */
+#define GAWK_PI_4_MED LDC(2.632161460363116600724708860070677474e-10)
/* 32-bit 2nd part */
+#define GAWK_PI_4_LOW LDC(1.500889318411548350422246024296643642e-19)
/* variable precision 3rd part */
static AWKLDBL taylor_exp(AWKLDBL x);
+static AWKLDBL taylor_cos(AWKLDBL x);
+static AWKLDBL taylor_sin(AWKLDBL x);
+static AWKLDBL euler_atan_1(AWKLDBL x);
static AWKLDBL gawk_frexpl(AWKLDBL x, int *exponent);
+/* gawk_sinl --- Compute sin(x) */
+
static AWKLDBL
gawk_sinl(AWKLDBL x)
{
- return sin( (double) x);
+ AWKLDBL sinval, dq;
+ int sign = 1;
+
+ if (isinf(x) || isnan(x))
+ return GAWK_NAN;
+
+ if (x == _0L) /* XXX: sin(-0.0) = -0.0 */
+ return x;
+
+ if (x < _0L) {
+ /* sin(-x) = - sin(x) */
+ sign = -1;
+ x = -x;
+ }
+
+ /* range reduction -- 0 <= x < pi / 4 */
+
+ dq = double_to_int(x / GAWK_PI_4);
+
+ if (dq < pow2ld(32)) {
+ gawk_uint_t q = (gawk_uint_t) dq;
+
+ sign *= ( (q / 4) % 2 ? -1 : 1 );
+ switch (q % 4) {
+ case 0:
+ x -= q * GAWK_PI_4_HIGH;
+ x -= q * GAWK_PI_4_MED;
+ x -= q * GAWK_PI_4_LOW;
+ sinval = taylor_sin(x);
+ break;
+ case 1:
+ q++;
+ x -= q * GAWK_PI_4_HIGH;
+ x -= q * GAWK_PI_4_MED;
+ x -= q * GAWK_PI_4_LOW;
+ sinval = taylor_cos(-x);
+ break;
+ case 2:
+ x -= q * GAWK_PI_4_HIGH;
+ x -= q * GAWK_PI_4_MED;
+ x -= q * GAWK_PI_4_LOW;
+ sinval = taylor_cos(x);
+ break;
+ case 3:
+ q++;
+ x -= q * GAWK_PI_4_HIGH;
+ x -= q * GAWK_PI_4_MED;
+ x -= q * GAWK_PI_4_LOW;
+ sinval = taylor_sin(-x);
+ break;
+ default:
+ break;
+ }
+ } else {
+ /* FIXME -- Payne and Hanek Reduction Algorithm ? */
+
+ sinval = (AWKLDBL) sin( (double) x);
+ }
+
+ return sign > 0 ? sinval : -sinval;
}
+/* gawk_cosl --- Compute cos(x) */
+
static AWKLDBL
gawk_cosl(AWKLDBL x)
{
- return cos( (double) x);
+ AWKLDBL cosval, dq;
+ int sign = 1;
+
+ if (isinf(x) || isnan(x))
+ return GAWK_NAN;
+ if (x < _0L) {
+ /* cos(-x) = cos(x) */
+ x = -x;
+ }
+ if (x == _0L)
+ return _1L;
+
+ /* range reduction -- 0 <= x < pi / 4 */
+
+ dq = double_to_int(x / GAWK_PI_4);
+
+ if (dq < pow2ld(32)) {
+ gawk_uint_t q = (gawk_uint_t) dq;
+
+ sign *= ( (q / 4) % 2 ? -1 : 1 );
+ switch (q % 4) {
+ case 0:
+ x -= q * GAWK_PI_4_HIGH;
+ x -= q * GAWK_PI_4_MED;
+ x -= q * GAWK_PI_4_LOW;
+ cosval = taylor_cos(x);
+ break;
+ case 1:
+ q++;
+ x -= q * GAWK_PI_4_HIGH;
+ x -= q * GAWK_PI_4_MED;
+ x -= q * GAWK_PI_4_LOW;
+ cosval = taylor_sin(-x);
+ break;
+ case 2:
+ x -= q * GAWK_PI_4_HIGH;
+ x -= q * GAWK_PI_4_MED;
+ x -= q * GAWK_PI_4_LOW;
+ cosval = -taylor_sin(x);
+ break;
+ case 3:
+ q++;
+ x -= q * GAWK_PI_4_HIGH;
+ x -= q * GAWK_PI_4_MED;
+ x -= q * GAWK_PI_4_LOW;
+ cosval = -taylor_cos(-x);
+ break;
+ default:
+ break;
+ }
+ } else {
+ /* FIXME Payne and Hanek Reduction Algorithm ? */
+
+ cosval = (AWKLDBL) cos( (double) x);
+ }
+
+ return sign > 0 ? cosval : -cosval;
}
+
+/* gawk_atan2l(x) --- Compute atan2(y, x) */
+
static AWKLDBL
gawk_atan2l(AWKLDBL y, AWKLDBL x)
{
- return atan2( (double) y, (double) x);
+ int sign;
+
+ /*
+ * Using Euler atan(1/x) and the identity:
+ * atan(x) = - pi / 2 - atan(1/x), x < 0
+ * = pi / 2 - atan(1/x), x > 0
+ */
+
+ if (isnan(x) || isnan(y))
+ return GAWK_NAN;
+
+ if (isinf(y)) {
+ if (y > _0L) {
+ if (isinf(x))
+ return x < _0L ? _3L * GAWK_PI_4 : GAWK_PI_4;
+ return GAWK_PI_2;
+ }
+ if (isinf(x))
+ return x < _0L ? -_3L * GAWK_PI_4 : -GAWK_PI_4;
+ return -GAWK_PI_2;
+ }
+
+ if (isinf(x)) {
+ if (x > _0L)
+ return y < _0L ? -_0L : _0L;
+ /* x = -Infinity */
+ return (y >= _0L) ? GAWK_PI : -GAWK_PI;
+ }
+
+ if (x > _0L) {
+ if (y == _0L)
+ return y; /* +0 or -0 */
+ sign = 1;
+ if (y < _0L) {
+ sign = -1;
+ y = -y;
+ }
+ if (y > x)
+ return sign * (GAWK_PI_2 - euler_atan_1(y / x));
+ return sign * euler_atan_1(x / y);
+ }
+
+ if (x < _0L) {
+ if (y == _0L)
+ return atan2( (double) y, -_1L) < 0.0 ? -GAWK_PI :
GAWK_PI;
+
+ x = -x;
+ if (y < _0L) {
+ y = -y;
+ if (y > x)
+ return -GAWK_PI_2 - euler_atan_1(y / x);
+ return euler_atan_1(x / y) - GAWK_PI;
+ }
+ /* y > 0 */
+ if (y > x)
+ return GAWK_PI_2 + euler_atan_1(y / x);
+ return - euler_atan_1(x / y) + GAWK_PI;
+ }
+
+ /*
+ * XXX -- using atan2 instead of trying to detect +0 or -0. No need for
signbit(x)!
+ * We need to make sure the value of y is in the double range;
+ * Replace positive (negative) y with 1 (-1), the actual value isn't
important.
+ */
+
+ if ( y > _0L)
+ y = _1L;
+ else if (y < _0L)
+ y = -_1L;
+
+ /* x = +0 or -0 */
+ if (atan2( (double) y, (double) x) < 0.0)
+ return -GAWK_PI;
+ if (atan2( (double) y, (double) x) > 0.0)
+ return GAWK_PI;
+ /* x = +0 and y = +0 or -0 */
+ return _0L;
}
+
/* gawk_logl --- Compute log(x) */
static AWKLDBL
@@ -459,7 +669,7 @@ taylor_exp(AWKLDBL x)
return _1L;
k = 1;
- while (x > 0.001) {
+ while (x > LDC(0.001)) {
/* XXX: For x <= 0.001, max(k) = 10, and max # of terms 6
(80-bit) / 10 (128-bit) */
x /= _2L;
@@ -482,3 +692,117 @@ taylor_exp(AWKLDBL x)
y = 2 * y + y * y;
return y + _1L;
}
+
+/* taylor_sin --- Compute sin(x) using Taylor series */
+
+static AWKLDBL
+taylor_sin(AWKLDBL x)
+{
+ AWKLDBL xpow, xpow_odd, sinval;
+ AWKLDBL err, term, fact;
+ unsigned int i;
+
+ assert(x >= _0L);
+
+ if (x == _0L)
+ return x;
+
+ i = 3;
+ fact = LDC(6.0); /* 3! */
+ xpow = x * x;
+ xpow_odd = xpow * x;
+ sinval = x - xpow_odd / fact;
+
+ do {
+ fact *= (AWKLDBL) ((i + 1) * (i + 2));
+ i += 2;
+ xpow_odd *= xpow;
+ term = xpow_odd / fact;
+
+ fact *= (AWKLDBL) ((i + 1) * (i + 2));
+ i += 2;
+ xpow_odd *= xpow;
+ term -= xpow_odd / fact;
+
+ sinval += term;
+ err = term / sinval;
+ } while (err > REL_ERROR);
+ return sinval;
+}
+
+/* taylor_cos --- Compute cos(x) using Taylor series */
+
+static AWKLDBL
+taylor_cos(AWKLDBL x)
+{
+ AWKLDBL xpow, xpow_even, cosval;
+ AWKLDBL err, term, fact;
+ unsigned int i;
+
+ if (x == _0L)
+ return _1L;
+
+ i = 2;
+ fact = _2L; /* 2! */
+ xpow = x * x;
+ xpow_even = xpow;
+ cosval = _1L - xpow / fact;
+
+ do {
+ fact *= (AWKLDBL) ((i + 1) * (i + 2));
+ i += 2;
+ xpow_even *= xpow;
+ term = xpow_even / fact;
+
+ fact *= (AWKLDBL) ((i + 1) * (i + 2));
+ i += 2;
+ xpow_even *= xpow;
+ term -= xpow_even / fact;
+
+ cosval += term;
+ err = term / cosval;
+ } while (err > REL_ERROR);
+ return cosval;
+}
+
+/* euler_atan_1 --- Compute Euler arctan(1/x) approximation */
+
+static AWKLDBL
+euler_atan_1(AWKLDBL x)
+{
+ AWKLDBL xpow2_plus_one, term, sum, err;
+ int sign = 1;
+ unsigned int i;
+
+ /*
+ * Euler atan formula
(http://mathworld.wolfram.com/InverseTangent.html):
+ * atan(x) = (y/x) (1 + 2/3 y + (2·4)/(3·5) y^2 +
(2·4·6)/(3·5·7) y^3 + ...)
+ * where
+ * y = (x^2) / (1 + x^2) and -1 <= x <= 1
+ *
+ * Substituting x = 1/x, for x >= 1
+ * atan(1/x) = (x / (1 + x^2)) + (2/3) * (x / (1 + x^2)^2)
+ * + (2*4/(3*5)) * (x / (1 + x^2)^3)
+ * + (2*4*6/(3*5*7)) * (x / (1 +
x^2)^4) + ...
+ */
+
+ if (x < _0L) {
+ sign = -1;
+ x = -x;
+ }
+
+ xpow2_plus_one = x * x + _1L;
+ term = x / xpow2_plus_one;
+ sum = term;
+ i = 0;
+
+ do {
+ term *= (AWKLDBL) (i + 2);
+ term /= ((AWKLDBL) (i + 3)) * xpow2_plus_one;
+ i += 2;
+ sum += term;
+ err = term / sum;
+ } while (err > REL_ERROR);
+
+ return sign > 0 ? sum : -sum;
+}
diff --git a/misc/ldbl-tests/pow.awk b/misc/ldbl-tests/pow.awk
index 0e75b56..73cf492 100644
--- a/misc/ldbl-tests/pow.awk
+++ b/misc/ldbl-tests/pow.awk
@@ -16,14 +16,8 @@ BEGIN {
PREC = "quad"
y1 += 0; y2 += 0; y3 += 0
- if (y1 <= 1.0e750) {
- printf("%*.*e\n", 0, DIG, y1)
- }
- if (y2 <= 1.0e750) {
- printf("%*.*e\n", 0, DIG, y2)
- }
- if (y3 <= 1.0e750) {
- printf("%*.*e\n", 0, DIG, y3)
- }
+ printf("%*.*e\n", 0, DIG, y1)
+ printf("%*.*e\n", 0, DIG, y2)
+ printf("%*.*e\n", 0, DIG, y3)
PREC = save_PREC
}
-----------------------------------------------------------------------
Summary of changes:
misc/ChangeLog | 7 +
misc/Makefile | 69 +++++++---
misc/floatcmp.awk | 27 +++-
misc/fp_math.awk | 2 +-
misc/gawk_math.c | 340 +++++++++++++++++++++++++++++++++++++++++++++-
misc/ldbl-tests/pow.awk | 12 +--
6 files changed, 415 insertions(+), 42 deletions(-)
hooks/post-receive
--
gawk
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gawk-diffs] [SCM] gawk branch, long-double, updated. 3e6f4db8c8724164a52b12216746c6bb3e050b69,
John Haque <=