[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gawk-diffs] [SCM] gawk branch, long-double, updated. 88304b52a0dc79b4b3
|
From: |
John Haque |
|
Subject: |
[gawk-diffs] [SCM] gawk branch, long-double, updated. 88304b52a0dc79b4b33b038acb6f20aaacf01092 |
|
Date: |
Wed, 16 Jan 2013 11:39:39 +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 88304b52a0dc79b4b33b038acb6f20aaacf01092 (commit)
from 5f4dbd4819a7dae080c6ff4833a06cf1ca8d502c (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=88304b52a0dc79b4b33b038acb6f20aaacf01092
commit 88304b52a0dc79b4b33b038acb6f20aaacf01092
Author: John Haque <address@hidden>
Date: Wed Jan 16 05:38:11 2013 -0600
Add files sprintf*_0Lf.awk, bin2dec.awk.
diff --git a/TODO.LDBL b/TODO.LDBL
index 3e2732f..f81fe53 100644
--- a/TODO.LDBL
+++ b/TODO.LDBL
@@ -9,10 +9,6 @@ Also, append suffix L to all long double constants.
* Convert gawk math routines to C for platforms lacking any. Without ceil()
and floor()
for long doubles, these can't even get 64 bit integers!
-- finish fp_fmod() in misc/fp_math.awk.
- -- write test script(s) to verify (some) output, but DON'T add it to
the dist.
- 80 bit vs 128 bit, and the C math library CANNOT be trusted to produce
- correctly rounded results. Double check any discrepencies using
arbitrary-precision
- math.
* Don't use adjust_uint(uintmax_t n) from floatcomp.c, it is for AWKNUM. What
is the point
of floor() and ceil() wrappers? Don't have a clue.
diff --git a/configure b/configure
index 4bf106b..0bb8f5f 100755
--- a/configure
+++ b/configure
@@ -5554,7 +5554,7 @@ $as_echo "yes" >&6; }
else
{ $as_echo "$as_me:${as_lineno-$LINENO}: result: no" >&5
$as_echo "no" >&6; }
- CFLAGS="$CFLAGS -DNDEBUG" # turn off assertions
+ CFLAGS="$CFLAGS" # DONOT TURN OF ASSERTIONS
fi
diff --git a/configure.ac b/configure.ac
index d3780bd..e8686e0 100644
--- a/configure.ac
+++ b/configure.ac
@@ -95,7 +95,7 @@ then
AC_MSG_RESULT([yes])
else
AC_MSG_RESULT([no])
- CFLAGS="$CFLAGS -DNDEBUG" # turn off assertions
+ CFLAGS="$CFLAGS" # DONOT TURN OF ASSERTIONS
fi
AC_SUBST(CFLAGS)
diff --git a/misc/bin2dec.awk b/misc/bin2dec.awk
new file mode 100644
index 0000000..04a82e1
--- /dev/null
+++ b/misc/bin2dec.awk
@@ -0,0 +1,83 @@
+# bcd.awk --- Calculate the magic numbers in:
+# Binary to Decimal Conversion in Limited Precision - Douglas B.
Jones
+# URL: homepage.cs.uiowa.edu/~jones/bcd/decimal.html
+#
+# Usage: ./gawk -M -f bin2dec.awk
+#
+
+BEGIN {
+ UINT32_MAX = 4294967296
+ UINT64_MAX = 18446744073709551616
+
+ N = 4 # number of chunks
+
+## 64-bit integers
+# CHUNK_SIZE = CHUNK3_SIZE = 16
+# DEC_BASE = 10000 # 4 decimal digits per chunk, 20 total
+
+## 96-bit integers
+# CHUNK_SIZE = CHUNK3_SIZE = 24
+# DEC_BASE = 10000000 # 7 decimal digits per chunk, 35 total
+
+# 128-bit integers
+ CHUNK_SIZE = 32
+ CHUNK3_SIZE = 17 # the useful size of the first (MSB) chunk <=
CHUNK_SIZE
+ # 17 + 32 + 32 + 32 = 113, precision of 128-bit
double
+ DEC_BASE = 10000000 # 7 decimal digits per chunk, 35 total
+
+######################################################
+
+ CHUNK_MAX = 2^CHUNK_SIZE
+ CHUNK3_MAX = 2^CHUNK3_SIZE
+
+ print "CHUNK_SIZE = ", CHUNK_SIZE
+ print "N = ", N
+ print "DEC_BASE = ", DEC_BASE
+
+ # powers of 2^CHUNK_SIZE
+ POWER_2_C[0] = 1
+ for (i = 1; i < N; i++)
+ POWER_2_C[i] = 2 ^(CHUNK_SIZE * i)
+
+ # powers of DEC_BASE
+ POWER_BASE[0] = 1
+ for (i = 1; i < N; i++)
+ POWER_BASE[i] = DEC_BASE^i
+
+ coeff[N-1][N-1] = int(POWER_2_C[N-1] / POWER_BASE[N-1])
+ for (i = N-2; i >= 0; i--) {
+ for (k = N-1; k >= i; k--) {
+ coeff[i][k] = int(POWER_2_C[k]/POWER_BASE[i])
+ for (m = 1; m <= k - i; m++)
+ coeff[i][k] -= POWER_BASE[m] * coeff[m + i][k]
+ }
+ }
+
+ for (i = N - 1; i >= 0; i--) {
+ for (j = N - 1; j >= i; j--) {
+ printf("%d ", coeff[i][j])
+ }
+ print "";
+ }
+
+ print "----- a0, a1, a2, a4 upper bounds -----"
+ print "UINT32_MAX =", UINT32_MAX
+ print "UINT64_MAX =", UINT64_MAX
+ for (i = N - 1; i >= 0; i--) {
+ a[i] = 0;
+ for (j = N - 1; j >= i; j--) {
+ if (j == 3)
+ a[i] += coeff[i][j] * CHUNK3_MAX
+ else
+ a[i] += coeff[i][j] * CHUNK_MAX
+ }
+ printf("a%d <= %d\n", i, a[i])
+ }
+
+ print "----- c1, c2, c3 upper bounds -----"
+ c[0] = 0
+ for (i = 1; i < N; i++) {
+ c[i] = (a[i - 1] + c[i - 1])/ DEC_BASE
+ printf("c%d <= %d\n", i, c[i])
+ }
+}
diff --git a/misc/sprintf128_0Lf.awk b/misc/sprintf128_0Lf.awk
new file mode 100644
index 0000000..e3e61a9
--- /dev/null
+++ b/misc/sprintf128_0Lf.awk
@@ -0,0 +1,308 @@
+#
+# sprintf128_0Lf.awk -- format integers stored as long doubles and
+# don't have %Lf in (s)printf.
+# This test script handles integers < 2^113 i.e. 128-bit long doubles.
+# Requires -M -vPREC=113 (or -B if have 128-bit long doubles)
+#
+# Usage:
+#
+# $ echo "2^113-1" | ./gawk -M -vPREC=113 -f sprintf128_0Lf.awk
+# 2^113-1 (1.03845937170696552570609926584401910e+34) is an integer
+# sprintf_int128 :10384593717069655257060992658440191
+# arbitrary-prec :10384593717069655257060992658440191
+# floor = 10384593717069655257060992658440191
+# ceil = 10384593717069655257060992658440191
+# $ echo "2^113" | ./gawk -M -vPREC=113 -f sprintf128_0Lf.awk
+# 2^113 (1.03845937170696552570609926584401920e+34) is an integer
+# sprintf_int128 :10384593717069655257060992658440192
+# arbitrary-prec :10384593717069655257060992658440192
+# floor = 10384593717069655257060992658440192
+# ceil = 10384593717069655257060992658440192
+# $ echo "2^113+1" | ./gawk -M -vPREC=113 -f sprintf128_0Lf.awk
+# 2^113+1 (1.03845937170696552570609926584401920e+34) is an integer
+# sprintf_int128 :10384593717069655257060992658440192 ***
+# arbitrary-prec :10384593717069655257060992658440193
+# floor = 10384593717069655257060992658440192 ***
+# ceil = 10384593717069655257060992658440192 ***
+#
+
+BEGIN {
+ # Only need 113 for 128-bit long double.
+
+ LDBL_MANT_DIG = 128
+ POW2[0] = 1
+ for (i = 1; i <= LDBL_MANT_DIG; i++)
+ POW2[i] = POW2[i - 1] * 2
+}
+
+{
+ #
+ # Compare apples to apples, setting PREC=128;
+ #
+
+ eval = sprintf("./gawk -M -vPREC=128 'BEGIN { printf(\"%%0.40e\", %s);
print \"\";\
+ printf(\"%%035d\", %s); print \"\";\
+ printf(\"%%d\", int(%s) == %s + 0); print \"\"}'", $0, $0, $0,
$0)
+ eval | getline num
+ eval | getline ap_intval
+ eval | getline ap_isint
+ close(eval)
+
+ pnum = num + 0
+ sign = ""
+ if (pnum < 0) {
+ pnum = -pnum
+ sign = "-"
+ }
+
+ isint = isinteger(pnum, W)
+ if (isint) {
+ if (! ap_isint)
+ print "!!!!!!!!!MPFR says it is not an integer!!!!!!!!!"
+ printf("%s (%s) is an integer\n", $0, sprintf("%0.35e", num))
+ printf("sprintf_int128 :%c%s\n", sign, sprintf_int128(W))
+ printf("arbitrary-prec :%s\n", ap_intval)
+ } else {
+ if (ap_isint)
+ print "!!!!!!!!!MPFR says it is an integer!!!!!!!!!"
+ printf("%s (%s) is not an integer\n", $0, sprintf("%0.35e",
num))
+ }
+
+#
+# Don't have %Lf !
+# printf("floor = %d\n", floor(num))
+# printf("ceil = %d\n", ceil(num))
+#
+ # use MPFR to format floor(), ceil() returned value
+ eval = sprintf("./gawk -M -vPREC=128 'BEGIN { printf(\"%%d\", %s);
print \"\";\
+ printf(\"%%d\", %s); print \"\"; }'", floor(num), ceil(num))
+ eval | getline floor_val
+ eval | getline ceil_val
+ close(eval)
+
+ print "floor =", floor_val
+ print "ceil =", ceil_val
+}
+
+
+function sprintf_int128(chunk, d0, d1, d2, d3, d4, q, k, dec_base)
+{
+ #
+ # Binary to Decimal Conversion in Limited Precision - Douglas B. Jones
+ # URL: homepage.cs.uiowa.edu/~jones/bcd/decimal.html
+ #
+
+ #
+ # Results from the application of the idea in the article to the
128-bit conversion
+ # problem:
+ #
+ # Divide a 128-bit number into 32-bit chunks, and use 10,000,000 as the
+ # working decimal base. This gives us 7 digits per chunk, and 35 total
decimal digits.
+ # Using the notations used in the article,
+ #
+ # a3 = 79228162 n3
+ # a2 = 5142643 n3 + 184467 n2
+ # a1 = 3400832 n3 + 4407370 n2 + 429 n1
+ # a0 = 3585223950336 n3 + 9551616 n2 + 4967296 n1 + 1 n0
+ #
+ # n3, n2, n1 and n0 are the 32-bit chunks, n3 represents the most
+ # significant 32 bits. Assuming these chunks are at their maximum
possible
+ # value, which is 2^32:
+ #
+ # UINT64_MAX = 18446744073709551616
+ # a3 <= 340282364712189952
+ # a2 <= 22879763232194560
+ # a1 <= 33537814771531776
+ # a0 <= 15398481973785556680704
+ # c1 <= 1539848197378555
+ # c2 <= 3507766296
+ # c3 <= 2287976673
+ #
+ # c1, c2 and c3 are carries from one digit position to the next
+ # (results from integer divisions by the decimal base; see below).
+ #
+ # The maximum value of a0 is > UINT64_MAX. If we restrict n3 to be
+ # in the range 0 < n3 < 2^17, corresponding to the 113-bit
precision available
+ # from a 128-bit long double, we instead get these upper bounds:
+ #
+ # UINT64_MAX = 18446744073709551616
+ # a3 <= 10384593649664
+ # a2 <= 792953788694528
+ # a1 <= 18931798306193408
+ # a0 <= 532280730126909440
+ # c1 <= 53228073012
+ # c2 <= 1893185153
+ # c3 <= 79295568
+ #
+ # We can use 64-bit signed or unsigned long type! QED.
+ #
+ # All these numbers calculated using the script misc/bin2dec.awk.
+ #
+
+
+ dec_base = 10000000
+
+ if (DEBUG) {
+ for (k = 0; k < 4; k++) {
+ if (chunk[k] > 4294967296)
+ printf("chunk[%d] > 4294967296, this ain't
right\n", k)
+ }
+ print "(32 MSB bits first):", chunk[3], "|", chunk[2], "|",
chunk[1], "|", chunk[0]
+ }
+
+ d0 = 3585223950336 * chunk[3] + 9551616 * chunk[2] + 4967296 *
chunk[1] + chunk[0];
+
+ q = int(d0 / dec_base); # int() may call floorl(), but double
version is good enough
+ # in this case. This is only a consideration in
the awk code,
+ # not in the C translation.
+
+ d0 = d0 - q * dec_base; # not using % operator, may not have
fmodl()
+
+ d1 = q + 3400832 * chunk[3] + 4407370 * chunk[2] + 429 * chunk[1];
+ q = int(d1 / dec_base);
+ d1 = d1 - q * dec_base;
+
+ d2 = q + 5142643 * chunk[3] + 184467 * chunk[2];
+ q = int(d2 / dec_base);
+ d2 = d2 - q * dec_base;
+
+ d3 = q + 79228162 * chunk[3];
+ q = int(d3 / dec_base);
+ d3 = d3 - q * dec_base;
+
+ d4 = q;
+ return sprintf("%7.7u%7.7u%7.7u%7.7u%7.7u", \
+ d4, d3, d2, d1, d0); # don't need %Lf support here
+}
+
+function floor(x, d, sign, chunk, isint)
+{
+ x = +x
+
+#
+# if (isnan(x) || isinf(x) || x == 0)
+# return x
+
+ sign = 1
+ if (x < 0) {
+ sign = -1
+ x = -x
+ }
+
+ #
+ # For x >= 2^LDBL_MANT_DIG, use floor for C-double. Bailing out here.
+ #
+
+ if (x >= POW2[LDBL_MANT_DIG]) {
+ printf("number %0.35e is too big (> %0.35e), bailing out\n", x,
POW2[LDBL_MANT_DIG] - 1)
+ exit(1);
+ }
+
+ isint = isinteger(x, chunk)
+
+ # d is long double type
+ d = chunk[3] * POW2[96] + chunk[2] * POW2[64] + chunk[1] * POW2[32] +
chunk[0]
+
+ if (isint) # d == x
+ return sign * d
+ return (sign > 0) ? d : -d - 1
+}
+
+
+function ceil(x, d, sign, chunk, isint)
+{
+ x = +x
+
+#
+# if (isnan(x) || isinf(x) || x == 0)
+# return x
+
+ sign = 1
+ if (x < 0) {
+ sign = -1
+ x = -x
+ }
+
+ #
+ # For x >= 2^LDBL_MANT_DIG, use ceil for C-double. Bailing out here.
+ #
+
+ if (x >= POW2[LDBL_MANT_DIG]) {
+ printf("number %0.35e is too big (> %0.35e), bailing out\n", x,
POW2[LDBL_MANT_DIG] - 1)
+ exit(1);
+ }
+
+ isint = isinteger(x, chunk)
+ # d is long double type
+ d = chunk[3] * POW2[96] + chunk[2] * POW2[64] + chunk[1] * POW2[32] +
chunk[0]
+
+ if (isint) # d == x
+ return sign * d
+ return (sign > 0) ? d + 1 : -d
+}
+
+
+function isinteger(x, D, j, low, high)
+{
+ # x >= 0
+
+ if (x >= POW2[LDBL_MANT_DIG]) {
+ printf("number %0.40e is too big (> %0.40e), bailing out\n", x,
2^LDBL_MANT_DIG - 1)
+ exit(1);
+ }
+
+ D[0] = D[1] = D[2] = D[3] = 0
+
+ # binary search
+ high = LDBL_MANT_DIG - 1
+ while (x >= 2) {
+ low = 0
+ while (low <= high) {
+ mid = int((low + high) / 2)
+ if (x < POW2[mid])
+ high = mid - 1
+ else
+ low = mid + 1
+ }
+
+ if (x < POW2[low]) {
+ x -= POW2[low - 1]
+
+ #
+ # Easier to seperate the 128 bits into 32-bit chunks now
+ # then bit fiddling later.
+ # | D[3] | D[2] | D[1] | D[0] |
+ # |<------- x (128 bits) ----->|
+ #
+
+ if (low <= 32) { D[0] += POW2[low - 1] }
+ else if (low <= 64) { D[1] += POW2[low - 33] }
+ else if (low <= 96) { D[2] += POW2[low - 65] }
+ else { D[3] += POW2[low - 97] }
+ }
+
+ high = low
+ }
+
+#
+# # XXX -- O(n^2) vs. O(n * log(n)) -- how much do we really gain?
+# # Consider extra additions and divisions, but lower # of <
comparisons.
+#
+# while (x >= 2) {
+# for (j = 1; j <= LDBL_MANT_DIG; j++) {
+# if (x < POW2[j]) {
+# x -= POW2[j - 1]
+# break
+# }
+# }
+# }
+#
+
+ if (x >= 1) {
+ D[0]++
+ x -= 1
+ }
+
+ return (x == 0)
+}
+
diff --git a/misc/sprintf64_0Lf.awk b/misc/sprintf64_0Lf.awk
new file mode 100644
index 0000000..2ec38ac
--- /dev/null
+++ b/misc/sprintf64_0Lf.awk
@@ -0,0 +1,247 @@
+#
+# sprintf64_0Lf.awk -- format integers stored as long doubles and
+# don't have %Lf in (s)printf. This test script handles integers < 2^64.
+# Requires both -M (for verification) and -B
+# Usage:
+# echo "2^63 + 2^54 + 2^13" | ./gawk -B -f misc/sprintf_0Lf.awk
+#
+
+BEGIN {
+ # 113 for 128-bit long double, but can only get 64-bit ints.
+ LDBL_MANT_DIG = 64
+ POW2[0] = 1
+ for (i = 1; i <= LDBL_MANT_DIG; i++)
+ POW2[i] = POW2[i - 1] * 2
+}
+
+{
+ #
+ # Compare apples to apples, setting PREC=64; With real 128-bit long
doubles,
+ # one probably should use PREC=quad or 113.
+ # %0.34e ain't the right format but then ..
+ #
+
+ eval = sprintf("./gawk -M -vPREC=64 'BEGIN { printf(\"%%0.34e\", %s);
print \"\";\
+ printf(\"%%020d\", %s); print \"\";\
+ printf(\"%%d\", int(%s) == %s + 0); print \"\"}'", $0, $0, $0,
$0)
+ eval | getline num
+ eval | getline ap_intval
+ eval | getline ap_isint
+ close(eval)
+
+ pnum = num + 0
+ sign = ""
+ if (pnum < 0) {
+ pnum = -pnum
+ sign = "-"
+ }
+
+ isint = isinteger(pnum, chunk16)
+ if (isint) {
+ if (! ap_isint)
+ print "!!!!!!!!!MPFR says it is not an integer!!!!!!!!!"
+ printf("%s (%s) is an integer\n", $0, sprintf("%0.34e", num))
+ printf("sprintf_int64 :%c%s\n", sign, sprintf_int64(chunk16))
+ printf("arbitrary-prec :%s\n", ap_intval)
+ printf("LDBL %%0.Lf :%c%020d\n", sign, pnum) # or
garbage if no %Lf
+ } else {
+ if (ap_isint)
+ print "!!!!!!!!!MPFR says it is an integer!!!!!!!!!"
+ printf("%s (%s) is not an integer\n", $0, sprintf("%0.34e",
num))
+ }
+
+#
+# Don't have %Lf !
+# printf("floor = %d\n", floor(num))
+# printf("ceil = %d\n", ceil(num))
+#
+ # use MPFR to format floor(), ceil() returned values
+ eval = sprintf("./gawk -M -vPREC=64 'BEGIN { printf(\"%%d\", %s); print
\"\";\
+ printf(\"%%d\", %s); print \"\"; }'", floor(num), ceil(num))
+ eval | getline floor_val
+ eval | getline ceil_val
+ close(eval)
+
+ print "floor =", floor_val
+ print "ceil =", ceil_val
+}
+
+
+function sprintf_int64(chunk, d0, d1, d2, d3, d4, q, k, dec_base)
+{
+ #
+ # Binary to Decimal Conversion in Limited Precision - Douglas B. Jones
+ # URL: homepage.cs.uiowa.edu/~jones/bcd/decimal.html
+ #
+
+ dec_base = 10000
+
+ if (DEBUG) {
+ for (k = 0; k < 4; k++) {
+ if (chunk[k] > 65535)
+ printf("chunk[%d] > 65535, this ain't right\n",
k)
+ }
+ print "(16 MSB bits first):", chunk[3], "|", chunk[2], "|",
chunk[1], "|", chunk[0]
+ }
+
+ d0 = 656 * chunk[3] + 7296 * chunk[2] + 5536 * chunk[1] + chunk[0];
+
+ q = int(d0 / dec_base); # int() may call floorl(), but double
version is good enough
+ # in this case. This is only a consideration in
the awk code,
+ # not in the C translation.
+
+ d0 = d0 - q * dec_base; # not using % operator, may not have
fmodl()
+
+ d1 = q + 7671 * chunk[3] + 9496 * chunk[2] + 6 * chunk[1];
+ q = int(d1 / dec_base);
+ d1 = d1 - q * dec_base;
+
+ d2 = q + 4749 * chunk[3] + 42 * chunk[2];
+ q = int(d2 / dec_base);
+ d2 = d2 - q * dec_base;
+
+ d3 = q + 281 * chunk[3];
+ q = int(d3 / dec_base);
+ d3 = d3 - q * dec_base;
+
+ d4 = q;
+ return sprintf("%4.4u%4.4u%4.4u%4.4u%4.4u", \
+ d4, d3, d2, d1, d0); # don't need %Lf support here
+}
+
+
+function floor(x, d, sign, chunk)
+{
+ x = +x
+
+#
+# if (isnan(x) || isinf(x) || x == 0)
+# return x
+
+ sign = 1
+ if (x < 0) {
+ sign = -1
+ x = -x
+ }
+
+ #
+ # For x >= 2^LDBL_MANT_DIG, use floor for C-double. Bailing out here.
+ #
+
+ if (x >= POW2[LDBL_MANT_DIG]) {
+ printf("number %0.34e is too big (> %0.34e), bailing out\n", x,
2^64 - 1)
+ exit(1);
+ }
+
+ isinteger(x, chunk)
+
+ # d is long double type
+ d = chunk[3] * POW2[48] + chunk[2] * POW2[32] + chunk[1] * POW2[16] +
chunk[0]
+
+ if (d == x)
+ return sign * d
+ return (sign > 0) ? d : -d - 1
+}
+
+
+function ceil(x, d, sign, chunk)
+{
+ x = +x
+
+#
+# if (isnan(x) || isinf(x) || x == 0)
+# return x
+
+ sign = 1
+ if (x < 0) {
+ sign = -1
+ x = -x
+ }
+
+ #
+ # For x >= 2^LDBL_MANT_DIG, use ceil for C-double. Bailing out here.
+ #
+
+ if (x >= POW2[LDBL_MANT_DIG]) {
+ printf("number %0.34e is too big (> %0.34e), bailing out\n", x,
2^64 - 1)
+ exit(1);
+ }
+
+ isinteger(x, chunk)
+ # d is long double type
+ d = chunk[3] * POW2[48] + chunk[2] * POW2[32] + chunk[1] * POW2[16] +
chunk[0]
+
+ if (d == x)
+ return sign * d
+ return (sign > 0) ? d + 1 : -d
+}
+
+
+function isinteger(x, D, j, low, high)
+{
+ # Assume x >= 0
+
+ #
+ # FIXME -- x >= 2^LDBL_MANT_DIG; C-double test is adequate, these
integers have same representations.
+ # output maybe incorrect anyway.
+ #
+
+ if (x >= POW2[LDBL_MANT_DIG]) {
+ printf("number %0.34e is too big (> %0.34e), bailing out\n", x,
2^64 - 1)
+ exit(1);
+ }
+
+ D[0] = D[1] = D[2] = D[3] = 0
+
+ # binary search
+ high = LDBL_MANT_DIG - 1
+ while (x >= 2) {
+ low = 0
+ while (low <= high) {
+ mid = int((low + high) / 2)
+ if (x < POW2[mid])
+ high = mid - 1
+ else
+ low = mid + 1
+ }
+
+ if (x < POW2[low]) {
+ x -= POW2[low - 1]
+
+ #
+ # Easier to seperate the 64 bits into 16-bit chunks now
+ # then bit fiddling later.
+ # | D[3] | D[2] | D[1] | D[0] |
+ # |<------- x (64 bits) ----->|
+ #
+
+ if (low <= 16) { D[0] += POW2[low - 1] }
+ else if (low <= 32) { D[1] += POW2[low - 17] }
+ else if (low <= 48) { D[2] += POW2[low - 33] }
+ else { D[3] += POW2[low - 49] }
+ }
+
+ high = low
+ }
+
+#
+# # XXX -- O(n^2) vs. O(n * log(n)) -- how much do we really gain?
+# # Consider extra additions and divisions, but lower # of <
comparisons.
+#
+# while (x >= 2) {
+# for (j = 1; j <= LDBL_MANT_DIG; j++) {
+# if (x < POW2[j]) {
+# x -= POW2[j - 1]
+# break
+# }
+# }
+# }
+#
+
+ if (x >= 1) {
+ D[0]++
+ x -= 1
+ }
+ return (x == 0)
+}
+
diff --git a/misc/sprintf_0Lf.awk b/misc/sprintf_0Lf.awk
deleted file mode 100644
index d504604..0000000
--- a/misc/sprintf_0Lf.awk
+++ /dev/null
@@ -1,172 +0,0 @@
-#
-# sprintf_0Lf.awk -- format integers stored as long doubles and
-# don't have %Lf in (s)printf.
-# This test script handles integers < 2^64 even with 128-bit long doubles.
-# Requires both -M (for verification) and -B
-# Usage: echo "2^63 + 2^54 + 2^13" | ./gawk -B -f misc/sprintf_0Lf.awk
-#
-
-BEGIN {
- # 113 for 128-bit long double, but can only get 64-bit ints.
- LDBL_MANT_DIG = 64
- POW2[0] = 1
- for (i = 1; i <= LDBL_MANT_DIG; i++)
- POW2[i] = POW2[i - 1] * 2
-}
-
-{
- #
- # Compare apples to apples, setting PREC=64; With real 128-bit long
doubles,
- # one probably should use PREC=quad or 113.
- #
-
- eval = sprintf("./gawk -M -vPREC=64 'BEGIN { printf(\"%%0.34e\", %s);
print \"\";\
- printf(\"%%020d\", %s); print \"\";\
- printf(\"%%d\", int(%s) == %s + 0); print \"\"}'", $0, $0, $0,
$0)
- eval | getline num
- eval | getline ap_intval
- eval | getline ap_isint
- close(eval)
-
- pnum = num + 0
- sign = ""
- if (pnum < 0) {
- pnum = -pnum
- sign = "-"
- }
-
- isint = isinteger(pnum, W)
- if (isint) {
- if (! ap_isint)
- print "!!!!!!!!!MPFR says it is not an integer!!!!!!!!!"
- printf("%s (%s) is an integer\n", $0, sprintf("%0.34e", num))
- printf("sprintf_int64 :%c%s\n", sign, sprintf_int64(W))
- printf("arbitrary-prec :%s\n", ap_intval)
- printf("LDBL %%0.Lf :%c%020d\n", sign, pnum) # or
garbage if no %Lf
- } else {
- if (ap_isint)
- print "!!!!!!!!!MPFR says it is an integer!!!!!!!!!"
- printf("%s (%s) is not an integer\n", $0, sprintf("%0.34e",
num))
- }
-}
-
-
-function sprintf_int64(d, d0, d1, d2, d3, d4, q, k)
-{
- #
- # Binary to Decimal Conversion in Limited Precision - Douglas B. Jones
- # URL: homepage.cs.uiowa.edu/~jones/bcd/decimal.html
- #
-
- if (DEBUG) {
- for (k = 0; k < 4; k++) {
- if (d[k] > 65535)
- printf("d[%d] > 65535, this ain't right\n", k)
- }
- print "(16 MSB bits first):", d[3], "|", d[2], "|", d[1], "|",
d[0]
- }
-
- d0 = 656 * d[3] + 7296 * d[2] + 5536 * d[1] + d[0];
-
- q = int(d0 / 10000); # int() may call floorl(), but double version
is good enough
- # in this case. This is only a consideration in
the awk code,
- # not in the C translation.
-
- d0 = d0 - q * 10000; # not using % operator, may not have fmodl()
-
- d1 = q + 7671 * d[3] + 9496 * d[2] + 6 * d[1];
- q = int(d1 / 10000);
- d1 = d1 - q * 10000;
-
- d2 = q + 4749 * d[3] + 42 * d[2];
- q = int(d2 / 10000);
- d2 = d2 - q * 10000;
-
- d3 = q + 281 * d[3];
- q = int(d3 / 10000);
- d3 = d3 - q * 10000;
-
- d4 = q;
- return sprintf("%4.4u%4.4u%4.4u%4.4u%4.4u", \
- d4, d3, d2, d1, d0); # don't need %Lf support here
-}
-
-
-function isinteger(x, D, j, low, high)
-{
- # Assume x >= 0
-
- #
- # XXX: x >= 2^LDBL_MANT_DIG
- # test using C-double is adequate at least in the range
- # 2^LDBL_MANT_DIG <= x <= DBL_MAX; the integers have same
- # representations. Here in awk code, we simply crap out.
- #
-
- if (x >= POW2[LDBL_MANT_DIG]) {
- printf("number %0.34e is too big (> %0.34e), bailing out\n", x,
2^64 - 1)
- exit(1);
- }
-
- D[0] = D[1] = D[2] = D[3] = 0
-
- # binary search
- high = LDBL_MANT_DIG - 1
- while (x >= 2) {
- low = 0
- while (low <= high) {
- mid = int((low + high) / 2)
- if (x < POW2[mid])
- high = mid - 1
- else
- low = mid + 1
- }
-
- if (x < POW2[low]) {
- x -= POW2[low - 1]
-
- #
- # Easier to seperate the 64 bits into 16-bit chunks now
- # then doing bit fiddling later.
- # | D[3] | D[2] | D[1] | D[0] |
- # |<------- x (64 bits) ----->|
- #
-
- if (low <= 16) { D[0] += POW2[low - 1] }
- else if (low <= 32) { D[1] += POW2[low - 17] }
- else if (low <= 48) { D[2] += POW2[low - 33] }
- else { D[3] += POW2[low - 49] }
- }
-
- high = low
- }
-
-#
-# # XXX -- O(n^2) vs. O(n * log(n)) -- how much do we really gain?
-# # Consider extra additions and divisions, but less # of <
comparisons.
-#
-# while (x >= 2) {
-# for (j = 1; j <= LDBL_MANT_DIG; j++) {
-# if (x < POW2[j]) {
-# x -= POW2[j - 1]
-# break
-# }
-# }
-# }
-#
-
- if (x == 0)
- return 1;
- if (x == 1) {
- if (and(D[0], 1)) {
- ######## XXX: MAYBE ADDING A SECOND TIME ??? #######
- print "Opps! there may be a bug in code! D[0] = ", D[0]
> "/dev/stderr"
- exit(1)
- }
-
- D[0]++
- return 1
- }
-
- return 0
-}
-----------------------------------------------------------------------
Summary of changes:
TODO.LDBL | 4 -
configure | 2 +-
configure.ac | 2 +-
misc/bin2dec.awk | 83 +++++++++++++
misc/sprintf128_0Lf.awk | 308 +++++++++++++++++++++++++++++++++++++++++++++++
misc/sprintf64_0Lf.awk | 247 +++++++++++++++++++++++++++++++++++++
misc/sprintf_0Lf.awk | 172 --------------------------
7 files changed, 640 insertions(+), 178 deletions(-)
create mode 100644 misc/bin2dec.awk
create mode 100644 misc/sprintf128_0Lf.awk
create mode 100644 misc/sprintf64_0Lf.awk
delete mode 100644 misc/sprintf_0Lf.awk
hooks/post-receive
--
gawk
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gawk-diffs] [SCM] gawk branch, long-double, updated. 88304b52a0dc79b4b33b038acb6f20aaacf01092,
John Haque <=