Pell佩尔方程漫谈

posted at 2022.7.1 13:04 by administrator

Pell佩尔方程有着悠久迷人的历史,它首次被记载在“阿基米德牛群问题”中。该问题涉及8种花色的牛,各种线性关系已经给出,要求读者确定各种花色的牛的数量。经过一系列计算,这一问题最终简化成求解佩尔方程:

x2-4729494y2=1

最小解由Amthor最先于1880年确定,其中的y是一个41位数,看起来阿基米德及其同时代的人是不可能确定这个解的,但能想到这样一个问题是令人神魂颠倒的。

现代欧洲关于佩尔方程的历史始于1657年,当时费马向他的数学家朋友提出求解方程x2-61y2=1的挑战,其中几位数学家找到了最小解(x,y)=(1766319049,226153980)1657年,布朗克尔(William Brouncker)描述了求解佩尔方程的一般方法,布朗克尔为了说明他的方法的有效性,仅用二个小时就求出方程:x2-313y2=1的一个解(32188120829134849,1819380158564160)

沃利斯(J.Wallis)在一本关于代数与数论的书中描述了布朗克尔的方法,并且沃利斯和费马都断言佩尔方程总是有解的。欧拉错误地认为沃利斯书中的方法属于另一位英国数学家佩尔,并且正是欧拉将这个方程称为人们目前所熟知的“佩尔方程”。这个误解使佩尔获得了不朽的数学名声!

那么,如何求解佩尔方程

           X2-Dy2=1

的一个正整数解xy?因式分解

( x – y*D0.5)( x + y*D0.5) = 1

1表成两个数的乘积,其中一个相当大,更精确地说,x + y*D0.5相当大,尤其在xy很大的时候更是如此,因此另一个因子

 x - y*D0.5= 1 / ( x + y*D0.5)

必定相当小。

那么x-y*D0.5能小到何种程度?即|x-y*D0.5|<?。我们使用方法--鸽笼原理来求出结果。

首先认识“鸽笼原理”(抽屉原理或Schubfachschluβ

这个奇妙的原理是说,如果鸽子比鸽笼多,那么至少有一只笼子里有两只以上的鸽子!尽管该原理看上去是明显和平凡的,但是恰当应用能产生丰富的数学成果。

我们将要做的是寻找两个不同的倍数x*D0.5y*D0.5使它们的差非常接近一个整数。为此,选择某个大数Y并考虑如下所有的倍数:

0*D0.51*D0.52*D0.53*D0.5Y*D0.5

将每一个倍数写成一个整数与一个介于0和他之间的小数之和。

0*D0.5=N0+F0   其中N0=0,F0=0

1*D0.5=N1+F1   其中N1是整数,0F1 <1

     2*D0.5=N2+F2   其中N2是整数,0F2<1

3*D0.5=N3+F3   其中N3是整数,0F13<1

       

Y*D0.5=Ny+Fy   其中Ny是整数,0F1y<1

我们这里的鸽子就是Y+1个数F0F1Fy。所有的鸽子都在01之间,所以它们全部落在区间0t<1中。我们将这个区间等分成Y个鸽笼,也就是说,将如下的区间取作鸽笼:

     鸽笼10/Yt <1/Y

鸽笼21/Yt <2/Y

鸽笼32/Yt <3/Y

鸽笼Y(Y-1)/Yt <Y/Y

每只鸽子都栖息在一个鸽笼中,且鸽子的数目多于鸽笼的数目,因此,鸽笼原理确保了某个鸽笼中至少有两只鸽子。比如,鸽子FmFn栖息在同一个鸽笼中。将这两只鸽子做上标记使得m<n。注意到鸽笼相当狭窄,长度仅有1/Y,所以FmFn之间的距离小于1/Y。用数学术语表示出来就是:

| Fm - Fn | < 1/Y

下面利用

m*D0.5=Nm+Fm  n*D0.5=Nn+Fn 

将上述不等式改写成

|m*D0.5- Nm-n*D0.5- Nn| < 1/Y

即:

|Nn - Nm-n-m*D0.5 | < 1/Y

注意到Nn - Nmn-m都是(正)整数。如果分别将它们记作xy那么就完成了使|x-y*D0.5|相当小的目标。

通过试算x y*D0.5无限逼近1/y值,我们就可以求出x,y值。 

C#程序代码:

         private static string Pell(double D)
        {
            DateTime dt = DateTime.Now;
            double m = Math.Pow(D, 0.5);
            string result = "1";
            double x=1,y=1,z;
            for (x = 1; x < 100000; x++)
            {
                for (y = 1; y < 100000; y++)
                {
                    if (Math.Abs( x - m * y) < 1 / y)
                    {    
                        z=Math.Pow(x,2)-Math.Pow(y,2)*D;
                        result += "("+x.ToString()+","+y.ToString()+","+z.ToString()+")";   
                    }
                    else
                    {  }                  
                }
            }
            DateTime dt2 = DateTime.Now;
            TimeSpan timespan = dt2 - dt;
            Console.WriteLine("耗时{0}秒", timespan.TotalSeconds);
            return result;
        }
 

我们解决了佩尔方程的基本求解问题,但是,大家有没有注意到,运用上述C#程序求解佩尔方程解时,其运行时间为ON2非常缓慢。例如:求解x2-61y2=1,耗时54秒之多。

有没有快速求解程序呢?回答有。它是采用连分数数学特性求解佩尔方程。

首先我们理解连分数与佩尔方程定理:

D0.5的连分数为

     D0.5= [ab1b2…,bm-1bm]

并令p/q=[ ab1b2…,bm-1],则佩尔方程

         x2 - Dy2 = 1

的最小正整数解为

       (x1y1)等于

                     1、若m为偶数,为pq

                     2m为奇数p2+q2D2pq

其他所有的解均由公式

                 xk+yk*D0.5=(x1+y1*D0.5)kk=123

给出。

下面举例说明求解x2-313y2=1

3130.5=[17,1,2,4,11,1,1,3,2,2,3,1,1,11,4,2,1,34]

按定理所述,去掉周期部分的最后一个数字,即34,并计算分数

126862368/7170685 = [17,1,2,4,11,1,1,3,2,2,3,1,1,11,4,2,1]

周期m等于17,故p,q=(1268623687170685)满足

       1268623682-313·71706852=-1

为求出佩尔方程的最小解,根据定理,计算

         P2+q2D=1268623682+313·71706852=32188120829134849

         2pq=2·126862368·7170685=1819380158564160

于是x2-313y2=1的最小解为

   (x,y)=( 32188120829134849, 1819380158564160)

如果要求下一个最小解,只需对

  32188120829134849+1819380158564160·3130.5

平方即可得答案

      x,y=(2072150245021969438104715652505601,

             117124856755987405647781716823680)

 

C#程序代码:

             private static string Lianfenshu(int D)  //x^2-D*y^2=1中的D

          {
            int[] d = LianfenshuQuery(D);
            //int[] d= LianfenshuQueryDictionary(D);
            double[] zum=new double[d.Length+1];
            zum[d.Length] = 0;
            zum[d.Length-1]=1;
            double sum =1;         
            string result = "";
            int x ;
            for (x = d.Length-2; x >=0; x--)
            {                  
                    sum=d[x]*sum+zum[x+2];
                    zum[x] = sum;
 
                result += "(" + sum.ToString() + "," + x.ToString() + ")";
 
            }
            if (d.Length%2!=0)
            {
                Console.WriteLine("最小解(" + zum[0].ToString() + ","+ zum[1].ToString() + ")") ;
 
            }
            else
            {
                double s= Math.Pow(zum[0], 2)+ Math.Pow(zum[1], 2)*D;
                double t = 2* zum[0]* zum[1];
                Console.WriteLine("最小解("+s.ToString()+","+ t.ToString()+")");
            }
           
            return result;
        }
        private static int[] LianfenshuQuery(int D)
        {
            int[][] lianfenshu=new int[360][];
            lianfenshu[2]=new int[]{ 1, 2 };
            lianfenshu[3] = new int[] { 1, 1, 2 };
            lianfenshu[5] = new int[] { 2, 4 };
            lianfenshu[6] = new int[] { 2, 2, 4 };
            lianfenshu[7] = new int[] { 2, 1, 1, 1, 4 };
            lianfenshu[8] = new int[] { 2, 1, 4 };
            lianfenshu[10] = new int[] { 3, 6};
            lianfenshu[11] = new int[] { 3, 3, 6 };
            lianfenshu[12] = new int[] { 3, 2, 6 };
            lianfenshu[13] = new int[] { 3, 1, 1, 1, 1, 6 };
            lianfenshu[14] = new int[] { 3, 1, 2, 1, 6 };
            lianfenshu[15] = new int[] { 3, 1, 6 };
            lianfenshu[17] = new int[] { 4, 8 };
            lianfenshu[18] = new int[] { 4, 4, 8 };
            lianfenshu[19] = new int[] { 4, 2, 1, 3, 1, 2, 8 };
            lianfenshu[20] = new int[] { 4, 2, 8 };
            lianfenshu[55] = new int[] { 7, 2, 2, 2, 14 };
            lianfenshu[61] = new int[] { 7, 1, 4, 3, 1, 2, 2, 1, 3, 4, 1, 14 };
            lianfenshu[71] = new int[] { 8, 2, 2, 1, 7, 1, 2, 2, 16 };
            lianfenshu[313] = new int[] { 17, 1, 2, 4, 11, 1, 1, 3, 2, 2, 3, 1, 1, 11, 4, 2, 1, 34 };  
            lianfenshu[4] = new int[] { 3,7,15,1,292,1,1,1,2,1,3,1,14,2,1,1,2 };                                    //Math.PI
            lianfenshu[9] = new int[] { 2,1,2,1,1,4,1,1,6,1,1,8,1,1,10,1,1,12,1,1,14 };                               //e
            int[] result = lianfenshu[D];
            return result;
         }

          效果图:

Tags: , , ,

IT技术

Add comment

  Country flag

biuquote
  • Comment
  • Preview
Loading