sgu159：Self-Replicating Numbers

题意：
给出B(2<=B<=36)，N(1<=N<=1000)。
求出B进制下长度为N的所有幂自守数(一个K位的幂自守数T就是T^2的后K位和T的K位一样)。
分析：
首先，这样的数肯定不会有很多，我们从最后一位开始构造。
引理：一个K位的幂自守数a1a2a3...ak，满足ak,ak-1ak,...,...,a2a3...ak，也为幂自守数。
根据这个引理，我们一步一步枚举上去。但是高精要注意常数优化。
具体来说：
①：假设已经枚举出了k位幂自守数T,我们枚举第k+1位；
②：设第k+1位为p，则(p*B^k+T)^2的后k+1位与p*B^k+T的k+1位相等；
③：化简上述式子，可得：(p*B^2k+2p*B^k*T+T^2)%B^(k+1)=2*pT%B+T^2%B^(k+1)=p*B^k+T；
④：T是已知的，p枚举就很好判断是否成立，关键是T^2要在时间限制内求出来，所以可以yy一下，用一下平方和公式乱搞什么的...

#include <cstdio>
#include <algorithm>
#include <ctime>
using namespace std;
const char c[40] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
const int MAXN = 2005;
int base, n;
struct bi
{
	int a[MAXN<<1];
	bi operator + (const bi x)
	{
		int lim = max(a[0], x.a[0]);
		bi tmp = {{1}};
		for(int i = 1; i <= lim; ++i)
		{
			tmp.a[i] += a[i]+x.a[i];
			tmp.a[i+1] += tmp.a[i]/base;
			tmp.a[i] %= base;	
		}
		if(tmp.a[lim+1]) tmp.a[0] = lim+1;
		else tmp.a[0] = lim;
		return tmp;
	}
	bi operator * (const int x)
	{
		bi tmp = {{1}};
		for(int i = 1; i <= a[0]; ++i)
		{
			tmp.a[i] += a[i]*x;
			tmp.a[i+1] += tmp.a[i]/base;
			tmp.a[i] %= base;
		}
		if(tmp.a[a[0]+1]) tmp.a[0] = a[0]+1;
		else tmp.a[0] = a[0];
		return tmp;
	}
}ans[10] = {{{1}}}, tmp = {{1}},  as = {{1}}, O = {{1}};
int top = 0;

bi ch(int p)
{
	bi re = O;
	while(p)
	{
		re.a[re.a[0]++] = p%base;
		p /= base;
	}
	if(re.a[0] > 1) re.a[0]--;
	return re;
}

bi mag(bi g, int q)
{
	for(int j = g.a[0]; j >= 1; --j)
	{
		g.a[j+q] = g.a[j];
		g.a[j] = 0;
	}	
	g.a[0] += q;
	return g;
}

bi T, T2, T3;
int t;
void solve()
{
	for(int o = 2; o <= n; ++o)
	{	
		if(o == 3 && tmp.a[1] == 6)
			o = 3;
		int flag = -1;
		for(int i = 0; i < base; ++i)
		{
			t = 2*i*tmp.a[1]+as.a[o];
			if(t%base == i)
			{
				T = ch(i*i);
				T2 = tmp*(2*i);
				T3 = as;
				for(int p = 1; p <= T.a[0]; ++p)
				{
					T3.a[p+2*(o-1)] += T.a[p];
					T3.a[p+2*(o-1)+1] += T3.a[p+2*(o-1)]/base;
					T3.a[p+2*(o-1)] %= base;
				}
				for(int p = 1; p <= T2.a[0]; ++p)
				{
					T3.a[p+(o-1)] += T2.a[p];
					T3.a[p+(o-1)+1] += T3.a[p+(o-1)]/base;
					T3.a[p+(o-1)] %= base;
				}
				if(T3.a[T.a[0]+2*(o-1)+1]) T3.a[0] = T.a[0]+2*(o-1)+1;
				else T3.a[0] = T.a[0]+2*(o-1);
				as = T3;
				tmp.a[o] = i;
				tmp.a[0] = o;
				flag = i;
				break;
			}
		}
		if(flag == -1) return ;
	}
	if(tmp.a[n] || n == 1)
		ans[++top] = tmp;
}

int main()
{
	scanf("%d%d", &base, &n);
	for(int i = 0; i < base; ++i)
		if(i*i%base == i)
		{
			as = tmp = O;
			tmp.a[1] = i;
			as = ch(i*i);
			solve();	
		}
	printf("%d\n", top);
	for(int i = 1; i <= top; ++i)
	{
		for(int j = n; j >= 1; --j)
			printf("%c", c[ans[i].a[j]]);
		printf("\n");	
	}
	return 0;	
}