数论模板 - 佩尔方程

佩尔方程

输入一个数k,求方程ans^2 = k*n*n+1 中ans的最小整数解(n是大于等于一的整数,无上限)

C语言模板：

#include <iostream>
#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <math.h>
using namespace std;

int can[1005] = {0};
int a[10005][605]= {0};
int x[6005], y[6005], t[6005];
int h1,h2;
int bb, ee, xx, yy, c, n;

void Pell(int ji,int many,int ma,int kk)
{
    if (ji < kk)
        Pell(ji + 1, a[ma][(ji-1)%a[ma][600]+1], ma, kk);
    else
    {
        h1 = 1;
        h2 = 1;
        x[1] = many;
        y[1] = 1;
        return;
    }
    for (int i = 1; i <= h1; i++)
        t[i] = x[i];
    for (int i = 1; i <= h2; i++)
        x[i] = y[i];
    for (int i = 1; i <= h1; i++)
        y[i] = t[i];
    c = h1;
    h1 = h2;
    h2 = c;
    for (int i = 1; i <= h2; i++)
        if (i <= h1)
            x[i] += many * y[i];
        else
            x[i] = many * y[i];
    if (h2 > h1)
        h1 = h2;
    for (int i = 1; i < h1; i++)
        if (x[i] >= 10)
        {
            x[i+1] += x[i] / 10;
            x[i] %= 10;
        }
    while(x[h1] >= 10)
    {
        x[h1+1] = x[h1] / 10;
        x[h1] %= 10;
        h1++;
    }
    x[0] = h1;
}

void solve()
{
    int i, j;
    for (j = 1; j <= 31; j++)
        can[j*j] = true;
    for(i = 1; i <= 1000; i++)
    {
        if(!can[i])
        {
            a[i][600]=1;
            bb = 1;
            ee = (int)sqrt((double)i);
            a[i][0] = ee;
            ee =- ee;
            xx = bb;
            yy = ee;
            xx =- yy;
            yy = i - yy * yy;
            n=0;
            while((xx - yy) * (xx - yy) < i || xx >= 0)
            {
                xx -= yy;
                n++;
            }
            a[i][1] = n;
            c = xx, xx = yy, yy = c;
            while(xx != bb || yy != ee)
            {
                a[i][600]++;
                c = xx;
                xx =- yy;
                yy = i - yy * yy;
                yy = yy / c;
                n = 0;
                while ((xx - yy) * (xx - yy) < i || xx >= 0)
                {
                    xx -= yy;
                    n++;
                }
                a[i][a[i][600]] = n;
                c = xx, xx = yy, yy = c;
            }
        }
    }
}

int main()
{
    int k;
    solve();
    while(scanf("%d",&k)!=EOF)///ans^2 = k * n ^ 2 + 1,求ans的最小值(k <= 1000)
    {
        if(!can[k])
        {
            if(a[k][600] % 2)
                Pell(1, a[k][0], k, a[k][600]*2);
            else
                Pell(1, a[k][0], k, a[k][600]);
            for (int j = x[0]; j >= 1; j--)
                printf("%d",x[j]);
            printf("\n");
        }
        else
            printf("no solution\n");
    }
    return 0;
}

Java模板（Java处理数据更大）：

import java.io.*;
import java.util.*;
import java.math.*;

public class Main{

    public static BigInteger[] pell(BigInteger n)
    {///求解pell方程，返回两个大数，即一个两个元素的大数数组，分别为x和y
        BigInteger[]p = new BigInteger[10];
        BigInteger[]q = new BigInteger[10];
        BigInteger[]h = new BigInteger[10];
        BigInteger[]g = new BigInteger[10];
        BigInteger[]a = new BigInteger[10];
        p[0] = BigInteger.ZERO;q[0] = BigInteger.ONE;
        p[1] = BigInteger.ONE;q[1] = BigInteger.ZERO;
        a[0] = sqrt(n);
        a[2] = a[0];
        g[1] = BigInteger.ZERO;h[1] = BigInteger.ONE;
        while(true){
            g[2] = a[2].multiply(h[1]).subtract(g[1]);
            h[2] = n.subtract(g[2].multiply(g[2])).divide(h[1]);
            a[3] = g[2].add(a[0]).divide(h[2]);
            p[2] = a[2].multiply(p[1]).add(p[0]);
            q[2] = a[2].multiply(q[1]).add(q[0]);
            if(p[2].multiply(p[2]).subtract(n.multiply(q[2].multiply(q[2]))).equals(BigInteger.ONE)==true
               &&p[2].equals(BigInteger.ZERO)==false&&q[2].equals(BigInteger.ZERO)==false)
            return new BigInteger[] {p[2],q[2]};
            g[0] = g[1];h[0] = h[1];
            g[1] = g[2];h[1] = h[2];
            a[2] = a[3];
            p[0] = p[1];p[1] = p[2];
            q[0] = q[1];q[1] = q[2];
        }

    }

    static BigInteger sqrt(BigInteger n)
    {///对于大数开根号
        BigInteger low = BigInteger.ZERO,high = n,ans = null;
        while(low.compareTo(high)<=0){
            BigInteger mid = low.add(high).shiftRight(1);
            if(mid.multiply(mid).compareTo(n)<=0){
                ans = mid;
                low = mid.add(BigInteger.ONE);
            }
            else high = mid.subtract(BigInteger.ONE);
        }
        return ans;
    }
    public static void main(String[]args) throws IOException
    {
        /*
        System.setIn(new BufferedInputStream(new FileInputStream("in.txt")));///文件重定向输入
        System.setOut(new PrintStream(new FileOutputStream("out.txt")));   ///文件重定向输出
        */
        Scanner cin = new Scanner(new BufferedInputStream(System.in));
        while(cin.hasNext()){
            int n = cin.nextInt(),ok = 1;
            for(int i = 1;i<=31;++i)
                if(i*i==n){ok = 0;break;}
            if(ok==0){System.out.println(-1);continue;}
            BigInteger r[] = pell(BigInteger.valueOf(n));
            System.out.println(r[0]);
        }
    }
}