/* Copyright 2006, 2007 Free Software Foundation, Inc.

This program is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 3 of the License, or (at your option) any later
version.

This program 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 General Public License for more details.

You should have received a copy of the GNU General Public License along with
this program.  If not, see http://www.gnu.org/licenses/.  */


#include <stdlib.h>		/* for exit, strtoul */
#include <string.h>		/* for strlen */
#include <stdio.h>		/* for printf */
#include <math.h>		/* for log10, fmod, pow */
#include <unistd.h>		/* for isatty */

#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"

#include "new.h"


static mp_size_t
sizearg (char *arg)
{
  if (arg[strlen (arg) - 1] == 'b')
    return strtoul (arg, 0, 0);
  else
    return strtoul (arg, 0, 0) * GMP_LIMB_BITS;
}


#ifdef CHECK
void
dumpy (mp_srcptr p, mp_size_t n)
{
  mp_size_t i;
  if (n > 20)
    {
      for (i = n - 1; i >= n - 4; i--)
	{
	  printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]);
	  printf (" ");
	}
      printf ("... ");
      for (i = 3; i >= 0; i--)
	{
	  printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]);
	  printf (" " + (i == 0));
	}
    }
  else
    {
      for (i = n - 1; i >= 0; i--)
	{
	  printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]);
	  printf (" " + (i == 0));
	}
    }
  puts ("");
}

unsigned long test, testh;
int err;

mp_limb_t rran0, rran1, qran0, qran1;

void
check_one (mp_srcptr qrefp, mp_srcptr rrefp, mp_ptr qp, mp_srcptr rp,
	   mp_srcptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn,
	   char *fname, mp_limb_t q_allowed_err)
{
  mp_size_t qn = nn - dn + 1;
  int qcmp, rcmp;
  mp_ptr tp;
  mp_size_t i;
  TMP_DECL;
  TMP_MARK;

  qcmp = mpn_cmp (qp, qrefp, qn);
  if (qcmp != 0 && rp == rrefp)	/* FIXME: rp=rrefp indicates we check a approx-q fn */
    {
      mp_limb_t qerr;

      tp = TMP_ALLOC_LIMBS (qn);
      qerr = mpn_sub_n (tp, qp, qrefp, qn);

      for (i = 1; i < qn; i++)
	qerr |= tp[i];

      qcmp = qerr | (tp[0] > q_allowed_err);
    }

  rcmp = mpn_cmp (rp, rrefp, dn);

  if (qcmp != 0 || rcmp != 0)
    {
      mp_size_t lo, hi;
      printf ("\r*******************************************************************************\n");
      printf ("%s and mpn_sb_divrem_mn disagree in test %lu,,%lu\n",
	      fname, testh, test);
      printf ("N=   "); dumpy (np, nn);
      printf ("D=   "); dumpy (dp, dn);
      printf ("Q=   "); dumpy (qp, qn);
      printf ("refQ="); dumpy (qrefp, qn);
      printf ("R=   "); dumpy (rp, dn);
      printf ("refR="); dumpy (rrefp, dn);
      printf ("nn = %ld, dn = %ld, qn = %ld", nn, dn, nn-dn+1);
      if (qcmp != 0)
	{
	  for (lo = 0; qp[lo] == qrefp[lo]; lo++);
	  for (hi = nn - dn; qp[hi] == qrefp[hi]; hi--);
	  printf (", quotient bad %ld--%ld", hi, lo);
	}
      if (rcmp != 0)
	{
	  for (lo = 0; rp[lo] == rrefp[lo]; lo++);
	  for (hi = dn - 1; rp[hi] == rrefp[hi]; hi--);
	  printf (", remainder bad %ld--%ld", hi, lo);
	}
      printf ("\n*******************************************************************************\n");
      if (++err >= 3)
	abort ();
    }

  TMP_FREE;
}

#include <getopt.h>

struct nameflag {char *name; int *flagp;};

int flag_mpn_sb_div_qr, flag_mpn_sb_divappr_q, flag_mpn_sb_div_q;
int flag_mpn_dc_div_qr, flag_mpn_dc_divappr_q, flag_mpn_dc_div_q;
int flag_mpn_mu_div_qr, flag_mpn_mu_divappr_q, flag_mpn_mu_div_q;

static struct nameflag foo[] = {
  {"mpn_sb_div_qr",     &flag_mpn_sb_div_qr},
  {"mpn_sb_divappr_q",  &flag_mpn_sb_divappr_q},
  {"mpn_sb_div_q",      &flag_mpn_sb_div_q},
  {"mpn_dc_div_qr",     &flag_mpn_dc_div_qr},
  {"mpn_dc_divappr_q",  &flag_mpn_dc_divappr_q},
  {"mpn_dc_div_q",      &flag_mpn_dc_div_q},
  {"mpn_mu_div_qr",     &flag_mpn_mu_div_qr},
  {"mpn_mu_divappr_q",  &flag_mpn_mu_divappr_q},
  {"mpn_mu_div_q",      &flag_mpn_mu_div_q},
  {NULL,  NULL}
};

int testname;

static struct option longopts[] = {
  {"test",   required_argument,    &testname, 1},
  {NULL,     0,                    NULL,      0}
};

int
main (int argc, char **argv)
{
  gmp_randstate_t rs;
  unsigned long maxnbits, maxdbits, nbits, dbits;
  mpz_t n, d;
  mp_size_t maxnn, maxdn, nn, dn, clearn, i;
  mp_ptr np, dp, qp, rp, qrefp, rrefp;
  int work = 0;
  mp_limb_t di_32, di_21, dip[2], *dixp;
  mp_limb_t cy;
  int specific_test_opt_seen = 0;
  char *progname = argv[0];
  TMP_DECL;

  /* Default to testing all functions.  */
  for (i = 0; foo[i].name != 0; i++)
    *(foo[i].flagp) = 1;

  for (;;)
    {
      int ch;
      testname = 0;
      ch = getopt_long (argc, argv, "", longopts, NULL);
      if (ch == -1)
	break;
      if (testname)
	{
	  if (! specific_test_opt_seen)
	    for (i = 0; foo[i].name != 0; i++)
	      *(foo[i].flagp) = 0;

	  specific_test_opt_seen = 1;

	  for (i = 0; foo[i].name != 0; i++)
	    {
	      if (strcmp (foo[i].name, optarg) == 0)
		{
		  *(foo[i].flagp) = 1;
		  break;
		}
	    }
	  if (foo[i].name == 0)
	    {
	      printf ("%s: no function matching argument %s to -test option\n",
		      progname, optarg);
	      exit (1);
	    }
	}
    }
  argc -= optind;
  argv += optind;

  if (argc == 1)
    {
      maxdbits = sizearg (argv[0]);
      maxnbits = 2 * maxdbits;
    }
  else if (argc == 2)
    {
      maxnbits = sizearg (argv[0]);
      maxdbits = sizearg (argv[1]);
    }
  else
    {
      printf ("usage: %s nsize dsize\n", progname);
      printf ("   or: %s size\n", progname);
      exit (1);
    }

  if (maxnbits <= maxdbits)
    {
      printf ("the dividend needs to be larger than the divisor\n");
      exit (1);
    }

  if (maxdbits <= GMP_NUMB_BITS)
    {
      printf ("the divisor needs to be at least 2 limbs (%d bits)\n",
	      1 + GMP_NUMB_BITS);
      exit (1);
    }

  gmp_randinit_default (rs);

  mpz_init (n);
  mpz_init (d);

  maxnn = maxnbits / GMP_NUMB_BITS + 1;
  maxdn = maxdbits / GMP_NUMB_BITS + 1;

  TMP_MARK;

  qp = TMP_ALLOC_LIMBS (maxnn + 2) + 1;
  rp = TMP_ALLOC_LIMBS (maxnn + 2) + 1;
  qrefp = TMP_ALLOC_LIMBS (maxnn);
  rrefp = TMP_ALLOC_LIMBS (maxnn);
  dixp = TMP_ALLOC_LIMBS (maxdn + 1);

  testh = 0;
  for (test = 0;;)
    {
#ifndef FIXED_SIZE
      nbits = random () % maxnbits + 1;
      if (maxdbits > nbits)
	dbits = random () % nbits + 1;
      else
	dbits = random () % maxdbits + 1;
#else
      nbits = maxnbits;
      dbits = maxdbits;
#endif

#if RAND_UNIFORM
	  mpz_urandomb (n, rs, nbits);
#else
	  mpz_rrandomb (n, rs, nbits);
#endif
      do
	{
#if RAND_UNIFORM
	  mpz_urandomb (d, rs, dbits);
#else
	  mpz_rrandomb (d, rs, dbits);
#endif
	}
      while (mpz_sgn (d) == 0);

      np = PTR (n);
      dp = PTR (d);
      nn = SIZ (n);
      dn = SIZ (d);

      dp[dn - 1] |= GMP_NUMB_HIGHBIT;

#if LEAD_D_MAX
      dp[dn - 1] = GMP_NUMB_MAX;
#endif

      if (random() % 2 == 0)
	{
	  clearn = random () % (nn + 1);
	  for (i = clearn; i < nn; i++)
	    np[i] = 0;
	}

      if (mpz_cmp (n, d) < 0)
	continue;

      if (dn < 2)
	continue;

      test++;
      testh += test == 0;

      work += nbits - dbits;
      if (work >= 3123451)
	{
	  if (isatty (fileno (stdout)))
	    printf ("\r%lu,,%lu", testh, test);
#ifdef DEBUG
	  printf ("\n");
#endif
	  fflush (stdout);
	  work = 0;
	}

#ifdef DEBUG
      printf ("n="); mpn_dump (np, nn);
      printf ("d="); mpn_dump (dp, dn);
#endif

      mpn_tdiv_qr (qrefp, rrefp, 0, np, nn, dp, dn);

      mp_limb_t xp[2];
      cy = mpn_add_1 (xp, dp + dn - 2, 2, 1);
      if (cy != 0)
	{
	  dip[0] = dip[1] = 0;
	  di_32 = 0;
	}
      else
	{
	  mpn_invert (dip, xp, 2, NULL);
	  di_32 = dip[1];
	}
      invert_limb (di_21, dp[dn - 1]);

      MPN_COPY (rp + 1, dp, dn);
      rp[0] = 0;
      mpn_invert (dixp, rp, dn + 1, NULL);

      mpn_random (&rran0, 1);
      mpn_random (&rran1, 1);
      mpn_random (&qran0, 1);
      mpn_random (&qran1, 1);

      qp[-1] = qran0;
      qp[nn - dn + 1] = qran1;
      rp[-1] = rran0;

      if ((double) (nn - dn) * dn < 1e8)
	{
#if 1
	  if (1)
	    {
	      MPN_COPY (rp, np, nn);
	      MPN_ZERO (qp, nn - dn);
	      qp[nn - dn] = mpn_div_n_qr_pi1_32 (qp, rp, nn, dp, dn, di_32);
	      check_one (qrefp, rrefp, qp, rp, np, nn, dp, dn, "mpn_div_n_qr_pi1_32", 0);
	    }
	  if (dn >= 3)
	    {
	      MPN_COPY (rp, np, nn);
	      MPN_ZERO (qp, nn - dn);
	      qp[nn - dn] = mpn_div_n_qr_pi1_21 (qp, rp, nn, dp, dn, di_21);
	      check_one (qrefp, rrefp, qp, rp, np, nn, dp, dn, "mpn_div_n_qr_pi1_21", 0);
	    }
	  if (0 && nn - dn <= dn)
	    {
	      MPN_COPY (rp, np, nn);
	      MPN_ZERO (qp, nn - dn);
	      qp[nn - dn] = mpn_div_n_qr_pin (qp, rp, nn, dp, dn, dixp);
	      check_one (qrefp, rrefp, qp, rp, np, nn, dp, dn, "mpn_div_n_qr_pin", 0);
	    }
#endif
	  if (flag_mpn_sb_div_qr)
	    {
	      MPN_COPY (rp, np, nn);
	      MPN_ZERO (qp, nn - dn);
	      qp[nn - dn] = mpn_sb_div_qr (qp, rp, nn, dp, dn, dip);
	      check_one (qrefp, rrefp, qp, rp, np, nn, dp, dn, "mpn_sb_div_qr", 0);
	    }

	  if (flag_mpn_sb_divappr_q)
	    {
	      MPN_COPY (rp, np, nn);
	      MPN_ZERO (qp, nn - dn);
	      qp[nn - dn] = mpn_sb_divappr_q (qp, rp, nn, dp, dn, dip);
	      check_one (qrefp, rrefp, qp, rrefp, np, nn, dp, dn, "mpn_sb_divappr_q", 1);
	    }

	  if (flag_mpn_sb_div_q)
	    {
	      MPN_COPY (rp, np, nn);
	      MPN_ZERO (qp, nn - dn);
	      qp[nn - dn] = mpn_sb_div_q (qp, rp, nn, dp, dn, dip);
	      check_one (qrefp, rrefp, qp, rrefp, np, nn, dp, dn, "mpn_sb_div_q", 0);
	    }
	}

      if (flag_mpn_dc_div_qr)
	{
	  MPN_COPY (rp, np, nn);
	  MPN_ZERO (qp, nn - dn);
	  qp[nn - dn] = mpn_dc_div_qr (qp, rp, nn, dp, dn);
	  ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
	  ASSERT_ALWAYS (rp[-1] == rran0);
	  check_one (qrefp, rrefp, qp, rp, np, nn, dp, dn, "mpn_dc_div_qr", 0);
	}

      if (flag_mpn_dc_divappr_q)
	{
	  MPN_COPY (rp, np, nn);
	  MPN_ZERO (qp, nn - dn);
	  qp[nn - dn] = mpn_dc_divappr_q (qp, rp, nn, dp, dn);
	  ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
	  ASSERT_ALWAYS (rp[-1] == rran0);
	  check_one (qrefp, rrefp, qp, rrefp, np, nn, dp, dn, "mpn_dc_divappr_q", 1);
	}

      if (flag_mpn_dc_div_q)
	{
	  MPN_COPY (rp, np, nn);
	  MPN_ZERO (qp, nn - dn);
	  qp[nn - dn] = mpn_dc_div_q (qp, rp, nn, dp, dn);
	  ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
	  ASSERT_ALWAYS (rp[-1] == rran0);
	  check_one (qrefp, rrefp, qp, rrefp, np, nn, dp, dn, "mpn_dc_div_q", 0);
	}

      if (nn - dn <= 5 || mpn_cmp (np + nn - dn, dp, dn) >= 0)
	continue;

      if (flag_mpn_mu_div_qr)
	{
	  mp_size_t itch;
	  mp_ptr scratch;
	  mp_limb_t ran;
	  itch = mpn_mu_div_qr_itch (nn, dn, 0);
	  scratch = __GMP_ALLOCATE_FUNC_LIMBS (itch + 1);
	  mpn_random (scratch + itch, 1);  ran = scratch[itch];
	  MPN_ZERO (qp, nn - dn);
	  MPN_ZERO (rp, dn);
	  rp[dn] = rran1;
	  qp[nn - dn] = mpn_mu_div_qr (qp, rp, np, nn, dp, dn, scratch);
	  ASSERT_ALWAYS (ran == scratch[itch]);
	  ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
	  ASSERT_ALWAYS (rp[-1] == rran0);  ASSERT_ALWAYS (rp[dn] == rran1);
	  check_one (qrefp, rrefp, qp, rp, np, nn, dp, dn, "mpn_mu_div_qr", 0);
	  __GMP_FREE_FUNC_LIMBS (scratch, itch);
	}

      if (flag_mpn_mu_divappr_q)
	{
	  mp_size_t itch;
	  mp_ptr scratch;
	  mp_limb_t ran;
	  itch = mpn_mu_divappr_q_itch (nn, dn, 0);
	  scratch = __GMP_ALLOCATE_FUNC_LIMBS (itch + 1);
	  mpn_random (scratch + itch, 1);  ran = scratch[itch];
	  MPN_ZERO (qp, nn - dn);
	  qp[nn - dn] = mpn_mu_divappr_q (qp, np, nn, dp, dn, scratch);
	  ASSERT_ALWAYS (ran == scratch[itch]);
	  ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
	  check_one (qrefp, rrefp, qp, rrefp, np, nn, dp, dn, "mpn_mu_divappr_q", 4);
	  __GMP_FREE_FUNC_LIMBS (scratch, itch);
	}

      if (flag_mpn_mu_div_q)
	{
	  mp_size_t itch;
	  mp_ptr scratch;
	  mp_limb_t ran;
	  itch = mpn_mu_divappr_q_itch (nn, dn, 0); /* FIXME: wrong itch function */
	  scratch = __GMP_ALLOCATE_FUNC_LIMBS (itch + 1);
	  mpn_random (scratch + itch, 1);  ran = scratch[itch];
	  MPN_ZERO (qp, nn - dn);
	  qp[nn - dn] = mpn_mu_div_q (qp, np, nn, dp, dn, scratch);
	  ASSERT_ALWAYS (ran == scratch[itch]);
	  ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
	  check_one (qrefp, rrefp, qp, rrefp, np, nn, dp, dn, "mpn_mu_div_q", 0);
	  __GMP_FREE_FUNC_LIMBS (scratch, itch);
	}
    }

  TMP_FREE;
}
#endif

#ifdef TIMING

#include "cputime.h"

#define TIME(t,func)							\
  do { long int __t0, __times, __t, __tmp;				\
    __times = 1;							\
    for (;;)								\
      {									\
	__t0 = cputime ();						\
	for (__t = 0; __t < __times; __t++)				\
	  {func;}							\
	__tmp = cputime () - __t0;					\
	if (__tmp > 1000) break;					\
	__times <<= 1;							\
      }									\
    (t) = (double) __tmp / __times;					\
  } while (0)

void
printres (double t, int width)
{
  int lg;

  if (t <= 0)
    {
      printf ("%*s", width, "badm");
      return;
    }

  lg = floor (log10 (t));

  if (lg < 0)
    printf (" %*.*fns", width-3, -1 - lg, t * 1000);
  else if (lg < 3)		/* 1-999 */
    printf (" %*.*fµs", width-3, 2 - lg, t);
  else if (lg < 6)		/* 1000-999999 */
    printf (" %*.*fms", width-3, 5 - lg, t * 0.001);
  else if (lg < 9)			/* 1000000-oo */
    printf (" %*.*fs", width-2, 8 - lg, t * 0.000001);
  else
    {
      t = t - fmod (t, pow(10.0, lg - 9));
      printf (" %*.*fs", width-2, 0, t * 0.000001);
    }
}

int
main (int argc, char **argv)
{
  gmp_randstate_t rs;
  unsigned long nbits, dbits;
  mpz_t n, d, q;
  mp_size_t nn, dn;
  mp_ptr np, dp, qp, rp;
  double t, tc;
  mp_limb_t di_32, di_21, dip[2];
  mp_limb_t cy;
  mp_size_t itch;
  mp_ptr scratch;
  TMP_DECL;

  TMP_MARK;

#ifdef POWERD
  unsigned long i;
  for (i = 4000000000u; i != 0; i--)
    ;
#endif

  if (argc == 2)
    {
      dbits = sizearg (argv[1]);
      nbits = 2 * dbits;
    }
  else if (argc == 3)
    {
      nbits = sizearg (argv[1]);
      dbits = sizearg (argv[2]);
    }
  else
    {
      printf ("usage: %s nsize dsize\n", argv[0]);
      printf ("   or: %s size\n", argv[0]);
      exit (1);
    }

  if (nbits <= dbits)
    {
      printf ("the dividend needs to be larger than the divisor\n");
      exit (1);
    }

  if (dbits <= GMP_NUMB_BITS)
    {
      printf ("the divisor needs to be at least 2 limbs (%d bits)\n",
	      1 + GMP_NUMB_BITS);
      exit (1);
    }

  printf ("DC_DIV_QR_THRESHOLD=%u  ", DC_DIV_QR_THRESHOLD);
  printf ("DC_DIV_Q_THRESHOLD=%u  ", DC_DIV_Q_THRESHOLD);
  printf ("DC_DIVAPPR_Q_THRESHOLD=%u\n", DC_DIVAPPR_Q_THRESHOLD);
  printf ("INV_NEWTON_THRESHOLD=%u\n", INV_NEWTON_THRESHOLD);

  gmp_randinit_default (rs);

  mpz_init (n);
  mpz_init (d);
  mpz_init (q);

  nn = nbits / GMP_NUMB_BITS + 1;
  dn = dbits / GMP_NUMB_BITS + 1;

  qp = TMP_ALLOC_LIMBS (nn);
  rp = TMP_ALLOC_LIMBS (nn);

  mpz_urandomb (n, rs, nbits);

  do
    mpz_urandomb (d, rs, dbits);
  while (mpz_sgn (d) == 0);

  np = PTR (n);
  dp = PTR (d);
  nn = SIZ (n);
  dn = SIZ (d);

  dp[dn - 1] |= GMP_NUMB_HIGHBIT;

  /* For timing purposes, avoid the somewhat atypical situation where the
     quotient is nn-dn+1 limbs.  Also, the MU functions don't handle that.  */
  if (mpn_cmp (np + nn - dn, dp, dn) >= 0)
    mpn_sub_n (np + nn - dn, np + nn - dn, dp, dn);

  printf ("                    SB                      DC                      MU\n");
  printf ("%8ld  div_qr divappr_q div_q  div_qr divappr_q div_q  div_qr divappr_q div_q\n", (long) nn);
  printf ("%8ld", (long) dn);

  mp_limb_t xp[2];
  cy = mpn_add_1 (xp, dp + dn - 2, 2, 1);
  if (cy != 0)
    {
      dip[0] = dip[1] = 0;
      di_32 = 0;
    }
  else
    {
      mpn_invert (dip, xp, 2, NULL);
      di_32 = dip[1];
    }
  invert_limb (di_21, dp[dn - 1]);

  int cnt;
  mp_limb_t tp[dn + 2];
  mp_limb_t ip[dn + 2];

      count_leading_zeros (cnt, dp[dn - 1]);
      if (cnt != 0)
	{
	  mpn_lshift (tp + 1, dp, dn, cnt);
	  tp[0] = 0;
	  mpn_invert (ip, tp, dn + 1, NULL);
	  ip[dn + 1] = 1;
	  mpn_lshift (ip, ip, dn + 2, cnt);
	}
      else
	{
	  MPN_COPY (tp + 1, dp, dn);
	  tp[0] = 0;
	  mpn_invert (ip, tp, dn + 1, NULL);
	  ip[dn + 1] = 1;
	}


  TIME (tc, MPN_COPY (rp, np, nn););

  /* SB functions */
  if ((double) (nn - dn) * dn < 1e9)
    {
      TIME (t, MPN_COPY (rp, np, nn); mpn_sb_div_qr (qp, rp, nn, dp, dn, dip));
      printres ((t - tc) * 1000, 8);
      fflush (stdout);
    }
  else printf ("%8s", "");

  if ((double) (nn - dn) * MIN(dn,nn-dn) < 1e9)
    {
      TIME (t, MPN_COPY (rp, np, nn); mpn_sb_divappr_q (qp, rp, nn, dp, dn, dip));
      printres ((t - tc) * 1000, 8);
      fflush (stdout);

      TIME (t, MPN_COPY (rp, np, nn); mpn_sb_div_q (qp, rp, nn, dp, dn, dip));
      printres ((t - tc) * 1000, 8);
      fflush (stdout);
    }
  else printf ("%16s", "");

  /* DC functions */
  TIME (t, MPN_COPY (rp, np, nn); mpn_dc_div_qr (qp, rp, nn, dp, dn));
  printres ((t - tc) * 1000, 8);
  fflush (stdout);

  TIME (t, MPN_COPY (rp, np, nn); mpn_dc_divappr_q (qp, rp, nn, dp, dn));
  printres ((t - tc) * 1000, 8);
  fflush (stdout);

  TIME (t, MPN_COPY (rp, np, nn); mpn_dc_div_q (qp, rp, nn, dp, dn));
  printres ((t - tc) * 1000, 8);
  fflush (stdout);

  /* MU functions */
  itch = mpn_mu_div_qr_itch (nn, dn, 0);
  scratch = __GMP_ALLOCATE_FUNC_LIMBS (itch);
  TIME (t, mpn_mu_div_qr (qp, rp, np, nn, dp, dn, scratch));
  printres (t * 1000, 8);
  fflush (stdout);
  __GMP_FREE_FUNC_LIMBS (scratch, itch);

  itch = mpn_mu_divappr_q_itch (nn, dn, 0);
  scratch = __GMP_ALLOCATE_FUNC_LIMBS (itch);
  TIME (t, mpn_mu_divappr_q (qp, np, nn, dp, dn, scratch));
  printres (t * 1000, 8);
  fflush (stdout);
  __GMP_FREE_FUNC_LIMBS (scratch, itch);

  itch = mpn_mu_divappr_q_itch (nn, dn, 0); /* FIXME: wrong itch function */
  scratch = __GMP_ALLOCATE_FUNC_LIMBS (itch);
  TIME (t, mpn_mu_div_q (qp, np, nn, dp, dn, scratch));
  printres (t * 1000, 8);
  fflush (stdout);
  __GMP_FREE_FUNC_LIMBS (scratch, itch);

  puts ("");

  printf ("              Legacy and Experimental Routines\n");
  printf ("%8ld    div_n_qr_pi1_21  div_n_qr_pi1_32  mpn_sb_pi1_div_qr  mpn_sb_divrem_mn  mpn_tdiv_qr\n", (long) nn);
  printf ("%8ld", (long) dn);

  if ((double) (nn - dn) * dn < 1e9)
    {
      TIME (t, MPN_COPY (rp, np, nn); mpn_div_n_qr_pi1_21 (qp, rp, nn, dp, dn, di_21));
      printres ((t - tc) * 1000, 14);
      fflush (stdout);
    }
  else printf ("%14s", "");
  if ((double) (nn - dn) * dn < 1e9)
    {
      TIME (t, MPN_COPY (rp, np, nn); mpn_div_n_qr_pi1_32 (qp, rp, nn, dp, dn, di_32));
      printres ((t - tc) * 1000, 17);
      fflush (stdout);
    }
  else printf ("%17s", "");
  if ((double) (nn - dn) * dn < 1e9)
    {
#if 0
      TIME (t, MPN_COPY (rp, np, nn); mpn_sb_pi1_div_qr (qp, rp, nn, dp, dn, di));
      printres ((t - tc) * 1000, 18);
      fflush (stdout);
#else
      printf ("%18s", "");
#endif

      if (dn > 2)
	{
	  TIME (t, MPN_COPY (rp, np, nn); mpn_sb_divrem_mn (qp, rp, nn, dp, dn));
	  printres ((t - tc) * 1000, 18);
	  fflush (stdout);
	}
      else printf ("%18s", "");
    }
  else printf ("%36s", "");

  TIME (t, mpn_tdiv_qr (qp, rp, 0, np, nn, dp, dn));
  printres (t * 1000, 16);
  fflush (stdout);

  puts ("");

  printf ("%8ld    div_n_qr_pin\n", (long) nn);
  printf ("%8ld", (long) dn);

  TIME (t, MPN_COPY (rp, np, nn); mpn_div_n_qr_pin (qp, np, nn, dp, dn, ip));
      printres ((t - tc) * 1000, 18);
  fflush (stdout);

  puts ("");

  TMP_FREE;
  exit (0);
}
#endif
