{
\\ This program splits the calculation of the Collatz sequence into smaller batches for faster computation. The time complexity is in practice about O(log(n)^1.5) for an input number n, compared to approximately O(log(n)^2) for a simple sequential algorithm, meaning if the number of digits of n increases 10-fold, the computation time increases about 30-fold (instead of 100-fold). We are also saving about 15% by starting the computation at 3^g[h]-1 instead of 2^g[h]-1, as any Mersenne number 2^p-1 reaches 3^p-1 after 2p steps.

\\ Known Mersenne prime exponents:
g=[2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689, 9941, 11213, 19937, 21701, 23209, 44497, 86243, 110503, 132049, 216091, 756839, 859433, 1257787, 1398269, 2976221, 3021377, 6972593, 13466917, 20996011, 24036583, 25964951, 30402457, 32582657, 37156667, 42643801, 43112609, 57885161, 74207281, 77232917, 82589933];

\\ The data for computing 16 steps in a row is stored in a vector. Here, one step always includes the halving step (so each n/2 and (3n+1)/2 is counted as one step respectively).
s=16;
t=2^s;
a=vector(t);
b=a;
for(n=0, t-1,
	k=n;
	j=0;
	for(i=1, s,
		l=bittest(k,0);
		k=(k*3^l+l)/2;
		j+=l
	);
	a[n+1]=k;
	b[n+1]=j
);

\\ Now computing the Collatz sequence for the given Mersenne prime exponents g[h]. Relevant for output: i = number of all steps, j = number of odd steps.
for(h=1, #g,
	n=(3^g[h]-1);
	j=g[h];
	i=j;
	while(n>1,
		u=max(1024,2^floor(1+log(log(n)))); \\ The Collatz sequence for the rightmost u bits of n is computed seperately, and then applied to the whole number.
		v=2^u;
		if(n>v,
			m=n%v;
			n>>=u;
			y=0;
			for(k=1, u/s,
				l=m%t;
				m>>=s;
				m=m*3^b[l+1]+a[l+1];
				y+=b[l+1]
			);
			n=n*3^y+m;
			i+=u;
			j+=y;
			print1("2^"g[h]"-1... step: "i,"  #digits: "floor(1+log(n)/log(10))"   ",Strchr(13)),
			if(n>t,
				l=n%t;
				n>>=s;
				n=n*3^b[l+1]+a[l+1];
				i+=s;
				j+=b[l+1]
			,
				l=bittest(n,0);
				n=(n*3^l+l)/2;
				i++;
				j+=l
			)
		)
	);
	print("2^"g[h]"-1: "j" odd steps, "i+j" steps total") \\ For the output, the odd steps don't include the halving step and are counted seperately.
);
}