/* 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 "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;
  for (i = n - 1; i >= 0; i--)
    {
      printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]);
      printf (" " + (i == 0));
    }
  puts ("");
}

int
main (int argc, char **argv)
{
  gmp_randstate_t rs;
  unsigned long maxnbits, maxdbits, nbits, dbits;
  mpz_t d;
  mp_size_t maxdn, dn;
  mp_ptr ip;
  mp_ptr scratch;
  int test, work = 0;
  mp_size_t itch;
  mp_size_t i;
  mp_limb_t ran, sum;
  TMP_DECL;

  TMP_MARK;

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

  gmp_randinit_default (rs);

  mpz_init (d);

  maxdn = maxdbits / GMP_NUMB_BITS + 1;

  ip = TMP_ALLOC_LIMBS (maxdn);

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

#if RAND_UNIFORM
      mpz_urandomb (d, rs, dbits);
#else
      mpz_rrandomb (d, rs, dbits);
#endif
      mpz_setbit (d, 0);		/* mpn_binvert requirement */

      dn = SIZ (d);

      work += dbits;
      if (work >= 3123451)
	{
	  printf ("\r%d", test);
#ifdef DEBUG
	  printf ("\n");
#endif
	  fflush (stdout);
	  work = 0;
	}

#ifdef DEBUG
      printf ("d="); mpn_dump (PTR(d), dn);
#endif

      itch = mpn_binvert_itch (dn);
      scratch = __GMP_ALLOCATE_FUNC_LIMBS (itch + 1);

      mpn_random (&ran, 1);
      scratch[itch] = ran;

      mpn_binvert (ip, PTR(d), dn, scratch);

      if (scratch[itch] != ran)
	{
	  printf ("clobbered end of scratch for dn=%ld itch=%ld\n", dn, itch);
	}

      mpn_mullow_n (scratch, ip, PTR(d), dn);

      sum = 0;
      for (i = 1; i < dn; i++)
	sum |= scratch[i];

      if (sum != 0 || scratch[0] != 1)
	{
	  printf ("**********************************************************\n");
	  printf ("mpn_binvert wrong in %d\n", test);
	  printf ("d=    "); dumpy (PTR(d), dn);
	  printf ("ip=   "); dumpy (ip, dn);
	  printf ("**********************************************************\n");
	  abort ();
	}

      __GMP_FREE_FUNC_LIMBS (scratch, itch + 1);
    }

  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 > 250) break;						\
	__times <<= 1;							\
      }									\
    (t) = (double) __tmp / __times;					\
  } while (0)

void
printres (char *fmt, double t)
{
  int prc;
  t = 1000.0 * t;
  prc = 2 - floor (log10 (t));
  if (prc < 0)
    {
      t = t - fmod (t, pow (10.0, -prc));
      prc = 0;
    }
  printf (fmt, prc, t);
}

int
main (int argc, char **argv)
{
  gmp_randstate_t rs;
  unsigned long nbits, dbits;
  mpz_t d;
  mp_size_t dn;
  mp_ptr ip;
  mp_ptr scratch;
  double t;
  mp_size_t itch;
  char *progname = argv[0];
  int nh;
  TMP_DECL;

  TMP_MARK;

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

  nh = argc >= 2 && strcmp (argv[1], "-nh") == 0;
  if (nh)
    {
      argc--;
      argv++;
    }

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

  gmp_randinit_default (rs);

  mpz_init (d);

  dn = dbits / GMP_NUMB_BITS + 1;
  ip = TMP_ALLOC_LIMBS (dn);

  mpz_urandomb (d, rs, dbits);
  mpz_setbit (d, 0);			/* mpn_binvert requirement */

  dn = SIZ (d);

  itch = mpn_binvert_itch (dn);
  scratch = TMP_ALLOC_LIMBS (itch);

  if (! nh)
    printf ("        n    mpn_binvert\n");

  printf ("%9ld", dn);

  TIME (t, mpn_binvert (ip, PTR(d), dn, scratch));
  printres ("%12.*fµs\n", t);

  TMP_FREE;
  exit (0);
}
#endif