/*
 * $Id: pb-newton.c 59279 2013-03-26 07:40:41Z vinc17/xvii $
 *
 * Compute phi with (1 + sqrt(5)) / 2 and with Newton's iteration.
 * Arguments: mode, precision n
 *   mode = 0: compute everything in precision n.
 *   mode = 1: double the precision at each iteration.
 *   mode = 2: ditto with better initial precision.
 *   mode = 3: ditto with the traditional iteration (avoid
 *             the division in the full target precision).
 *
 * f(x) = x^2 - x - 1, f'(x) = 2x - 1
 * x - f(x)/f'(x) = (x^2 + 1) / (2x - 1)
 */

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

int main (int argc, char *argv[])
{
  mpfr_t phi, x, y, num, den, s;
  int mode, n, i, p = 0, bp = 0, rp;
  clock_t c;
  double t;

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

  mode = atoi (argv[1]);
  if (mode < 0 || mode > 3)
    {
      fprintf (stderr, "pb-newton: incorrect mode (must be 0 or 1)\n");
      exit (EXIT_FAILURE);
    }

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

  c = clock ();
  mpfr_init2 (phi, n + 8);
  mpfr_init2 (s, 3);  /* s is a 3-bit precision number */
  mpfr_set_ui (s, 5, MPFR_RNDN);
  mpfr_sqrt (phi, s, MPFR_RNDN);
  mpfr_add_ui (phi, phi, 1, MPFR_RNDN);
  mpfr_div_2ui (phi, phi, 1, MPFR_RNDN);
  t = (double) (clock () - c) / CLOCKS_PER_SEC;
  printf ("Time: %.2f s for phi = (1 + sqrt(5)) / 2\n", t);

  if (mode >= 2)
    {
      bp = n;
      while (bp > 100)
        bp = (bp + 1) / 2;
    }

  c = clock ();
  mpfr_inits2 (mode ? 16 : n, x, y, num, den, (void *) 0);
  mpfr_set_ui (x, 1, MPFR_RNDN);
  for (i = 1; p < n; i++)
    {
      p = i == 3 ? 11 : 2 * p;
      if (bp && p >= bp)
        {
          p = bp;
          bp = 0;
        }
      if (p > n)
        p = n;
      if (mode && p > 16)
        mpfr_set_prec (y, p);
      mpfr_sqr (y, x, MPFR_RNDN);
      if (mode < 3)
        {
          mpfr_add_ui (y, y, 1, MPFR_RNDN);  /* x^2 + 1 */
          mpfr_mul_2ui (x, x, 1, MPFR_RNDN);
          mpfr_sub_ui (x, x, 1, MPFR_RNDN);  /* 2x - 1 */
          if (mode && p > 16)
            mpfr_prec_round (x, p, MPFR_RNDN);
          mpfr_div (x, y, x, MPFR_RNDN);     /* (x^2 + 1) / (2x - 1) */
        }
      else
        {
          mpfr_sub (y, y, x, MPFR_RNDN);
          mpfr_sub_ui (num, y, 1, MPFR_RNDN);
          mpfr_set_ui_2exp (s, 1, -1, MPFR_RNDN);
          mpfr_sub (den, x, s, MPFR_RNDN);
          mpfr_mul_2ui (den, den, 1, MPFR_RNDN);
          mpfr_div (num, num, den, MPFR_RNDN);  /* f(x)/f'(x) */
          if (p > 16)
            mpfr_prec_round (x, p, MPFR_RNDN);
          mpfr_sub (x, x, num, MPFR_RNDN);
          if (p > 16)
            {
              mpfr_set_prec (num, p);
              mpfr_set_prec (den, p);
            }
        }
      mpfr_sub (s, x, phi, MPFR_RNDN);
      rp = mpfr_get_exp (s) <= 0 ? (int) (1 - mpfr_get_exp (s)) : 0;
      printf ("i = %d, prec estimate = %d, real prec = %d\n", i, p, rp);
    }
  t = (double) (clock () - c) / CLOCKS_PER_SEC;
  printf ("Time: %.2f s\n", t);

  mpfr_clears (phi, x, y, num, den, s, (void *) 0);
  return 0;
}

/*
$ ./pb-newton 0 6000000
Time: 1.54 s for phi = (1 + sqrt(5)) / 2
i = 1, prec estimate = 0, real prec = 2
i = 2, prec estimate = 0, real prec = 5
i = 3, prec estimate = 11, real prec = 10
i = 4, prec estimate = 22, real prec = 21
i = 5, prec estimate = 44, real prec = 44
i = 6, prec estimate = 88, real prec = 88
i = 7, prec estimate = 176, real prec = 177
i = 8, prec estimate = 352, real prec = 355
i = 9, prec estimate = 704, real prec = 710
i = 10, prec estimate = 1408, real prec = 1421
i = 11, prec estimate = 2816, real prec = 2843
i = 12, prec estimate = 5632, real prec = 5686
i = 13, prec estimate = 11264, real prec = 11374
i = 14, prec estimate = 22528, real prec = 22748
i = 15, prec estimate = 45056, real prec = 45497
i = 16, prec estimate = 90112, real prec = 90995
i = 17, prec estimate = 180224, real prec = 181991
i = 18, prec estimate = 360448, real prec = 363982
i = 19, prec estimate = 720896, real prec = 727965
i = 20, prec estimate = 1441792, real prec = 1455930
i = 21, prec estimate = 2883584, real prec = 2911861
i = 22, prec estimate = 5767168, real prec = 5823723
i = 23, prec estimate = 6000000, real prec = 6000001
Time: 49.96 s
$ ./pb-newton 1 6000000
Time: 1.54 s for phi = (1 + sqrt(5)) / 2
i = 1, prec estimate = 0, real prec = 2
i = 2, prec estimate = 0, real prec = 5
i = 3, prec estimate = 11, real prec = 10
i = 4, prec estimate = 22, real prec = 22
i = 5, prec estimate = 44, real prec = 44
i = 6, prec estimate = 88, real prec = 88
i = 7, prec estimate = 176, real prec = 176
i = 8, prec estimate = 352, real prec = 353
i = 9, prec estimate = 704, real prec = 705
i = 10, prec estimate = 1408, real prec = 1410
i = 11, prec estimate = 2816, real prec = 2824
i = 12, prec estimate = 5632, real prec = 5635
i = 13, prec estimate = 11264, real prec = 11265
i = 14, prec estimate = 22528, real prec = 22529
i = 15, prec estimate = 45056, real prec = 45058
i = 16, prec estimate = 90112, real prec = 90114
i = 17, prec estimate = 180224, real prec = 180226
i = 18, prec estimate = 360448, real prec = 360450
i = 19, prec estimate = 720896, real prec = 720898
i = 20, prec estimate = 1441792, real prec = 1441793
i = 21, prec estimate = 2883584, real prec = 2883585
i = 22, prec estimate = 5767168, real prec = 5767168
i = 23, prec estimate = 6000000, real prec = 6000001
Time: 5.36 s
$ ./pb-newton 2 6000000
Time: 1.54 s for phi = (1 + sqrt(5)) / 2
i = 1, prec estimate = 0, real prec = 2
i = 2, prec estimate = 0, real prec = 5
i = 3, prec estimate = 11, real prec = 10
i = 4, prec estimate = 22, real prec = 22
i = 5, prec estimate = 44, real prec = 44
i = 6, prec estimate = 88, real prec = 88
i = 7, prec estimate = 92, real prec = 98
i = 8, prec estimate = 184, real prec = 186
i = 9, prec estimate = 368, real prec = 369
i = 10, prec estimate = 736, real prec = 736
i = 11, prec estimate = 1472, real prec = 1477
i = 12, prec estimate = 2944, real prec = 2946
i = 13, prec estimate = 5888, real prec = 5891
i = 14, prec estimate = 11776, real prec = 11778
i = 15, prec estimate = 23552, real prec = 23553
i = 16, prec estimate = 47104, real prec = 47107
i = 17, prec estimate = 94208, real prec = 94212
i = 18, prec estimate = 188416, real prec = 188418
i = 19, prec estimate = 376832, real prec = 376833
i = 20, prec estimate = 753664, real prec = 753666
i = 21, prec estimate = 1507328, real prec = 1507331
i = 22, prec estimate = 3014656, real prec = 3014656
i = 23, prec estimate = 6000000, real prec = 6000001
Time: 3.48 s
$ ./pb-newton 3 6000000
Time: 1.54 s for phi = (1 + sqrt(5)) / 2
i = 1, prec estimate = 0, real prec = 2
i = 2, prec estimate = 0, real prec = 5
i = 3, prec estimate = 11, real prec = 10
i = 4, prec estimate = 22, real prec = 22
i = 5, prec estimate = 44, real prec = 44
i = 6, prec estimate = 88, real prec = 88
i = 7, prec estimate = 92, real prec = 98
i = 8, prec estimate = 184, real prec = 186
i = 9, prec estimate = 368, real prec = 369
i = 10, prec estimate = 736, real prec = 736
i = 11, prec estimate = 1472, real prec = 1477
i = 12, prec estimate = 2944, real prec = 2946
i = 13, prec estimate = 5888, real prec = 5891
i = 14, prec estimate = 11776, real prec = 11778
i = 15, prec estimate = 23552, real prec = 23553
i = 16, prec estimate = 47104, real prec = 47107
i = 17, prec estimate = 94208, real prec = 94212
i = 18, prec estimate = 188416, real prec = 188418
i = 19, prec estimate = 376832, real prec = 376833
i = 20, prec estimate = 753664, real prec = 753666
i = 21, prec estimate = 1507328, real prec = 1507331
i = 22, prec estimate = 3014656, real prec = 3014656
i = 23, prec estimate = 6000000, real prec = 6000000
Time: 1.58 s
*/