Added mul_si

git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@26 211d60ee-9f03-0410-a15a-8952a2c7a4e4

5 files changed, 182 insertions, 6 deletions

diff --git a/makefile b/makefile
index f2f0bf7..03bf9a7 100644
--- a/makefile
+++ b/makefile
@@ -41,8 +41,8 @@ VERSION=0.4.5
 INCLUDES=-I. -I$(MPFR)/include -I$(GMP)/include
 LIBS=$(MPFR)/lib/libmpfr.a $(GMP)/lib/libgmp.a
 
-OBJECTS= abs.o add.o add_fr.o add_ui.o clear.o cmp.o conj.o div.o div_2exp.o div_fr.o div_ui.o exp.o init.o init2.o init3.o inp_str.o mul.o mul_2exp.o mul_fr.o mul_ui.o neg.o norm.o out_str.o random.o random2.o set.o set_d_d.o set_dfl_prec.o set_prec.o set_ui_fr.o set_si_si.o set_ui_ui.o sqr.o sqrt.o sub.o sub_ui.o ui_div.o uceil_log2.o
-SOURCES= abs.c add.c add_fr.c add_ui.c clear.c cmp.c conj.c div.c div_2exp.c div_fr.c div_ui.c exp.c init.c init2.c init3.c inp_str.c mul.c mul_2exp.c mul_fr.c mul_ui.c neg.c norm.c out_str.c random.c random2.c set.c set_d_d.c set_dfl_prec.c set_prec.c set_ui_fr.c set_si_si.c set_ui_ui.c sqr.c sqrt.c sub.c sub_ui.c ui_div.c uceil_log2.c
+OBJECTS= abs.o add.o add_fr.o add_ui.o clear.o cmp.o conj.o div.o div_2exp.o div_fr.o div_ui.o exp.o init.o init2.o init3.o inp_str.o mul.o mul_2exp.o mul_fr.o mul_ui.o mul_si.o neg.o norm.o out_str.o random.o random2.o set.o set_d_d.o set_dfl_prec.o set_prec.o set_ui_fr.o set_si_si.o set_ui_ui.o sqr.o sqrt.o sub.o sub_ui.o ui_div.o uceil_log2.o
+SOURCES= abs.c add.c add_fr.c add_ui.c clear.c cmp.c conj.c div.c div_2exp.c div_fr.c div_ui.c exp.c init.c init2.c init3.c inp_str.c mul.c mul_2exp.c mul_fr.c mul_ui.c mul_si.c neg.c norm.c out_str.c random.c random2.c set.c set_d_d.c set_dfl_prec.c set_prec.c set_ui_fr.c set_si_si.c set_ui_ui.c sqr.c sqrt.c sub.c sub_ui.c ui_div.c uceil_log2.c
 TESTS= test.c tmul.c tsqr.c tdiv.c texp.c
 DIST= $(SOURCES) $(TESTS) makefile mpc.h mpc-impl.h COPYING.LIB mpc.texi INSTALL
 
diff --git a/mpc.h b/mpc.h
index b677edd..b99ac56 100644
--- a/mpc.h
+++ b/mpc.h
@@ -1,6 +1,6 @@
 /* mpc.h -- Include file for mpc.
 
-Copyright (C) 2002, 2003, 2004 Andreas Enge, Paul Zimmermann
+Copyright (C) 2002, 2003, 2004, 2005 Andreas Enge, Paul Zimmermann
 
 This file is part of the MPC Library.
 
@@ -108,6 +108,7 @@ int  mpc_sub_ui __MPC_PROTO ((mpc_ptr, mpc_srcptr, unsigned long int, mpc_rnd_t)
 int  mpc_mul    __MPC_PROTO ((mpc_ptr, mpc_srcptr, mpc_srcptr, mpc_rnd_t));
 int  mpc_mul_fr __MPC_PROTO ((mpc_ptr, mpc_srcptr, mpfr_srcptr, mpc_rnd_t));
 int  mpc_mul_ui __MPC_PROTO ((mpc_ptr, mpc_srcptr, unsigned long int, mpc_rnd_t));
+int  mpc_mul_si __MPC_PROTO ((mpc_ptr, mpc_srcptr, long int, mpc_rnd_t));
 int  mpc_sqr    __MPC_PROTO ((mpc_ptr, mpc_srcptr, mpc_rnd_t));
 int  mpc_div    __MPC_PROTO ((mpc_ptr, mpc_srcptr, mpc_srcptr, mpc_rnd_t));
 int  mpc_div_fr __MPC_PROTO ((mpc_ptr, mpc_srcptr, mpfr_srcptr, mpc_rnd_t));
diff --git a/mpc.texi b/mpc.texi
index 59147fe..4b1c0ec 100644
--- a/mpc.texi
+++ b/mpc.texi
@@ -30,7 +30,7 @@ END-INFO-DIR-ENTRY
 @ifinfo
 This file documents MPC, a library for multiple precision complex arithmetic
 
-Copyright (C) 2002, 2003, 2004 Andreas Enge, Paul Zimmermann
+Copyright (C) 2002, 2003, 2004, 2005 Andreas Enge, Paul Zimmermann
 
 Permission is granted to make and distribute verbatim copies of
 this manual provided the copyright notice and this permission notice
@@ -76,7 +76,7 @@ by the Foundation.
 
 @page
 @vskip 0pt plus 1filll
-Copyright @copyright{} 2002, 2003, 2004 Andreas Enge, Paul Zimmermann
+Copyright @copyright{} 2002, 2003, 2004, 2005 Andreas Enge, Paul Zimmermann
 
 @sp 2
 
@@ -596,6 +596,7 @@ Set @var{rop} to @var{op1} @minus{} @var{op2} rounded according to @var{rnd}.
 
 @deftypefun int mpc_mul (mpc_t @var{rop}, mpc_t @var{op1}, mpc_t @var{op2}, mpc_rnd_t @var{rnd})
 @deftypefunx int mpc_mul_ui (mpc_t @var{rop}, mpc_t @var{op1}, unsigned long int @var{op2}, mpc_rnd_t @var{rnd})
+@deftypefunx int mpc_mul_si (mpc_t @var{rop}, mpc_t @var{op1}, long int @var{op2}, mpc_rnd_t @var{rnd})
 @deftypefunx int mpc_mul_fr (mpc_t @var{rop}, mpc_t @var{op1}, mpfr_t @var{op2}, mpc_rnd_t @var{rnd})
 Set @var{rop} to @var{op1} times @var{op2} rounded according to @var{rnd}.
 @end deftypefun
diff --git a/mul_si.c b/mul_si.c
new file mode 100644
index 0000000..81f30e5
--- /dev/null
+++ b/mul_si.c
@@ -0,0 +1,54 @@
+/* mpc_mul_si -- Multiply a complex number by a signed integer.
+
+Copyright (C) 2005 Andreas Enge
+
+This file is part of the MPC Library.
+
+The MPC 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 MPC 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 MPC 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 "mpfr.h"
+#include "mpc.h"
+#include "mpc-impl.h"
+
+int
+mpc_mul_si (mpc_ptr a, mpc_srcptr b, long int c, mpc_rnd_t rnd)
+{
+  int inex_re, inex_im;
+
+  /* We could use mpfr_mul_si with newer versions of mpfr (2.1.0), but 
+     not with the version contained in gmp-4.1.4. So copy and paste for 
+     the time being                                                     */
+  if (c >= 0)
+  /* inlining by copy and paste */
+  {
+     inex_re = mpfr_mul_ui (MPC_RE(a), MPC_RE(b), (unsigned long int) c,
+           MPC_RND_RE(rnd));
+     inex_im = mpfr_mul_ui (MPC_IM(a), MPC_IM(b), (unsigned long int) c,
+           MPC_RND_IM(rnd));
+  }
+  else
+  {
+     inex_re = -mpfr_mul_ui (MPC_RE(a), MPC_RE(b), (unsigned long int) (-c),
+           INV_RND (MPC_RND_RE(rnd)));
+     MPFR_CHANGE_SIGN (MPC_RE (a));
+     inex_im = -mpfr_mul_ui (MPC_IM(a), MPC_IM(b), (unsigned long int) (-c),
+           INV_RND (MPC_RND_IM(rnd)));
+     MPFR_CHANGE_SIGN (MPC_IM (a));
+  }
+  
+  return MPC_INEX(inex_re, inex_im);
+}
diff --git a/tmul.c b/tmul.c
index f541cab..d35c5f0 100644
--- a/tmul.c
+++ b/tmul.c
@@ -1,6 +1,6 @@
 /* tmul -- test file for mpc_mul.
 
-Copyright (C) 2002 Andreas Enge, Paul Zimmermann
+Copyright (C) 2002, 2005 Andreas Enge, Paul Zimmermann
 
 This file is part of the MPC Library.
 
@@ -28,6 +28,8 @@ MA 02111-1307, USA. */
 #include "mpc-impl.h"
 
 void cmpmul (mpc_srcptr, mpc_srcptr, mpc_rnd_t);
+void cmpmului (mpc_srcptr, unsigned long int, mpc_rnd_t);
+void cmpmulsi (mpc_srcptr, long int, mpc_rnd_t);
 void testmul (long, long, long, long, mp_prec_t, mpc_rnd_t);
 void special (void);
 void timemul (void);
@@ -184,6 +186,117 @@ void cmpmul (mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
 }
 
 
+void cmpmului (mpc_srcptr x, unsigned long int y, mpc_rnd_t rnd)
+   /* computes the product of x and y using mpc_mul_fr and mpc_mul_ui      */
+   /* using the rounding mode rnd and compares the results and return      */
+   /* values.                                                              */
+   /* In our current test suite, the real and imaginary parts of x have    */
+   /* the same precision, and we use this precision also for the result.   */
+
+{
+   mpc_t z, t;
+   mpfr_t yf;
+   int   inexact_z, inexact_t;
+
+   mpc_init2 (z, MPC_MAX_PREC (x));
+   mpc_init2 (t, MPC_MAX_PREC (x));
+   mpfr_init2 (yf, 8 * sizeof (long int));
+   mpfr_set_si (yf, y, GMP_RNDN);
+   
+   inexact_z = mpc_mul_fr (z, x, yf, rnd);
+   inexact_t = mpc_mul_ui (t, x, y, rnd);
+
+   if (mpc_cmp (z, t))
+   {
+      fprintf (stderr, "mul_fr and mul_ui differ for rnd=(%s,%s) \nx=",
+               mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
+               mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
+      mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
+      fprintf (stderr, "\nand y=%li", y);
+      fprintf (stderr, "\nmpc_mul_fr gives ");
+      mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);
+      fprintf (stderr, "\nmpc_mul_ui gives ");
+      mpc_out_str (stderr, 2, 0, t, MPC_RNDNN);
+      fprintf (stderr, "\n");
+      exit (1);
+   }
+   if (inexact_z != inexact_t)
+   {
+      fprintf (stderr, "The return values of mul_fr and mul_ui differ for ");
+      fprintf (stderr, "rnd=(%s,%s) \nx=",
+               mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
+               mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
+      mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
+      fprintf (stderr, "\nand y=%li", y);
+      fprintf (stderr, "\nand x*y=");
+      mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);
+      fprintf (stderr, "\nmpc_mul_fr gives %i", inexact_z);
+      fprintf (stderr, "\nmpc_mul_ui gives %i", inexact_t);
+      fprintf (stderr, "\n");
+      exit (1);
+   }
+   
+   mpc_clear (z);
+   mpc_clear (t);
+   mpfr_clear (yf);
+}
+
+void cmpmulsi (mpc_srcptr x, long int y, mpc_rnd_t rnd)
+   /* computes the product of x and y using mpc_mul_fr and mpc_mul_si      */
+   /* using the rounding mode rnd and compares the results and return      */
+   /* values.                                                              */
+   /* In our current test suite, the real and imaginary parts of x have    */
+   /* the same precision, and we use this precision also for the result.   */
+
+{
+   mpc_t z, t;
+   mpfr_t yf;
+   int   inexact_z, inexact_t;
+
+   mpc_init2 (z, MPC_MAX_PREC (x));
+   mpc_init2 (t, MPC_MAX_PREC (x));
+   mpfr_init2 (yf, 8 * sizeof (long int));
+   mpfr_set_si (yf, y, GMP_RNDN);
+   
+   inexact_z = mpc_mul_fr (z, x, yf, rnd);
+   inexact_t = mpc_mul_si (t, x, y, rnd);
+
+   if (mpc_cmp (z, t))
+   {
+      fprintf (stderr, "mul_fr and mul_si differ for rnd=(%s,%s) \nx=",
+               mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
+               mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
+      mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
+      fprintf (stderr, "\nand y=%li", y);
+      fprintf (stderr, "\nmpc_mul_fr gives ");
+      mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);
+      fprintf (stderr, "\nmpc_mul_si gives ");
+      mpc_out_str (stderr, 2, 0, t, MPC_RNDNN);
+      fprintf (stderr, "\n");
+      exit (1);
+   }
+   if (inexact_z != inexact_t)
+   {
+      fprintf (stderr, "The return values of mul_fr and mul_si differ for ");
+      fprintf (stderr, "rnd=(%s,%s) \nx=",
+               mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
+               mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
+      mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
+      fprintf (stderr, "\nand y=%li", y);
+      fprintf (stderr, "\nand x*y=");
+      mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);
+      fprintf (stderr, "\nmpc_mul_fr gives %i", inexact_z);
+      fprintf (stderr, "\nmpc_mul_si gives %i", inexact_t);
+      fprintf (stderr, "\n");
+      exit (1);
+   }
+   
+   mpc_clear (z);
+   mpc_clear (t);
+   mpfr_clear (yf);
+}
+
+
 void
 testmul (long a, long b, long c, long d, mp_prec_t prec, mpc_rnd_t rnd)
 {
@@ -291,6 +404,7 @@ main()
   mpc_rnd_t rnd_re, rnd_im;
   mp_prec_t prec;
   int i;
+  long int ysi;
 
 /*
   timemul ();
@@ -323,10 +437,16 @@ main()
 
       mpc_random (x);
       mpc_random (y);
+      ysi = rand ();
 
       for (rnd_re = 0; rnd_re < 4; rnd_re ++)
         for (rnd_im = 0; rnd_im < 4; rnd_im ++)
+         {
             cmpmul (x, y, RNDC(rnd_re, rnd_im));
+            cmpmului (x, (unsigned long int) ysi, RNDC(rnd_re, rnd_im));
+            cmpmulsi (x, ysi, RNDC(rnd_re, rnd_im));
+            cmpmulsi (x, -ysi, RNDC(rnd_re, rnd_im));
+         }
     }
   }