New or modified function

[Paul Zimmermann]

       Hi Wojciech,

it's been a long time!

> From: nisse at lysator.liu.se (Niels =?iso-8859-1?Q?M=F6ller?=)
> Date: Thu, 19 Nov 2009 11:34:41 +0100
> 
> Wojciech Florek <Wojciech.Florek at amu.edu.pl> writes:
> 
> > exmul(a,b)=(a/g)*(b/g).
> 
> I think it would be simplest to implement this on top of the mpz
> interface, not messing with internals and things like TMP_DECL.
> 
> > BTW, I wonder if exmul=(a*b)/(g^2) is faster (2 divexact and 1 mul
> > versus 1 mul, 1 square and 1 divexact, but for larger numbers).
> 
> I'd expect (a/g) * (b/g) to be a lot faster when g is large (since the
> running time for both the divisions and the product essentially depend
> on the size of the quotients, which are small when g is large), while
> (a*b)/g^2 may be sligtly faster when g is small (since in this case
> the quotients are large, and the cost of each division will be more
> significant).
> 
> But real benchmarking is the only way to find out for sure...
> 
> Regards,
> /Niels

real benchmarking says that (a/g)*(b/g) is always faster:

patate% ./a.out 100000 100000
a has 200000 limbs, b has 200000 limbs, g has 100000 limbs
(a/g)*(b/g) took 1214ms
(a*b)/(g^2) took 1593ms
patate% ./a.out 1000000 1
a has 1000001 limbs, b has 1000001 limbs, g has 1 limbs
(a/g)*(b/g) took 1036ms
(a*b)/(g^2) took 1072ms
patate% ./a.out 1 1000000
a has 1000001 limbs, b has 1000001 limbs, g has 1000000 limbs
(a/g)*(b/g) took 0ms
(a*b)/(g^2) took 1843ms

Paul Zimmermann

#include "stdio.h"
#include "gmp.h"
#include <sys/types.h>
#include <sys/resource.h>

int
cputime ()
{
  struct rusage rus;

  getrusage (0, &rus);
  return rus.ru_utime.tv_sec * 1000 + rus.ru_utime.tv_usec / 1000;
}

int
main (int argc, char *argv[])
{
  mpz_t g, a, b, c;
  size_t n = atoi (argv[1]);
  size_t m = atoi (argv[2]);
  int st;
  mpz_init (a);
  mpz_init (b);
  mpz_init (c);
  mpz_init (g);
  mpz_random (a, n);
  mpz_random (b, n);
  mpz_random (g, m);
  mpz_mul (a, a, g);
  mpz_mul (b, b, g);
  printf ("a has %lu limbs, b has %lu limbs, g has %lu limbs\n",
          mpz_size (a), mpz_size (b), mpz_size (g));
  st = cputime ();
  mpz_divexact (a, a, g);
  mpz_divexact (b, b, g);
  mpz_mul (c, a, b);
  printf ("(a/g)*(b/g) took %dms\n", cputime () - st);
  mpz_mul (a, a, g);
  mpz_mul (b, b, g);
  st = cputime ();
  mpz_mul (c, a, b);
  mpz_mul (g, g, g);
  mpz_divexact (c, c, g);
  printf ("(a*b)/(g^2) took %dms\n", cputime () - st);
}