/* ******************************************
Reduce to 1 and find max steps.
This is for the A=0 case (to be run before you can use a 2^k sieve).
Note that the n=2 record of 1 step won't be found due to my 3^2 sieve.
I count (3*n+1)/2 as 1 step.
If you want to count it as 2 steps, there are 4 lines of code
to uncomment or switch to.
Collatz conjecture is the following...
Repeated application of following f eventually reduces any n>1 integer to 1.
f(n) = (3*n+1)/2 if n%2 = 1
f(n) = n/2 if n%2 = 0
For example, let's start at 9...
9 -> 14 -> 7 -> 11 -> 17 -> 26 -> 13 -> 20 -> 10 -> 5 -> 8 -> 4 -> 2 -> 1
You'll need a 64-bit computer with a not-ancient version of gcc
in order for __uint128_t to work. This gives a 128-bit integer!
Compile using something like...
gcc -O3 collatzSieveless_reduceTo1_Aequals0.c
I use __builtin_ctzll(), which should be at least for 64-bit integers.
Note that I use the "long long" function strtoull() when reading in the arguments.
Storage drive not used.
Minimal RAM used to store the 2^k2 lookup table.
Here is a great read for some basic optimizations...
https://en.wikipedia.org/wiki/Collatz_conjecture#Optimizations
Read it, or else the following won't make sense!
I am using the idea of this paper...
http://www.ijnc.org/index.php/ijnc/article/download/135/144
Here is a brilliant algorithm that I use to make k2 lookup table...
https://github.com/xbarin02/collatz/blob/master/doc/ALGORITHM.md
https://rdcu.be/b5nn1
In binary form, break n into A and B, first choosing k>1...
B is the k least significant bits.
A is the rest of the bits.
For k=7...
n = 0bAAAABBBBBBB
For example...
643 = 0b1010000011
gives, for k=7...
A = 0b101 = 5
B = 0b0000011 = 3
5*2^7 + 3 = 643
Breaking numbers into this form is useful for doing multiple steps at a time!
A nice optimization is to not run n%3==2 numbers because n = 3*N + 2 always
follows the already tested 2*N + 1 number (note that 2*N+1 is all odd numbers).
It is important that 2*N + 1 is less than 3*N + 2. This is called a 3^1 sieve.
A better sieve is to also exclude n%9==4 numbers. This is a 3^2 sieve.
Any n = 9*N + 4 number always follows the already tested 8*N + 3.
There are tons of these rules, but almost all exclusions occur with the above rules.
For example, all mod 3^9 rules will block only a few percent more n than a 3^2 sieve,
but you won't gain because checking against these rules takes time.
To have a well defined right bit shift (>>), use an unsigned type for n.
__uint128_t has a max of 2^128 - 1.
I do lots of bitwise operations in this code!
For this paragraph, I will use a double slash, //, for integer division.
(3*n+1) / 2 = 3*(n//2) + 2
since n//2 = (n-1)/2 since n is odd.
It is better to use 3*(n//2) + 2 because you calculate n//2 THEN multiply by 3
allowing n to be twice as large without causing an overflow.
This is a minor improvement worth doing!
(Though David Barina's algorithm no longer requires this calculation!)
If overflow is reached, you may have found a number that
could disprove the conjecture! Overflow must be carefully checked.
If the code reaches the step limit,
you may also have disproved the conjecture!
This code requires an argument...
./a.out task_id
Here is a neat way to run this code...
seq -f %1.0f 0 7 | parallel -P 2 ./a.out |tee -a log.txt &
To pick up where you left off, change the start (and stop) value of seq.
TASK_SIZE should not change between runs.
Change the -P argument of parallel to run more CPU threads at a time!
274133054632352106267 will overflow 128 bits and, for not-crazy k2, should be
noticed by this code!
This number is over 64-bits, so it cannot be made into an integer literal.
To test 128-bit overflow...
task_id = 274133054632352106267 // (1 << TASK_SIZE)
(c) 2021 Bradley Knockel
****************************************** */
#include
#include
#include
#include
struct timeval tv1, tv2;
#define min(a,b) (((a)<(b))?(a):(b))
/*
The following variables are in log2
TASK_SIZE should at least be 10 to give each process something to run!
*/
const int TASK_SIZE = 33;
/*
lookup table size
k2 < 37
*/
const int k2 = 15;
/*
Where to start the search for max steps.
Put your old result here
*/
int stepsMax = 762;
// Prints __uint128_t numbers since printf("%llu\n", x) doesn't work
// since "long long" is only 64-bit in gcc.
// This function works for any non-negative integer less than 128 bits.
void print128(__uint128_t n) {
char a[40] = { '\0' };
char *p = a + 39;
if (n==0) { *--p = (char)('0'); }
else { for (; n != 0; n /= 10) *--p = (char)('0' + n % 10); }
printf("%s\n", p);
fflush(stdout);
}
int main(int argc, char *argv[]) {
if( argc < 2 ) {
printf("Too few arguments. Aborting.\n");
return 0;
}
uint64_t task_id = (uint64_t)strtoull(argv[1], NULL, 10);
printf("task_id = ");
print128(task_id);
printf("TASK_SIZE = ");
print128(TASK_SIZE);
printf(" k2 = %i\n", k2);
fflush(stdout);
int j;
const __uint128_t UINTmax = -1; // trick to get all bits to be 1
// 3^c = 2^(log2(3)*c) = 2^(1.585*c),
// so c=80 is the max to fit in 128-bit numbers.
// Note that c3[0] = 3^0
//const int lenC3 = 81; // for 128-bit numbers
const int lenC3 = 65;
__uint128_t n, nStart;
int nMod;
// calculate lookup table for c3, which is 3^c
__uint128_t* c3 = (__uint128_t*)malloc(lenC3*sizeof(__uint128_t));
c3[0] = 1;
for (j=1; j>= alpha;
L *= c3[alpha];
Salpha += alpha;
if (R == 0) {
L--;
goto next;
}
} while (!(L & 1));
L--;
do {
if ((uint64_t)L == 0) beta = 64; // __builtin_ctzll(0) is undefined
else beta = __builtin_ctzll(L);
beta = min(beta, (size_t)R);
R -= beta;
L >>= beta;
if (L == 1 && reducedTo1 == 0) reducedTo1 = k2 - R;
//if (L == 1 && reducedTo1 == 0) reducedTo1 = k2 - R + Salpha; // if (3*n + 1)/2 is 2 steps
if (R == 0) goto next;
} while (!(L & 1));
} while (1);
next:
/* stores both L and Salpha */
arrayk2[index] = L + ((uint64_t)Salpha << 58);
size_t Ssteps = 0;
if (L == 0) goto finish;
if (reducedTo1 > 0) {
if (index == 1) delayk2[1] = 0;
else delayk2[index] = reducedTo1;
continue;
}
do {
L++;
do {
if ((uint64_t)L == 0) alpha = 64; // __builtin_ctzll(0) is undefined
else alpha = __builtin_ctzll(L);
L >>= alpha;
L *= c3[alpha];
//Salpha += alpha; // if (3*n + 1)/2 is 2 steps
Ssteps += alpha;
} while (!(L & 1));
L--;
do {
if ((uint64_t)L == 0) beta = 64; // __builtin_ctzll(0) is undefined
else beta = __builtin_ctzll(L);
L >>= beta;
Ssteps += beta;
if (L == 1) goto finish;
} while (!(L & 1));
} while (1);
finish:
delayk2[index] = k2 + Ssteps;
//delayk2[index] = k2 + Ssteps + Salpha; // if (3*n + 1)/2 is 2 steps
}
////////////////////////////////////////////////////////////////
//////// test integers that aren't excluded by 3^2 sieve
////////////////////////////////////////////////////////////////
__uint128_t bStart = ( (__uint128_t)1 << TASK_SIZE )*task_id;
__uint128_t bEnd = bStart + ( (__uint128_t)1 << TASK_SIZE );
if (bStart == 0) bStart = 1;
for (nStart = bStart; nStart < bEnd; nStart++) {
// a trick for modNum == 9 to get nStart%9 faster
uint64_t r = 0;
r += (uint64_t)(nStart) & 0xfffffffffffffff;
r += (uint64_t)(nStart >> 60) & 0xfffffffffffffff;
r += (uint64_t)(nStart >> 120);
nMod = r%9;
if (nMod == 2 || nMod == 4 || nMod == 5 || nMod == 8) { // modNum = 9
continue;
}
n = nStart;
int steps = 0;
/* do k2 steps at a time */
while ( (n >> k2) > 0 ) {
size_t index = n & ( ((uint64_t)1<> 58; // just 6 bits gives c
newB &= 0x3ffffffffffffff; // rest of bits gives b
/* find the new n */
//n = (n >> k2)*c3[newC] + newB;
n >>= k2;
if (n > maxNs[newC]) {
printf("Overflow! nStart = ");
print128(nStart);
break;
}
n *= c3[newC];
if (n > UINTmax - newB) {
printf("Overflow! nStart = ");
print128(nStart);
break;
}
n += newB;
steps += k2;
//steps += k2 + (int)newC; // if (3*n + 1)/2 is 2 steps
if (steps > 1000000) { // must be checked so that steps doesn't overflow
printf("Steps limit reached! nStart = ");
print128(nStart);
break;
}
}
// will crash if steps limit was reached or if overflow, but who cares?
steps += delayk2[n];
if (steps > stepsMax) {
stepsMax = steps;
printf(" steps = %i found. nStart = ", steps);
print128(nStart);
}
}
gettimeofday(&tv2, NULL);
printf(" %e seconds\n\n",
(double)(tv2.tv_usec - tv1.tv_usec) / 1000000.0 + (double)(tv2.tv_sec - tv1.tv_sec));
// free memory (cuz why not?)
free(maxNs);
free(c3);
free(arrayk2);
free(delayk2);
return 0;
}