RosettaCodeData/Task/Pell-numbers/PARI-GP/pell-numbers.parigp
2025-08-11 18:05:26 -07:00

59 lines
1.7 KiB
Text
Raw Permalink Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

\\--- pell.gp ------------------------------------------------------------
\\ 1) Define Pell numbers and PellLucas numbers by simple recursion
pell(n) =
{
if (n==0, return(0));
if (n==1, return(1));
return(2*pell(n-1) + pell(n-2));
}
pellL(n) =
{
if (n==0, return(2));
if (n==1, return(2));
return(2*pellL(n-1) + pellL(n-2));
}
\\ 2) Compute the first 10 Pell and PellLucas numbers
pns = vector(10, i, pell(i-1)); \\ pell(0) .. pell(9)
plns = vector(10, i, pellL(i-1)); \\ pellL(0) .. pellL(9)
print("Pell numbers (n=0..9): ", pns);
print("PellLucas numbers (n=0..9): ", plns);
\\ 3) Form the ratios (PellLucas)/2 over Pell, skipping n=0
\\ i.e. for n=1..9
denom = vector(9, i, pns[i+1]); \\ pns[2..10]
numr = vector(9, i, plns[i+1]/2); \\ plns[2..10]/2
approx = vector(9, i, numr[i]/denom[i]);
print("Ratios (pellL(n)/2) / pell(n) for n=1..9:");
for (i=1,9, print(i, ": ", approx[i]));
print("Numeric floats:");
for (i=1,9, print(i, ": ", real(approx[i])));
\\ 4) List (n, pell(n)) for n=0..100 and pick the first 10 where pell(n) is prime
plist = vector(101, i, [i-1, pell(i-1)]);
primespell = select(x->isprime(x[2]), plist);
print("First 10 Pell(n) that are prime:");
for (i=1,10, print(primespell[i]));
\\ 5) Define PellS(n) = If[n==0,1, pell(2n) + pell(2n+1)]
pellS(n) = if(n==0, 1, pell(2*n) + pell(2*n+1));
print("PellS(n) for n=0..19:");
print(vector(20, n, pellS(n)));
\\ 6) Generate Pythagorean triples from Pell numbers:
\\ short = sum_{k=1..2n} pell(k), long = short+1, hypo = pell(2n+1)
pythag(n) =
{
short = sum(k=1, 2*n, pell(k));
long = short + 1;
hypo = pell(2*n + 1);
return([short, long, hypo]);
}
print("Pythagorean triples for n=0..9:");
for (n=0,9, print(n, ": ", pythag(n)));