求Pell方程解

#include <iostream>#include <cstdio>#include <cmath>using namespace std;struct PellAns{    int p,q;};struct Node{    int g,h;};PellAns Solve(int n){    PellAns s[4];       //迭代过程中的答案    Node w[4];    int a[4];    s[0].p=0;s[0].q=1;    s[1].p=1;s[1].q=0;    a[0]=(int)floor(sqrt((double)n));       //拿到 √d    a[2]=a[0];    w[1].g=0;w[1].h=1;    while(1)        //求解迭代的过程    {        w[2].g=-w[1].g+a[2]*w[1].h;        w[2].h=(n-w[2].g*w[2].g)/w[1].h;        a[3]=(int)floor((double)(w[2].g+a[0])/w[2].h);        s[2].p=a[2]*s[1].p+s[0].p;        s[2].q=a[2]*s[1].q+s[0].q;        if((s[2].p*s[2].p-n*s[2].q*s[2].q)==1&&s[2].p>0&&s[2].q>0)      //找到了正整数的解            return s[2];        //保留结果        w[0]=w[1];w[1]=w[2];        a[2]=a[3];        s[0]=s[1];s[1]=s[2];    }}int main(){    int n;    PellAns ans;    while(scanf("%d",&n)!=EOF)       //n为完全平方数时没有非平凡解    {        ans = Solve(n);        printf("%d\t%d\n",ans.p,ans.q);    }    return 0;}