```
def pr_sol(maxval,N,ee,ff):
count=0
aa,bb=0,1
gg=N*ff
while count<maxval:
if bb**2 - 1 ==N*aa**2:
print(bb,'^2-1=',N,'*',aa,'^2')
else:
print(bb,'^2+1=',N,'*',aa,'^2')
anew=ee*aa+ff*bb
bnew=gg*aa+ee*bb
aa=anew
bb=bnew
count+=1
return
```

pr_sol(20,92,1151,120) #print first 20 solutions when N=92

""" the smallest solution is 1151,120 where
1151 ^2-1= 92 * 120 ^2 (from wikipedia).
(I'm not sure if it's the same as wikipedia recursion,
but in this case the following recursion applies:

```
a_{t+1}=1151*a_{t} +120*b_{t}
b_{t+1}=120*92*a_{t} +1151*b_{t}
(in general you can generate such a recursion, which gives
values for x^2-1=Ny^2 AND x^2+1=Ny^2
(but x^2 ==91 mod 92 has no solutions so
there are no sol to x^2+1=Ny^2 ))
The first few solutions are:
1 ^2-1= 92 * 0 ^2 (a,b)=( 0 , 1 )
1151 ^2-1= 92 * 120 ^2 (a,b)=( 120 , 1151 )
2649601 ^2-1= 92 * 276240 ^2 (a,b)=( 276240 , 2649601 )
6099380351 ^2-1= 92 * 635904360 ^2
14040770918401 ^2-1= 92 * 1463851560480 ^2
32321848554778751 ^2-1= 92 * 3369785656320600 ^2
I've just been doing some investigation into these equations
and found this recursion. (which has some advantage. For
example N=313 has very large first solution to x^2-1=Ny^2 but
you can generate the numbers from the much smaller solution:
126862368^2 +1 =313 * 7170685^2
and the recursion:
a_{t+1}= 126862368 a_t+7170685 b_t
b_{t+1}=313*7170685a_t+126862368 b_t
so
pr_sol(20,313,126862368,7170685)
prints first 20 solutions for N=313
(note they alternate between x^2+1=313y^2 and x^2-1=313y^2)
```

"""