[gawk-diffs] [SCM] gawk branch, num-handler, updated. d898d83434007253f314c8f3fabcb2686820026a

[John Haque]

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, num-handler has been updated
       via  d898d83434007253f314c8f3fabcb2686820026a (commit)
       via  5089b03ccc67b6a5399a852b6d9acee8f20349a2 (commit)
      from  9d7c03c2c55d28d9664881eff0a5e11b030cc67c (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=d898d83434007253f314c8f3fabcb2686820026a

commit d898d83434007253f314c8f3fabcb2686820026a
Author: John Haque <address@hidden>
Date:   Sun Jan 6 09:01:49 2013 -0600

    Add repl_math.awk file.

diff --git a/awklib/ChangeLog b/awklib/ChangeLog
index 98bdc3d..178ce8b 100644
--- a/awklib/ChangeLog
+++ b/awklib/ChangeLog
@@ -1,3 +1,7 @@
+2013-01-06         John Haque            <address@hidden>
+
+       * eg/lib/repl_math.awk: New file.
+
 2012-12-24         Arnold D. Robbins     <address@hidden>
 
        * 4.0.2: Release tar ball made.
diff --git a/awklib/eg/lib/repl_math.awk b/awklib/eg/lib/repl_math.awk
new file mode 100644
index 0000000..a57be27
--- /dev/null
+++ b/awklib/eg/lib/repl_math.awk
@@ -0,0 +1,294 @@
+# repl_math.awk --- arbitrary-precision math functions
+
+#
+#      repl_sin        -- Compute sin(x)
+#      repl_cos        -- Compute cos(x)
+#      repl_atan2      -- Compute atan2(y, x)
+#
+
+#
+# Machin's formula to compute pi 
(http://mathworld.wolfram.com/MachinsFormula.html):
+#      pi / 4 = 4 acot(5) - acot(239)
+#            = 4 atan(1/5) - atan(1/239)
+#
+# 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) = (1 / (1 + x^2))    + (2/3) * (1 / (1 + x^2)^2)
+#                                     + (2*4/(3*5)) * (1 / (1 + x^2)^3)
+#                                     + (2*4*6/(3*5*7)) * (1 / (1 + x^2)^4)  + 
... 
+#
+
+function euler_atan_one_over(x,        \
+               xpow2_plus_one, term, sum, i, err)
+{
+       xpow2_plus_one = x * x + 1
+       term = x / xpow2_plus_one
+       sum = term
+       i = 0
+       err = 1.0
+
+       while (err > __EPSILON__) {
+               term *= (i + 2) / (i + 3)
+               err = term /= xpow2_plus_one
+               i += 2
+               sum += term
+               if (err < 0)
+                       err = -err
+       }
+       return sum
+}
+
+function setup_repl_math( \
+               prec, digits, extra_prec)
+{
+       switch (PREC) {
+       case "half":    prec = 11; break;
+       case "single":  prec = 24; break;
+       case "double":  prec = 53; break;
+       case "quad":    prec = 113; break;
+       case "oct":     prec = 237; break;
+       default:        prec = PREC + 0;
+       }
+
+       if (prec <= 0) {
+               # double or long double ?
+               print "PREC value not specified" > "/dev/stderr"
+               exit(1)
+       }
+       __SAVE_PREC__ = PREC
+       extra_prec = 10
+       prec += extra_prec      # temporarily raise precision
+       if (prec != __PI_PREC__) {
+               # compute PI only once for a given precision
+               digits = int ((prec - extra_prec) / 3.32193)
+               __EPSILON__ = sprintf("1.0e-%d", digits + 2) + 0
+               __PI_OVER_4__ = 4 * euler_atan_one_over(5) - 
euler_atan_one_over(239)
+               __PI_PREC__ = prec
+       }
+}
+
+#
+# atan2(y, x) =        atan(y/x),              x > 0
+#             =        atan(y/x) + pi,         x < 0, y >= 0
+#             = atan(y/x) - pi,                x < 0, y < 0
+#             = pi/2,                  x = 0, y > 0
+#             = -pi/2,                 x = 0, y < 0
+#             = ?                      x = 0, y = 0
+#
+
+function euler_atan2(y, x,     \
+               sign, plus_inf, minus_inf)
+{
+       # Using Euler atan(1/x) and the identity:
+       #       atan(x) = - pi / 2 - atan(1/x),         x < 0
+       #               = pi / 2 - atan(1/x),           x > 0
+
+       plus_inf = "+inf" + 0
+       minus_inf = "-inf" + 0
+
+       # detect all the "abnormal" cases first
+       x = + x                 # or x += 0.0 or x = x + 0.0
+       y = + y
+       if (x == "+nan" + 0 || x == "-nan" + 0 ||
+                       y == "+nan" + 0 || y == "-nan" + 0)
+               return "nan" + 0
+
+       if (y == plus_inf) {
+               if (x == minus_inf)
+                       return 3 * __PI_OVER_4__
+               else if (x == plus_inf)
+                       return __PI_OVER_4__
+               else
+                       return 2 * __PI_OVER_4__
+       } else if (y == minus_inf) {
+               if (x == minus_inf)
+                       return - 3 * __PI_OVER_4__
+               else if (x == plus_inf)
+                       return - __PI_OVER_4__
+               else
+                       return - 2 * __PI_OVER_4__
+       }
+
+       if (x == plus_inf)
+               return atan2(y, x)      # use builtin, -0 or -0
+       if (x == minus_inf) {
+               if (y >= 0)
+                       return 4 * __PI_OVER_4__
+               # y < 0
+               return - 4 * __PI_OVER_4__
+       }
+
+       if (x > 0) {
+               if (y == 0)
+                       return atan2(y, x);     # use builtin; returns +0 or -0
+               sign = 1
+               if (y < 0) {
+                       sign = -1
+                       y = -y
+               }
+               if (y > x)
+                       return sign * (2 * __PI_OVER_4__ - 
euler_atan_one_over(y / x))
+               return sign * euler_atan_one_over(x / y)
+       } else if (x < 0) {
+               if (y == 0) {
+                       if (atan2(y, x) < 0)    # use builtin to detect sign
+                               return - 4 * __PI_OVER_4__
+                       return 4 * __PI_OVER_4__
+               }
+
+               if (y < 0) {
+                       y = -y
+                       x = -x
+                       if (y > x)
+                               return - 2 * __PI_OVER_4__ - 
euler_atan_one_over(y / x)
+                       return euler_atan_one_over(x / y) - 4 * __PI_OVER_4__
+               }
+               # y > 0
+               x = -x
+               if (y > x)
+                       return 2 * __PI_OVER_4__ + euler_atan_one_over(y / x)
+               return - euler_atan_one_over(x / y) + 4 * __PI_OVER_4__
+       }
+
+       # x == +0/-0 and y == +0/-0
+       return atan2(y, x);     # use builtin
+}
+
+
+function taylor_sin(x, \
+               i, fact, xpow2, xpow_odd, sum, term, sign, err)
+{
+       i = 1
+       fact = 1
+       xpow2 = x * x
+       xpow_odd = x
+       sum = x
+       sign = 1
+       err = 1.0
+
+       while (err > __EPSILON__) {
+               fact *= (i + 1) * (i + 2)
+               i += 2
+               xpow_odd *= xpow2
+               sign *= -1
+               err = term = xpow_odd / fact  * sign
+               sum += term
+               if (err < 0)
+                       err = -err
+       }
+       return sum
+}
+
+function taylor_cos(x, \
+               i, fact, xpow2, xpow_even, sum, term, sign, err)
+{
+       i = 0
+       fact = 1
+       xpow2 = x * x
+       xpow_even = 1
+       sum = 1
+       sign = 1
+       err = 1.0
+
+       while (err > __EPSILON__) {
+               fact *= (i + 1) * (i + 2)
+               i += 2
+               xpow_even *= xpow2
+               sign *= -1
+               err = term = xpow_even / fact * sign
+               sum += term
+               if (err < 0)
+                       err = -err
+       }
+       return sum
+}
+
+#
+# For 0 <= x <= PI/4, using Taylor series approximation for sin(x):
+#      x - x^3/5! + x^5/7! - ...
+#
+# for PI/4 < x <= PI/2, use identity sin(x) = cos(PI/2 - x).
+#
+#
+
+function repl_sin(x,   \
+               k, sign, y, sv)
+{
+       setup_repl_math()
+       x = + x         # or x += 0.0 or x = x + 0.0 
+       if (x == "+inf" + 0 || x == "-inf" + 0 ||
+                       x == "+nan" + 0 || x == "-nan" + 0)
+               return "nan" + 0
+
+       if (x < 0) {
+               # sin(x) = - sin(x)
+               sign = -1
+               x = -x
+       } else
+               sign = 1
+
+       # range reduction -- 0 <= y <= pi / 4
+
+       k = int(x / __PI_OVER_4__)
+       sign *= ( int(k / 4) % 2 ? -1 : 1 )
+       switch (k % 4) {
+       case 0: y = x - k * __PI_OVER_4__; sv = taylor_sin(y); break;
+       case 1: y = (k + 1) * __PI_OVER_4__ - x; sv = taylor_cos(y); break;
+       case 2: y = x - k * __PI_OVER_4__; sv = taylor_cos(y); break;
+       case 3: y = (k + 1) * __PI_OVER_4__ - x; sv = taylor_sin(y); break;
+       }
+       sv *= sign
+
+       PREC = __SAVE_PREC__
+       return +sv;     # unary plus returns a number with current precision 
+}
+
+#
+# Using Taylor series approximation for sin(x) for 0 <= x <= PI/4:
+#      1 - x^2/2! + x^4/4! - ...
+# for PI/4 < x < PI/2, use identity cos(x) = sin(PI/2 - x).
+#
+
+
+function repl_cos(x,   \
+               k, sign, y, cv)
+{
+       setup_repl_math()
+
+       x = + x         # or x += 0.0 or x = x + 0.0
+       if (x == "+inf" + 0 || x == "-inf" + 0 ||
+                       x == "+nan" + 0 || x == "-nan" + 0)
+               return "nan" + 0
+
+       if (x < 0)      # cos(x) = cos(-x)
+               x = -x
+
+       # range reduction -- 0 <= y <= pi / 4
+
+       k = int(x / __PI_OVER_4__)
+       sign = ( int(k / 4) % 2 ? -1 : 1 )
+       switch (k % 4) {
+       case 0: y = x - k * __PI_OVER_4__; cv = taylor_cos(y); break;
+       case 1: y = (k + 1) * __PI_OVER_4__ - x; cv = taylor_sin(y); break;
+       case 2: y = x - k * __PI_OVER_4__; cv = -taylor_sin(y); break;
+       case 3: y = (k + 1) * __PI_OVER_4__ - x; cv = -taylor_cos(y); break;
+       }
+       cv *= sign
+
+       PREC = __SAVE_PREC__
+       return +cv      # unary plus to apply current precision
+}
+
+function repl_atan2(y, x, \
+               tv)
+{
+       setup_repl_math()
+       tv = euler_atan2(y, x)
+
+       PREC = __SAVE_PREC__
+       return +tv      # unary plus to apply current precision
+}

http://git.sv.gnu.org/cgit/gawk.git/commit/?id=5089b03ccc67b6a5399a852b6d9acee8f20349a2

commit 5089b03ccc67b6a5399a852b6d9acee8f20349a2
Author: John Haque <address@hidden>
Date:   Sun Jan 6 08:51:19 2013 -0600

    Fix bug in GMP to MPFR number conversion.

diff --git a/mpfr.c b/mpfr.c
index e3215c9..edb8fd7 100644
--- a/mpfr.c
+++ b/mpfr.c
@@ -1737,7 +1737,7 @@ mpfp_sprintf(const char *mesg, ...)
 static mpfr_ptr
 mpz2mpfr(mpz_ptr mpz_val, mpfr_ptr mpfr_val)
 {
-       size_t prec;
+       long prec, prec1;
        int tval;
 
        /*
@@ -1755,10 +1755,10 @@ mpz2mpfr(mpz_ptr mpz_val, mpfr_ptr mpfr_val)
        /* estimate minimum precision for exact conversion */
        prec = mpz_sizeinbase(mpz_val, 2);      /* most significant 1 bit 
position starting at 1 */
 
-       if (mpfr_val != NULL && mpfr_get_prec(mpfr_val) >= prec)
+       if (mpfr_val != NULL && (prec1 = mpfr_get_prec(mpfr_val)) >= prec)
                goto finish;
 
-       prec -= (size_t) mpz_scan1(mpz_val, 0); /* least significant 1 bit 
index starting at 0 */
+       prec -= (long) mpz_scan1(mpz_val, 0);   /* least significant 1 bit 
index starting at 0 */
        if (prec < MPFR_PREC_MIN)
                prec = MPFR_PREC_MIN;
        else if (prec > MPFR_PREC_MAX)
@@ -1767,10 +1767,11 @@ mpz2mpfr(mpz_ptr mpz_val, mpfr_ptr mpfr_val)
        if (mpfr_val == NULL) {
                emalloc(mpfr_val, mpfr_ptr, sizeof (mpfr_t), "mpz2mpfr");
                mpfr_init2(mpfr_val, prec);
-       } else if (prec < mpfr_get_prec(mpfr_val))
+       } else if (prec > prec1)
                mpfr_set_prec(mpfr_val, prec);
 
 finish:
+
        tval = mpfr_set_z(mpfr_val, mpz_val, ROUND_MODE);
        IEEE_FMT(mpfr_val, tval);
        return mpfr_val;

-----------------------------------------------------------------------

Summary of changes:
 awklib/ChangeLog            |    4 +
 awklib/eg/lib/repl_math.awk |  294 +++++++++++++++++++++++++++++++++++++++++++
 mpfr.c                      |    9 +-
 3 files changed, 303 insertions(+), 4 deletions(-)
 create mode 100644 awklib/eg/lib/repl_math.awk


hooks/post-receive
-- 
gawk