\\ PARI/GP Code and Notes for A010342
\\ Kevin Ryde, March 2023

\\ A010342 is integers k for which the continued fraction expansion of
\\ sqrt(k) has periodic part of the form {1,...,1,E}, with one or more 1s.
\\ The code here is
\\
\\     * predicate test for such k
\\     * vector-making function for terms to given maximum value
\\     * number of terms <= given k (the inverse function)
\\
\\ Requires PARI/GP 2.13 or higher for sqrtint() "remainder" parameter in
\\ is_A010342().

default(recover,0);
default(strictargs,1);


\\-----------------------------------------------------------------------------
\\ is_A010342()  Predicate

\\ Fibonacci_to_Lucas() takes a possible Fibonacci number f >= 1.
\\ If f is Fibonacci number f = F(x) = A000045(x), then return
\\ corresponding Lucas number L(x) = A000032(x).
\\ If f is not a Fibonacci number, then return 0.
\\ f=1 is taken to be F(1)=1 and the return is L(1)=1.
Fibonacci_to_Lucas(f) =
{
  my(t=5*f^2, s);
  if(issquare(t-4,&s) || issquare(t+4,&s), s, 0);
}

\\ Return 1 if k is a term of A010342, or return 0 if not.
is_A010342(k) =
{
  my(b,q=sqrtint(k,&b));
  b>=2 || return(0);  \\ k cannot be square or square+1

  my(r=(b-1)/(2*q-1),
     Fb=numerator(r),
     Lb=Fibonacci_to_Lucas(Fb));  \\ L(m) if Fb is Fibonacci F(m)
  (Lb
   && Lb + Fb == 2*denominator(r));  \\ denominator must be F(m+1)
}

{
  print1("A010342 = ");
  for(k=0,100,
     if(is_A010342(k), print1(k,",")));
   print();
}

\\-----------------------------------------------------------------------------
\\ A010342_upto()  Terms to Given Maximum Value

\\ Return a vector of terms of A010342 which are <= "limit".
A010342_upto(limit) =
{
  my(l=List([]), m=0);
  while(1,
        m += m%3 + 1;   \\ m>=1 and == 0 or 1 mod 3
        my(Fb=fibonacci(m),
           Fq=fibonacci(m+1),  \\ Fq odd
           d=1);               \\ d odd >= 1
        while(1,
              my(k=((d*Fq+1)/2)^2 + d*Fb + 1);
              if(k>limit, break(if(d==1,2,1)));
              listput(l,k);
              d+=2));
  Set(l);
}

\\ check versus predicate
A010342_upto(1000) == select(is_A010342,[0..1000])  || error();


\\-----------------------------------------------------------------------------
\\ A010342_num_terms_le()  Inverse

\\ A010342_num_terms_le(k) returns the number of terms of A010342
\\ which are <= k.
\\ If k is a term A010342(n) = k then this is its index n,
\\ or if k not a term then the index n of the next smaller term.
A010342_num_terms_le(k) =
{
  my(ret=0,
     alt = 1, \\ (-1)^m starting from m=0
     f=1,     \\ F(m+1)
     g=0);    \\ F(m)
  while(1,
        [f,g] = [f+g,f]; alt = -alt;
        bitand(f,1) || next;   \\ must have F(m+1) odd
        my(ff = f^2,
           e = (ff - f - 2*g + sqrtint(4*(k*ff - alt)))
               \ (2*ff));
        if(e==0, return(ret)); \\ when m too big
        ret += e);
}

{
  \\ check against predicate
  my(want=0);
  for(k=0,1000,
     want += is_A010342(k);
     want == A010342_num_terms_le(k) || error());
}

\\-----------------------------------------------------------------------------
\\ Implementation Notes

\\ The key points for efficient calculation are in Thomas Ordowski's
\\ comments in the sequence.  They can be expanded on here.
\\
\\ The first term of the continued fraction expansion of sqrt(k) is simply
\\
\\     q = floor(sqrt(k)) = sqrtint(k)
\\
\\ A textbook result (by Lagrange's theorem) shows that the expansion is
\\ periodic after q, with form
\\
\\     sqrt(k) = continued fraction [q; {..., 2q} repeating}]
\\
\\ The "..." in the periodic part is zero or more terms before 2q.
\\ The sequence requires these are one or more 1s.
\\ Calculating convergents (which are recurrences) shows that
\\
\\     continued fraction [q; {1,1,...,1,1, 2q} repeating ]
\\                            with m>=1 many 1s
\\     = sqrt(k)
\\ where
\\     k = q^2 + (2*q-1)*F(m)/F(m+1) + 1
\\
\\ and F(x) = A000045(x) is the Fibonacci numbers in usual numbering
\\ F(0)=0, F(1)=1, etc.

F(n) = fibonacci(n);
k_from_q_m(q,m) = q^2 + (2*q-1)*F(m)/F(m+1) + 1;

\\ Let b be the "remainder" from q = floor(sqrt(k)) so that
\\
\\     k = q^2 + b        in code q=sqrtint(k,&b)
\\
\\ Then sequence terms have
\\
\\     b = (2*q-1)*F(m)/F(m+1) + 1
\\
\\ and so ratio
\\
\\     (b-1)/(2q-1) = F(m)/F(m+1)


\\ is_A010342()
\\ ------------
\\
\\ The ratio above allows a term k to be identified from its b-1 and 2q-1.
\\
\\ F(m)/F(m+1) is in least terms since consecutive Fibonacci F(m) and F(m+1)
\\ have no common factor (by induction for instance),
\\
\\     gcd(F(m), F(m+1)) = 1

for(m=1,100, gcd(F(m),F(m+1))==1 || error());

\\ So when (b-1)/(2q-1) is reduced to least terms by dividing out any
\\ common factor d, the resulting numerator and denominator are exactly
\\ F(m) and F(m+1).
\\
\\     Fb = (b-1)/d     must be F(m)
\\     Fq = (2q-1)/d    must be F(m+1)
\\     where d = gcd(b-1, 2q-1)  common factor
\\
\\ In the code, Fibonacci_to_Lucas(f) identifies Fibonacci number F(m).
\\ See below for notes on its implementation.  If Fb is indeed some
\\ Fibonacci number F(m), then it returns Lucas number L(m) = Lb.  
\\ Those together determine F(m+1) by identity
\\
\\     F(m) + L(m) = 2*F(m+1)
\\     where L(x) = A000032(x) = Lucas number
\\
\\ and hence must have Fb + Lb = 2*Fq.
\\
\\ If Fb = 1 then that could be either F(1)=1 or F(2)=1, ie. m=1 or m=2.
\\ But m=2 never occurs for integer k, since then F(m+1) = 2 is even whereas
\\ Fq = (2q-1)/d is always odd.  Fibonacci_to_Lucas() treats input f=1 as
\\ F(1) for that reason.
\\
\\ Finally, the requirement that m>=1 means F(m) >= 1 which is (b-1)/d >= 1
\\ and so b>=2.  This immediately excludes any k which is a perfect square
\\ or square + 1.  Their continued fractions are not the required form,
\\
\\     sqrt(q^2)   = continued fraction [q; ]       terminating
\\     sqrt(q^2+1) = continued fraction [q; {2q} repeating ]
\\       so q^2 and q^2+1 not terms


\\ Fibonacci_to_Lucas()
\\ --------------------
\\
\\ Fibonacci_to_Lucas() uses Pell equation tests as set out in the Fibonacci
\\ characteristic sequence A010056.  Usual theory of Pell equations shows
\\ Fibonacci numbers of even and odd indices are the sole solutions f in
\\
\\     5*f^2 + 4 = s^2, for some s,   iff f = some F(2n)
\\     5*f^2 - 4 = s^2, for some s,   iff f = some F(2n+1)
\\
\\ The s value in both cases is the corresponding Lucas number, from
\\ identities
\\
\\     5*F(2n)^2   + 4 = L(2n)^2
\\     5*F(2n+1)^2 - 4 = L(2n+1)^2
\\
\\ Input f=1 satisfies both equations,
\\
\\     5*1^2 + 4 = 9 is square
\\     5*1^2 - 4 = 1 is square
\\
\\ f=1 is taken to be F(1)=1 (not F(2)=1) by checking the odd case 5*f^2 - 4
\\ first.
\\
\\ In a general purpose Fibonacci to Lucas, might have to worry about what
\\ f=1 should do.  But here only f=1 as F(1)=1 returning L(1)=1 is wanted.


\\ A010342_upto()
\\ --------------
\\
\\ For enumerating terms, it's convenient to go by each m>=1 and odd d>=1,
\\ making b-1 = d*F(m) and 2q-1 = d*F(m+1).  All terms k are then
\\
\\     k = ((d*F(m+1) + 1)/2)^2 + d*F(m) + 1
\\     where d>=1 odd
\\           m>=1 and == 0 or 1 mod 3
\\
\\ 2q-1 is odd so must have both d and F(m+1) odd.
\\ F(m+1) is odd iff m == 0 or 1 (mod 3).

for(m=0,100, (F(m+1)%2==1) == (m%3==0 || m%3==1) || error());

\\ k is an increasing function of d and of m.  Any convenient enumeration of
\\ those could be used.  The code has outer loop m and inner loop d.  When
\\ k>limit is reached, the d loop stops and if that happened at d=1 then the
\\ m loop stops too.  (Any bigger m will exceed limit at d=1 too.)
\\
\\ Fibonacci numbers F(m),F(m+1) could be stepped through by their
\\ recurrence, but the number of terms <= limit is roughly sqrt(limit) so a
\\ very large number of terms before any large Fibonacci number is needed.


\\ A010342_num_terms_le()
\\ ----------------------
\\
\\ The above k in terms of d and m is a quadratic in d for given m.
\\ So at given k, the number of terms <= k can be calculated by taking the
\\ inverse of the quadratic for each m.
\\
\\ If k is a sequence term, then the number of terms <= k is its "inverse"
\\ going from k=a(n) to the index n.
\\
\\ It's convenient to change variables from d to instead an e count,
\\
\\     e = (d+1)/2     ie. d = 2e-1, and so starting at e=1
\\
\\ The quadratic to solve is
\\
\\     a*e^2 + b*e + c = 0       quadratic in e
\\     where a = F(m+1)^2
\\           b = -F(m+1)^2 + F(m+1) + 2*F(m)
\\           c = F(m+1)^2/4 - F(m+1)/2 - F(m) + 5/4 - k
\\
\\ The discriminant multiplied out is
\\ 
\\     b^2 - 4ac = (4*k-4)*F(m+1)^2 + 4*F(m+1)*F(m) + 4*F(m)^2
\\
\\ but identity 
\\ 
\\     F(m+1)*F(m) = F(m+1)^2 - F(m)^2 - (-1)^m
\\     
\\ simplifies to
\\ 
\\     b^2 - 4ac = 4*(k-(-1)^m) * F(m+1)^2
\\     
\\ The positive solution to the quadratic is then
\\
\\               F(m+1)^2 - F(m+1) - 2*F(m)
\\                + sqrt( 4*(k-(-1)^m) * F(m+1)^2 )
\\     e(k,m) =  ----------------------------------
\\                         2*F(m+1)^2
\\
\\ floor(e(k,m)) is the number of terms due to this m which are <= k.
\\ Summing over m is
\\
\\     num_le(k) =   Sum          floor(e(k,m))
\\                m == 0,1 mod 3
\\                and >= 1
\\
\\ e(k,m) form shown above is used in the code.
\\ In the numerator, F(m+1)+2*F(m) is Lucas number L(m+1) if preferred.
\\
\\ The sqrt() in e(k,m) can be sqrtint() in the code since the rest of the
\\ numerator and denominator are integers so any fractional part of the
\\ sqrt() cannot change the integer part of the final e(k,m).


\\ Limit Growth
\\ ------------
\\
\\ e(k,m) grows as sqrt(k) for given m.  Dividing through sqrt(k) is
\\
\\     e(k,m)          / 1 - (-1)^m/k  \     1/2 + L(m+1)/F(m+1)^2
\\     -------  =  sqrt| ------------  |  +  ---------------------
\\     sqrt(k)         \  F(m+1)^2 )   /            sqrt(k)
\\
\\              -> 1 / F(m+1)^2   as k -> infinity
\\              
\\ floor(e(k,m)) is at most a constant -1 smaller than e(k,m) and so the
\\ same limit as e(k,m)/sqrt(k).
\\
\\ Consequently,
\\
\\               n
\\       Lim ----------  =  C
\\           sqrt(a(n))
\\
\\ where
\\                        
\\     C  =    Sum           1 / F(m+1)      odd Fibonaccis
\\           m == 0,1 mod 3   
\\           and >= 1
\\
\\        =  1.696383...  =  A360957 - 1
\\
\\ Or equivalently
\\
\\     a(n) ~ (1/C^2) * n^2    where C^2 = 2.877717...
\\                             or 1/C^2 = 0.347497...
\\


\\-----------------------------------------------------------------------------
\\ Other Ways -- Predicate with Lookup Table

\\ is_A010342() could use a lookup table to identify Fibonacci numbers.
\\ This limits to 2q-1 <= table maximum and trades memory use for some
\\ arithmetic and issquare().
\\
\\ The following measures a little faster, but not enough to worry about.
\\ Remember issquare() often identifies a non-square quickly.


\\ Fibonacci_table[m] = F(m+1)
\\ Can test all x <= Fibonacci_table_max.
Fibonacci_table = vector(20,m, F(m+1));
Fibonacci_table_max = F(#Fibonacci_table+2) - 1;

is_A010342_by_table(k) =
{
  my(b,q=sqrtint(k,&b));
  b>=2 || return(0);
  my(r=(b-1)/(2*q-1),
     Fb=numerator(r),
     Fq=denominator(r));
  if(Fq > Fibonacci_table_max,
     error("Fibonacci_table[] too small"));

  my(m=setsearch(Fibonacci_table,Fq));
  m>=1 && numerator(r) == if(m==1,1,Fibonacci_table[m-1]);
}

for(k=0,1000, is_A010342_by_table(k) == is_A010342(k) || error());


\\-----------------------------------------------------------------------------
\\ Other Ways -- Predicate by Continued Fraction

\\ Actually calculating the continued fraction expansion of sqrt(k) would be
\\ possible.  Its terms must be 1 until some final term E > 1 goes back to
\\ the start of the periodic part, ie. the complete quotient after first
\\ term q.
\\
\\ If there's many 1s then that might be a lot of one by one steps so
\\ efficiency might depend on the expected inputs (ie. how many 1s most of
\\ the time).


\\-----------------------------------------------------------------------------
print("end");