mini-gmp and mpq

[Marco Bodrato]

Ciao,

Il Lun, 26 Febbraio 2018 9:03 pm, Bradley Lucier ha scritto:
> 2.  I stumbled on this issue when approximating the golden ratio $\phi$,
> the root of $\phi^2-\phi-1=0$, using the naive iteration
> $\phi_{n+1}=1+1/\phi_n$ with $\phi_1=1$, or, in Scheme:

Ok, I tested the following code with the naïve implementation I sent to
this list.

{
  mpq_t  x, u;
  int    i;

  mpq_init (x);
  mpq_init (u);

  mpq_set_ui (x, 1, 1);
  mpq_set (u, x);
  for (i = 1000; --i >= 0; )
    {
      mpq_inv (x, x);
      mpq_add (x, u, x);
    }

  printf ("%lu / %lu -> %g\n",
	  mpz_sizeinbase (mpq_numref(x),2),
	  mpz_sizeinbase (mpq_denref(x),2),
	  mpq_get_d (x));

  mpq_clear (x);
  mpq_clear (u);
}

> Mzscheme:
>
>  > (time (fib 1000))
> cpu time: 379490 real time: 454623 gc time: 450
>
> Gambit:
>
> (time (fib 1000))
>      8 ms real time
>      10 ms cpu time (10 user, 0 system)
>
> So these numbers weren't really that big (about 700 bits in numerator
> and denominator) and there was about a factor of 40,000 difference in
> performance.

On my computer (a slow one!) I obtain:

$ time mini-gmp/tests/t-mpq_gr
695 / 694 -> 1.61803
real	0m0.076s
user	0m0.070s

I'll be happy to know if your measures differs.

Of course, if we handle larger fraction, I mean
  for (i = 15000; --i >= 0; )
    {
      mpq_inv (x, x);
      mpq_add (x, u, x);
    }
the timings are longer.

$ time mini-gmp/tests/t-mpq_gr
10414 / 10414 -> 1.61803
real	0m47.719s
user	0m46.240s

With a not-so-naïve implementation of mini-mpq the timings are reduced to
$ time mini-gmp/tests/t-mpq_gr10414 / 10414 -> 1.61803
real	0m0.212s
user	0m0.200s

A factor 250 for sizes larger than 10'000 bits, by substituting the simple
8-lines mpq_add with the 24 lines I attach... Is it worth doing?

The realway to speed-up that loop is to use the suggestion in gmpxx.h to
compute a rational+1, both for mini-gmp and the full GMP library:

struct __gmp_unary_increment
{
  static void eval(mpz_ptr z) { mpz_add_ui(z, z, 1); }
  static void eval(mpq_ptr q)
  { mpz_add(mpq_numref(q), mpq_numref(q), mpq_denref(q)); }
  static void eval(mpf_ptr f) { mpf_add_ui(f, f, 1); }
};


Ĝis,
m

---------------------------
void
mpq_add (mpq_t r, const mpq_t a, const mpq_t b)
{
  mpz_t t;

  mpz_init (t);
  mpz_gcd (t, DEN(a), DEN(b));
  if (mpz_cmp_ui (t, 1) == 0)
    {
      mpz_mul (t, NUM(a), DEN(b));
      mpz_addmul (t, NUM(b), DEN(a));
      mpz_mul (DEN(r), DEN(a), DEN(b));
      mpz_swap (NUM(r), t);
    }
  else
    {
      mpz_t x, y;
      mpz_init (x);
      mpz_init (y);

      mpz_div_qr (x, NULL, DEN(b), t, GMP_DIV_TRUNC);
      mpz_div_qr (y, NULL, DEN(a), t, GMP_DIV_TRUNC);
      mpz_mul (x, NUM(a), x);
      mpz_addmul (x, NUM (b), y);

      mpz_gcd (t, x, t);
      mpz_div_qr (NUM (r), NULL, x, t, GMP_DIV_TRUNC);
      mpz_div_qr (x, NULL, DEN(b), t, GMP_DIV_TRUNC);
      mpz_mul (DEN (r), x, y);

      mpz_clear (x);
      mpz_clear (y);
    }
  mpz_clear (t);
}


-- 
http://bodrato.it/