[Gmp-commit] /var/hg/gmp: 5 new changesets

mercurial at gmplib.org mercurial at gmplib.org
Mon Dec 31 01:02:29 CET 2012


details:   /var/hg/gmp/rev/561f67240dd8
changeset: 15222:561f67240dd8
user:      Torbjorn Granlund <tege at gmplib.org>
date:      Mon Dec 31 00:03:33 2012 +0100
description:
Whitespace cleanup.

details:   /var/hg/gmp/rev/bf2ba4e63461
changeset: 15223:bf2ba4e63461
user:      Torbjorn Granlund <tege at gmplib.org>
date:      Mon Dec 31 00:55:47 2012 +0100
description:
(check_random): Limit exponents on vax.

details:   /var/hg/gmp/rev/d07282c2fa0f
changeset: 15224:d07282c2fa0f
user:      Torbjorn Granlund <tege at gmplib.org>
date:      Mon Dec 31 01:00:56 2012 +0100
description:
(double_extract): New union type for vax D floats.

details:   /var/hg/gmp/rev/808dff309bfd
changeset: 15225:808dff309bfd
user:      Torbjorn Granlund <tege at gmplib.org>
date:      Mon Dec 31 01:02:09 2012 +0100
description:
Minor reorg, add vax D code.

details:   /var/hg/gmp/rev/eb121660202e
changeset: 15226:eb121660202e
user:      Torbjorn Granlund <tege at gmplib.org>
date:      Mon Dec 31 01:02:22 2012 +0100
description:
ChangeLog

diffstat:

 AUTHORS             |   26 +++---
 ChangeLog           |    8 ++
 gmp-impl.h          |   15 ++++
 mpn/generic/get_d.c |  181 +++++++++++++++++++++++++++++++++------------------
 tests/mpq/t-get_d.c |   13 ++-
 5 files changed, 161 insertions(+), 82 deletions(-)

diffs (truncated from 369 to 300 lines):

diff -r 6c283c477e58 -r eb121660202e AUTHORS
--- a/AUTHORS	Sun Dec 30 09:06:09 2012 +0100
+++ b/AUTHORS	Mon Dec 31 01:02:22 2012 +0100
@@ -4,20 +4,21 @@
 
 John Amanatides		Original version of mpz/pprime_p.c
 
-Paul Zimmermann		mpn/generic/mul_fft.c, dc_divrem_n.c, rootrem.c,
-			old mpz/powm.c, old toom3 code.
+Paul Zimmermann		mpn/generic/mul_fft.c, now defunct dc_divrem_n.c,
+			rootrem.c, old mpz/powm.c, old toom3 code.
 
-Ken Weber		mpn/generic/bdivmod.c, old mpn/generic/gcd.c
+Ken Weber		Now defunct mpn/generic/bdivmod.c, old mpn/generic/gcd.c
 
-Bennet Yee		mpz/jacobi.c mpz/legendre.c
+Bennet Yee		Previous versions of mpz/jacobi.c mpz/legendre.c
 
 Andreas Schwab		mpn/m68k/lshift.asm, mpn/m68k/rshift.asm
 
-Robert Harley		Old mpn/generic/mul_n.c, many files in mpn/arm
+Robert Harley		Old mpn/generic/mul_n.c, previous versions of files in
+			mpn/arm
 
 Linus Nordberg		Random number framework, original autoconfery
 
-Kent Boortz		MacOS 9 port
+Kent Boortz		MacOS 9 port, now defunct.
 
 Kevin Ryde		Most x86 assembly, new autoconfery, and countless other
 			things (please see the GMP manual for complete list)
@@ -27,12 +28,12 @@
 Pedro Gimeno		Mersenne Twister random generator, other random number
 			revisions
 
-Jason Moxham		mpz/fac_ui.c and gen-fac_ui.c
+Jason Moxham		Previous versions of mpz/fac_ui.c and gen-fac_ui.c
 
 Niels Möller		gen-jacobitab.c,
-      			mpn/generic/hgcd2.c, hgcd.c, hgcd_step.c,
+			mpn/generic/hgcd2.c, hgcd.c, hgcd_step.c,
 			hgcd_appr.c, hgcd_matrix.c, hgcd_reduce.c,
- 			gcd.c, gcdext.c, matrix22_mul.c,
+			gcd.c, gcdext.c, matrix22_mul.c,
 			gcdext_1.c, gcd_subdiv_step.c, gcd_lehmer.c,
 			gcdext_subdiv_step.c, gcdext_lehmer.c,
 			jacobi_2.c, jacbase.c, hgcd_jacobi.c, hgcd2_jacobi.c
@@ -43,15 +44,14 @@
 			toom_eval_pm1.c, toom_eval_pm2.c, toom_eval_pm2exp.c,
 			divexact.c, mod_1_1.c, div_qr_2.c,
 			div_qr_2n_pi1.c, div_qr_2u_pi1.c, broot.c,
-      			brootinv.c,
+			brootinv.c,
 			mpn/x86/k7/invert_limb.asm, mod_1_1.asm,
 			mpn/x86_64/invert_limb.asm,
-      			invert_limb_table.asm, mod_1_1.asm,
+			invert_limb_table.asm, mod_1_1.asm,
 			div_qr_2n_pi1.asm, div_qr_2u_pi1.asm,
 			mpn/x86_64/core2/aorsmul_1.asm,
 			mpz/nextprime.c, divexact.c, gcd.c, gcdext.c,
-      			jacobi.c, combit.c, mini-gmp/mini-gmp.c.
-				
+			jacobi.c, combit.c, mini-gmp/mini-gmp.c.
 
 Marco Bodrato		mpn/generic/toom44_mul.c, toom4_sqr.c, toom53_mul.c,
 			toom62_mul.c, toom43_mul.c, toom52_mul.c, toom54_mul.c,
diff -r 6c283c477e58 -r eb121660202e ChangeLog
--- a/ChangeLog	Sun Dec 30 09:06:09 2012 +0100
+++ b/ChangeLog	Mon Dec 31 01:02:22 2012 +0100
@@ -1,3 +1,11 @@
+2012-12-31  Torbjorn Granlund  <tege at gmplib.org>
+
+	* mpn/generic/get_d.c: Minor reorg, add vax D code.
+
+	* gmp-impl.h (double_extract): New union type for vax D floats.
+
+	* tests/mpq/t-get_d.c (check_random): Limit exponents on vax.
+
 2012-12-30 Marco Bodrato <bodrato at mail.dm.unipi.it>
 
 	* tests/mpz/bit.c (check_clr_extend): Check _set shrink.
diff -r 6c283c477e58 -r eb121660202e gmp-impl.h
--- a/gmp-impl.h	Sun Dec 30 09:06:09 2012 +0100
+++ b/gmp-impl.h	Mon Dec 31 01:02:22 2012 +0100
@@ -3722,6 +3722,21 @@
 };
 #endif
 
+#if HAVE_DOUBLE_VAX_D
+union double_extract
+{
+  struct
+    {
+      gmp_uint_least32_t man3:7;	/* highest 7 bits */
+      gmp_uint_least32_t exp:8;		/* excess-128 exponent */
+      gmp_uint_least32_t sig:1;
+      gmp_uint_least32_t man2:16;
+      gmp_uint_least32_t man1:16;
+      gmp_uint_least32_t man0:16;	/* lowest 16 bits */
+    } s;
+  double d;
+};
+#endif
 
 /* Use (4.0 * ...) instead of (2.0 * ...) to work around buggy compilers
    that don't convert ulong->double correctly (eg. SunOS 4 native cc).  */
diff -r 6c283c477e58 -r eb121660202e mpn/generic/get_d.c
--- a/mpn/generic/get_d.c	Sun Dec 30 09:06:09 2012 +0100
+++ b/mpn/generic/get_d.c	Mon Dec 31 01:02:22 2012 +0100
@@ -60,16 +60,12 @@
 /* Return the value {ptr,size}*2^exp, and negative if sign<0.  Must have
    size>=1, and a non-zero high limb ptr[size-1].
 
-   {ptr,size} is truncated towards zero.  This is consistent with other gmp
-   conversions, like mpz_set_f or mpz_set_q, and is easy to implement and test.
+   When we know the fp format, the result is truncated towards zero.  This is
+   consistent with other gmp conversions, like mpz_set_f or mpz_set_q, and is
+   easy to implement and test.
 
-   In the past conversions had attempted (imperfectly) to let the hardware
-   float rounding mode take effect, but that gets tricky since multiple
-   roundings need to be avoided, or taken into account, and denorms mean the
-   effective precision of the mantissa is not constant.  (For reference,
-   mpz_get_d on IEEE systems was ok, except it operated on the absolute value.
-   mpf_get_d and mpq_get_d suffered from multiple roundings and from not always
-   using enough bits to get the rounding right.)
+   When we do not know the format, such truncation seems much harder.  One
+   would need to defeat any rounding mode, including round-up.
 
    It's felt that GMP is not primarily concerned with hardware floats, and
    really isn't enhanced by getting involved with hardware rounding modes
@@ -81,7 +77,7 @@
    limb and is done with shifts and masks.  The 64-bit case in particular
    should come out nice and compact.
 
-   The generic code used to work one bit at a time, which was ont only slow,
+   The generic code used to work one bit at a time, which was not only slow,
    but implicitly relied upon denoms for intermediates, since the lowest bits'
    weight of a perfectly valid fp number underflows in non-denorm.  Therefore,
    the generic code now works limb-per-limb, initially creating a number x such
@@ -120,6 +116,8 @@
 
 
 
+#undef FORMAT_RECOGNIZED
+
 double
 mpn_get_d (mp_srcptr up, mp_size_t size, mp_size_t sign, long exp)
 {
@@ -151,61 +149,7 @@
       exp += GMP_NUMB_BITS * size;
     }
 
-#if ! _GMP_IEEE_FLOATS
-    {      /* Non-IEEE or strange limb size, do something generic. */
-      mp_size_t i;
-      double d, weight;
-      unsigned long uexp;
-
-      /* First generate an fp number disregarding exp, instead keeping things
-	 within the numb base factor from 1, which should prevent overflow and
-	 underflow even for the most exponent limited fp formats.  The
-	 termination criteria should be refined, since we now include too many
-	 limbs.  */
-      weight = 1/MP_BASE_AS_DOUBLE;
-      d = up[size - 1];
-      for (i = size - 2; i >= 0; i--)
-	{
-	  d += up[i] * weight;
-	  weight /= MP_BASE_AS_DOUBLE;
-	  if (weight == 0)
-	    break;
-	}
-
-      /* Now apply exp.  */
-      exp -= GMP_NUMB_BITS;
-      if (exp > 0)
-	{
-	  weight = 2.0;
-	  uexp = exp;
-	}
-      else
-	{
-	  weight = 0.5;
-	  uexp = 1 - (unsigned long) (exp + 1);
-	}
-#if 1
-      /* Square-and-multiply exponentiation.  */
-      if (uexp & 1)
-	d *= weight;
-      while (uexp >>= 1)
-	{
-	  weight *= weight;
-	  if (uexp & 1)
-	    d *= weight;
-	}
-#else
-      /* Plain exponentiation.  */
-      while (uexp > 0)
-	{
-	  d *= weight;
-	  uexp--;
-	}
-#endif
-
-      return sign >= 0 ? d : -d;
-    }
-#else /* _GMP_IEEE_FLOATS */
+#if _GMP_IEEE_FLOATS
     {
       union ieee_double_extract u;
 
@@ -346,5 +290,112 @@
       u.s.sig = (sign < 0);
       return u.d;
     }
+#define FORMAT_RECOGNIZED 1
+#endif
+
+#if HAVE_DOUBLE_VAX_D
+    {
+      union double_extract u;
+
+      up += size;
+
+      mhi = up[-1];
+
+      count_leading_zeros (lshift, mhi);
+      exp -= lshift;
+      mhi <<= lshift;
+
+      mlo = 0;
+      if (size > 1)
+	{
+	  mlo = up[-2];
+	  if (lshift != 0)
+	    mhi += mlo >> (GMP_LIMB_BITS - lshift);
+	  mlo <<= lshift;
+
+	  if (size > 2 && lshift > 8)
+	    {
+	      x = up[-3];
+	      mlo += x >> (GMP_LIMB_BITS - lshift);
+	    }
+	}
+
+      if (UNLIKELY (exp >= 128))
+	{
+	  /* overflow, return maximum number */
+	  mhi = 0xffffffff;
+	  mlo = 0xffffffff;
+	  exp = 127;
+	}
+      else if (UNLIKELY (exp < -128))
+	{
+	  return 0.0;	 /* underflows to zero */
+	}
+
+      u.s.man3 = mhi >> 24;	/* drop msb, since implicit */
+      u.s.man2 = mhi >> 8;
+      u.s.man1 = (mhi << 8) + (mlo >> 24);
+      u.s.man0 = mlo >> 8;
+      u.s.exp = exp + 128;
+      u.s.sig = sign < 0;
+      return u.d;
+    }
+#define FORMAT_RECOGNIZED 1
+#endif
+
+#if ! FORMAT_RECOGNIZED
+    {      /* Non-IEEE or strange limb size, do something generic. */
+      mp_size_t i;
+      double d, weight;
+      unsigned long uexp;
+
+      /* First generate an fp number disregarding exp, instead keeping things
+	 within the numb base factor from 1, which should prevent overflow and
+	 underflow even for the most exponent limited fp formats.  The
+	 termination criteria should be refined, since we now include too many
+	 limbs.  */
+      weight = 1/MP_BASE_AS_DOUBLE;
+      d = up[size - 1];
+      for (i = size - 2; i >= 0; i--)
+	{
+	  d += up[i] * weight;
+	  weight /= MP_BASE_AS_DOUBLE;
+	  if (weight == 0)
+	    break;
+	}
+
+      /* Now apply exp.  */
+      exp -= GMP_NUMB_BITS;
+      if (exp > 0)
+	{
+	  weight = 2.0;
+	  uexp = exp;
+	}
+      else
+	{
+	  weight = 0.5;


More information about the gmp-commit mailing list