/* solu.c -- (C) Tapio Rajala, Dec 2010

   Searches for integers n that are products of small primes
   for which the equation
              x*y*(y-x-1) = n
   has many positive integer solutions (x,y).

     Notice that there are serious restrictions
     due to the use of 64-bit integers !!   

   Compile:
   gcc -Wall -O2 -o solu solu.c

   Usage:
   solu

*/

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

#define REPORT 22

/* Return 1 if p is prime. */
static int isprime(uint32_t p)
{
  uint32_t i, r;

  if (p <= 2 || p%2==0)
    return (p == 2);

  for (i = 3, r = sqrt(p); i <= r; i += 2)
    if (p % i == 0)
      return 0;

  return 1;
}

uint32_t primes[100];
uint32_t n_factors[100];
uint32_t x_factors[100];
uint32_t y_factors[100];
uint32_t top[100];
uint32_t start[100];
uint32_t total_primes, temppi, m, k, record;
uint64_t record_n;
uint32_t solutions[1000];

void test(uint64_t n, uint64_t x)
{
 uint64_t y, t;
 for (m = 1; m < 2; m++)
 {
  t = (long long int)sqrt((long double)((x+m)*(x+m)+4*n/x)+0.5); // This should find at least roughly the right root..
  while (t*t < (x+m)*(x+m)+4*n/x) // but let's find it anyways.
   t++;
  while (t*t > (x+m)*(x+m)+4*n/x)
   t--;
  y = (x+m+t)/2;
  if ((y > 0) && (x*y*(y-x-m) == n))
   {
    solutions[m]++;
   }
 }
}

void set_x(int index, uint64_t n, uint64_t x)
{
  uint64_t x_new;
  uint32_t temp;
  x_new = x;
  for (temp = 0; (temp < n_factors[index]+1) && (x_new < 2147483648); temp++, x_new = x_new*primes[index]) // x < 2^31
  {
   x_factors[index] = temp;
   if (index > 0)
    set_x(index-1, n, x_new);
   else
    test(n, x_new);
  }
}

void set_n(int index, uint64_t n)
{
  uint64_t n_new;
  uint32_t temp;
  n_new = n;   
  for (temp = 0; temp < start[index]; temp++)
   n_new = n_new*primes[index];
  for (temp = start[index]; (temp < top[index]+1) && (n_new < 288230376151711744); temp++, n_new = n_new*primes[index]) // n < 2^58
  {
   n_factors[index] = temp;
   if (index > 0)
   {
     set_n(index-1, n_new);
   }
   else
   {
    for (m = 1; m < 2; m++)
     solutions[m] = 0;
    for (temppi = 0; temppi < total_primes; temppi++)
    {
     printf("%i ",n_factors[temppi]);
     if (n_factors[temppi] < 10)
      printf(" ");
    }
    printf("\r");
   set_x(total_primes, n_new, 1);
    for (m = 1; m < 2; m++)
    {
     if (solutions[m] > REPORT-1) //Print out n:s with at least REPORT solutions.
      printf("Found %i positive solutions for n = %li                                 \n",solutions[m],n_new);   
     if (solutions[m] > record)
     {
      if (solutions[m] > REPORT-1)
       printf(" A new record for this run ! \n");
      record = solutions[m];
      record_n = n_new;
     }
    }
   }
  }
}

int main(int argc, char **argv)
{
  uint64_t temp1, temp2;

  for (temp1 = 0; temp1 < 100; temp1++)
  {
   top[temp1] = 1;
   start[temp1] = 0;
  }

  total_primes = 19;

  //manually set starting powers and maximum powers (if you want)
  start[0] = 3;  top[0] = 12; //2
  start[1] = 2;  top[1] = 8; //3
  start[2] = 1;  top[2] = 5;  //5
  start[3] = 1;  top[3] = 4;  //7
  start[4] = 0;  top[4] = 3;  //11
  start[5] = 0;  top[5] = 2;  //13
  start[6] = 0;  top[6] = 2;  //17
  start[7] = 0;  top[7] = 2;  //19
  start[8] = 0;  top[8] = 2;  //23
  start[9] = 0;  top[9] = 2;  //29
  start[10] = 0; top[10] = 2; //31
  start[11] = 0; top[11] = 2; //37
  start[12] = 0; top[12] = 2; //41


  temp1 = 2;  temp2 = 0;
  while (temp2 < total_primes)
  {
   if (isprime(temp1))
   {
    primes[temp2] = temp1;
    temp2++;
   }
   temp1++;
  }

  record = 0;
  record_n = 0;

  printf("\n");
  printf("          -----   T. Rajala - Dec. 2010   -----\n");
  printf("      Searching solutions for equations of the form\n");
  printf("                     n = x*y*(x-y-1),\n");
  printf("       where n < 2^64 is a product of small primes.\n");
  printf("(Reporting equations with at least %i positive solutions.)\n\n",REPORT);

  set_n(total_primes-1,1);

  printf("The record was %i with n = %li\n",record,record_n);
  return 0;
}