求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;
}