Newcoder 143 B.div（数论+高精度）

Description

定义$n$是好的当且仅当存在$x\in [n^2+1,n^2+2n]$满足$x|n^4$，给出一整数$m$，求不小于$m$的最小的好数$n$

Input

一个整数$m(1\le m\le 10^{1000})$

Output

输出不小于$m$的最小的好数$n$

Sample Input

4

Sample Output

6

Solution

假设$n^2+a|n^4$，那么有$n^2+a|(n^4-(n^2+a)(n^2-a))=a^2$，但注意到$a^2\le 4n^2<4(n^2+a)$，故有$a^2=t(n^2+a),t=1,2,3$

$1.t=1$时，则$a(a-1)=n^2$，无解

$2.t=2$时，则$a^2=2n^2+2a$，也即$(a-1{)}^2-2n^2=1$，考虑椭圆曲线$x^2-2y^2=1$，初始解为$(1,0),(3,2)$，故通解为$a_0=0,a_1=2,a_{k+2}=6a_{k+1}-a_k$

$3.t=3$时，则$a^2=3n^2+3a$，由于$3|a(a-3)$，故$3|a$，进而$3|n$，假设$a=3b,n=3m$，那么有$9b^2=27m^2+9b$，也即$(2b-1{)}^2-12n^2=1$，考虑椭圆曲线$x^2-12y^2=1$，初始解为$(1,0),(7,2)$，故通解为$b_0=0,b_1=6,b_{k+2}=14b_{k+1}-b_k$

直接求出两个序列即得到所有的解，$m$比较大用高精度即可

Code

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int maxn=300;
struct BigInt
{
	const static int mod=10000;
	const static int LEN=4;
	int a[maxn],len;
	BigInt() 
	{
		memset(a,0,sizeof(a));
		len=1;
	}
	void init(int x)
	{
		memset(a,0,sizeof(a));
		len=0;
		do
		{
			a[len++]=x%mod;
			x/=mod;
		}while(x);
	}
	void Init(const char s[])
	{
		memset(a,0,sizeof(a));
		int l=strlen(s),res=0;
		len=l/LEN;
		if(l%LEN)len++;
		for(int i=l-1;i>=0;i-=LEN)
		{
			int t=0,k=max(i-LEN+1,0);
			for(int j=k;j<=i;j++)t=t*10+s[j]-'0';
			a[res++]=t;
		}
	}
	int Compare(const BigInt &b)
	{
		if(len<b.len)return -1;
		if(len>b.len)return 1;
		for(int i=len-1;i>=0;i--)
			if(a[i]<b.a[i])return -1;
			else if(a[i]>b.a[i])return 1;
		return 0;
	}
	BigInt operator +(const BigInt &b)const
	{
		BigInt ans;
		ans.len=max(len,b.len);
		for(int i=0;i<=ans.len;i++)ans.a[i]=0;
		for(int i=0;i<ans.len;i++)
		{
			ans.a[i]+=((i<len)?a[i]:0)+((i<b.len)?b.a[i]:0);
			ans.a[i+1]+=ans.a[i]/mod;
			ans.a[i]%=mod;
		}
		if(ans.a[ans.len]>0)ans.len++;
		return ans;
	}
	BigInt operator -(const BigInt &b)const
	{
		BigInt ans;
		ans.len=len;
		int k=0;
		for(int i=0;i<ans.len;i++)
		{
			ans.a[i]=a[i]+k-b.a[i];
			if(ans.a[i]<0)ans.a[i]+=mod,k=-1;
			else k=0;			
		}
		while(ans.a[ans.len-1]==0&&ans.len>1)ans.len--;
		return ans;
	}
	BigInt operator *(const BigInt &b)const
	{
		BigInt ans;
		for(int i=0;i<len;i++)
		{
			int k=0;
			for(int j=0;j<b.len;j++)
			{
				int temp=a[i]*b.a[j]+ans.a[i+j]+k;
				ans.a[i+j]=temp%mod;
				k=temp/mod;
			}
			if(k!=0)ans.a[i+b.len]=k;
		}
		ans.len=len+b.len;
		while(ans.a[ans.len-1]==0&&ans.len>1)ans.len--;
		return ans;
	}
	BigInt operator /(const int &n)const
	{
		BigInt ans;
		ans.len=len;
		int k=0;
		for(int i=ans.len-1;i>=0;i--)
		{
			k=k*mod+a[i];
			ans.a[i]=k/n;
			k=k%n;
		}
		while(ans.a[ans.len-1]==0&&ans.len>1)ans.len--;
		return ans;
	}
	void output()
	{
		printf("%d",a[len-1]);
		for(int i=len-2;i>=0;i--)
			printf("%04d",a[i]);
		printf("\n");
	}
};
char s[1005];
BigInt t6,t14,a[3],b[3],c,d;
int main()
{
    scanf("%s",s);
    c.Init(s);
   	a[0].init(0),a[1].init(2);
   	b[0].init(0),b[1].init(6);
   	t6.init(6);
   	t14.init(14);
   	for(int i=1;;i++)
   	{
   		if(i>1)a[i%3]=a[(i+2)%3]*t6-a[(i+1)%3];
   		if(a[i%3].Compare(c)>=0)
   		{
   			d=a[i%3];
   			break;
   		}
   	}
   	for(int i=1;;i++)
   	{
   		if(i>1)b[i%3]=b[(i+2)%3]*t14-b[(i+1)%3];
   		if(b[i%3].Compare(c)>=0)
   		{
   			if(d.Compare(b[i%3])>0)d=b[i%3];
   			break;
   		}
   	}
   	d.output();
   	return 0;
}