佩尔方程求解问题

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

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 2typedef 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());           }    }  }