/* 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 #include #include #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; }