[Gmp-commit] /home/hgfiles/gmp: Provide mpn_toom63_mul with unit testing prog...
mercurial at gmplib.org
mercurial at gmplib.org
Fri Dec 18 03:50:53 CET 2009
details: /home/hgfiles/gmp/rev/effd5e60bbfa
changeset: 13118:effd5e60bbfa
user: Torbjorn Granlund <tege at gmplib.org>
date: Fri Dec 18 03:50:43 2009 +0100
description:
Provide mpn_toom63_mul with unit testing program.
diffstat:
ChangeLog | 12 +
configure.in | 3 +-
gmp-impl.h | 15 +-
mpn/generic/toom63_mul.c | 378 ++++++++++++++++++++++++++++++++++++
mpn/generic/toom_interpolate_8pts.c | 170 ++++++++++++++++
tests/mpn/Makefile.am | 2 +-
tests/mpn/t-toom63.c | 8 +
7 files changed, 585 insertions(+), 3 deletions(-)
diffs (truncated from 663 to 300 lines):
diff -r c58defebaeee -r effd5e60bbfa ChangeLog
--- a/ChangeLog Thu Dec 17 17:35:44 2009 +0100
+++ b/ChangeLog Fri Dec 18 03:50:43 2009 +0100
@@ -1,3 +1,15 @@
+2009-12-18 Torbjorn Granlund <tege at gmplib.org>
+
+ * tests/mpn/t-toom63.c: New test program.
+
+2009-12-18 Marco Bodrato <bodrato at mail.dm.unipi.it>
+
+ * mpn/generic/toom63_mul.c: New file.
+ * mpn/generic/toom_interpolate_8pts.c: New file.
+ * gmp-impl.h: Provide corresponding declarations.
+ * configure.in (gmp_mpn_functions): List toom_interpolate_8pts and
+ toom63_mul.
+
2009-12-17 Torbjorn Granlund <tege at gmplib.org>
* mpn/generic/mul.c: Move allocation of ws to where it is used.
diff -r c58defebaeee -r effd5e60bbfa configure.in
--- a/configure.in Thu Dec 17 17:35:44 2009 +0100
+++ b/configure.in Fri Dec 18 03:50:43 2009 +0100
@@ -2490,12 +2490,13 @@
tdiv_qr jacbase get_d \
matrix22_mul hgcd2 hgcd mullo_n mullo_basecase \
toom22_mul toom32_mul toom42_mul toom52_mul toom62_mul \
- toom33_mul toom43_mul toom53_mul \
+ toom33_mul toom43_mul toom53_mul toom63_mul \
toom44_mul \
toom2_sqr toom3_sqr toom4_sqr \
toom_eval_dgr3_pm1 toom_eval_dgr3_pm2 \
toom_eval_pm1 toom_eval_pm2 toom_eval_pm2exp \
toom_interpolate_5pts toom_interpolate_6pts toom_interpolate_7pts \
+ toom_interpolate_8pts \
invertappr invert binvert mulmod_bnm1 sqrmod_bnm1 \
sbpi1_div_q sbpi1_div_qr sbpi1_divappr_q \
dcpi1_div_q dcpi1_div_qr dcpi1_divappr_q \
diff -r c58defebaeee -r effd5e60bbfa gmp-impl.h
--- a/gmp-impl.h Thu Dec 17 17:35:44 2009 +0100
+++ b/gmp-impl.h Fri Dec 18 03:50:43 2009 +0100
@@ -1061,6 +1061,9 @@
#define mpn_toom_interpolate_7pts __MPN(toom_interpolate_7pts)
__GMP_DECLSPEC void mpn_toom_interpolate_7pts __GMP_PROTO ((mp_ptr, mp_size_t, enum toom7_flags, mp_ptr, mp_ptr, mp_ptr, mp_ptr, mp_size_t, mp_ptr));
+#define mpn_toom_interpolate_8pts __MPN(toom_interpolate_8pts)
+__GMP_DECLSPEC void mpn_toom_interpolate_8pts __GMP_PROTO ((mp_ptr, mp_size_t, mp_ptr, mp_ptr, mp_size_t, mp_ptr));
+
#define mpn_toom_eval_dgr3_pm1 __MPN(toom_eval_dgr3_pm1)
__GMP_DECLSPEC int mpn_toom_eval_dgr3_pm1 __GMP_PROTO ((mp_ptr, mp_ptr, mp_srcptr, mp_size_t, mp_size_t, mp_ptr));
@@ -1103,6 +1106,9 @@
#define mpn_toom53_mul __MPN(toom53_mul)
__GMP_DECLSPEC void mpn_toom53_mul __GMP_PROTO ((mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t, mp_ptr));
+#define mpn_toom63_mul __MPN(toom63_mul)
+__GMP_DECLSPEC void mpn_toom63_mul __GMP_PROTO ((mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t, mp_ptr));
+
#define mpn_toom3_sqr __MPN(toom3_sqr)
__GMP_DECLSPEC void mpn_toom3_sqr __GMP_PROTO ((mp_ptr, mp_srcptr, mp_size_t, mp_ptr));
@@ -1739,7 +1745,7 @@
#endif
#ifndef SQRMOD_BNM1_THRESHOLD
-#define SQRMOD_BNM1_THRESHOLD 16
+#define SQRMOD_BNM1_THRESHOLD 16
#endif
#if HAVE_NATIVE_mpn_addmul_2 || HAVE_NATIVE_mpn_redc_2
@@ -4427,6 +4433,13 @@
return 10 * n + 10;
}
+static inline mp_size_t
+mpn_toom63_mul_itch (mp_size_t an, mp_size_t bn)
+{
+ mp_size_t n = 1 + (3 * an >= 6 * bn ? (an - 1) / (size_t) 6 : (bn - 1) / (size_t) 3);
+ return 9 * n + 3;
+}
+
#ifdef __cplusplus
/* A little helper for a null-terminated __gmp_allocate_func string.
diff -r c58defebaeee -r effd5e60bbfa mpn/generic/toom63_mul.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/mpn/generic/toom63_mul.c Fri Dec 18 03:50:43 2009 +0100
@@ -0,0 +1,378 @@
+/* Implementation of the algorithm for Toom-Cook 4.5-way.
+
+ Contributed to the GNU project by Marco Bodrato.
+
+ THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE. IT IS ONLY
+ SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
+ GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
+
+Copyright 2009 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+#if HAVE_NATIVE_mpn_addlsh_n
+#define DO_mpn_addlsh_n(dst,src,n,s,ws) mpn_addlsh_n (dst,dst,src,n,s)
+#else
+static mp_limb_t
+DO_mpn_addlsh_n (mp_ptr dst, mp_srcptr src, mp_size_t n, unsigned int s, mp_ptr ws)
+{
+#if USE_MUL_1
+ return mpn_addmul_1(dst,src,n,CNST_LIMB(1) <<(s));
+#else
+ mp_limb_t __cy = mpn_lshift (ws,src,n,s);
+ return __cy + mpn_add_n (dst,dst,ws,n);
+#endif
+}
+#endif
+
+
+/* Stores |{ap,n}-{bp,n}| in {rp,n}, returns the sign. */
+static int
+abs_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
+{
+ mp_limb_t x, y;
+ while (--n >= 0)
+ {
+ x = ap[n];
+ y = bp[n];
+ if (x != y)
+ {
+ n++;
+ if (x > y)
+ {
+ mpn_sub_n (rp, ap, bp, n);
+ return 0;
+ }
+ else
+ {
+ mpn_sub_n (rp, bp, ap, n);
+ return 1;
+ }
+ }
+ rp[n] = 0;
+ }
+ return 0;
+}
+
+static void
+toom_couple_handling (mp_ptr pp, mp_size_t n, mp_ptr np, int nsign, mp_size_t off, int ps, int ns)
+{
+ if (nsign) {
+#ifdef HAVE_NATIVE_mpn_rsh1sub_n
+ mpn_rsh1sub_n (np, pp, np, n);
+#else
+ mpn_sub_n (np, pp, np, n);
+ mpn_rshift (np, np, n, 1);
+#endif
+ } else {
+#ifdef HAVE_NATIVE_mpn_rsh1add_n
+ mpn_rsh1add_n (np, pp, np, n);
+#else
+ mpn_add_n (np, pp, np, n);
+ mpn_rshift (np, np, n, 1);
+#endif
+ }
+
+#ifdef HAVE_NATIVE_mpn_rsh1sub_n
+ if (ps == 1)
+ mpn_rsh1sub_n (pp, pp, np, n);
+ else
+#endif
+ {
+ mpn_sub_n (pp, pp, np, n);
+ if (ps > 0)
+ mpn_rshift (pp, pp, n, ps);
+ }
+ if (ns > 0)
+ mpn_rshift (np, np, n, ns);
+ pp[n-1] += mpn_add_n (pp+off, pp+off, np, n>>1);
+ ASSERT_NOCARRY( mpn_add_1(pp+n-1, np+off, n-off, pp[n-1]) );
+}
+
+static int
+abs_sub_add_n (mp_ptr rm, mp_ptr rp, mp_srcptr rs, mp_size_t n) {
+ int result = abs_sub_n (rm, rp, rs, n);
+ ASSERT_NOCARRY(mpn_add_n (rp, rp, rs, n));
+ return result;
+}
+
+static int
+mpn_toom_ev_pm1(mp_ptr rp, mp_ptr rm,
+ mp_srcptr ap, unsigned int q, mp_size_t n, mp_size_t t,
+ mp_ptr ws)
+{
+ /* {ap,q*n+t} -> {rp,n+1} {rm,n+1} , with {ws, n+1}*/
+ ASSERT( n >= t );
+ ASSERT( q > 2 );
+ if ( (q & 1) == 0) {
+ rp[n] = mpn_add (rp, ap+n*(q-2), n, ap+n*q, t);
+ q--;
+ ws[n] = mpn_add_n (ws, ap+n*(q-2), ap+n*q, n);
+ q-=3;
+ rp[n]+= mpn_add_n (rp, rp, ap+n*q, n);
+ } else {
+ ws[n] = mpn_add (ws, ap+n*(q-2), n, ap+n*q, t);
+ q-=3;
+ rp[n] = mpn_add_n (rp, ap+n*q, ap+n*(q+2), n);
+ }
+ while (q) {
+ q--;
+ ws[n] += mpn_add_n (ws, ws, ap+n*q, n);
+ q--;
+ rp[n] += mpn_add_n (rp, rp, ap+n*q, n);
+ }
+ return abs_sub_add_n (rm, rp, ws, n + 1);
+}
+
+static int
+mpn_toom_ev_lsh (mp_ptr rp, mp_ptr rm,
+ mp_srcptr ap, unsigned int q, mp_size_t n, mp_size_t t,
+ unsigned int s, mp_ptr ws, mp_ptr scratch)
+{
+ /* {ap,q*n+t} -> {rp,n+1} {rm,n+1} , with {ws, n+1}*/
+ ASSERT( n >= t );
+ ASSERT( s != 0 ); /* or _ev_pm1 should be used */
+ ASSERT( q > 2 );
+ ASSERT( s*q < GMP_NUMB_BITS );
+ if ((q & 1) == 0) {
+ rp[t] = mpn_lshift (rp, ap+n*q, t, s*q);
+ MPN_ZERO (rp + t + 1, n - t);
+ q--;
+ ws[n] = mpn_lshift (ws, ap+n*q, n, s*q);
+ q--;
+ rp[n] += DO_mpn_addlsh_n (rp, ap+n*q, n, s*q, scratch);
+ } else {
+ ws[t] = mpn_lshift (ws, ap+n*q, t, s*q);
+ MPN_ZERO(ws + t + 1, n - t);
+ q--;
+ rp[n] = mpn_lshift (rp, ap+n*q, n, s*q);
+ }
+ do {
+ q--;
+ ws[n] += DO_mpn_addlsh_n (ws, ap+n*q, n, s*q, scratch);
+ q--;
+ if (q != 0)
+ rp[n] += DO_mpn_addlsh_n (rp, ap+n*q, n, s*q, scratch);
+ else {
+ rp[n] += mpn_add_n (rp, rp, ap, n);
+ break;
+ }
+ } while (1);
+ return abs_sub_add_n (rm, rp, ws, n + 1);
+}
+
+#ifdef WANT_TOOM_EV_MUL
+static int
+mpn_toom_ev_mul (mp_ptr rp, mp_ptr rm,
+ mp_srcptr ap, unsigned int q, mp_size_t n, mp_size_t t,
+ unsigned int s, mp_ptr ws)
+{
+ unsigned int m = s;
+ ASSERT( n >= t );
+ ASSERT( s > 2 );
+ ASSERT( q > 2 );
+ /* {ap,q*n+t} -> {rp,n+1} {rm,n+1} , with {ws, n+1}*/
+ ap += n*q;
+ {unsigned int i; for (i=q; --i; ) {m *= s;}; };
+ ASSERT( (m*s-1)/(s-1) <= GMP_NUMB_MAX );
+ if ((q & 1) == 0) {
+ rp[t] = mpn_mul_1(rp, ap, t, m);
+ MPN_ZERO (rp + t + 1, n - t);
+ ap-=n;m/=s;
+ ws[n] = mpn_mul_1(ws, ap, n, m);
+ ap-=n;m/=s;
+ rp[n] += mpn_addmul_1 (rp, ap, n, m);
+ } else {
+ ws[t] = mpn_mul_1(ws, ap, t, m);
+ MPN_ZERO(ws + t + 1, n - t);
+ ap-=n;m/=s;
+ rp[n] = mpn_mul_1(rp, ap, n, m);
+ }
+ do {
+ ap-=n;m/=s;
+ ws[n] += mpn_addmul_1 (ws, ap, n, m);
+ ap-=n;m/=s;
+ if (m != 1)
+ rp[n] += mpn_addmul_1 (rp, ap, n, m);
+ else {
More information about the gmp-commit
mailing list