佩尔方程求解问题

关于佩尔方程，下面推荐几个不错的博客：

http://www.matrix67.com/blog/archives/5556

http://blog.sina.com.cn/s/blog_5d06e2390100ll92.html

http://m.blog.csdn.net/blog/wh2124335/8871535

HDU3292  No more tricks, Mr Nanguo    

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=3292

题目大意：讲述的是滥竽充数的故事。有一个乐队，可以站成一个边长为x的正方形，也可以减去一个人，然后站成若干个边长为y的正方形，很明显乐队人数不是唯一的，现在让你找出第k大的乐队人数。

分析：根据题意可以提炼出二元不定方程x^2-D*y^2=1，典型的佩尔方程。首先我们知道，若D是完全平方数，那么方程无解；否则，可以暴力（也可以用连分数）求出方程的特解，然后矩阵快速幂求出方程第k大的解。

这题我用暴力求的特解，实现代码如下：

#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
#define mod 8191
#define N 2
typedef struct
{
    int m[N][N];
}matrix;
matrix per,d;
int n;
int x,y,D;
void search()
{
    y=1;
    while(true)
    {
        x=(int)sqrt(D*y*y+1.0);
        if(x*x-D*y*y==1) break;
        y++;
    }
}
void init()
{
    d.m[0][0]=x%mod;
    d.m[0][1]=D*y%mod;
    d.m[1][0]=y%mod;
    d.m[1][1]=x%mod;
    for(int i=0;i<N;i++)
      for(int j=0;j<N;j++)
        per.m[i][j]=(i==j);
}
matrix multi(matrix a,matrix b)
{
    matrix c;
    for(int i=0;i<N;i++)
      for(int j=0;j<N;j++)
      {
          c.m[i][j]=0;
          for(int k=0;k<N;k++)
            c.m[i][j]+=a.m[i][k]*b.m[k][j];
        c.m[i][j]%=mod;
      }
      return c;
}
matrix quick_mod(int k)
{
    matrix p=d,ans=per;
    while(k)
    {
        if(k&1)
        {
            ans=multi(ans,p);
            k--;
        }
        k>>=1;
        p=multi(p,p);
    }
    return ans;
}
int main()
{
    int K;
    while(cin>>D>>K)
    {
        int ad=(int)sqrt(D+0.0);
        if(ad*ad==D)
        {
            puts("No answers can meet such conditions");
            continue;
        }
        search();
        init();
        d=quick_mod(K-1);
        printf("%d\n",(d.m[0][0]*x%mod+d.m[0][1]*y%mod)%mod);
    }
    return 0;
}

POJ1320  Street Numbers

题目连接：http://poj.org/problem?id=1320

题目大意：一个程序员住在一条街道上，街道上的房屋按1到m编号，程序员出家门（假设程序员家的编号为n）左转，把每个房屋的编号（不包括自己家）加起来的值，等于出家门右转，把每个房屋（同样不包括自己家）的编号加起来的值，让你找出前10组符合这样n和m的值。

分析：题目描述可以转换为，求解两个不相等的正整数n和m（m>n），使得1+2+...+(n-1)=(n+1)+...+m，两边都为等差数列，化简一下可以得到：(n)×(n-1)=（m-n）×(m+n+1)，进一步化简得：(2m+1)^2-8*n^2=1，即变为了佩尔方程的求解问题。对于方程x^2-8*y^2=1，很容易看出特解为x=3,y=1，那么由递推关系式xn=xn-1*x1+d*yn-1*y1;yn=xn-1*y1+yn-1*x1，得到xn+1=3*xn+8*yn,yn+1=xn+3*yn；

实现代码如下：

#include <cstdio>
using namespace std;
int main()
{
    int x,y;
    int x0=3,y0=1;
    for(int i=1;i<=10;i++)
    {
        x=3*x0+8*y0;
        y=x0+3*y0;
        printf("%10d%10d\n",y,(x-1)/2);
        x0=x,y0=y;
    }
    return 0;
}

poj2427  Smith's Problem

题目连接：http://poj.org/problem?id=2427

题目大意：求方程x^2-n*y^2=1的最小正整数解。

分析：连分数求佩尔方程特解。由于要用高精度，用java。

实现代码如下：

import java.math.BigInteger;  
import java.util.Scanner;  

public class Main
{  
    public static void solve(int n) 
    {  
        BigInteger N, p1, p2, q1, q2, a0, a1, a2, g1, g2, h1, h2,p,q;  
        g1 = q2 = p1 = BigInteger.ZERO;  
        h1 = q1 = p2 = BigInteger.ONE;  
        a0 = a1 = BigInteger.valueOf((int)Math.sqrt(1.0*n));
        BigInteger ans=a0.multiply(a0);
        if(ans.equals(BigInteger.valueOf(n)))
        {
        	System.out.println("No solution!");
        	return;
        }
        N = BigInteger.valueOf(n);  
        while (true) 
        {  
            g2 = a1.multiply(h1).subtract(g1);     
            h2 = N.subtract(g2.pow(2)).divide(h1);  
            a2 = g2.add(a0).divide(h2);          
            p = a1.multiply(p2).add(p1);           
            q = a1.multiply(q2).add(q1);          
            if (p.pow(2).subtract(N.multiply(q.pow(2))).compareTo(BigInteger.ONE) == 0) break;
            g1 = g2;h1 = h2;a1 = a2;  
            p1 = p2;p2 = p;  
            q1 = q2;q2 = q;  
        }
        System.out.println(p+" "+q);
    }  
   
    public static void main(String[] args) 
    {  
        Scanner cin = new Scanner(System.in);  
        while(cin.hasNextInt())
        {
            solve(cin.nextInt());   
        }
    }  
}