/*
 * $Id: pb-muller1.c 47490 2011-11-10 17:07:08Z vinc17/ypig $
 *
 * Link with -lmpfr -lgmp *in that order*
 *
 * Example by Jean-Michel Muller:
 *   u_0 = 2
 *   u_1 = -4
 *   u_{n+1} = 111 - 1130 / u_n + 3000 / (u_n u_{n-1})
 * --> 6 (not 100)
 *
 * Arguments: precision, n
 */

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <mpfr.h>

static const mp_rnd_t rnd[3] = { MPFR_RNDD, MPFR_RNDN, MPFR_RNDU };
static int nep = 0;

static void check_sign (mpfr_t x[3])
{
  if (!nep && MPFR_SIGN (x[0]) != MPFR_SIGN (x[2]))
    nep = 1;
}

static void mdiv (mpfr_t tmp, unsigned long c, mpfr_t x[3], int r)
{
  int neg;

  neg = MPFR_SIGN (x[0]) < 0;
  mpfr_set_ui (tmp, c, rnd[neg ? 2-r : r]);
  mpfr_div (tmp, tmp, x[2-r], rnd[r]);
}

static void display (mpfr_t x[3], long k)
{
  int i;

  printf ("%6ld   ", k);
  for (i = 0; i < 3; i++)
    {
      size_t len;
      len = (i | nep) == i ? mpfr_out_str (stdout, 10, 17, x[i], rnd[i]) : 0;
      if (i < 2)
        while (len++ < 24)
          putchar (' ');
    }
  putchar ('\n');
}

int main (int argc, char *argv[])
{
  int prec, n, i, k;
  mpfr_t u[3], v[3], p[3], tmp;

  if (argc != 3)
    {
      fprintf (stderr, "Usage: pb-muller1 <prec> <n>\n");
      exit (EXIT_FAILURE);
    }

  prec = atoi (argv[1]);
  if (prec < MPFR_PREC_MIN || prec > MPFR_PREC_MAX)
    {
      fprintf (stderr, "pb-muller1: incorrect precision %d\n", prec);
      exit (EXIT_FAILURE);
    }
  mpfr_set_default_prec (prec);

  n = atoi (argv[2]);
  if (n <= 0)
    {
      fprintf (stderr, "pb-muller1: n must be positive\n");
      exit (EXIT_FAILURE);
    }

  for (i = 0; i < 3; i++)
    {
      mpfr_init (u[i]);
      mpfr_init (v[i]);
      mpfr_init (p[i]);
      mpfr_set_si (u[i], 2, rnd[i]);
      mpfr_set_si (v[i], -4, rnd[i]);
    }
  mpfr_init (tmp);

  display (u, 0);
  display (v, 1);

  for (k = 2; k <= n; k++)
    {
      int s;

      check_sign (v);
      s = MPFR_SIGN(u[1]) * MPFR_SIGN(v[1]);
      for (i = 0; i < 3; i++)
        if ((i | nep) == i)
          {
            mpfr_mul (p[s < 0 ? 2-i : i], u[i], v[i], rnd[i]);
            mpfr_set (u[i], v[i], rnd[i]);
          }
      check_sign (p);
      for (i = 0; i < 3; i++)
        if ((i | nep) == i)
          {
            mpfr_set_si (v[i], 111, rnd[i]);
            mdiv (tmp, 1130, u, 2-i);
            mpfr_sub (v[i], v[i], tmp, rnd[2-i]);
            mdiv (tmp, 3000, p, i);
            mpfr_add (v[i], v[i], tmp, rnd[i]);
          }
      display (v, k);
    }

  for (i = 0; i < 3; i++)
    {
      mpfr_clear (u[i]);
      mpfr_clear (v[i]);
      mpfr_clear (p[i]);
    }
  mpfr_clear (tmp);

  return 0;
}

/*
First 11 decimals of u_{100}: 6.0000000160, with prec >= 437,
and >= 577 to guarantee the result in interval arithmetic.
*/