2.1 Example

• One of Pell’s equations is to find pairs \((p,q)\) satisfying \(p^2-8q^2=1\). More details about Pell’s equation can be found at https://en.wikipedia.org/wiki/Pell%27s_equation.

3.1 Method 1: Avoid square

• Define the difference as a fuction \(f(p,q)=p^2-8q^2-1\).
• Codes are provided below.

int pf8int1(int pair1[], int pair2[], int max_iter){    
    int p=1,q=1,count=0,times=0,f0=p*p-8*q*q-1;
    while(count<max_iter){
        times++;
        if(f0<0){f0=f0+4*(p+1);p+=2;}
        else if(f0>0){f0=f0-8*(2*q+1);q++;}
        else{
            pair1[count]=p;
            pair2[count]=q;
            count++;            
            f0=4*(p+1)-8*(2*q+1);
            p+=2;q++;
        }       
    }
    return times;
}

3.2 Method 2: Transfer problem to low order \(\sqrt{p}\).

• Let \(x=(p-1)/2\). Pell’s equation becomes \(x(x+1)/2=q^2\). It is obvious that \(x/2\) and \(x+1\) are square numbers if \(x\) is even, and \(x\) and \((x+1)/2\) are square numbers if \(x\) is odd. Thus, we can reduce the iteration to approximate \(\sqrt{p}\). This idea is provided by a student.
• Codes are provided below.

int ifsquare(int x) {
    // if x is a square number, return its square root; else return 0;
    int i = (int) sqrt((float) x);
    int isquare = i * i;
    while (isquare < x) { isquare = isquare + 2 * i + 1; i++; }
    if (isquare == x) return i;
    else return 0;
}

int pf8new(int pair1[], int pair2[], int max_iter) {
    int t, p, q, tsquare, i, indicator,times=0;
    i = 0, t = 1, tsquare = t * t;
    while (i < max_iter) {
        times++;
        // tsquare = t*t; test if (tsquare+1)/2 is a square number;
        indicator = ifsquare( (tsquare + 1) / 2 );
        if (indicator)
        {
            p = tsquare * 2 + 1;
            q = t * indicator;
            pair1[i]=p;
            pair2[i]=q;
            i++;
        }
        else {
            // if we define x = tsquare - 1, then x + 1 is a square number.
            // test if x/2 (which equal to (tsquare - 1)/2) is a square number;
            indicator = ifsquare( (tsquare - 1) / 2 );
            if (indicator)
            {
                p = tsquare * 2 - 1; // p = x * 2 + 1;
                q = t * indicator;
                pair1[i]=p;
                pair2[i]=q;
                i++;
            }
        }
        tsquare = tsquare + t * 4 + 4; // set tsquare = (t + 2)*(t + 2);
        t = t + 2;
    }
    return times;
}

3.3 Method 3: Recursion formula

• We consider a general Pell’s equation \[\begin{align} p^2-Nq^2=1, \end{align}\] where \(p\) and \(q\) are positive integers, and \(N\) is positive nonsquare integer. Let pair \((p_0,q_0)\) be initial. The recursion formula can be given as following. \[\begin{align*} \binom{p_{n+1}}{q_{n+1}}=\left( \begin{array}{cc} p_0 & Nq_0 \\ q_0 & p_0 \\ \end{array} \right)\binom{p_n}{q_n}, \end{align*}\] where \((p_n,q_n)\) is the \(n^{th}\) pair solution to Pell’s equation.
• Here \((p_0,q_0)=(3,1)\) if \(N=8\). Codes are provided below.

void pf8iter(int pair1[], int pair2[], int max_iter){
    int p=3,q=1,count=0,tmp1;
    pair1[count]=p;
    pair2[count++]=q;
    while(count<max_iter){
        tmp1 = 3*p+8*q;
        q = 3*q+p;
        p=tmp1;
        pair1[count]=p;
        pair2[count++]=q;
    }
}

3.4 Main function calling above subfunctions

//run on VS 2010
#include<iostream>
int main(){
    int i,j;
    int pair1[10],pair2[10];
    for(i=0;i<10;i++)pair1[i]=pair2[i]=0;
    
    int times = pf8new(pair1,pair2,10);
    //int times = pf8int1(pair1,pair2,10);
    
    for(int i=0;i<10;i++)
        printf("(p,q)=(%d,%d)\n",(int)pair1[i],(int)pair2[i]);
    printf("times=%d\n",times);

    pf8iter(pair1,pair2,10);
    for(int i=0;i<13;i++)
    printf("(p,q)=(%d,%d)\n",pair1[i],pair2[i]);
    
    return 0;
}

3.5 Drawback

• The drawback of above three methods is that we only can output the first 10 pairs. Outputting more pairs needs arrays to save big numbers. We will bring C codes with array to compute more pairs next time.

4.1 Approximate \(\sqrt{N}\)

• One of usages of Pell’s equation is to approximate \(\sqrt{N}\) with very high-precision.
• We can transfer Pell’s equation to \(\sqrt{N}=\sqrt{p^2/q^2-1/q^2}\). Let \(f(x)=\sqrt{p^2/q^2-x}\). It is obvious that \[\begin{align} f(0)=\frac{p}{q}, \quad\mbox{and}\quad f'(0)=-\frac{q}{2p}. \end{align}\] Then it follows by Taylor series to the first devirative \[\begin{align} \sqrt{N}\approx \frac{p}{q}-\frac{1}{2pq}, \end{align}\] which gives an approximation of \(\sqrt{N}\) with high-precision as following \[\begin{align} \sqrt{N}\approx \frac{p}{q}, \end{align}\] when \(pq\) is large enough. For example, when \(p\approx 10^{50}\) and \(q\approx 10^{50}\), we have \(pq\approx 10^{100}\), and consequently \(\sqrt{N}\) can be approximated to 100 digits if we give the 100 digits of \(p/q\). For the problem of approximation of \(\sqrt{N}\), it remains to calculate \(p/q\) with very high-precision.

5.1 A web-based calculator is available, including any pair \((p,q)\) satisfying Pell’s equation \(p^2-Nq^2=1\) with \(N=2,3,5,6,7,8\).

5.2 The first 50th pairs \((p,q)\) satisfying Pell’s equation \(p^2-2q^2=1\) are below.

# N=2

p= 3 q= 2 

p= 17 q= 12 

p= 99 q= 70 

p= 577 q= 408 

p= 3363 q= 2378 

p= 19601 q= 13860 

p= 114243 q= 80782 

p= 665857 q= 470832 

p= 3880899 q= 2744210 

p= 22619537 q= 15994428 

p= 131836323 q= 93222358 

p= 768398401 q= 543339720 

p= 4478554083 q= 3166815962 

p= 26102926097 q= 18457556052 

p= 152139002499 q= 107578520350 

p= 886731088897 q= 627013566048 

p= 5168247530883 q= 3654502875938 

p= 30122754096401 q= 21300003689580 

p= 175568277047523 q= 124145519261542 

p= 1023286908188737 q= 723573111879672 

p= 5964153172084899 q= 4217293152016490 

p= 34761632124320657 q= 24580185800219268 

p= 202605639573839043 q= 143263821649299118 

p= 1180872205318713601 q= 835002744095575440 

p= 6882627592338442563 q= 4866752642924153522 

p= 40114893348711941777 q= 28365513113449345692 

p= 233806732499933208099 q= 165326326037771920630 

p= 1362725501650887306817 q= 963592443113182178088 

p= 7942546277405390632803 q= 5616228332641321147898 

p= 46292552162781456490001 q= 32733777552734744709300 

p= 269812766699283348307203 q= 190786436983767147107902 

p= 1572584048032918633353217 q= 1111984844349868137938112 

p= 9165691521498228451812099 q= 6481122629115441680520770 

p= 53421565080956452077519377 q= 37774750930342781945186508 

p= 311363698964240484013304163 q= 220167382952941249990598278 

p= 1814760628704486452002305601 q= 1283229546787304717998403160 

p= 10577200073262678228000529443 q= 7479209897770887057999820682 

p= 61648439810871582916000871057 q= 43592029839838017630000520932 

p= 359313438791966819268004696899 q= 254072969141257218722003304910 

p= 2094232192940929332692027310337 q= 1480845785007705294702019308528 

p= 12206079718853609176884159165123 q= 8631001740904974549490112546258 

p= 71142246120180725728612927680401 q= 50305164660422142002238655969020 

p= 414647397002230745194793406917283 q= 293199986221627877463941823267862 

p= 2416742135893203745440147513823297 q= 1708894752669345122781412283638152 

p= 14085805418356991727446091676022499 q= 9960168529794442859224531878561050 

p= 82098090374248746619236402542311697 q= 58052116426097312032565778987728148 

p= 478502736827135487987972323577847683 q= 338352530026789429336170142047807838 

p= 2788918330588564181308597538924774401 q= 1972063063734639263984455073299118880 

p= 16255007246704249599863612909970798723 q= 11494025852381046154570560297746905442 

p= 94741125149636933417873079920900017937 q= 66992092050551637663438906713182313772