COME ON CODE ON

A blog about programming and more programming.

Pell Equation

leave a comment »

Pell’s equation is any Diophantine equation of the form
                                     x2 – Dy2 = 1
where D is a nonsquare integer and x and y are integers. Lagrange proved that for any natural number D that is not a perfect square there are x and y > 0 that satisfy Pell’s equation. Moreover, infinitely many such solutions of this equation exist.

There are many ways to get the solutiont to these equations some of which are given on Wikipedia, MathWorld and probably a good explanation here.

Though at first read it may seem difficult to understand the solution. The basic aim here is to get the convergents of Sqrt(D).  Sqrt(D) = [a0,a1,a2,….,ar] such that ar+1=2a0.
Very important point is that if r is even then solution is p2r+1 where pn is the numerator in nth convergent.

The following function is a direct application of the logic given on Mathworld and works for all D<=1000.

double PellSolution(double D)
{
double x1,y1;
double Pn[100],Qn[100],a[100],p[100],q[100];
long long n,r;
if (sqrt(D)!=floor(sqrt(D)))
{
//Initialization
Pn[0]=0;Qn[0]=1;
a[0]=floor(sqrt(D));
p[0]=a[0];q[0]=1;

n=1;
Pn[1]=a[0];Qn[1]=D-a[0]*a[0];
a[1]=floor((a[0]+Pn[1])/Qn[1]);
p[1]=a[0]*a[1]+1.0;
q[1]=a[1];

while (a[n]!=2.0*a[0])
{
n=n+1;
Pn[n]=a[n-1]*Qn[n-1]-Pn[n-1];
Qn[n]=(D-Pn[n]*Pn[n])/Qn[n-1];
a[n]=floor((a[0]+Pn[n])/Qn[n]);
p[n]=a[n]*p[n-1]+p[n-2];
q[n]=a[n]*q[n-1]+q[n-2];
}
r=n-1;
if (r%2==0)
{
for (n=r+2;n<=2*r+1;n++) { Pn[n]=a[n-1]*Qn[n-1]-Pn[n-1]; Qn[n]=(D-Pn[n]*Pn[n])/Qn[n-1]; a[n]=floor((a[0]+Pn[n])/Qn[n]); p[n]=a[n]*p[n-1]+p[n-2]; q[n]=a[n]*q[n-1]+q[n-2]; } x1=p[2*r+1]; y1=q[2*r+1]; } else { x1=p[r]; y1=q[r]; } return x1; } else return 0; }[/sourcecode] Note that the above program only returns the value of x for fundamental solution. To get fundamental y you can return y1. Also to get the nth solution u can use the following recurrence relation x_nth=((x1+y1*Sqrt_D)^nth + (x1-y1*Sqrt_D)^nth)/2 y_nth=((x1+y1*Sqrt_D)^nth - (x1-y1*Sqrt_D)^nth)/(2*Sqrt_D) -fR0D

Advertisements

Written by fR0DDY

April 5, 2009 at 9:36 AM

Posted in Maths, Programming

Tagged with , ,

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: