/*
---- calculation of PI(= 3.14159...) using FFT and AGM ----
    by T.Ooura, Nov. 1999.

Example compilation:
    GNU      : gcc -O -funroll-loops -fomit-frame-pointer pi_fftcs.c fftsg_h.c -lm -o pi_css5
    SUN      : cc -fast pi_fftcs.c fftsg_h.c -lm -o pi_css5
    HP       : aCC -fast pi_fftcs.c fftsg_h.c -lm -o pi_css5
    Microsoft: cl -O2 pi_fftcs.c fftsg_h.c -o pi_css5
    ...
    etc.
*/

#define PI_FFTC_VER "ver. LG1.1.2-MP1.5.2a.memsave"

/* Please check the following macros before compiling */
#ifndef DBL_ERROR_MARGIN
#define DBL_ERROR_MARGIN 0.4	/* must be < 0.5 */
#endif

#define DGTINT short int	/* sizeof(DGTINT) == 2 */
#define DGTINT_MAX SHRT_MAX

#define DGT_PACK 10
#define DGT_PACK_LINE 5
#define DGT_LINE_BLOCK 20

#include <math.h>
#include <limits.h>
#include <float.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <ctype.h>

void mp_load_0 (int n, int radix, int out[]);
void mp_load_1 (int n, int radix, int out[]);
void mp_round (int n, int radix, int m, int inout[]);
int mp_cmp (int n, int radix, int in1[], int in2[]);
void mp_add (int n, int radix, int in1[], int in2[], int out[]);
void mp_sub (int n, int radix, int in1[], int in2[], int out[]);
void mp_imul (int n, int radix, int in1[], int in2, int out[]);
int mp_idiv (int n, int radix, int in1[], int in2, int out[]);
void mp_idiv_2 (int n, int radix, int in[], int out[]);
double mp_mul_radix_test (int n, int radix, int nfft, double tmpfft[]);
void mp_mul (int n, int radix, int in1[], int in2[], int out[],
	     int tmp[], int nfft, double tmp1fft[], double tmp2fft[],
	     double tmp3fft[]);
void mp_squ (int n, int radix, int in[], int out[], int tmp[],
	     int nfft, double tmp1fft[], double tmp2fft[]);
void mp_mulhf (int n, int radix, int in1[], int in2[], int out[],
	       int tmp[], int nfft, double in1fft[], double tmpfft[]);
void mp_mulhf_use_in1fft (int n, int radix, double in1fft[], int in2[],
			  int out[], int tmp[], int nfft, double tmpfft[]);
void mp_squhf_use_infft (int n, int radix, double infft[], int in[],
			 int out[], int tmp[], int nfft, double tmpfft[]);
void mp_mulh (int n, int radix, int in1[], int in2[], int out[],
	      int nfft, double in1fft[], double outfft[]);
void mp_squh (int n, int radix, int in[], int out[],
	      int nfft, double outfft[]);
int mp_inv (int n, int radix, int in[], int out[],
	    int tmp1[], int tmp2[], int nfft,
	    double tmp1fft[], double tmp2fft[]);
int mp_sqrt (int n, int radix, int in[], int out[],
	     int tmp1[], int tmp2[], int nfft,
	     double tmp1fft[], double tmp2fft[]);
int mp_invisqrt (int n, int radix, int in, int out[],
		 int tmp1[], int tmp2[], int nfft,
		 double tmp1fft[], double tmp2fft[]);
void mp_sprintf (int n, int log10_radix, int in[], char out[]);
void mp_sscanf (int n, int log10_radix, char in[], int out[]);


int
main (int argc, char *argv[])
{
  int nfft, log2_nfft, radix, log10_radix, n, j = 0, k = 0, l = 0, npow, nprc;
  double err;
  int *a, *b, *c, *e, *i1, *i2;
  double *d1, *d2, *d3;
  char *dgt, filename[100];
  clock_t start_time;
  double elap_time, loop_time;
  FILE *f_out;
#ifndef QUIET_OUT
  fprintf (stdout, "Calculation of PI using FFT and AGM, %s\n", PI_FFTC_VER);
#endif

  /* if not run with the proper parameters from the command line */
  if (argc != 2)
    {
      printf ("\nUsage: %s digits\n", argv[0]);
      printf ("\nNumber of digits of pi to calculate?\n");
      scanf ("%d", &nfft);
    }
  else
    nfft = atoi (argv[1]);

#ifndef QUIET_OUT
  fprintf (stdout, "initializing...\n");
#endif
  nfft /= 4;
  start_time = clock ();
  for (log2_nfft = 1; (1 << log2_nfft) < nfft; log2_nfft++);
  nfft = 1 << log2_nfft;
  n = nfft + 2;
  a = (int *) malloc (2 * sizeof (int) + n * sizeof (DGTINT));
  b = (int *) malloc (2 * sizeof (int) + n * sizeof (DGTINT));
  c = (int *) malloc (2 * sizeof (int) + n * sizeof (DGTINT));
  e = (int *) malloc (2 * sizeof (int) + n * sizeof (DGTINT));
  i1 = (int *) malloc (2 * sizeof (int) + n * sizeof (DGTINT));
  i2 = (int *) malloc (2 * sizeof (int) + n * sizeof (DGTINT));
  d1 = (double *) malloc ((nfft + 2) * sizeof (double));
  d2 = (double *) malloc ((nfft + 2) * sizeof (double));
  d3 = (double *) malloc ((nfft + 2) * sizeof (double));
  if (d3 == NULL)
    {
      printf ("Allocation Failure!\n");
      exit (1);
    }
  /* ---- radix test ---- */
  log10_radix = 1;
  radix = 10;
  err = mp_mul_radix_test (n, radix, nfft, d1);
  err += DBL_EPSILON * (n * radix * radix / 4);
  while (100 * err < DBL_ERROR_MARGIN && radix <= DGTINT_MAX / 20)
    {
      err *= 100;
      log10_radix++;
      radix *= 10;
    }
#ifndef QUIET_OUT
  fprintf (stdout, "nfft= %d\nradix= %d\nerror_margin= %g\n", nfft, radix,
	   err);
  fprintf (stdout, "calculating %d digits of PI...\n", log10_radix * (n - 2));
#endif
  /*
   * ---- a formula based on the AGM (Arithmetic-Geometric Mean) ----
   *   c = sqrt(0.125);
   *   a = 1 + 3 * c;
   *   b = sqrt(a);
   *   e = b - 0.625;
   *   b = 2 * b;
   *   c = e - c;
   *   a = a + e;
   *   npow = 4;
   *   do {
   *       npow = 2 * npow;
   *       e = (a + b) / 2;
   *       b = sqrt(a * b);
   *       e = e - b;
   *       b = 2 * b;
   *       c = c - e;
   *       a = e + b;
   *   } while (e > SQRT_SQRT_EPSILON);
   *   e = e * e / 4;
   *   a = a + b;
   *   pi = (a * a - e - e / 2) / (a * c - e) / npow;
   * ---- modification ----
   *   This is a modified version of Gauss-Legendre formula
   *   (by T.Ooura). It is faster than original version.
   * ---- reference ----
   *   1. E.Salamin, 
   *      Computation of PI Using Arithmetic-Geometric Mean, 
   *      Mathematics of Computation, Vol.30 1976.
   *   2. R.P.Brent, 
   *      Fast Multiple-Precision Evaluation of Elementary Functions, 
   *      J. ACM 23 1976.
   *   3. D.Takahasi, Y.Kanada, 
   *      Calculation of PI to 51.5 Billion Decimal Digits on 
   *      Distributed Memoriy Parallel Processors, 
   *      Transactions of Information Processing Society of Japan, 
   *      Vol.39 No.7 1998.
   *   4. T.Ooura, 
   *      Improvement of the PI Calculation Algorithm and 
   *      Implementation of Fast Multiple-Precision Computation, 
   *      Information Processing Society of Japan SIG Notes, 
   *      98-HPC-74, 1998.
   */
  /* ---- c = 1 / sqrt(8) ---- */
  mp_invisqrt (n, radix, 8, c, i1, i2, nfft, d1, d2);
  /* ---- a = 1 + 3 * c ---- */
  mp_imul (n, radix, c, 3, e);
  mp_sscanf (n, log10_radix, "1", a);
  mp_add (n, radix, a, e, a);
  /* ---- b = sqrt(a) ---- */
  mp_sqrt (n, radix, a, b, i1, i2, nfft, d1, d2);
  /* ---- e = b - 0.625 ---- */
  mp_sscanf (n, log10_radix, "0.625", e);
  mp_sub (n, radix, b, e, e);
  /* ---- b = 2 * b ---- */
  mp_add (n, radix, b, b, b);
  /* ---- c = e - c ---- */
  mp_sub (n, radix, e, c, c);
  /* ---- a = a + e ---- */
  mp_add (n, radix, a, e, a);
#ifndef QUIET_OUT
  fprintf (stdout, "AGM iteration\n");
#endif
  npow = 4;
  elap_time = ((double) clock () - (double) start_time) / CLOCKS_PER_SEC;

  do
    {
      clock_t start_loop_time = clock ();
      npow *= 2;
      /* ---- e = (a + b) / 2 ---- */
      mp_add (n, radix, a, b, e);
      mp_idiv_2 (n, radix, e, e);
      /* ---- b = sqrt(a * b) ---- */
      mp_mul (n, radix, a, b, a, i1, nfft, d1, d2, d3);
      mp_sqrt (n, radix, a, b, i1, i2, nfft, d1, d2);
      /* ---- e = e - b ---- */
      mp_sub (n, radix, e, b, e);
      /* ---- b = 2 * b ---- */
      mp_add (n, radix, b, b, b);
      /* ---- c = c - e ---- */
      mp_sub (n, radix, c, e, c);
      /* ---- a = e + b ---- */
      mp_add (n, radix, e, b, a);
      /* ---- convergence check ---- */
      nprc = -e[1];
      if (e[0] == 0)
	{
	  nprc = n;
	}
      loop_time =
	((double) clock () - (double) start_loop_time) / CLOCKS_PER_SEC;
      elap_time += loop_time;
#ifndef QUIET_OUT
      fprintf (stdout, "precision= %d: %0.2f sec\n", 4 * nprc * log10_radix,
	       loop_time);
#endif
    }
  while (4 * nprc <= n);
  start_time = clock ();
  /* ---- e = e * e / 4 (half precision) ---- */
  mp_idiv_2 (n, radix, e, e);
  mp_squh (n, radix, e, e, nfft, d1);
  /* ---- a = a + b ---- */
  mp_add (n, radix, a, b, a);
  /* ---- a = (a * a - e - e / 2) / (a * c - e) / npow ---- */
  mp_mulhf (n, radix, a, c, c, i1, nfft, d1, d2);
  mp_sub (n, radix, c, e, c);
  mp_inv (n, radix, c, b, i1, i2, nfft, d2, d3);
  mp_squhf_use_infft (n, radix, d1, a, a, i1, nfft, d2);
  mp_sub (n, radix, a, e, a);
  mp_idiv_2 (n, radix, e, e);
  mp_sub (n, radix, a, e, a);
  mp_mul (n, radix, a, b, a, i1, nfft, d1, d2, d3);
  mp_idiv (n, radix, a, npow, a);
  /* ---- output ---- */
  dgt = (char *) d1;
  mp_sprintf (n - 1, log10_radix, a, dgt);
  elap_time += ((double) clock () - (double) start_time) / CLOCKS_PER_SEC;

  sprintf (filename, "pi%i.txt", log10_radix * (n - 2));

  f_out = fopen (filename, "w");
#ifndef QUIET_OUT
  fprintf (stdout, "writing %s...\n", filename);
#endif
  do
    {
      if (!isdigit (*dgt))
	{
	  if (isalpha (*dgt) != 0)
	    {
	      fputc ('\n', f_out);
	      fputc ('\n', f_out);
	    }
	  fputc (*dgt, f_out);
	  fputc ('\n', f_out);
	  fputc ('\n', f_out);
	  j = 0;
	  k = 0;
	  l = 0;
	  continue;
	}
      fputc (*dgt, f_out);
      if (++j >= DGT_PACK)
	{
	  j = 0;
	  if (++k >= DGT_PACK_LINE)
	    {
	      k = 0;
	      fputc ('\n', f_out);
	      if (++l >= DGT_LINE_BLOCK)
		{
		  l = 0;
		  fputc ('\n', f_out);
		}
	    }
	  else
	    {
	      fputc (' ', f_out);
	    }
	}
    }
  while (*dgt++ && *dgt!='e');
  fputc ('\n', f_out);
  fprintf (f_out, "%s\n", dgt);
  fclose (f_out);
  free (d3);
  free (d2);
  free (d1);
  free (i2);
  free (i1);
  free (e);
  free (c);
  free (b);
  free (a);
  /* ---- difftime ---- */
  fprintf (stdout, "%0.2f sec. (real time)\n", elap_time);
  
  /* don't quite before allowing the user to see all output */
  if(argc!=2)
  {
  	fgets(filename,99,stdin);
	fprintf(stdout,"Hit RETURN to exit.\n");
	fgets(filename,99,stdin);
  }

  return (0);
}


/* -------- multiple precision routines -------- */


#include <math.h>
#include <float.h>
#include <stdio.h>

/* ---- floating point format ----
    data := data[0] * pow(radix, data[1]) * 
            (data[2] + data[3]/radix + data[4]/radix/radix + ...), 
    data[0]       : sign (1;data>0, -1;data<0, 0;data==0)
    data[1]       : exponent (0;data==0)
    data[2...n+1] : digits
   ---- function prototypes ----
    void mp_load_0(int n, int radix, int out[]);
    void mp_load_1(int n, int radix, int out[]);
    void mp_round(int n, int radix, int m, int inout[]);
    int mp_cmp(int n, int radix, int in1[], int in2[]);
    void mp_add(int n, int radix, int in1[], int in2[], int out[]);
    void mp_sub(int n, int radix, int in1[], int in2[], int out[]);
    void mp_imul(int n, int radix, int in1[], int in2, int out[]);
    int mp_idiv(int n, int radix, int in1[], int in2, int out[]);
    void mp_idiv_2(int n, int radix, int in[], int out[]);
    double mp_mul_radix_test(int n, int radix, int nfft, 
            double tmpfft[]);
    void mp_mul(int n, int radix, int in1[], int in2[], int out[], 
            int tmp[], int nfft, double tmp1fft[], double tmp2fft[], 
            double tmp3fft[]);
    void mp_squ(int n, int radix, int in[], int out[], int tmp[], 
            int nfft, double tmp1fft[], double tmp2fft[]);
    void mp_mulhf(int n, int radix, int in1[], int in2[], int out[], 
            int tmp[], int nfft, double in1fft[], double tmpfft[]);
    void mp_mulhf_use_in1fft(int n, int radix, double in1fft[], int in2[], 
            int out[], int tmp[], int nfft, double tmpfft[]);
    void mp_squhf_use_infft(int n, int radix, double infft[], int in[], 
            int out[], int tmp[], int nfft, double tmpfft[]);
    void mp_mulh(int n, int radix, int in1[], int in2[], int out[], 
            int nfft, double in1fft[], double outfft[]);
    void mp_squh(int n, int radix, int in[], int out[], 
            int nfft, double outfft[]);
    int mp_inv(int n, int radix, int in[], int out[], 
            int tmp1[], int tmp2[], int nfft, 
            double tmp1fft[], double tmp2fft[]);
    int mp_sqrt(int n, int radix, int in[], int out[], 
            int tmp1[], int tmp2[], int nfft, 
            double tmp1fft[], double tmp2fft[]);
    int mp_invisqrt(int n, int radix, int in, int out[], 
            int tmp1[], int tmp2[], int nfft, 
            double tmp1fft[], double tmp2fft[]);
    void mp_sprintf(int n, int log10_radix, int in[], char out[]);
    void mp_sscanf(int n, int log10_radix, char in[], int out[]);
   ----
*/


/* -------- mp_load routines -------- */


void
mp_load_0 (int n, int radix, int out[])
{
  int j;
  DGTINT *outr;

  outr = ((DGTINT *) & out[2]) - 2;
  out[0] = 0;
  out[1] = 0;
  for (j = 2; j <= n + 1; j++)
    {
      outr[j] = 0;
    }
}


void
mp_load_1 (int n, int radix, int out[])
{
  int j;
  DGTINT *outr;

  outr = ((DGTINT *) & out[2]) - 2;
  out[0] = 1;
  out[1] = 0;
  outr[2] = 1;
  for (j = 3; j <= n + 1; j++)
    {
      outr[j] = 0;
    }
}


void
mp_round (int n, int radix, int m, int inout[])
{
  int j, x;
  DGTINT *inoutr;

  inoutr = ((DGTINT *) & inout[2]) - 2;
  if (m < n)
    {
      for (j = n + 1; j > m + 2; j--)
	{
	  inoutr[j] = 0;
	}
      x = 2 * inoutr[m + 2];
      inoutr[m + 2] = 0;
      if (x >= radix)
	{
	  for (j = m + 1; j >= 2; j--)
	    {
	      x = inoutr[j] + 1;
	      if (x < radix)
		{
		  inoutr[j] = (DGTINT) x;
		  break;
		}
	      inoutr[j] = 0;
	    }
	  if (x >= radix)
	    {
	      inoutr[2] = 1;
	      inout[1]++;
	    }
	}
    }
}


/* -------- mp_add routines -------- */


int
mp_cmp (int n, int radix, int in1[], int in2[])
{
  int mp_unsgn_cmp (int n, int in1[], int in2[]);

  if (in1[0] > in2[0])
    {
      return 1;
    }
  else if (in1[0] < in2[0])
    {
      return -1;
    }
  return in1[0] * mp_unsgn_cmp (n, &in1[1], &in2[1]);
}


void
mp_add (int n, int radix, int in1[], int in2[], int out[])
{
  int mp_unsgn_cmp (int n, int in1[], int in2[]);
  int mp_unexp_add (int n, int radix, int expdif,
		    DGTINT in1[], DGTINT in2[], DGTINT out[]);
  int mp_unexp_sub (int n, int radix, int expdif,
		    DGTINT in1[], DGTINT in2[], DGTINT out[]);
  int outsgn, outexp, expdif;

  expdif = in1[1] - in2[1];
  outexp = in1[1];
  if (expdif < 0)
    {
      outexp = in2[1];
    }
  outsgn = in1[0] * in2[0];
  if (outsgn >= 0)
    {
      if (outsgn > 0)
	{
	  outsgn = in1[0];
	}
      else
	{
	  outsgn = in1[0] + in2[0];
	  outexp = in1[1] + in2[1];
	  expdif = 0;
	}
      if (expdif >= 0)
	{
	  outexp += mp_unexp_add (n, radix, expdif,
				  (DGTINT *) & in1[2], (DGTINT *) & in2[2],
				  (DGTINT *) & out[2]);
	}
      else
	{
	  outexp += mp_unexp_add (n, radix, -expdif,
				  (DGTINT *) & in2[2], (DGTINT *) & in1[2],
				  (DGTINT *) & out[2]);
	}
    }
  else
    {
      outsgn = mp_unsgn_cmp (n, &in1[1], &in2[1]);
      if (outsgn >= 0)
	{
	  expdif = mp_unexp_sub (n, radix, expdif,
				 (DGTINT *) & in1[2], (DGTINT *) & in2[2],
				 (DGTINT *) & out[2]);
	}
      else
	{
	  expdif = mp_unexp_sub (n, radix, -expdif,
				 (DGTINT *) & in2[2], (DGTINT *) & in1[2],
				 (DGTINT *) & out[2]);
	}
      outexp -= expdif;
      outsgn *= in1[0];
      if (expdif == n)
	{
	  outsgn = 0;
	}
    }
  if (outsgn == 0)
    {
      outexp = 0;
    }
  out[0] = outsgn;
  out[1] = outexp;
}


void
mp_sub (int n, int radix, int in1[], int in2[], int out[])
{
  int mp_unsgn_cmp (int n, int in1[], int in2[]);
  int mp_unexp_add (int n, int radix, int expdif,
		    DGTINT in1[], DGTINT in2[], DGTINT out[]);
  int mp_unexp_sub (int n, int radix, int expdif,
		    DGTINT in1[], DGTINT in2[], DGTINT out[]);
  int outsgn, outexp, expdif;

  expdif = in1[1] - in2[1];
  outexp = in1[1];
  if (expdif < 0)
    {
      outexp = in2[1];
    }
  outsgn = in1[0] * in2[0];
  if (outsgn <= 0)
    {
      if (outsgn < 0)
	{
	  outsgn = in1[0];
	}
      else
	{
	  outsgn = in1[0] - in2[0];
	  outexp = in1[1] + in2[1];
	  expdif = 0;
	}
      if (expdif >= 0)
	{
	  outexp += mp_unexp_add (n, radix, expdif,
				  (DGTINT *) & in1[2], (DGTINT *) & in2[2],
				  (DGTINT *) & out[2]);
	}
      else
	{
	  outexp += mp_unexp_add (n, radix, -expdif,
				  (DGTINT *) & in2[2], (DGTINT *) & in1[2],
				  (DGTINT *) & out[2]);
	}
    }
  else
    {
      outsgn = mp_unsgn_cmp (n, &in1[1], &in2[1]);
      if (outsgn >= 0)
	{
	  expdif = mp_unexp_sub (n, radix, expdif,
				 (DGTINT *) & in1[2], (DGTINT *) & in2[2],
				 (DGTINT *) & out[2]);
	}
      else
	{
	  expdif = mp_unexp_sub (n, radix, -expdif,
				 (DGTINT *) & in2[2], (DGTINT *) & in1[2],
				 (DGTINT *) & out[2]);
	}
      outexp -= expdif;
      outsgn *= in1[0];
      if (expdif == n)
	{
	  outsgn = 0;
	}
    }
  if (outsgn == 0)
    {
      outexp = 0;
    }
  out[0] = outsgn;
  out[1] = outexp;
}


/* -------- mp_add child routines -------- */


int
mp_unsgn_cmp (int n, int in1[], int in2[])
{
  int j, cmp;
  DGTINT *in1r, *in2r;

  in1r = ((DGTINT *) & in1[1]) - 1;
  in2r = ((DGTINT *) & in2[1]) - 1;
  cmp = in1[0] - in2[0];
  for (j = 1; j <= n && cmp == 0; j++)
    {
      cmp = in1r[j] - in2r[j];
    }
  if (cmp > 0)
    {
      cmp = 1;
    }
  else if (cmp < 0)
    {
      cmp = -1;
    }
  return cmp;
}


int
mp_unexp_add (int n, int radix, int expdif,
	      DGTINT in1[], DGTINT in2[], DGTINT out[])
{
  int j, x, carry;

  carry = 0;
  if (expdif == 0 && in1[0] + in2[0] >= radix)
    {
      x = in1[n - 1] + in2[n - 1];
      carry = x >= radix ? -1 : 0;
      for (j = n - 1; j > 0; j--)
	{
	  x = in1[j - 1] + in2[j - 1] - carry;
	  carry = x >= radix ? -1 : 0;
	  out[j] = (DGTINT) (x - (radix & carry));
	}
      out[0] = (DGTINT) - carry;
    }
  else
    {
      if (expdif > n)
	{
	  expdif = n;
	}
      for (j = n - 1; j >= expdif; j--)
	{
	  x = in1[j] + in2[j - expdif] - carry;
	  carry = x >= radix ? -1 : 0;
	  out[j] = (DGTINT) (x - (radix & carry));
	}
      for (j = expdif - 1; j >= 0; j--)
	{
	  x = in1[j] - carry;
	  carry = x >= radix ? -1 : 0;
	  out[j] = (DGTINT) (x - (radix & carry));
	}
      if (carry != 0)
	{
	  for (j = n - 1; j > 0; j--)
	    {
	      out[j] = out[j - 1];
	    }
	  out[0] = (DGTINT) - carry;
	}
    }
  return -carry;
}


int
mp_unexp_sub (int n, int radix, int expdif,
	      DGTINT in1[], DGTINT in2[], DGTINT out[])
{
  int j, x, borrow, ncancel;

  if (expdif > n)
    {
      expdif = n;
    }
  borrow = 0;
  for (j = n - 1; j >= expdif; j--)
    {
      x = in1[j] - in2[j - expdif] + borrow;
      borrow = x < 0 ? -1 : 0;
      out[j] = (DGTINT) (x + (radix & borrow));
    }
  for (j = expdif - 1; j >= 0; j--)
    {
      x = in1[j] + borrow;
      borrow = x < 0 ? -1 : 0;
      out[j] = (DGTINT) (x + (radix & borrow));
    }
  ncancel = 0;
  for (j = 0; j < n && out[j] == 0; j++)
    {
      ncancel = j + 1;
    }
  if (ncancel > 0 && ncancel < n)
    {
      for (j = 0; j < n - ncancel; j++)
	{
	  out[j] = out[j + ncancel];
	}
      for (j = n - ncancel; j < n; j++)
	{
	  out[j] = 0;
	}
    }
  return ncancel;
}


/* -------- mp_imul routines -------- */


void
mp_imul (int n, int radix, int in1[], int in2, int out[])
{
  void mp_unsgn_imul (int n, double dradix, int in1[], double din2,
		      int out[]);

  if (in2 > 0)
    {
      out[0] = in1[0];
    }
  else if (in2 < 0)
    {
      out[0] = -in1[0];
      in2 = -in2;
    }
  else
    {
      out[0] = 0;
    }
  mp_unsgn_imul (n, radix, &in1[1], in2, &out[1]);
  if (out[0] == 0)
    {
      out[1] = 0;
    }
}


int
mp_idiv (int n, int radix, int in1[], int in2, int out[])
{
  void mp_load_0 (int n, int radix, int out[]);
  void mp_unsgn_idiv (int n, double dradix, int in1[], double din2,
		      int out[]);

  if (in2 == 0)
    {
      return -1;
    }
  if (in2 > 0)
    {
      out[0] = in1[0];
    }
  else
    {
      out[0] = -in1[0];
      in2 = -in2;
    }
  if (in1[0] == 0)
    {
      mp_load_0 (n, radix, out);
      return 0;
    }
  mp_unsgn_idiv (n, radix, &in1[1], in2, &out[1]);
  return 0;
}


void
mp_idiv_2 (int n, int radix, int in[], int out[])
{
  int j, ix, carry, shift;
  DGTINT *inr, *outr;

  inr = ((DGTINT *) & in[2]) - 2;
  outr = ((DGTINT *) & out[2]) - 2;
  out[0] = in[0];
  shift = 0;
  if (inr[2] == 1)
    {
      shift = 1;
    }
  out[1] = in[1] - shift;
  carry = -shift;
  for (j = 2; j <= n + 1 - shift; j++)
    {
      ix = inr[j + shift] + (radix & carry);
      carry = -(ix & 1);
      outr[j] = (DGTINT) (ix >> 1);
    }
  if (shift > 0)
    {
      outr[n + 1] = (DGTINT) ((radix & carry) >> 1);
    }
}


/* -------- mp_imul child routines -------- */


void
mp_unsgn_imul (int n, double dradix, int in1[], double din2, int out[])
{
  int j, carry, shift;
  double x, d1_radix;
  DGTINT *in1r, *outr;

  in1r = ((DGTINT *) & in1[1]) - 1;
  outr = ((DGTINT *) & out[1]) - 1;
  d1_radix = 1.0 / dradix;
  carry = 0;
  for (j = n; j >= 1; j--)
    {
      x = din2 * in1r[j] + carry + 0.5;
      carry = (int) (d1_radix * x);
      outr[j] = (DGTINT) (x - dradix * carry);
    }
  shift = 0;
  x = carry + 0.5;
  while (x > 1)
    {
      x *= d1_radix;
      shift++;
    }
  out[0] = in1[0] + shift;
  if (shift > 0)
    {
      while (shift > n)
	{
	  carry = (int) (d1_radix * carry + 0.5);
	  shift--;
	}
      for (j = n; j >= shift + 1; j--)
	{
	  outr[j] = outr[j - shift];
	}
      for (j = shift; j >= 1; j--)
	{
	  x = carry + 0.5;
	  carry = (int) (d1_radix * x);
	  outr[j] = (DGTINT) (x - dradix * carry);
	}
    }
}


void
mp_unsgn_idiv (int n, double dradix, int in1[], double din2, int out[])
{
  int j, ix, carry, shift;
  double x, d1_in2;
  DGTINT *in1r, *outr;

  in1r = ((DGTINT *) & in1[1]) - 1;
  outr = ((DGTINT *) & out[1]) - 1;
  d1_in2 = 1.0 / din2;
  shift = 0;
  x = 0;
  do
    {
      shift++;
      x *= dradix;
      if (shift <= n)
	{
	  x += in1r[shift];
	}
    }
  while (x < din2 - 0.5);
  x += 0.5;
  ix = (int) (d1_in2 * x);
  carry = (int) (x - din2 * ix);
  outr[1] = (DGTINT) ix;
  shift--;
  out[0] = in1[0] - shift;
  if (shift >= n)
    {
      shift = n - 1;
    }
  for (j = 2; j <= n - shift; j++)
    {
      x = in1r[j + shift] + dradix * carry + 0.5;
      ix = (int) (d1_in2 * x);
      carry = (int) (x - din2 * ix);
      outr[j] = (DGTINT) ix;
    }
  for (j = n - shift + 1; j <= n; j++)
    {
      x = dradix * carry + 0.5;
      ix = (int) (d1_in2 * x);
      carry = (int) (x - din2 * ix);
      outr[j] = (DGTINT) ix;
    }
}


/* -------- mp_mul routines -------- */


double
mp_mul_radix_test (int n, int radix, int nfft, double tmpfft[])
{
  void mp_mul_csqu (int nfft, double d1[]);
  double mp_mul_d2i_test (int radix, int nfft, double din[]);
  int j, ndata, radix_2;

  ndata = (nfft >> 1) + 1;
  if (ndata > n)
    {
      ndata = n;
    }
  tmpfft[nfft + 1] = radix - 1;
  for (j = nfft; j > ndata; j--)
    {
      tmpfft[j] = 0;
    }
  radix_2 = (radix + 1) / 2;
  for (j = ndata; j > 2; j--)
    {
      tmpfft[j] = radix_2;
    }
  tmpfft[2] = radix;
  tmpfft[1] = radix - 1;
  tmpfft[0] = 0;
  mp_mul_csqu (nfft, tmpfft);
  return 2 * mp_mul_d2i_test (radix, nfft, tmpfft);
}


void
mp_mul (int n, int radix, int in1[], int in2[], int out[],
	int tmp[], int nfft, double tmp1fft[], double tmp2fft[],
	double tmp3fft[])
{
  void mp_add (int n, int radix, int in1[], int in2[], int out[]);
  void mp_mul_i2d (int n, int radix, int nfft, int shift,
		   int in[], double dout[]);
  void mp_mul_cmul_nt_out (int nfft, double d1[], double d2[]);
  void mp_mul_cmul_nt_d2 (int nfft, double d1[], double d2[]);
  void mp_mul_cmul_nt_d1_add (int nfft, double d1[], double d2[],
			      double d3[]);
  void mp_mul_d2i (int n, int radix, int nfft, double din[], int out[]);
  int n_h, shift;
  DGTINT *in1r, *in2r;

  in1r = ((DGTINT *) & in1[2]) - 2;
  in2r = ((DGTINT *) & in2[2]) - 2;
  shift = (nfft >> 1) + 1;
  while (n > shift)
    {
      if (in1r[shift + 2] + in2r[shift + 2] != 0)
	{
	  break;
	}
      shift++;
    }
  n_h = n / 2 + 1;
  if (n_h < n - shift)
    {
      n_h = n - shift;
    }
  /* ---- tmp3fft = (upper) in1 * (lower) in2 ---- */
  mp_mul_i2d (n, radix, nfft, 0, in1, tmp1fft);
  mp_mul_i2d (n, radix, nfft, shift, in2, tmp3fft);
  mp_mul_cmul_nt_out (nfft, tmp1fft, tmp3fft);
  /* ---- tmp = (upper) in1 * (upper) in2 ---- */
  mp_mul_i2d (n, radix, nfft, 0, in2, tmp2fft);
  mp_mul_cmul_nt_d2 (nfft, tmp2fft, tmp1fft);
  mp_mul_d2i (n, radix, nfft, tmp1fft, tmp);
  /* ---- tmp3fft += (upper) in2 * (lower) in1 ---- */
  mp_mul_i2d (n, radix, nfft, shift, in1, tmp1fft);
  mp_mul_cmul_nt_d1_add (nfft, tmp2fft, tmp1fft, tmp3fft);
  /* ---- out = tmp + tmp3fft ---- */
  mp_mul_d2i (n_h, radix, nfft, tmp3fft, out);
  mp_add (n, radix, out, tmp, out);
}


void
mp_squ (int n, int radix, int in[], int out[], int tmp[],
	int nfft, double tmp1fft[], double tmp2fft[])
{
  void mp_add (int n, int radix, int in1[], int in2[], int out[]);
  void mp_mul_i2d (int n, int radix, int nfft, int shift,
		   int in[], double dout[]);
  void mp_mul_cmul (int nfft, double d1[], double d2[]);
  void mp_mul_csqu_nt_d1 (int nfft, double d1[]);
  void mp_mul_d2i (int n, int radix, int nfft, double din[], int out[]);
  int n_h, shift;
  DGTINT *inr;

  inr = ((DGTINT *) & in[2]) - 2;
  shift = (nfft >> 1) + 1;
  while (n > shift)
    {
      if (inr[shift + 2] != 0)
	{
	  break;
	}
      shift++;
    }
  n_h = n / 2 + 1;
  if (n_h < n - shift)
    {
      n_h = n - shift;
    }
  /* ---- tmp = 2 * (upper) in * (lower) in ---- */
  mp_mul_i2d (n, radix, nfft, 0, in, tmp1fft);
  mp_mul_i2d (n, radix, nfft, shift, in, tmp2fft);
  mp_mul_cmul (nfft, tmp1fft, tmp2fft);
  mp_mul_d2i (n_h, radix, nfft, tmp2fft, tmp);
  mp_add (n_h, radix, tmp, tmp, tmp);
  /* ---- out = tmp + ((upper) in)^2 ---- */
  mp_mul_csqu_nt_d1 (nfft, tmp1fft);
  mp_mul_d2i (n, radix, nfft, tmp1fft, out);
  mp_add (n, radix, out, tmp, out);
}


void
mp_mulhf (int n, int radix, int in1[], int in2[], int out[],
	  int tmp[], int nfft, double in1fft[], double tmpfft[])
{
  void mp_add (int n, int radix, int in1[], int in2[], int out[]);
  void mp_mul_i2d (int n, int radix, int nfft, int shift,
		   int in[], double dout[]);
  void mp_mul_cmul (int nfft, double d1[], double d2[]);
  void mp_mul_cmul_nt_d1 (int nfft, double d1[], double d2[]);
  void mp_mul_d2i (int n, int radix, int nfft, double din[], int out[]);
  int n_h, shift;
  DGTINT *in2r;

  in2r = ((DGTINT *) & in2[2]) - 2;
  shift = (nfft >> 1) + 1;
  while (n > shift)
    {
      if (in2r[shift + 2] != 0)
	{
	  break;
	}
      shift++;
    }
  n_h = n / 2 + 1;
  if (n_h < n - shift)
    {
      n_h = n - shift;
    }
  /* ---- tmp = (upper) in1 * (upper) in2 ---- */
  mp_mul_i2d (n, radix, nfft, 0, in1, in1fft);
  mp_mul_i2d (n, radix, nfft, 0, in2, tmpfft);
  mp_mul_cmul (nfft, in1fft, tmpfft);
  mp_mul_d2i (n, radix, nfft, tmpfft, tmp);
  /* ---- out = tmp + (upper) in1 * (lower) in2 ---- */
  mp_mul_i2d (n, radix, nfft, shift, in2, tmpfft);
  mp_mul_cmul_nt_d1 (nfft, in1fft, tmpfft);
  mp_mul_d2i (n_h, radix, nfft, tmpfft, out);
  mp_add (n, radix, out, tmp, out);
}


void
mp_mulhf_use_in1fft (int n, int radix, double in1fft[], int in2[],
		     int out[], int tmp[], int nfft, double tmpfft[])
{
  void mp_add (int n, int radix, int in1[], int in2[], int out[]);
  void mp_mul_i2d (int n, int radix, int nfft, int shift,
		   int in[], double dout[]);
  void mp_mul_cmul_nt_d1 (int nfft, double d1[], double d2[]);
  void mp_mul_d2i (int n, int radix, int nfft, double din[], int out[]);
  int n_h, shift;
  DGTINT *in2r;

  in2r = ((DGTINT *) & in2[2]) - 2;
  shift = (nfft >> 1) + 1;
  while (n > shift)
    {
      if (in2r[shift + 2] != 0)
	{
	  break;
	}
      shift++;
    }
  n_h = n / 2 + 1;
  if (n_h < n - shift)
    {
      n_h = n - shift;
    }
  /* ---- tmp = (upper) in1fft * (upper) in2 ---- */
  mp_mul_i2d (n, radix, nfft, 0, in2, tmpfft);
  mp_mul_cmul_nt_d1 (nfft, in1fft, tmpfft);
  mp_mul_d2i (n, radix, nfft, tmpfft, tmp);
  /* ---- out = tmp + (upper) in1 * (lower) in2 ---- */
  mp_mul_i2d (n, radix, nfft, shift, in2, tmpfft);
  mp_mul_cmul_nt_d1 (nfft, in1fft, tmpfft);
  mp_mul_d2i (n_h, radix, nfft, tmpfft, out);
  mp_add (n, radix, out, tmp, out);
}


void
mp_squhf_use_infft (int n, int radix, double infft[], int in[],
		    int out[], int tmp[], int nfft, double tmpfft[])
{
  void mp_add (int n, int radix, int in1[], int in2[], int out[]);
  void mp_mul_i2d (int n, int radix, int nfft, int shift,
		   int in[], double dout[]);
  void mp_mul_cmul_nt_d1 (int nfft, double d1[], double d2[]);
  void mp_mul_csqu_nt_d1 (int nfft, double d1[]);
  void mp_mul_d2i (int n, int radix, int nfft, double din[], int out[]);
  int n_h, shift;
  DGTINT *inr;

  inr = ((DGTINT *) & in[2]) - 2;
  shift = (nfft >> 1) + 1;
  while (n > shift)
    {
      if (inr[shift + 2] != 0)
	{
	  break;
	}
      shift++;
    }
  n_h = n / 2 + 1;
  if (n_h < n - shift)
    {
      n_h = n - shift;
    }
  /* ---- tmp = (upper) infft * (lower) in ---- */
  mp_mul_i2d (n, radix, nfft, shift, in, tmpfft);
  mp_mul_cmul_nt_d1 (nfft, infft, tmpfft);
  mp_mul_d2i (n_h, radix, nfft, tmpfft, tmp);
  /* ---- out = tmp + ((upper) infft)^2 ---- */
  mp_mul_csqu_nt_d1 (nfft, infft);
  mp_mul_d2i (n, radix, nfft, infft, out);
  mp_add (n, radix, out, tmp, out);
}


void
mp_mulh (int n, int radix, int in1[], int in2[], int out[],
	 int nfft, double in1fft[], double outfft[])
{
  void mp_mul_i2d (int n, int radix, int nfft, int shift,
		   int in[], double dout[]);
  void mp_mul_cmul (int nfft, double d1[], double d2[]);
  void mp_mul_d2i (int n, int radix, int nfft, double din[], int out[]);

  mp_mul_i2d (n, radix, nfft, 0, in1, in1fft);
  mp_mul_i2d (n, radix, nfft, 0, in2, outfft);
  mp_mul_cmul (nfft, in1fft, outfft);
  mp_mul_d2i (n, radix, nfft, outfft, out);
}


void
mp_mulh_use_in1fft (int n, int radix, double in1fft[],
		    int shift, int in2[], int out[], int nfft,
		    double outfft[])
{
  void mp_mul_i2d (int n, int radix, int nfft, int shift,
		   int in[], double dout[]);
  void mp_mul_cmul_nt_d1 (int nfft, double d1[], double d2[]);
  void mp_mul_d2i (int n, int radix, int nfft, double din[], int out[]);
  DGTINT *in2r;

  in2r = ((DGTINT *) & in2[2]) - 2;
  while (n > shift)
    {
      if (in2r[shift + 2] != 0)
	{
	  break;
	}
      shift++;
    }
  mp_mul_i2d (n, radix, nfft, shift, in2, outfft);
  mp_mul_cmul_nt_d1 (nfft, in1fft, outfft);
  mp_mul_d2i (n, radix, nfft, outfft, out);
}


void
mp_squh (int n, int radix, int in[], int out[], int nfft, double outfft[])
{
  void mp_mul_i2d (int n, int radix, int nfft, int shift,
		   int in[], double dout[]);
  void mp_mul_csqu (int nfft, double d1[]);
  void mp_mul_d2i (int n, int radix, int nfft, double din[], int out[]);

  mp_mul_i2d (n, radix, nfft, 0, in, outfft);
  mp_mul_csqu (nfft, outfft);
  mp_mul_d2i (n, radix, nfft, outfft, out);
}


void
mp_squh_save_infft (int n, int radix, int in[], int out[],
		    int nfft, double infft[], double outfft[])
{
  void mp_mul_i2d (int n, int radix, int nfft, int shift,
		   int in[], double dout[]);
  void mp_mul_csqu_save_d1 (int nfft, double d1[], double d2[]);
  void mp_mul_d2i (int n, int radix, int nfft, double din[], int out[]);

  mp_mul_i2d (n, radix, nfft, 0, in, infft);
  mp_mul_csqu_save_d1 (nfft, infft, outfft);
  mp_mul_d2i (n, radix, nfft, outfft, out);
}


void
mp_squh_use_in1fft (int n, int radix, double inoutfft[], int out[], int nfft)
{
  void mp_mul_csqu_nt_d1 (int nfft, double d1[]);
  void mp_mul_d2i (int n, int radix, int nfft, double din[], int out[]);

  mp_mul_csqu_nt_d1 (nfft, inoutfft);
  mp_mul_d2i (n, radix, nfft, inoutfft, out);
}


/* -------- mp_mul child routines -------- */


void
mp_mul_i2d (int n, int radix, int nfft, int shift, int in[], double dout[])
{
  int j, x, carry, ndata, radix_2, topdgt;
  DGTINT *inr;

  inr = ((DGTINT *) & in[2]) - 2;
  ndata = 0;
  topdgt = 0;
  if (n > shift)
    {
      topdgt = inr[shift + 2];
      ndata = (nfft >> 1) + 1;
      if (ndata > n - shift)
	{
	  ndata = n - shift;
	}
    }
  dout[nfft + 1] = in[0] * topdgt;
  for (j = nfft; j > ndata; j--)
    {
      dout[j] = 0;
    }
  /* ---- abs(dout[j]) <= radix/2 (to keep FFT precision) ---- */
  if (ndata > 1)
    {
      radix_2 = radix / 2;
      carry = 0;
      for (j = ndata + 1; j > 3; j--)
	{
	  x = inr[j + shift] - carry;
	  carry = x >= radix_2 ? -1 : 0;
	  dout[j - 1] = x - (radix & carry);
	}
      dout[2] = inr[shift + 3] - carry;
    }
  dout[1] = topdgt;
  dout[0] = in[1] - shift;
}


void
mp_mul_cmul (int nfft, double d1[], double d2[])
{
  void cdft (int n, int isgn, double *a);
  void mp_mul_rcmul (int n, double *a, double *b);
  double xr, xi;

  cdft (nfft, 1, &d1[1]);
  cdft (nfft, 1, &d2[1]);
  d2[0] += d1[0];
  xr = d1[1] * d2[1] + d1[2] * d2[2];
  xi = d1[1] * d2[2] + d1[2] * d2[1];
  d2[1] = xr;
  d2[2] = xi;
  if (nfft > 2)
    {
      mp_mul_rcmul (nfft, &d1[1], &d2[1]);
    }
  d2[nfft + 1] *= d1[nfft + 1];
  cdft (nfft, -1, &d2[1]);
}


void
mp_mul_cmul_nt_d1 (int nfft, double d1[], double d2[])
{
  void cdft (int n, int isgn, double *a);
  void mp_mul_rcmul_nt_in1 (int n, double *a, double *b);
  double xr, xi;

  cdft (nfft, 1, &d2[1]);
  d2[0] += d1[0];
  xr = d1[1] * d2[1] + d1[2] * d2[2];
  xi = d1[1] * d2[2] + d1[2] * d2[1];
  d2[1] = xr;
  d2[2] = xi;
  if (nfft > 2)
    {
      mp_mul_rcmul_nt_in1 (nfft, &d1[1], &d2[1]);
    }
  d2[nfft + 1] *= d1[nfft + 1];
  cdft (nfft, -1, &d2[1]);
}


void
mp_mul_cmul_nt_d2 (int nfft, double d1[], double d2[])
{
  void cdft (int n, int isgn, double *a);
  void mp_mul_rcmul_nt_in2 (int n, double *a, double *b);
  double xr, xi;

  cdft (nfft, 1, &d1[1]);
  d2[0] += d1[0];
  xr = d1[1] * d2[1] + d1[2] * d2[2];
  xi = d1[1] * d2[2] + d1[2] * d2[1];
  d2[1] = xr;
  d2[2] = xi;
  if (nfft > 2)
    {
      mp_mul_rcmul_nt_in2 (nfft, &d1[1], &d2[1]);
    }
  d2[nfft + 1] *= d1[nfft + 1];
  cdft (nfft, -1, &d2[1]);
}


void
mp_mul_cmul_nt_out (int nfft, double d1[], double d2[])
{
  void cdft (int n, int isgn, double *a);
  void mp_mul_rcmul_nt_out (int n, double *a, double *b);
  double xr, xi;

  cdft (nfft, 1, &d1[1]);
  cdft (nfft, 1, &d2[1]);
  d2[0] += d1[0];
  xr = d1[1] * d2[1] + d1[2] * d2[2];
  xi = d1[1] * d2[2] + d1[2] * d2[1];
  d2[1] = xr;
  d2[2] = xi;
  if (nfft > 2)
    {
      mp_mul_rcmul_nt_out (nfft, &d1[1], &d2[1]);
    }
  d2[nfft + 1] *= d1[nfft + 1];
}


void
mp_mul_cmul_nt_d1_add (int nfft, double d1[], double d2[], double d3[])
{
  void cdft (int n, int isgn, double *a);
  void mp_mul_rcmul_nt_in1_add (int n, double *a, double *b, double *badd);
  double xr, xi;

  cdft (nfft, 1, &d2[1]);
  xr = d1[1] * d2[1] + d1[2] * d2[2];
  xi = d1[1] * d2[2] + d1[2] * d2[1];
  d3[1] += xr;
  d3[2] += xi;
  if (nfft > 2)
    {
      mp_mul_rcmul_nt_in1_add (nfft, &d1[1], &d2[1], &d3[1]);
    }
  d3[nfft + 1] += d1[nfft + 1] * d2[nfft + 1];
  cdft (nfft, -1, &d3[1]);
}


void
mp_mul_csqu (int nfft, double d1[])
{
  void cdft (int n, int isgn, double *a);
  void mp_mul_rcsqu (int n, double *a);
  double xr, xi;

  cdft (nfft, 1, &d1[1]);
  d1[0] *= 2;
  xr = d1[1] * d1[1] + d1[2] * d1[2];
  xi = 2 * d1[1] * d1[2];
  d1[1] = xr;
  d1[2] = xi;
  if (nfft > 2)
    {
      mp_mul_rcsqu (nfft, &d1[1]);
    }
  d1[nfft + 1] *= d1[nfft + 1];
  cdft (nfft, -1, &d1[1]);
}


void
mp_mul_csqu_save_d1 (int nfft, double d1[], double d2[])
{
  void cdft (int n, int isgn, double *a);
  void mp_mul_rcsqu_save (int n, double *a, double *b);
  double xr, xi;

  cdft (nfft, 1, &d1[1]);
  d2[0] = 2 * d1[0];
  xr = d1[1] * d1[1] + d1[2] * d1[2];
  xi = 2 * d1[1] * d1[2];
  d2[1] = xr;
  d2[2] = xi;
  if (nfft > 2)
    {
      mp_mul_rcsqu_save (nfft, &d1[1], &d2[1]);
    }
  d2[nfft + 1] = d1[nfft + 1] * d1[nfft + 1];
  cdft (nfft, -1, &d2[1]);
}


void
mp_mul_csqu_nt_d1 (int nfft, double d1[])
{
  void cdft (int n, int isgn, double *a);
  void mp_mul_rcsqu_nt_in (int n, double *a);
  double xr, xi;

  d1[0] *= 2;
  xr = d1[1] * d1[1] + d1[2] * d1[2];
  xi = 2 * d1[1] * d1[2];
  d1[1] = xr;
  d1[2] = xi;
  if (nfft > 2)
    {
      mp_mul_rcsqu_nt_in (nfft, &d1[1]);
    }
  d1[nfft + 1] *= d1[nfft + 1];
  cdft (nfft, -1, &d1[1]);
}


void
mp_mul_d2i (int n, int radix, int nfft, double din[], int out[])
{
  int j, carry, carry1, carry2, shift, ndata;
  double x, scale, d1_radix, d1_radix2, pow_radix, topdgt;
  DGTINT *outr;

  outr = ((DGTINT *) & out[2]) - 2;
  scale = 2.0 / nfft;
  d1_radix = 1.0 / radix;
  d1_radix2 = d1_radix * d1_radix;
  topdgt = din[nfft + 1];
  x = topdgt < 0 ? -topdgt : topdgt;
  shift = x + 0.5 >= radix ? 1 : 0;
  /* ---- correction of cyclic convolution of din[1] ---- */
  x *= nfft * 0.5;
  din[nfft + 1] = din[1] - x;
  din[1] = x;
  /* ---- output of digits ---- */
  ndata = n;
  if (n > nfft + 1 + shift)
    {
      ndata = nfft + 1 + shift;
      for (j = n + 1; j > ndata + 1; j--)
	{
	  outr[j] = 0;
	}
    }
  x = 0;
  pow_radix = 1;
  for (j = ndata + 1 - shift; j <= nfft + 1; j++)
    {
      x += pow_radix * din[j];
      pow_radix *= d1_radix;
      if (pow_radix < DBL_EPSILON)
	{
	  break;
	}
    }
  x = d1_radix2 * (scale * x + 0.5);
  carry2 = ((int) x) - 1;
  carry = (int) (radix * (x - carry2) + 0.5);
  for (j = ndata; j > 1; j--)
    {
      x = d1_radix2 * (scale * din[j - shift] + carry + 0.5);
      carry = carry2;
      carry2 = ((int) x) - 1;
      x = radix * (x - carry2);
      carry1 = (int) x;
      outr[j + 1] = (DGTINT) (radix * (x - carry1));
      carry += carry1;
    }
  x = carry + ((double) radix) * carry2 + 0.5;
  if (shift == 0)
    {
      x += scale * din[1];
    }
  carry = (int) (d1_radix * x);
  outr[2] = (DGTINT) (x - ((double) radix) * carry);
  if (carry > 0)
    {
      for (j = n + 1; j > 2; j--)
	{
	  outr[j] = outr[j - 1];
	}
      outr[2] = (DGTINT) carry;
      shift++;
    }
  /* ---- output of exp, sgn ---- */
  x = din[0] + shift + 0.5;
  shift = ((int) x) - 1;
  out[1] = shift + ((int) (x - shift));
  out[0] = topdgt > 0.5 ? 1 : -1;
  if (outr[2] == 0)
    {
      out[0] = 0;
      out[1] = 0;
    }
}


double
mp_mul_d2i_test (int radix, int nfft, double din[])
{
  int j, carry, carry1, carry2;
  double x, scale, d1_radix, d1_radix2, err;

  scale = 2.0 / nfft;
  d1_radix = 1.0 / radix;
  d1_radix2 = d1_radix * d1_radix;
  /* ---- correction of cyclic convolution of din[1] ---- */
  x = din[nfft + 1] * nfft * 0.5;
  if (x < 0)
    {
      x = -x;
    }
  din[nfft + 1] = din[1] - x;
  /* ---- check of digits ---- */
  err = 0;
  carry = 0;
  carry2 = 0;
  for (j = nfft + 1; j > 1; j--)
    {
      x = d1_radix2 * (scale * din[j] + carry + 0.5);
      carry = carry2;
      carry2 = ((int) x) - 1;
      x = radix * (x - carry2);
      carry1 = (int) x;
      x = radix * (x - carry1);
      carry += carry1;
      x = x - 0.5 - ((int) x);
      if (x > err)
	{
	  err = x;
	}
      else if (-x > err)
	{
	  err = -x;
	}
    }
  return err;
}


/* -------- mp_mul child^2 routines (mix RFFT routines) -------- */


#ifndef M_PI_2
#define M_PI_2      1.570796326794896619231321691639751442098584699687
#endif


#ifndef RDFT_LOOP_DIV		/* control of the RDFT's speed & tolerance */
#define RDFT_LOOP_DIV 64
#endif


void
mp_mul_rcmul (int n, double *a, double *b)
{
  int i, i0, j, k;
  double ec, w1r, w1i, wkr, wki, wdr, wdi, ss;
  double xr, xi, yr, yi, ajr, aji, akr, aki, bjr, bji, bkr, bki;

  ec = 2 * M_PI_2 / n;
  wkr = 0;
  wki = 0;
  wdi = cos (ec);
  wdr = sin (ec);
  wdi *= wdr;
  wdr *= wdr;
  w1r = 1 - 2 * wdr;
  w1i = 2 * wdi;
  ss = 2 * w1i;
  i = n >> 1;
  xr = a[i];
  xi = a[i + 1];
  yr = b[i];
  yi = b[i + 1];
  b[i] = xr * yr - xi * yi;
  b[i + 1] = xr * yi + xi * yr;
  for (;;)
    {
      i0 = i - 4 * RDFT_LOOP_DIV;
      if (i0 < 2)
	{
	  i0 = 2;
	}
      for (j = i - 2; j >= i0; j -= 2)
	{
	  k = n - j;
	  xr = wkr + ss * wdi;
	  xi = wki + ss * (0.5 - wdr);
	  wkr = wdr;
	  wki = wdi;
	  wdr = xr;
	  wdi = xi;
	  /* ---- transform CFFT data a[] into RFFT data ---- */
	  xr = a[j] - a[k];
	  xi = a[j + 1] + a[k + 1];
	  yr = wkr * xr - wki * xi;
	  yi = wkr * xi + wki * xr;
	  ajr = a[j] - yr;
	  aji = a[j + 1] - yi;
	  akr = a[k] + yr;
	  aki = a[k + 1] - yi;
	  a[j] = ajr;
	  a[j + 1] = aji;
	  a[k] = akr;
	  a[k + 1] = aki;
	  /* ---- transform CFFT data b[] into RFFT data ---- */
	  xr = b[j] - b[k];
	  xi = b[j + 1] + b[k + 1];
	  yr = wkr * xr - wki * xi;
	  yi = wkr * xi + wki * xr;
	  xr = b[j] - yr;
	  xi = b[j + 1] - yi;
	  yr = b[k] + yr;
	  yi = b[k + 1] - yi;
	  /* ---- cmul ---- */
	  bjr = ajr * xr - aji * xi;
	  bji = ajr * xi + aji * xr;
	  bkr = akr * yr - aki * yi;
	  bki = akr * yi + aki * yr;
	  /* ---- transform RFFT data bxx into CFFT data ---- */
	  xr = bjr - bkr;
	  xi = bji + bki;
	  yr = wkr * xr + wki * xi;
	  yi = wkr * xi - wki * xr;
	  b[j] = bjr - yr;
	  b[j + 1] = bji - yi;
	  b[k] = bkr + yr;
	  b[k + 1] = bki - yi;
	}
      if (i0 == 2)
	{
	  break;
	}
      wkr = 0.5 * sin (ec * i0);
      wki = 0.5 * cos (ec * i0);
      wdr = 0.5 - (wkr * w1r - wki * w1i);
      wdi = wkr * w1i + wki * w1r;
      wkr = 0.5 - wkr;
      i = i0;
    }
}


void
mp_mul_rcmul_nt_in1 (int n, double *a, double *b)
{
  int i, i0, j, k;
  double ec, w1r, w1i, wkr, wki, wdr, wdi, ss;
  double xr, xi, yr, yi, bjr, bji, bkr, bki;

  ec = 2 * M_PI_2 / n;
  wkr = 0;
  wki = 0;
  wdi = cos (ec);
  wdr = sin (ec);
  wdi *= wdr;
  wdr *= wdr;
  w1r = 1 - 2 * wdr;
  w1i = 2 * wdi;
  ss = 2 * w1i;
  i = n >> 1;
  xr = a[i];
  xi = a[i + 1];
  yr = b[i];
  yi = b[i + 1];
  b[i] = xr * yr - xi * yi;
  b[i + 1] = xr * yi + xi * yr;
  for (;;)
    {
      i0 = i - 4 * RDFT_LOOP_DIV;
      if (i0 < 2)
	{
	  i0 = 2;
	}
      for (j = i - 2; j >= i0; j -= 2)
	{
	  k = n - j;
	  xr = wkr + ss * wdi;
	  xi = wki + ss * (0.5 - wdr);
	  wkr = wdr;
	  wki = wdi;
	  wdr = xr;
	  wdi = xi;
	  /* ---- transform CFFT data b[] into RFFT data ---- */
	  xr = b[j] - b[k];
	  xi = b[j + 1] + b[k + 1];
	  yr = wkr * xr - wki * xi;
	  yi = wkr * xi + wki * xr;
	  xr = b[j] - yr;
	  xi = b[j + 1] - yi;
	  yr = b[k] + yr;
	  yi = b[k + 1] - yi;
	  /* ---- cmul ---- */
	  bjr = a[j] * xr - a[j + 1] * xi;
	  bji = a[j] * xi + a[j + 1] * xr;
	  bkr = a[k] * yr - a[k + 1] * yi;
	  bki = a[k] * yi + a[k + 1] * yr;
	  /* ---- transform RFFT data bxx into CFFT data ---- */
	  xr = bjr - bkr;
	  xi = bji + bki;
	  yr = wkr * xr + wki * xi;
	  yi = wkr * xi - wki * xr;
	  b[j] = bjr - yr;
	  b[j + 1] = bji - yi;
	  b[k] = bkr + yr;
	  b[k + 1] = bki - yi;
	}
      if (i0 == 2)
	{
	  break;
	}
      wkr = 0.5 * sin (ec * i0);
      wki = 0.5 * cos (ec * i0);
      wdr = 0.5 - (wkr * w1r - wki * w1i);
      wdi = wkr * w1i + wki * w1r;
      wkr = 0.5 - wkr;
      i = i0;
    }
}


void
mp_mul_rcmul_nt_in2 (int n, double *a, double *b)
{
  int i, i0, j, k;
  double ec, w1r, w1i, wkr, wki, wdr, wdi, ss;
  double xr, xi, yr, yi, bjr, bji, bkr, bki;

  ec = 2 * M_PI_2 / n;
  wkr = 0;
  wki = 0;
  wdi = cos (ec);
  wdr = sin (ec);
  wdi *= wdr;
  wdr *= wdr;
  w1r = 1 - 2 * wdr;
  w1i = 2 * wdi;
  ss = 2 * w1i;
  i = n >> 1;
  xr = a[i];
  xi = a[i + 1];
  yr = b[i];
  yi = b[i + 1];
  b[i] = xr * yr - xi * yi;
  b[i + 1] = xr * yi + xi * yr;
  for (;;)
    {
      i0 = i - 4 * RDFT_LOOP_DIV;
      if (i0 < 2)
	{
	  i0 = 2;
	}
      for (j = i - 2; j >= i0; j -= 2)
	{
	  k = n - j;
	  xr = wkr + ss * wdi;
	  xi = wki + ss * (0.5 - wdr);
	  wkr = wdr;
	  wki = wdi;
	  wdr = xr;
	  wdi = xi;
	  /* ---- transform CFFT data a[] into RFFT data ---- */
	  xr = a[j] - a[k];
	  xi = a[j + 1] + a[k + 1];
	  yr = wkr * xr - wki * xi;
	  yi = wkr * xi + wki * xr;
	  xr = a[j] - yr;
	  xi = a[j + 1] - yi;
	  yr = a[k] + yr;
	  yi = a[k + 1] - yi;
	  a[j] = xr;
	  a[j + 1] = xi;
	  a[k] = yr;
	  a[k + 1] = yi;
	  /* ---- cmul ---- */
	  bjr = b[j] * xr - b[j + 1] * xi;
	  bji = b[j] * xi + b[j + 1] * xr;
	  bkr = b[k] * yr - b[k + 1] * yi;
	  bki = b[k] * yi + b[k + 1] * yr;
	  /* ---- transform RFFT data bxx into CFFT data ---- */
	  xr = bjr - bkr;
	  xi = bji + bki;
	  yr = wkr * xr + wki * xi;
	  yi = wkr * xi - wki * xr;
	  b[j] = bjr - yr;
	  b[j + 1] = bji - yi;
	  b[k] = bkr + yr;
	  b[k + 1] = bki - yi;
	}
      if (i0 == 2)
	{
	  break;
	}
      wkr = 0.5 * sin (ec * i0);
      wki = 0.5 * cos (ec * i0);
      wdr = 0.5 - (wkr * w1r - wki * w1i);
      wdi = wkr * w1i + wki * w1r;
      wkr = 0.5 - wkr;
      i = i0;
    }
}


void
mp_mul_rcmul_nt_out (int n, double *a, double *b)
{
  int i, i0, j, k;
  double ec, w1r, w1i, wkr, wki, wdr, wdi, ss;
  double xr, xi, yr, yi, ajr, aji, akr, aki;

  ec = 2 * M_PI_2 / n;
  wkr = 0;
  wki = 0;
  wdi = cos (ec);
  wdr = sin (ec);
  wdi *= wdr;
  wdr *= wdr;
  w1r = 1 - 2 * wdr;
  w1i = 2 * wdi;
  ss = 2 * w1i;
  i = n >> 1;
  xr = a[i];
  xi = a[i + 1];
  yr = b[i];
  yi = b[i + 1];
  b[i] = xr * yr - xi * yi;
  b[i + 1] = xr * yi + xi * yr;
  for (;;)
    {
      i0 = i - 4 * RDFT_LOOP_DIV;
      if (i0 < 2)
	{
	  i0 = 2;
	}
      for (j = i - 2; j >= i0; j -= 2)
	{
	  k = n - j;
	  xr = wkr + ss * wdi;
	  xi = wki + ss * (0.5 - wdr);
	  wkr = wdr;
	  wki = wdi;
	  wdr = xr;
	  wdi = xi;
	  /* ---- transform CFFT data a[] into RFFT data ---- */
	  xr = a[j] - a[k];
	  xi = a[j + 1] + a[k + 1];
	  yr = wkr * xr - wki * xi;
	  yi = wkr * xi + wki * xr;
	  ajr = a[j] - yr;
	  aji = a[j + 1] - yi;
	  akr = a[k] + yr;
	  aki = a[k + 1] - yi;
	  a[j] = ajr;
	  a[j + 1] = aji;
	  a[k] = akr;
	  a[k + 1] = aki;
	  /* ---- transform CFFT data b[] into RFFT data ---- */
	  xr = b[j] - b[k];
	  xi = b[j + 1] + b[k + 1];
	  yr = wkr * xr - wki * xi;
	  yi = wkr * xi + wki * xr;
	  xr = b[j] - yr;
	  xi = b[j + 1] - yi;
	  yr = b[k] + yr;
	  yi = b[k + 1] - yi;
	  /* ---- cmul ---- */
	  b[j] = ajr * xr - aji * xi;
	  b[j + 1] = ajr * xi + aji * xr;
	  b[k] = akr * yr - aki * yi;
	  b[k + 1] = akr * yi + aki * yr;
	}
      if (i0 == 2)
	{
	  break;
	}
      wkr = 0.5 * sin (ec * i0);
      wki = 0.5 * cos (ec * i0);
      wdr = 0.5 - (wkr * w1r - wki * w1i);
      wdi = wkr * w1i + wki * w1r;
      wkr = 0.5 - wkr;
      i = i0;
    }
}


void
mp_mul_rcmul_nt_in1_add (int n, double *a, double *b, double *badd)
{
  int i, i0, j, k;
  double ec, w1r, w1i, wkr, wki, wdr, wdi, ss;
  double xr, xi, yr, yi, bjr, bji, bkr, bki;

  ec = 2 * M_PI_2 / n;
  wkr = 0;
  wki = 0;
  wdi = cos (ec);
  wdr = sin (ec);
  wdi *= wdr;
  wdr *= wdr;
  w1r = 1 - 2 * wdr;
  w1i = 2 * wdi;
  ss = 2 * w1i;
  i = n >> 1;
  xr = a[i];
  xi = a[i + 1];
  yr = b[i];
  yi = b[i + 1];
  badd[i] += xr * yr - xi * yi;
  badd[i + 1] += xr * yi + xi * yr;
  for (;;)
    {
      i0 = i - 4 * RDFT_LOOP_DIV;
      if (i0 < 2)
	{
	  i0 = 2;
	}
      for (j = i - 2; j >= i0; j -= 2)
	{
	  k = n - j;
	  xr = wkr + ss * wdi;
	  xi = wki + ss * (0.5 - wdr);
	  wkr = wdr;
	  wki = wdi;
	  wdr = xr;
	  wdi = xi;
	  /* ---- transform CFFT data b[] into RFFT data ---- */
	  xr = b[j] - b[k];
	  xi = b[j + 1] + b[k + 1];
	  yr = wkr * xr - wki * xi;
	  yi = wkr * xi + wki * xr;
	  xr = b[j] - yr;
	  xi = b[j + 1] - yi;
	  yr = b[k] + yr;
	  yi = b[k + 1] - yi;
	  /* ---- cmul + add ---- */
	  bjr = badd[j] + (a[j] * xr - a[j + 1] * xi);
	  bji = badd[j + 1] + (a[j] * xi + a[j + 1] * xr);
	  bkr = badd[k] + (a[k] * yr - a[k + 1] * yi);
	  bki = badd[k + 1] + (a[k] * yi + a[k + 1] * yr);
	  /* ---- transform RFFT data bxx into CFFT data ---- */
	  xr = bjr - bkr;
	  xi = bji + bki;
	  yr = wkr * xr + wki * xi;
	  yi = wkr * xi - wki * xr;
	  badd[j] = bjr - yr;
	  badd[j + 1] = bji - yi;
	  badd[k] = bkr + yr;
	  badd[k + 1] = bki - yi;
	}
      if (i0 == 2)
	{
	  break;
	}
      wkr = 0.5 * sin (ec * i0);
      wki = 0.5 * cos (ec * i0);
      wdr = 0.5 - (wkr * w1r - wki * w1i);
      wdi = wkr * w1i + wki * w1r;
      wkr = 0.5 - wkr;
      i = i0;
    }
}


void
mp_mul_rcsqu (int n, double *a)
{
  int i, i0, j, k;
  double ec, w1r, w1i, wkr, wki, wdr, wdi, ss;
  double xr, xi, yr, yi, ajr, aji, akr, aki;

  ec = 2 * M_PI_2 / n;
  wkr = 0;
  wki = 0;
  wdi = cos (ec);
  wdr = sin (ec);
  wdi *= wdr;
  wdr *= wdr;
  w1r = 1 - 2 * wdr;
  w1i = 2 * wdi;
  ss = 2 * w1i;
  i = n >> 1;
  xr = a[i];
  xi = a[i + 1];
  a[i] = xr * xr - xi * xi;
  a[i + 1] = 2 * xr * xi;
  for (;;)
    {
      i0 = i - 4 * RDFT_LOOP_DIV;
      if (i0 < 2)
	{
	  i0 = 2;
	}
      for (j = i - 2; j >= i0; j -= 2)
	{
	  k = n - j;
	  xr = wkr + ss * wdi;
	  xi = wki + ss * (0.5 - wdr);
	  wkr = wdr;
	  wki = wdi;
	  wdr = xr;
	  wdi = xi;
	  /* ---- transform CFFT data a[] into RFFT data ---- */
	  xr = a[j] - a[k];
	  xi = a[j + 1] + a[k + 1];
	  yr = wkr * xr - wki * xi;
	  yi = wkr * xi + wki * xr;
	  xr = a[j] - yr;
	  xi = a[j + 1] - yi;
	  yr = a[k] + yr;
	  yi = a[k + 1] - yi;
	  /* ---- csqu ---- */
	  ajr = xr * xr - xi * xi;
	  aji = 2 * xr * xi;
	  akr = yr * yr - yi * yi;
	  aki = 2 * yr * yi;
	  /* ---- transform RFFT data axx into CFFT data ---- */
	  xr = ajr - akr;
	  xi = aji + aki;
	  yr = wkr * xr + wki * xi;
	  yi = wkr * xi - wki * xr;
	  a[j] = ajr - yr;
	  a[j + 1] = aji - yi;
	  a[k] = akr + yr;
	  a[k + 1] = aki - yi;
	}
      if (i0 == 2)
	{
	  break;
	}
      wkr = 0.5 * sin (ec * i0);
      wki = 0.5 * cos (ec * i0);
      wdr = 0.5 - (wkr * w1r - wki * w1i);
      wdi = wkr * w1i + wki * w1r;
      wkr = 0.5 - wkr;
      i = i0;
    }
}


void
mp_mul_rcsqu_save (int n, double *a, double *b)
{
  int i, i0, j, k;
  double ec, w1r, w1i, wkr, wki, wdr, wdi, ss;
  double xr, xi, yr, yi, ajr, aji, akr, aki;

  ec = 2 * M_PI_2 / n;
  wkr = 0;
  wki = 0;
  wdi = cos (ec);
  wdr = sin (ec);
  wdi *= wdr;
  wdr *= wdr;
  w1r = 1 - 2 * wdr;
  w1i = 2 * wdi;
  ss = 2 * w1i;
  i = n >> 1;
  xr = a[i];
  xi = a[i + 1];
  b[i] = xr * xr - xi * xi;
  b[i + 1] = 2 * xr * xi;
  for (;;)
    {
      i0 = i - 4 * RDFT_LOOP_DIV;
      if (i0 < 2)
	{
	  i0 = 2;
	}
      for (j = i - 2; j >= i0; j -= 2)
	{
	  k = n - j;
	  xr = wkr + ss * wdi;
	  xi = wki + ss * (0.5 - wdr);
	  wkr = wdr;
	  wki = wdi;
	  wdr = xr;
	  wdi = xi;
	  /* ---- transform CFFT data a[] into RFFT data ---- */
	  xr = a[j] - a[k];
	  xi = a[j + 1] + a[k + 1];
	  yr = wkr * xr - wki * xi;
	  yi = wkr * xi + wki * xr;
	  xr = a[j] - yr;
	  xi = a[j + 1] - yi;
	  yr = a[k] + yr;
	  yi = a[k + 1] - yi;
	  a[j] = xr;
	  a[j + 1] = xi;
	  a[k] = yr;
	  a[k + 1] = yi;
	  /* ---- csqu ---- */
	  ajr = xr * xr - xi * xi;
	  aji = 2 * xr * xi;
	  akr = yr * yr - yi * yi;
	  aki = 2 * yr * yi;
	  /* ---- transform RFFT data axx into CFFT data ---- */
	  xr = ajr - akr;
	  xi = aji + aki;
	  yr = wkr * xr + wki * xi;
	  yi = wkr * xi - wki * xr;
	  b[j] = ajr - yr;
	  b[j + 1] = aji - yi;
	  b[k] = akr + yr;
	  b[k + 1] = aki - yi;
	}
      if (i0 == 2)
	{
	  break;
	}
      wkr = 0.5 * sin (ec * i0);
      wki = 0.5 * cos (ec * i0);
      wdr = 0.5 - (wkr * w1r - wki * w1i);
      wdi = wkr * w1i + wki * w1r;
      wkr = 0.5 - wkr;
      i = i0;
    }
}


void
mp_mul_rcsqu_nt_in (int n, double *a)
{
  int i, i0, j, k;
  double ec, w1r, w1i, wkr, wki, wdr, wdi, ss;
  double xr, xi, yr, yi, ajr, aji, akr, aki;

  ec = 2 * M_PI_2 / n;
  wkr = 0;
  wki = 0;
  wdi = cos (ec);
  wdr = sin (ec);
  wdi *= wdr;
  wdr *= wdr;
  w1r = 1 - 2 * wdr;
  w1i = 2 * wdi;
  ss = 2 * w1i;
  i = n >> 1;
  xr = a[i];
  xi = a[i + 1];
  a[i] = xr * xr - xi * xi;
  a[i + 1] = 2 * xr * xi;
  for (;;)
    {
      i0 = i - 4 * RDFT_LOOP_DIV;
      if (i0 < 2)
	{
	  i0 = 2;
	}
      for (j = i - 2; j >= i0; j -= 2)
	{
	  k = n - j;
	  xr = wkr + ss * wdi;
	  xi = wki + ss * (0.5 - wdr);
	  wkr = wdr;
	  wki = wdi;
	  wdr = xr;
	  wdi = xi;
	  /* ---- csqu ---- */
	  xr = a[j];
	  xi = a[j + 1];
	  yr = a[k];
	  yi = a[k + 1];
	  ajr = xr * xr - xi * xi;
	  aji = 2 * xr * xi;
	  akr = yr * yr - yi * yi;
	  aki = 2 * yr * yi;
	  /* ---- transform RFFT data axx into CFFT data ---- */
	  xr = ajr - akr;
	  xi = aji + aki;
	  yr = wkr * xr + wki * xi;
	  yi = wkr * xi - wki * xr;
	  a[j] = ajr - yr;
	  a[j + 1] = aji - yi;
	  a[k] = akr + yr;
	  a[k + 1] = aki - yi;
	}
      if (i0 == 2)
	{
	  break;
	}
      wkr = 0.5 * sin (ec * i0);
      wki = 0.5 * cos (ec * i0);
      wdr = 0.5 - (wkr * w1r - wki * w1i);
      wdi = wkr * w1i + wki * w1r;
      wkr = 0.5 - wkr;
      i = i0;
    }
}


/* -------- mp_inv routines -------- */


int
mp_inv (int n, int radix, int in[], int out[],
	int tmp1[], int tmp2[], int nfft, double tmp1fft[], double tmp2fft[])
{
  int mp_get_nfft_init (int radix, int nfft_max);
  void mp_inv_init (int n, int radix, int in[], int out[]);
  int mp_inv_newton (int n, int radix, int in[], int inout[],
		     int tmp1[], int tmp2[], int nfft, double tmp1fft[],
		     double tmp2fft[]);
  int n_nwt, nfft_nwt, thr, prc;

  if (in[0] == 0)
    {
      return -1;
    }
  nfft_nwt = mp_get_nfft_init (radix, nfft);
  n_nwt = nfft_nwt + 2;
  if (n_nwt > n)
    {
      n_nwt = n;
    }
  mp_inv_init (n_nwt, radix, in, out);
  thr = 8;
  do
    {
      n_nwt = nfft_nwt + 2;
      if (n_nwt > n)
	{
	  n_nwt = n;
	}
      prc = mp_inv_newton (n_nwt, radix, in, out,
			   tmp1, tmp2, nfft_nwt, tmp1fft, tmp2fft);
#ifdef DEBUG
      printf ("n=%d, nfft=%d, prc=%d\n", n_nwt, nfft_nwt, prc);
#endif
      if (thr * nfft_nwt >= nfft)
	{
	  thr = 0;
	  if (2 * prc <= n_nwt - 2)
	    {
	      nfft_nwt >>= 1;
	    }
	}
      else
	{
	  if (3 * prc < n_nwt - 2)
	    {
	      nfft_nwt >>= 1;
	    }
	}
      nfft_nwt <<= 1;
    }
  while (nfft_nwt <= nfft);
  return 0;
}


int
mp_sqrt (int n, int radix, int in[], int out[],
	 int tmp1[], int tmp2[], int nfft, double tmp1fft[], double tmp2fft[])
{
  void mp_load_0 (int n, int radix, int out[]);
  int mp_get_nfft_init (int radix, int nfft_max);
  void mp_sqrt_init (int n, int radix, int in[], int out[], int out_rev[]);
  int mp_sqrt_newton (int n, int radix, int in[], int inout[],
		      int inout_rev[], int tmp[], int nfft, double tmp1fft[],
		      double tmp2fft[], int *n_tmp1fft);
  int n_nwt, nfft_nwt, thr, prc, n_tmp1fft;

  if (in[0] < 0)
    {
      return -1;
    }
  else if (in[0] == 0)
    {
      mp_load_0 (n, radix, out);
      return 0;
    }
  nfft_nwt = mp_get_nfft_init (radix, nfft);
  n_nwt = nfft_nwt + 2;
  if (n_nwt > n)
    {
      n_nwt = n;
    }
  mp_sqrt_init (n_nwt, radix, in, out, tmp1);
  n_tmp1fft = 0;
  thr = 8;
  do
    {
      n_nwt = nfft_nwt + 2;
      if (n_nwt > n)
	{
	  n_nwt = n;
	}
      prc = mp_sqrt_newton (n_nwt, radix, in, out,
			    tmp1, tmp2, nfft_nwt, tmp1fft, tmp2fft,
			    &n_tmp1fft);
#ifdef DEBUG
      printf ("n=%d, nfft=%d, prc=%d\n", n_nwt, nfft_nwt, prc);
#endif
      if (thr * nfft_nwt >= nfft)
	{
	  thr = 0;
	  if (2 * prc <= n_nwt - 2)
	    {
	      nfft_nwt >>= 1;
	    }
	}
      else
	{
	  if (3 * prc < n_nwt - 2)
	    {
	      nfft_nwt >>= 1;
	    }
	}
      nfft_nwt <<= 1;
    }
  while (nfft_nwt <= nfft);
  return 0;
}


int
mp_invisqrt (int n, int radix, int in, int out[],
	     int tmp1[], int tmp2[], int nfft,
	     double tmp1fft[], double tmp2fft[])
{
  int mp_get_nfft_init (int radix, int nfft_max);
  void mp_invisqrt_init (int n, int radix, int in, int out[]);
  int mp_invisqrt_newton (int n, int radix, int in, int inout[],
			  int tmp1[], int tmp2[], int nfft, double tmp1fft[],
			  double tmp2fft[]);
  int n_nwt, nfft_nwt, thr, prc;

  if (in <= 0)
    {
      return -1;
    }
  nfft_nwt = mp_get_nfft_init (radix, nfft);
  n_nwt = nfft_nwt + 2;
  if (n_nwt > n)
    {
      n_nwt = n;
    }
  mp_invisqrt_init (n_nwt, radix, in, out);
  thr = 8;
  do
    {
      n_nwt = nfft_nwt + 2;
      if (n_nwt > n)
	{
	  n_nwt = n;
	}
      prc = mp_invisqrt_newton (n_nwt, radix, in, out,
				tmp1, tmp2, nfft_nwt, tmp1fft, tmp2fft);
#ifdef DEBUG
      printf ("n=%d, nfft=%d, prc=%d\n", n_nwt, nfft_nwt, prc);
#endif
      if (thr * nfft_nwt >= nfft)
	{
	  thr = 0;
	  if (2 * prc <= n_nwt - 2)
	    {
	      nfft_nwt >>= 1;
	    }
	}
      else
	{
	  if (3 * prc < n_nwt - 2)
	    {
	      nfft_nwt >>= 1;
	    }
	}
      nfft_nwt <<= 1;
    }
  while (nfft_nwt <= nfft);
  return 0;
}


/* -------- mp_inv child routines -------- */


int
mp_get_nfft_init (int radix, int nfft_max)
{
  int nfft_init;
  double r;

  r = radix;
  nfft_init = 1;
  do
    {
      r *= r;
      nfft_init <<= 1;
    }
  while (DBL_EPSILON * r < 1 && nfft_init < nfft_max);
  return nfft_init;
}


void
mp_inv_init (int n, int radix, int in[], int out[])
{
  void mp_unexp_d2mp (int n, int radix, double din, DGTINT out[]);
  double mp_unexp_mp2d (int n, int radix, DGTINT in[]);
  int outexp;
  double din;

  out[0] = in[0];
  outexp = -in[1];
  din = 1.0 / mp_unexp_mp2d (n, radix, (DGTINT *) & in[2]);
  while (din < 1)
    {
      din *= radix;
      outexp--;
    }
  out[1] = outexp;
  mp_unexp_d2mp (n, radix, din, (DGTINT *) & out[2]);
}


void
mp_sqrt_init (int n, int radix, int in[], int out[], int out_rev[])
{
  void mp_unexp_d2mp (int n, int radix, double din, DGTINT out[]);
  double mp_unexp_mp2d (int n, int radix, DGTINT in[]);
  int outexp;
  double din;

  out[0] = 1;
  out_rev[0] = 1;
  outexp = in[1];
  din = mp_unexp_mp2d (n, radix, (DGTINT *) & in[2]);
  if (outexp % 2 != 0)
    {
      din *= radix;
      outexp--;
    }
  outexp /= 2;
  din = sqrt (din);
  if (din < 1)
    {
      din *= radix;
      outexp--;
    }
  out[1] = outexp;
  mp_unexp_d2mp (n, radix, din, (DGTINT *) & out[2]);
  outexp = -outexp;
  din = 1.0 / din;
  while (din < 1)
    {
      din *= radix;
      outexp--;
    }
  out_rev[1] = outexp;
  mp_unexp_d2mp (n, radix, din, (DGTINT *) & out_rev[2]);
}


void
mp_invisqrt_init (int n, int radix, int in, int out[])
{
  void mp_unexp_d2mp (int n, int radix, double din, DGTINT out[]);
  int outexp;
  double dout;

  out[0] = 1;
  outexp = 0;
  dout = sqrt (1.0 / in);
  while (dout < 1)
    {
      dout *= radix;
      outexp--;
    }
  out[1] = outexp;
  mp_unexp_d2mp (n, radix, dout, (DGTINT *) & out[2]);
}


void
mp_unexp_d2mp (int n, int radix, double din, DGTINT out[])
{
  int j, x;

  for (j = 0; j < n; j++)
    {
      x = (int) din;
      if (x >= radix)
	{
	  x = radix - 1;
	  din = radix;
	}
      din = radix * (din - x);
      out[j] = (DGTINT) x;
    }
}


double
mp_unexp_mp2d (int n, int radix, DGTINT in[])
{
  int j;
  double d1_radix, dout;

  d1_radix = 1.0 / radix;
  dout = 0;
  for (j = n - 1; j >= 0; j--)
    {
      dout = d1_radix * dout + in[j];
    }
  return dout;
}


int
mp_inv_newton (int n, int radix, int in[], int inout[],
	       int tmp1[], int tmp2[], int nfft, double tmp1fft[],
	       double tmp2fft[])
{
  void mp_load_1 (int n, int radix, int out[]);
  void mp_round (int n, int radix, int m, int inout[]);
  void mp_add (int n, int radix, int in1[], int in2[], int out[]);
  void mp_sub (int n, int radix, int in1[], int in2[], int out[]);
  void mp_mulh (int n, int radix, int in1[], int in2[], int out[],
		int nfft, double in1fft[], double outfft[]);
  void mp_mulh_use_in1fft (int n, int radix, double in1fft[],
			   int shift, int in2[], int out[], int nfft,
			   double outfft[]);
  int n_h, shift, prc;

  shift = (nfft >> 1) + 1;
  n_h = n / 2 + 1;
  if (n_h < n - shift)
    {
      n_h = n - shift;
    }
  /* ---- tmp1 = inout * (upper) in (half to normal precision) ---- */
  mp_round (n, radix, shift, inout);
  mp_mulh (n, radix, inout, in, tmp1, nfft, tmp1fft, tmp2fft);
  /* ---- tmp2 = 1 - tmp1 ---- */
  mp_load_1 (n, radix, tmp2);
  mp_sub (n, radix, tmp2, tmp1, tmp2);
  /* ---- tmp2 -= inout * (lower) in (half precision) ---- */
  mp_mulh_use_in1fft (n, radix, tmp1fft, shift, in, tmp1, nfft, tmp2fft);
  mp_sub (n_h, radix, tmp2, tmp1, tmp2);
  /* ---- get precision ---- */
  prc = -tmp2[1];
  if (tmp2[0] == 0)
    {
      prc = nfft + 1;
    }
  /* ---- tmp2 *= inout (half precision) ---- */
  mp_mulh_use_in1fft (n_h, radix, tmp1fft, 0, tmp2, tmp2, nfft, tmp2fft);
  /* ---- inout += tmp2 ---- */
  mp_add (n, radix, inout, tmp2, inout);
  return prc;
}


int
mp_sqrt_newton (int n, int radix, int in[], int inout[],
		int inout_rev[], int tmp[], int nfft, double tmp1fft[],
		double tmp2fft[], int *n_tmp1fft)
{
  void mp_round (int n, int radix, int m, int inout[]);
  void mp_add (int n, int radix, int in1[], int in2[], int out[]);
  void mp_sub (int n, int radix, int in1[], int in2[], int out[]);
  void mp_idiv_2 (int n, int radix, int in[], int out[]);
  void mp_mulh (int n, int radix, int in1[], int in2[], int out[],
		int nfft, double in1fft[], double outfft[]);
  void mp_squh (int n, int radix, int in[], int out[],
		int nfft, double outfft[]);
  void mp_squh_use_in1fft (int n, int radix, double inoutfft[], int out[],
			   int nfft);
  int n_h, nfft_h, shift, prc;

  nfft_h = nfft >> 1;
  shift = nfft_h + 1;
  if (nfft_h < 2)
    {
      nfft_h = 2;
    }
  n_h = n / 2 + 1;
  if (n_h < n - shift)
    {
      n_h = n - shift;
    }
  /* ---- tmp = inout_rev^2 (1/4 to half precision) ---- */
  mp_round (n_h, radix, (nfft_h >> 1) + 1, inout_rev);
  if (*n_tmp1fft != nfft_h)
    {
      mp_squh (n_h, radix, inout_rev, tmp, nfft_h, tmp1fft);
    }
  else
    {
      mp_squh_use_in1fft (n_h, radix, tmp1fft, tmp, nfft_h);
    }
  /* ---- tmp = inout_rev - inout * tmp (half precision) ---- */
  mp_round (n, radix, shift, inout);
  mp_mulh (n_h, radix, inout, tmp, tmp, nfft, tmp1fft, tmp2fft);
  mp_sub (n_h, radix, inout_rev, tmp, tmp);
  /* ---- inout_rev += tmp ---- */
  mp_add (n_h, radix, inout_rev, tmp, inout_rev);
  /* ---- tmp = in - inout^2 (half to normal precision) ---- */
  mp_squh_use_in1fft (n, radix, tmp1fft, tmp, nfft);
  mp_sub (n, radix, in, tmp, tmp);
  /* ---- get precision ---- */
  prc = in[1] - tmp[1];
  if (((DGTINT *) & in[2])[0] > ((DGTINT *) & tmp[2])[0])
    {
      prc++;
    }
  if (tmp[0] == 0)
    {
      prc = nfft + 1;
    }
  /* ---- tmp = tmp * inout_rev / 2 (half precision) ---- */
  mp_round (n_h, radix, shift, inout_rev);
  mp_mulh (n_h, radix, inout_rev, tmp, tmp, nfft, tmp1fft, tmp2fft);
  *n_tmp1fft = nfft;
  mp_idiv_2 (n_h, radix, tmp, tmp);
  /* ---- inout += tmp ---- */
  mp_add (n, radix, inout, tmp, inout);
  return prc;
}


int
mp_invisqrt_newton (int n, int radix, int in, int inout[],
		    int tmp1[], int tmp2[], int nfft, double tmp1fft[],
		    double tmp2fft[])
{
  void mp_load_1 (int n, int radix, int out[]);
  void mp_round (int n, int radix, int m, int inout[]);
  void mp_add (int n, int radix, int in1[], int in2[], int out[]);
  void mp_sub (int n, int radix, int in1[], int in2[], int out[]);
  void mp_imul (int n, int radix, int in1[], int in2, int out[]);
  void mp_idiv_2 (int n, int radix, int in[], int out[]);
  void mp_squh_save_infft (int n, int radix, int in[], int out[],
			   int nfft, double infft[], double outfft[]);
  void mp_mulh_use_in1fft (int n, int radix, double in1fft[],
			   int shift, int in2[], int out[], int nfft,
			   double outfft[]);
  int n_h, shift, prc;

  shift = (nfft >> 1) + 1;
  n_h = n / 2 + 1;
  if (n_h < n - shift)
    {
      n_h = n - shift;
    }
  /* ---- tmp1 = in * inout^2 (half to normal precision) ---- */
  mp_round (n, radix, shift, inout);
  mp_squh_save_infft (n, radix, inout, tmp1, nfft, tmp1fft, tmp2fft);
  mp_imul (n, radix, tmp1, in, tmp1);
  /* ---- tmp2 = 1 - tmp1 ---- */
  mp_load_1 (n, radix, tmp2);
  mp_sub (n, radix, tmp2, tmp1, tmp2);
  /* ---- get precision ---- */
  prc = -tmp2[1];
  if (tmp2[0] == 0)
    {
      prc = nfft + 1;
    }
  /* ---- tmp2 *= inout / 2 (half precision) ---- */
  mp_mulh_use_in1fft (n_h, radix, tmp1fft, 0, tmp2, tmp2, nfft, tmp2fft);
  mp_idiv_2 (n_h, radix, tmp2, tmp2);
  /* ---- inout += tmp2 ---- */
  mp_add (n, radix, inout, tmp2, inout);
  return prc;
}


/* -------- mp_io routines -------- */


void
mp_sprintf (int n, int log10_radix, int in[], char out[])
{
  int j, k, x, y, outexp, shift;
  DGTINT *inr;

  inr = ((DGTINT *) & in[2]) - 2;
  if (in[0] < 0)
    {
      *out++ = '-';
    }
  x = inr[2];
  shift = log10_radix;
  for (k = log10_radix; k > 0; k--)
    {
      y = x % 10;
      x /= 10;
      out[k] = '0' + y;
      if (y != 0)
	{
	  shift = k;
	}
    }
  out[0] = out[shift];
  out[1] = '.';
  for (k = 1; k <= log10_radix - shift; k++)
    {
      out[k + 1] = out[k + shift];
    }
  outexp = log10_radix - shift;
  out += outexp + 2;
  for (j = 3; j <= n + 1; j++)
    {
      x = inr[j];
      for (k = log10_radix - 1; k >= 0; k--)
	{
	  y = x % 10;
	  x /= 10;
	  out[k] = '0' + y;
	}
      out += log10_radix;
    }
  *out++ = 'e';
  outexp += log10_radix * in[1];
  sprintf (out, "%d", outexp);
}


void
mp_sscanf (int n, int log10_radix, char in[], int out[])
{
  char *s;
  int j, x, outexp, outexp_mod;
  DGTINT *outr;

  outr = ((DGTINT *) & out[2]) - 2;
  while (*in == ' ')
    {
      in++;
    }
  out[0] = 1;
  if (*in == '-')
    {
      out[0] = -1;
      in++;
    }
  else if (*in == '+')
    {
      in++;
    }
  while (*in == ' ' || *in == '0')
    {
      in++;
    }
  outexp = 0;
  for (s = in; *s != '\0'; s++)
    {
      if (*s == 'e' || *s == 'E' || *s == 'd' || *s == 'D')
	{
	  if (sscanf (++s, "%d", &outexp) != 1)
	    {
	      outexp = 0;
	    }
	  break;
	}
    }
  if (*in == '.')
    {
      do
	{
	  outexp--;
	  while (*++in == ' ');
	}
      while (*in == '0' && *in != '\0');
    }
  else if (*in != '\0')
    {
      s = in;
      while (*++s == ' ');
      while (*s >= '0' && *s <= '9' && *s != '\0')
	{
	  outexp++;
	  while (*++s == ' ');
	}
    }
  x = outexp / log10_radix;
  outexp_mod = outexp - log10_radix * x;
  if (outexp_mod < 0)
    {
      x--;
      outexp_mod += log10_radix;
    }
  out[1] = x;
  x = 0;
  j = 2;
  for (s = in; *s != '\0'; s++)
    {
      if (*s == '.' || *s == ' ')
	{
	  continue;
	}
      if (*s < '0' || *s > '9')
	{
	  break;
	}
      x = 10 * x + (*s - '0');
      if (--outexp_mod < 0)
	{
	  if (j > n + 1)
	    {
	      break;
	    }
	  outr[j++] = (DGTINT) x;
	  x = 0;
	  outexp_mod = log10_radix - 1;
	}
    }
  while (outexp_mod-- >= 0)
    {
      x *= 10;
    }
  while (j <= n + 1)
    {
      outr[j++] = (DGTINT) x;
      x = 0;
    }
  if (outr[2] == 0)
    {
      out[0] = 0;
      out[1] = 0;
    }
}