#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include <string.h>
#include <math.h>

#include "../bignum.h"

#ifndef lint
static char sccsid[] = "@(#)factor.c	4.1 (Wollongong) 6/13/83";
#endif

/*
 *              factor [ number ]
 *
 * Written to replace factor.s in Bell V7 distribution
 */
/* hacked by jrm to use my bcd bignum routines 5/1996 */

void printfactors(BigNum *n);
BigNum *factor(BigNum *n, BigNum *f);

#ifndef CGI_SCRIPT

int main (int argc, char *argv[])
{
int i, err;
BigNum	n;
char str[ISIZE+8];

if (argc >= 2)
  {
  err = ascii2bcd(argv[1],&n);
  if ( err )
    {
    fprintf(stderr,"%s: Could not get a number from command line\n",argv[0]);
    return 1;
    }
  else
    printfactors(&n);
  }
#ifndef DONT_USE_STDIN
else
  {
  while (fgets(str,sizeof(str),stdin) != NULL)
    {
    if (str[strlen(str)-1] == '\n' )
      str[strlen(str)-1] = '\0';
    err = ascii2bcd(str,&n);
    if ( err )
      {
      fprintf(stderr,"%s: Could not get a number %d\n",argv[0],err);
      return 1;
      }
    printfactors(&n);
    printf("\n");
    }
  }
#else
else
  fprintf(stderr,"Usage: %s integer\n",argv[0]);
#endif
}
#endif /* CGI_SCRIPT */

/*
 * Print all prime factors of integer n > 0, smallest first, one to a line
 */
void printfactors(BigNum *n)
{
BigNum *prime;
BigNum p, tmp, r, q;
char str[ISIZE+8];

if ( /* n == 1 */ n->c == 1 && n->x[0] == 1 )
  {
#ifndef CGI_SCRIPT
  printf("\t1\n");
#else
  printf("1<br>\n");
#endif
  }
else
while ( /* n != 1 */ n->c > 1 || n->x[0] != 1)
  {
  int2bcd(0,&p);
  prime = factor(n,&p);
  bcd2ascii(str, sizeof(str), prime);
#ifndef CGI_SCRIPT
  fprintf(stdout,"  %s\n",str);
#else
  fprintf(stdout,"  %s<br>\n",str);
  fflush(stdout);
#endif
  /* n /= prime; */
  int2bcd(0,&q);
  int2bcd(0,&r);
  bcddiv(n,prime,&r,&q);
  *n = q;
  }
}

/*
 * Return smallest prime factor of integer N > 0
 *
 * Algorithm from E.W. Dijkstra (A Discipline of Programming, Chapter 20)
 */

struct AR {
	int	hib;
	BigNum	val[ISIZE];
} ar;

BigNum *factor(BigNum *N, BigNum *retval)
{
register int k;
int i,j;
BigNum	*p;
BigNum  f;

BigNum x, y, r, q, t;

p = retval;

ar.hib = -1;
x = *N;
/* y = 2; */
int2bcd(2, &y);

switch ( N->x[0] )
  {
  case 0:
    if( N->c == 0 )
      {
      *p = *N;
      return (p);
      }
  case 2:
  case 4:
  case 6:
  case 8:
    *p = y;
    return (p); 	/* 2 is a factor */
  default:
    break;
  }
for(i=0, j=0; i < N->c; i++)
  {
  j += N->x[i];
  }
if ( (j % 3) == 0 )
  {
  int2bcd(3, p);
  return (p);		/* 3 is a factor */
  }
if ( N->x[0] == 5 )
  {
  int2bcd(5, p);
  return (p);		/* 5 is a factor */
  }

while ( /* x != 0 */ x.c >=1 || x.x[0] != 0 )
  {
  /* ar.val[++ar.hib] = x % y; */
  /* x /= y; */
  bcddiv(&x,&y,&r,&q);
  ar.hib++;
  ar.val[ar.hib] = r;
  if ( ar.hib >= ISIZE ) fprintf(stderr,"\n\n  array overflow  \n\n");
  x = q;
  /* y += 1; */
  bcdinc(y.x,y.c);
  for ( y.c = 0, k = ISIZE -1 ; k >= 0; k-- )
    if ( y.x[k] != 0 )
      {
      y.c = k+1;
      break;
      }
  }



/* f = 2 */
int2bcd(2, &f);

while ( (ar.val[0].c >= 1 || ar.val[0].x[0] != 0) && ar.hib > 1)
  {
  /* f += 1; */
  bcdinc(f.x,f.c);
  for ( f.c = 0, k = ISIZE -1 ; k >= 0; k-- )
    {
    if ( f.x[k] != 0 )
      {
      f.c = k+1;
      break;
      }
    }
/*  printf("f.x[0] = %d, f.c=%d. f.isign= %d\n",f.x[0],f.c,f.isign); */
  i = 0;
/*  while (i != ar.hib) */
  while (i < ar.hib)
    {
    j = i + 1;
    /* ar.val[i] -= j * ar.val[j]; */
    int2bcd(j, &r);
    int2bcd(0, &q);
    if(ar.val[j].c > 0 && r.c > 0)
      {
      bcdmul(ar.val[j].x, r.x, q.x, ar.val[j].c, r.c); /* no multiply if 0 */
/*      printf("ar.val[j].isign = %d, r.isign = %d \n",ar.val[j].isign,r.isign); */
      if ( ar.val[j].isign == r.isign )
        q.isign = 1;
      else
        q.isign = -1;
      for ( q.c = 0, k = ISIZE -1 ; k >= 0; k-- )
        {
        if ( q.x[k] != 0 )
          {
          q.c = k+1;
          break;
          }
        }
      }
    r = ar.val[i];
    q.isign *= -1;
    addbignum(&r,&q,&(ar.val[i]));

    while (ar.val[i].isign < 0)
      {
      int2bcd(0,&t);
      /* ar.val[i] += f + i; */
      r = ar.val[i];
      addbignum(&r,&f,&t);
      int2bcd(i,&r);
      addbignum(&t,&r,&(ar.val[i]));
      int2bcd(1,&r);
      r.isign = -1;
      /* ar.val[j] -= 1; */
      t = ar.val[j];
      addbignum(&t,&r,&(ar.val[j]));
      }
    i = j;
    }
/*  printf(" ar.val[ar.hib].c=%d, ar.val[%d].x[0]=%d\n",ar.val[ar.hib].c,ar.hib,ar.val[ar.hib].x[0]); */
  while ( /* ar.val[ar.hib] == 0 */
	ar.val[ar.hib].c <= 1 && ar.val[ar.hib].x[0] == 0)
    {
    ar.hib--;
/*    printf("ar.hib=%d\n",ar.hib); */
    }
  }

if ( /* ar.val[0] == 0 */ ar.val[0].c <= 1 && ar.val[0].x[0] == 0 )
  *p = f;
else
  *p = *N;

return(p);
}