\\ describe a fraction:
\\  - integer part
\\  - transient part
\\  - periodic part
fracdigits(v, base=10) = {
	\\ integer part
	my (int=digits(floor(v), base));

	\\ period length
	my (den=denominator(v), pr=factor(base)[,1]~, p=znorder(Mod(base, den / prod(i=1, #pr, pr[i]^valuation(den,pr[i])))));

	\\ transient part
	my (tra=[], fra=frac(v));
	while (fra != frac(fra*base^p),
		tra = concat(tra, floor(fra*base));
		fra = frac(fra*base);
	);

	\\ periodic part
	my (per = vector(p, x, my (d=floor(fra*base)); fra = frac(fra*base); d));

	return ([int,tra,per]);
}

\\ rebuild a fraction
fromfracdigits(f, base=10) = fromdigits(f[1],base) + ( fromdigits(f[2],base) + fromdigits(f[3],base)/base^#f[3] / (1-1/base^#f[3])) / base^#f[2];

\\ 9-based n-th digits after radix point in fraction f
nth(n, f) = {
	if (n < #f[2],
			f[2][1+n],
			f[3][1+(n-#f[2])%#f[3]]
	);
}

\\ unique rational q such that for any m
\\ floor(2^m/n) AND floor(2^m/k) = floor(2^m*q)
f(n,k) = {
	my (fn=fracdigits(1/n, 2),
	    fk=fracdigits(1/k, 2),
		fq=[0,0,0]);

	fq[1] = binary(bitand(fromdigits(fn[1], 2),
				          fromdigits(fk[1], 2)));
	fq[2] = vector(max(#fn[2], #fk[2]), k, nth(k-1, fn)
										 * nth(k-1, fk));
	fq[3] = vector(lcm(#fn[3], #fk[3]), k, nth(#fq[2]+k-1, fn)
										 * nth(#fq[2]+k-1, fk));

	fromfracdigits(fq, 2);
}

A(n,k) = numerator(f(n,k))

{
	my (m=0);
	for (d=1, 141,
		for (k=1, d,
			print (m++ " " A(d+1-k,k));
		);
	);
}

quit