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

/* In the following I assume this file has been
   renamed olson.c 

   Only tested under linux (gcc), but should work with
   slightly modification on any system.
*/

/* This is a very simple program for A218030
   (for k=2) and similar sequences.

   Saved as olson.c, it make use of the GMP library for
   large numbers which must be installed on your system.
   I compiled it with: gcc -O6 olson.c -lgmp -o olson

   Invoked as : olson L f2 f3 f5 f7
   it finds solutions  N < 2^L to :
   k*N = product of nonzero digits of N^2,
   where k = 2^f2 * 3^f3 * 5^f5 * 7^f7 and the f may be negative.
   For example, for k = 35/8, with N<2^100, the command
   is: olsen 100 -3 0 1 1 

   The values found are saved in a file named 
   Ols-k-L.txt, where k can be of the form x:y.
   Note the values of N are usually not in order,
   so they must be sorted using some other program.

   NOTE: It could be optimized. 
   It is not to be taken as an example of good programming,
   because it is not.
   Since I'm lazy, I did not implement any check on the 
   consistency of the input.

   Giovanni Resta, Oct 20 2012.
 */

void Usage()
{
  printf("\n---------------------------------------\n");
  printf("USAGE: olson L f2 f3 f5 f7\n\n");
  printf("it finds solutions  N < 2^L to :\n");
  printf("k*N = product of nonzero digits of N^2,\n");
  printf("where k = 2^f2 * 3^f3 * 5^f5 * 7^f7\n");
  printf("Example, for k = 35/8 and N<2^100, it is\n");
  printf("olson 100 -3 0 1 1\n\n");
  printf("N are saved in a file names Ols-k-L.txt\n");
  printf("and are not sorted\n");
  printf("---------------------------------------\n\n");
  exit(0);
}

// used to create the file name
void setK(int v, int m, int *num, int *den)
{
  int i;
  if (v>0) for (i=1; i<=v; i++) *num*=m;
  else for (i=-1; i>=v; i--) *den*=m;
}

int main(int argc, char *argv[])
{
  // the search is limited to numbers below 2^Lim
  unsigned Lim;

  // exponents in the factorization of current N
  int e2, e3, e5, e7;

  // exponents in the product of the digit of N^2
  int c2, c3, c5, c7;

  // exponents of 2,3,5 and 7 in the factor K 
  // in the solution of 
  // K*N = product of nonzero digits of N^2
  int f2, f3, f5, f7;

  // z2,z3,z5,zN will contain partial products
  // zLim = 2^Lim is the search limit
  // zQ holds the square of N
  mpz_t z2,z3,z5,zN, zLim, zQ;

  // it holds the actual digits of N^2
  char *S;
  
  // name of the output file and output file
  char fname[64]; 
  FILE *f;

  // working vars for generating the output file name
  int num,den;

  // solutions found until now
  int nsol=0;

  // ================== START ==========

  // not the input I expected, print usage and exit
  if (argc!=6) Usage();

  // input the input
  Lim = atoi(argv[1]);
  f2 = atoi(argv[2]);
  f3 = atoi(argv[3]);
  f5 = atoi(argv[4]);
  f7 = atoi(argv[5]);

  // TODO: (not by me... )sanity check on the input values

  // I build the name of the file
  num=den=1;
  setK(f2,2,&num,&den);  setK(f3,3,&num,&den);
  setK(f5,5,&num,&den);  setK(f7,7,&num,&den);
  if (den==1)
    sprintf(fname,"Ols-%d-%d.txt",num,Lim);
  else
    sprintf(fname,"Ols-%d:%d-%d.txt",num,den,Lim);

  // I reserve a space to write the square digits
  S = malloc(2+(size_t)(0.301029995664*Lim*2));

  // init the GMP variables
  mpz_init(z2); mpz_init(z3); mpz_init(z5);
  mpz_init(zLim); mpz_init(zN);
  mpz_init(zQ);

  // set the search limit as zLim=2^Lim
  mpz_ui_pow_ui(zLim,2,Lim);

  // search all the numbers 2^e2*3^e3*5^e5*7^e7 < Lim
  mpz_set_ui(z2,1);
  e2=0;
  while (e2<Lim)
    {
      // since it could be slow, I print something
      printf("--> 2^%u [found = %d]\n",e2,nsol);
      mpz_set(z3,z2); e3=0;
      while (mpz_cmp(z3,zLim)<0)
	{
	  mpz_set(z5,z3); e5=0;
	  while (mpz_cmp(z5,zLim)<0)
	    {
	      mpz_set(zN,z5); e7=0;
	      while (mpz_cmp(zN,zLim)<0)
		{
		  char *p;
		  c2=c3=c5=c7=0;
		  // Inner loop, now I have my N and I check it
		  //--------------------
		  // zQ = zN^2
		  mpz_mul(zQ,zN,zN);
		  // S = string corresponding to zQ=zN^2
		  mpz_get_str(S,10,zQ);
		  /* I scan the digits of N^2 and
		     I update the exponents of the 
                     factorization of their product */
		  for (p=S; *p; p++)
		    switch (*p)
		      {
		      case '2': c2++; break;
		      case '3': c3++; break;
		      case '4': c2+=2; break;
		      case '5': c5++; break;
		      case '6': c2++; c3++; break;
		      case '7': c7++; break;
		      case '8': c2+=3; break;
		      case '9': c3+=2; break;
		      }
		  // I check if, based of their factors
		  // the product of digits is k*N
		  if (c2==f2+e2 && c3==f3+e3 &&
                      c5==f5+e5 && c7==f7+e7)
		    {
		      // I output the solution found
		      nsol++;
		      mpz_get_str(S,10,zN);
		      printf("%s\n",S);
		      f=fopen(fname,"a");
		      fprintf(f,"%s\n",S);
		      fclose(f);
		    }
		  //--------------------
		  mpz_mul_ui(zN,zN,7); e7++;
		}
	      mpz_mul_ui(z5,z5,5); e5++;
	    }
	  mpz_mul_ui(z3,z3,3); e3++;
	}
      mpz_mul_ui(z2,z2,2); e2++;
    }
  // I add a line to output file to certify I
  // ended the search
  f=fopen(fname,"a");
  fprintf(f,"# Search for n < 2^%d completed\n",Lim);
  fclose(f);
  return 0;
}