* mpz/aorsmul.c: New file.

1 files changed, 156 insertions, 0 deletions

diff --git a/mpz/aorsmul.c b/mpz/aorsmul.c
new file mode 100644
index 000000000..5695d7e56
--- /dev/null
+++ b/mpz/aorsmul.c
@@ -0,0 +1,156 @@
+/* mpz_addmul, mpz_submul -- add or subtract multiple.
+
+Copyright 2001 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 2.1 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; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+
+/* expecting x and y both with non-zero high limbs */
+#define mpn_cmp_twosizes_lt(xp,xsize, yp,ysize)                 \
+  ((xsize) < (ysize)                                            \
+   || ((xsize) == (ysize) && mpn_cmp (xp, yp, xsize) < 0))
+
+
+/* sub>=0 means an addmul w += x*y, sub<0 means a submul w -= x*y.
+
+   The signs of w, x and y are fully accounted for by each flipping "sub".
+
+   The sign of w is retained for the result, unless the absolute value
+   submul underflows, in which case it flips.  */
+
+static void __gmpz_aorsmul _PROTO ((REGPARM_3_1 (mpz_ptr w, mpz_srcptr x, mpz_srcptr y, mp_size_t sub))) REGPARM_ATTR (1);
+#define mpz_aorsmul(w,x,y,sub)  __gmpz_aorsmul (REGPARM_3_1 (w, x, y, sub))
+
+static void
+mpz_aorsmul (mpz_ptr w, mpz_srcptr x, mpz_srcptr y, mp_size_t sub)
+{
+  mp_size_t  xsize, ysize, tsize, wsize, wsize_signed;
+  mp_ptr     wp, tp;
+  mp_limb_t  c, high;
+  TMP_DECL (marker);
+
+  /* w unaffected if x==0 or y==0 */
+  xsize = SIZ(x);
+  ysize = SIZ(y);
+  if (xsize == 0 || ysize == 0)
+    return;
+
+  /* make x the bigger of the two */
+  if (ABS(ysize) > ABS(xsize))
+    {
+      MPZ_SRCPTR_SWAP (x, y);
+      MP_SIZE_T_SWAP (xsize, ysize);
+    }
+
+  sub ^= ysize;
+  ysize = ABS(ysize);
+
+  /* use mpn_addmul_1/mpn_submul_1 if possible */
+  if (ysize == 1)
+    {
+      mpz_aorsmul_1 (w, x, PTR(y)[0], sub);
+      return;
+    }
+
+  sub ^= xsize;
+  xsize = ABS(xsize);
+
+  wsize_signed = SIZ(w);
+  sub ^= wsize_signed;
+  wsize = ABS(wsize_signed);
+
+  tsize = xsize + ysize;
+  MPZ_REALLOC (w, MAX (wsize, tsize) + 1);
+  wp = PTR(w);
+
+  if (wsize_signed == 0)
+    {
+      /* Nothing to add to, just set w=x*y.  No w==x or w==y overlap here,
+         since we know x,y!=0 but w==0.  */
+      high = mpn_mul (wp, PTR(x),xsize, PTR(y),ysize);
+      tsize -= (high == 0);
+      SIZ(w) = (sub >= 0 ? tsize : -tsize);
+      return;
+    }
+
+  TMP_MARK (marker);
+  tp = TMP_ALLOC_LIMBS (tsize);
+
+  high = mpn_mul (tp, PTR(x),xsize, PTR(y),ysize);
+  tsize -= (high == 0);
+  ASSERT (tp[tsize-1] != 0);
+  if (sub >= 0)
+    {
+      mp_srcptr up    = wp;
+      mp_size_t usize = wsize;
+
+      if (usize < tsize)
+        {
+          up    = tp;
+          usize = tsize;
+          tp    = wp;
+          tsize = wsize;
+
+          wsize = usize;
+        }
+
+      c = mpn_add (wp, up,usize, tp,tsize);
+      wp[wsize] = c;
+      wsize += (c != 0);
+    }
+  else
+    {
+      mp_srcptr up    = wp;
+      mp_size_t usize = wsize;
+
+      if (mpn_cmp_twosizes_lt (up,usize, tp,tsize))
+        {
+          up    = tp;
+          usize = tsize;
+          tp    = wp;
+          tsize = wsize;
+
+          wsize = usize;
+          wsize_signed = -wsize_signed;
+        }
+
+      ASSERT_NOCARRY (mpn_sub (wp, up,usize, tp,tsize));
+      wsize = usize;
+      MPN_NORMALIZE (wp, wsize);
+    }
+
+  SIZ(w) = (wsize_signed >= 0 ? wsize : -wsize);
+
+  TMP_FREE (marker);
+}
+
+
+void
+mpz_addmul (mpz_ptr w, mpz_srcptr u, mpz_srcptr v)
+{
+  mpz_aorsmul (w, u, v, (mp_size_t) 0);
+}
+
+void
+mpz_submul (mpz_ptr w, mpz_srcptr u, mpz_srcptr v)
+{
+  mpz_aorsmul (w, u, v, (mp_size_t) -1);
+}