%%
%%  This is AmSLaTeX    m3011973.tex
%%
%% Here is a character listing to check to be sure that no
%% unwanted translations took place within the bowels of the net.
%% Upper case letters: ABCDEFGHIJKLMNOPQRSTUVWXYZ
%% Lower case letters: abcdefghijklmnopqrstuvwxyz
%% round parentheses, square brackets, curly braces: ()  []  {}
%% Exclaim, at, sharp, dollar, percent: ! @ # $ %
%% Caret, ampersand, star, underscore, hyphen: ^ & * _ -
%% vertical bar, backslash, tilde, backprime, plus: | \ ~ ` +
%% plus, equal, prime, quote, colon: + = ' " :
%% less than, greater than, slash, question, comma: < > / ? ,
%% period, semicolon: . ;
%%
%%
 \documentclass[11pt,reqno,twoside,final]{amsart}

%page size
\hsize=7.7677truein\vsize=10.5truein

% LaTeX's page layout parameters   Pedro Quaresma de Almeida
%                                   1995

\textwidth=6.4in
\textheight=9.25in   %9.5in
\oddsidemargin=.2in
\evensidemargin=.2in
\headheight=0.5in
\topmargin=-.05in
\headsep=0.25in
\voffset=-0.5in

\input epsf.tex    \epsfverbosetrue


\makeatletter
%\def\ddots{\mathinner{\mkern1mu\raise7\p@\vbox{\kern7\p@\hbox{.}}\mkern2mu
%    \raise4\p@\hbox{.}\mkern2mu\raise\p@\hbox{.}\mkern1mu}}
\def\secddots{\mathinner{\mkern1mu\raise\p@\vbox{\kern7\p@\hbox{.}}\mkern2mu
    \raise4\p@\hbox{.}\mkern2mu\raise7\p@\hbox{.}\mkern1mu}}
\makeatother

\usepackage{amssymb}
\usepackage{mathrsfs}
\usepackage{array}
\usepackage{cite}
\usepackage{epic}
\usepackage{eepic}
\usepackage{fancyhdr}

\def\mylmark{\ mylmark initial}
\markboth{}{} \markright{}

\pagestyle{fancy}
\lhead{\rmfamily \textsf{MAPA3011A\ 1997-98\ }-- New difference calculus \ \textbf{\thesection}\ --\ \mylmark\  --} \chead{} \rhead{\rmfamily\thepage}
\lfoot{}\cfoot{}\rfoot{}


\usepackage{pst-plot}

% tables des matieres avec ... : extrait de amsbook.cls avec dotfill au
%   de hfil
\makeatletter
\def\@tocline#1#2#3#4#5#6#7{\relax
  \ifnum #1>\c@tocdepth % then omit
  \else
    \par \addpenalty\@secpenalty\addvspace{#2}%
    \begingroup \hyphenpenalty\@M
    \@ifempty{#4}{%
      \@tempdima\csname r@tocindent\number#1\endcsname\relax
    }{%
      \@tempdima#4\relax
    }%
    \parindent\z@ \leftskip#3\relax \advance\leftskip\@tempdima\relax
    \rightskip\@pnumwidth plus1em \parfillskip-\@pnumwidth
    #5\leavevmode\hskip-\@tempdima #6\relax
    \ \ \dotfill\hbox to\@pnumwidth{\@tocpagenum{#7}}\par
    \nobreak
    \endgroup
  \fi}

% police de chapitre
\def\@makechapterhead#1{\global\topskip 7.5pc\relax
  \begingroup
  \fontsize{\@xivpt}{18}\bfseries\centering
    \ifnum\c@secnumdepth>\m@ne
      \leavevmode \hskip-\leftskip
      \rlap{\vbox to\z@{\vss
%         \centerline{\normalsize\bfseries
          \centerline{\Large    \bfseries
                        \@xp{\chaptername}\enspace\thechapter}
          \vskip 3pc}}\hskip\leftskip\fi
     #1\par \endgroup
  \skip@34\p@ \advance\skip@-\normalbaselineskip
  \vskip\skip@ }
\newcommand{\@chapapp}{\chaptername}

\def\section{\@startsection{section}{1}%
  \z@{.7\linespacing\@plus\linespacing}{.5\linespacing}%
% {\normalfont\bfseries\centering}}
  {\Large\normalfont\bfseries\centering}}



\makeatother

\begin{document}

\thispagestyle{empty}

\begin{center} \textbf{\Large
  MAPA 3011A Special topics in approximation theory.}

\smallskip

 \textbf{1997-1998: New difference calculus \& orthogonal polynomials.}

\bigskip
  Alphonse Magnus,    \\
  Institut de  Math\'ematique Pure et Appliqu\'ee,  \\
 Universit\'e Catholique de Louvain, \\
 Chemin du Cyclotron,2, \\
  B-1348 Louvain-la-Neuve \\
      (Belgium)\\% Research address for author one
        (0)(10)473157 , \texttt{magnus@anma.ucl.ac.be,
  http://www.math.ucl.ac.be/\~{}magnus/ }
%Paul Van Dooren, EULER, a119,  478040 , vandooren@anma.ucl.ac.be
\end{center}
\bigskip

This version: \today

\bigskip

\tableofcontents

\bigskip

We investigate a difference operator
which extends the usual derivative and the classical difference
operators, and appears as a kind of ultimate possibility within a certain
class, especially as far as mappings of polynomials are considered.

\smallskip

The main bulk of the study is on special orthogonal polynomials
satisfying difference equations. Although the new orthogonal
families do not seem (yet) to have many applications in approximation
problems\footnote{One of the first applications deals with a new kind
of discrete spherical harmonics \cite{25}. Interesting, as plain
spherical harmonics involve the first orthogonal polynomials
encountered in history, the Legendre polynomials!}, a unified theory
is always liable to shed more light on a given field.



\section{Difference calculus through the ages.}
\renewcommand{\mylmark}{Diff. calculus through the ages.}
\pagestyle{fancy}

\subsection{Arithmetic progressions} \  \\

Formula for, say, $1+2+\dotsb +n$ can be constructed by looking at
$f(n)$ such that $f(n+1)-f(n) = n+1$. $f(n)= n(n+1)/2$ is readily
found by trial and error, or by amusing (or deep??) links with
geometry (\textbf{\emph{polygonal numbers}}\footnote{See also the appendix.}):

\begin{center}\begin{tabular}{lcccc}
  triangular   &\quad  & square   &\quad  & pentagonal   \\
     &  &    &  & {\Small see \cite[pp.~24-27, 499]{Smi}}  \\
               &  &         &    &             \\
   $\bullet$ o o o   & &$\bullet\bullet\bullet$ & &    $\bullet$ \\
   $\bullet\ \!\bullet$ o o &  & $\bullet\bullet\bullet$ & & $\bullet\ \bullet$ \\
   $\bullet \bullet \bullet$ o &  & $\bullet\bullet\bullet$ & & $\bullet\ \bullet\ \bullet$ \\
                     &   &  & & $\bullet\ \bullet\ \bullet$ \\
                     &   &  & & $\bullet\ \bullet\ \bullet$
\end{tabular}\end{center}


\subsection{Summation formulas} \  \\


Not much more was found before the $17^{\text{th}}$ century, when
Jakob Bernoulli\footnote{(1654-1705) in his posthumus \textsl{Ars Conjectandi}
 (1713), see \cite{lehmer}}
 built a unified theory of the sums
\begin{equation*} %\label{skn}
   S_k(n) = 1 + 2^k +\dotsb + n^k .
\end{equation*}
Of course, we see that $S_k(n+1)-S_k(n)= (n+1)^k$, so that we look for
a discrete integral of $x^k$ through a function whose discrete
derivative (\textbf{\emph{difference}}) is given. The fundamental
theorem of discrete calculus is utterly trivial. The development of
such an elementary difference calculus had to wait for the birth
of the plain calculus, probably from competition.

\textbf{\emph{Moreover}}, there is an extremely interesting
interplay between difference and differential calculus: write \\
$S_k(x+1) -S_k(x) = (x+1)^k$ and\dots differentiate in $x$! \\
$S'_k(x+1) - S'_k(x) = k (x+1)^{k-1}$, showing that \\ $S'_k = k
S_{k-1}$, plus a constant. One finds, from $S_k(0)=0$ and
$S_k(1)=1$\footnote{The actual Bernoulli polynomials are not the
$S_k$ polynomials but $B_{k+1}(x) = (k+1)S_k(x-1) +$ a further constant such
that $\displaystyle{\int_0^1 B_{k+1}(x)\, dx =0}$:  1,\ $x-1/2$,\
$x^2-x+1/6$,\  $x^3- 3x^2/2 +x/2$,\  $x^4-2x^3+x^2-1/30$\dots }:

\begin{center}\begin{tabular}{ccccccc}
$S_0(x)$  & \quad  &  $S_1(x)$  & \quad & $S_2(x)$  & \quad & $S_3(x)$ \\
 $x$ &             & $x(x+1)/2$ &      & $x^3/3 +x^2/2+x/6$& & $x^4/4+x^3/2 +x^2/4 =
 S_1^2(x)$ \\
         &         &            &       &           &      &
     (incidentally (but see \cite{beardon}).)
 \end{tabular}\end{center}

A much simpler set of polynomials with respect to difference calculus
is the set of \textbf{\emph{factorial polynomials}}, for instance the
lower\footnote{written $(x-a)_n$: the subscript is the \emph{umbra}
 (shadow) of the exponent in the usual power \cite{dibuc}.} factorial polynomial
\\
$p_n(x-a) = (x-a)(x-a-1)\cdots (x-a-n+1)$. We immediately find
\begin{equation}    \label{deltapn}
  \Delta p_n(x-a) = p_n(x-a+1) -p_n(x-a) = n p_{n-1}(x-a).
\end{equation}

The simplest way to find the sum of a given polynomial is therefore
to expand it in the $p_n$'s  (Newton's interpolation expansion),
and to sum each of the $p_n$'s from \eqref{deltapn}.

\medskip
Of course, \eqref{deltapn} is the discrete extension of the derivative of
a power. \textbf{\emph{Hypergeometric expansions}} contain as well
factorial polynomials of some arguments and plain powers of other
arguments, allowing for a rich set of difference and differential
formulas.

\subsection{Classical difference operators}  \  \\

\cite{Mi}:

\begin{center}
\begin{tabular}{|c|c|c|c|c|}
 \hline
 $\Delta$  &  $\nabla$   &   $\delta$  &   $\mathsf{E}$  &   $\boldsymbol{\mu}$ \\
 \hline
 $f(x+h)-f(x)$ & $f(x)-f(x-h)$ & $f(x+h/2)-f(x-h/2)$ &  $f(x+h)$ &
                                (f(x+h/2)+f(x-h/2))/2  \\
 \hline
\end{tabular}
\end{center}


\subsection{Geometric progressions and the
            first $q-$difference operator} \  \\

Solving simple instances of $f(qx)- f(x) = g(x)$ amounts to a
sum of $g(x)$, $g(qx)$, $g(q^2x)$,etc.:
\begin{equation*} \begin{split}
           f(qx) -f(x) &= g(x) \\
           f(q^2x) -f(qx) &= g(qx) \\
           f(q^3x) -f(q^2x) &= g(q^2x) \\
              \cdots \ \    & \ \   \cdots  \\
           f(q^N x) -f(q^{N-1} x) &= g(q^{N-1} x) \\
         \rule[3pt]{106pt}{0.1pt}\  & \  \  \rule[3pt]{86pt}{0.1pt} \\
    f(q^N x) - f(x) &= g(x) +g(qx) +\dotsb + g(q^{N-1} x)
\end{split}\end{equation*}

Classical geometric progression, $g(x)=x$, asks for
$\displaystyle{ f(x) = \dfrac{x}{q-1} }$.

Remark that the difference of a polynomial of degree $n$ is no more
a polynomial of degree $n-1$. This is settled by using the
\textbf{\emph{divided difference}}
$\displaystyle{ D_q f(x) = \dfrac{f(qx) - f(x)}{(q-1)x} }$.

\medskip

Here again, shall we find simple factorial-like polynomials recovered
by the action of $D_q$? $p_n(x)$ should be a product such that
$p_n(qx)$ contains almost all the factors of $p_n(x)$. Le me try
$p_n(x) = (x-a)(x-qa)\dots (x-q^{n-1}a)$, then
$p_n(qx) = (qx-a)(qx-qa)\dots (qx-q^{n-1}a) =$, hmm, well, let's
go further,
$p_n(qx) = q^n (x-a/q)(x-a)\dots (x-q^{n-2}a)$, so,   \\
$p_n(qx) -p_n(x)= (x-a)\dots (x-q^{n-2}a) [q^n x -aq^{n-1} -x +q^{n-1}a]$,
$$ D_q (x-a)(x-qa)\dots (x-q^{n-1}a) =
   \dfrac{q^n -1}{q-1} (x-a) \dots (x-q^{n-2}a).
$$

\section{\emph{\textbf{The}} difference operator.}
\renewcommand{\mylmark}{\emph{\textbf{The}} diff. operator.}
\pagestyle{fancy}

\subsection{Introducing $\mathscr{D}$}  \   \\

\begin{minipage}[t]{100pt}
\  \par


\setlength{\unitlength}{1cm}
\begin{picture}(3,2)
 \put(0,0){\vector(1,0){2.5}}  \put(2.4,0.2){$x$}
 \put(0,0){\vector(0,1){2.5}}  \put(0.2,2.4){$y$}
%\put(2.5,0){\line(1,0){0.5}}
 \dashline{0.1}(1,0)(1,2) \put(1.1,0.15){$\eta_1(x)$}
                        \put(1.1,2.25){$\eta_2(x)$}
 \Thicklines
 \spline(0.5,0.5)(1,0.2)(2,1.2)
 \spline(0.5,1.7)(1,2.3)(2,1.6)
\end{picture}
\end{minipage}
 \begin{minipage}[t]{1mm}

 \ \par %\rule{0.2mm}{250pt}
 \end{minipage}
\begin{minipage}[t]{5.06in}
\  \par

\vspace*{-10pt}

We look at the most general (divided) difference operator
of the form
\begin{equation} \label{opdiff}
    (\mathscr{D}f)(x) = \dfrac{ f(\eta_2(x)) -f(\eta_1(x))}
                              { \eta_2(x) -\eta_1(x)}
\end{equation}
leaving a polynomial of degree $n-1$ when $f$ is a polynomial of
degree $n$. The first non trivial conditions are
\end{minipage}

\begin{equation*}\begin{split}
 \eta_1(x)+\eta_2(x) &= \text{ a polynomial of degree } 1;\\
 \eta_1^2(x)+\eta_1(x)\eta_2(x) +\eta_2^2(x) &=
 \text{ a polynomial of degree } 2.
\end{split}\end{equation*}
This is enough to show that $\eta_1(x)$ and $\eta_2(x)$ MUST be
the two roots $y$ of some quadratic equation
\begin{equation} \label{conic}
F(x,y) := \mathfrak{a} y^2 +2\mathfrak{b} xy +\mathfrak{c} x^2
  +2\mathfrak{d} y +2\mathfrak{e} x +\mathfrak{f}=0,
\end{equation}
and that
$$ \eta_1(x) = \mathfrak{p}(x) +\sqrt{\mathfrak{q}(x)},\quad
   \eta_2(x) = \mathfrak{p}(x) -\sqrt{\mathfrak{q}(x)},
$$
with $\mathfrak{p}$ and $\mathfrak{q}$ polynomials of degrees
$\leqslant 1$ and 2,
$$ \mathfrak{p}(x) = p_0x+p_1 = - \dfrac{ \mathfrak{b} x+\mathfrak{d}  }{\mathfrak{a}},\quad
 \mathfrak{q}(x) = q_0x^2 +q_1x +q_2 =\dfrac{ (\mathfrak{b}^2 - \mathfrak{ac}) x^2+
 (\mathfrak{bd}-\mathfrak{ae})x +\mathfrak{d}^2 -\mathfrak{af}  }{\mathfrak{a}^2}.
$$

 So, we have, up to now,

\begin{center}\begin{tabular}{ccc}
   $f(x)$  &  \quad &  $\mathscr{D}f(x)$ \\
    1      &        &     0           \\
   $x$      &        &     1         \\
   $x^2$    &        & $2\mathfrak{p}(x)$  \\
   $x^3$   &      & $3\mathfrak{p}^2(x) +\mathfrak{q}(x)$
 \end{tabular}\end{center}
the action of $\mathscr{D}$ on polynomials up to degree 3. Great!
We can check that $\mathscr{D}$ applied to a
polynomial of any degree yields another polynomial by expanding
powers of $\mathfrak{p}\pm \sqrt{\mathfrak{q}}$, a not very
elegant procedure. We can try to begin a list of algebraic properties
of $\mathscr{D}$. We need the companion operator
$$ \mathscr{M}f(x) = \dfrac{ f(\eta_1(x)) +f(\eta_2(x))}{2}$$
then
\begin{equation}\label{Dfg}\begin{split}
 (\mathscr{D}fg)(x) &=
    \dfrac{ f(\eta_2(x))g(\eta_2(x)) -f(\eta_1(x))g(\eta_1(x))}
                              { \eta_2(x) -\eta_1(x)}   \\
 &= \dfrac{
 \left[\mathscr{M}f  + \sqrt{\mathfrak{q}}\mathscr{D}f\right]
 \left[\mathscr{M}g  + \sqrt{\mathfrak{q}}\mathscr{D}g\right]
 -
 \left[\mathscr{M}f  - \sqrt{\mathfrak{q}}\mathscr{D}f\right]
 \left[\mathscr{M}g  - \sqrt{\mathfrak{q}}\mathscr{D}g\right]
  }{2\sqrt{q}}               \\
&=    \mathscr{M}f\; \mathscr{D}g +
      \mathscr{M}g\; \mathscr{D}f.
\end{split}  \end{equation}

\begin{equation}\label{Mfg}\raisetag{-92pt}\begin{split}
 (\mathscr{M}fg)(x) &=
    \dfrac{ f(\eta_2(x))g(\eta_2(x)) +f(\eta_1(x))g(\eta_1(x))}
                              {2}   \\
 &= \dfrac{
 \left[\mathscr{M}f  + \sqrt{\mathfrak{q}}\mathscr{D}f\right]
 \left[\mathscr{M}g  + \sqrt{\mathfrak{q}}\mathscr{D}g\right]
 +
 \left[\mathscr{M}f  - \sqrt{\mathfrak{q}}\mathscr{D}f\right]
 \left[\mathscr{M}g  - \sqrt{\mathfrak{q}}\mathscr{D}g\right]
  }{2}   \\
&=    \mathscr{M}f\; \mathscr{M}g +
     \mathfrak{q} \mathscr{D}f\; \mathscr{D}g.
\end{split} \end{equation}

So, the operators $\mathscr{D}$ and $\mathscr{M}$ map polynomials of degree $n$ to
polynomials of degrees $\leqslant n-1$ and $n$:
\begin{equation*}
   \begin{bmatrix}  \mathscr{D} xf \\
                    \mathscr{M} xf     \end{bmatrix} =
   \begin{bmatrix}
             \mathfrak{p}      &         1      \\
             \mathfrak{q}      &  \mathfrak{p}
   \end{bmatrix}
   \begin{bmatrix}  \mathscr{D} f \\
                    \mathscr{M} f     \end{bmatrix},
\end{equation*}
and also $\mathscr{D}x^2 f =
          2\mathfrak{p}\mathscr{D}x f +
          (\mathfrak{q}-\mathfrak{p}^2)\mathscr{D} f $.
 Is this $(\mathscr{D}, \mathscr{M})$
treatment so efficient? Painful experiences lead me to prefer often the plain
$f(\eta_1(x))$, $f(\eta_2(x))$ form. For instance, what is the coefficient of
$x^{n-1}$ of $\mathscr{D} x^n$? A little bit of brute force does no harm:
 $$  \mathscr{D} x^n = \dfrac{\eta_2^n(x) -\eta_1^n(x)}
                              {\eta_2(x) -\eta_1(x)}  =
     \dfrac{(\mathfrak{p}+\sqrt{\mathfrak{q}})^n -
            (\mathfrak{p}-\sqrt{\mathfrak{q}})^n }
            { 2\sqrt{\mathfrak{q}} } =
 \dfrac{  (p_0+\sqrt{q_0})^n -
          (p_0-\sqrt{q_0})^n }{2\sqrt{q_0}} x^{n-1} +\dotsb
$$

Remember that we had $\displaystyle{D_q x^n = \dfrac{q^n-1}{q-1} x^{n-1}+\dotsb }$. Will
$\displaystyle{  \dfrac{ p_0 +\sqrt{q_0}}{p_0-\sqrt{q_0}} }$ act as $q$? Suspense\dots

\medskip
Difference of a ratio:
\begin{equation*}\begin{split}
 \left(\mathscr{D}\dfrac{f}{g}\right)(x) &=
    \dfrac{ f(\eta_2(x))/g(\eta_2(x)) -f(\eta_1(x))/g(\eta_1(x))}
                              { \eta_2(x) -\eta_1(x)}    \\
 &= \dfrac{
 \left[\mathscr{M}f  + \sqrt{\mathfrak{q}}\mathscr{D}f\right]
 \left[\mathscr{M}g  - \sqrt{\mathfrak{q}}\mathscr{D}g\right]
 -
 \left[\mathscr{M}f  - \sqrt{\mathfrak{q}}\mathscr{D}f\right]
 \left[\mathscr{M}g  + \sqrt{\mathfrak{q}}\mathscr{D}g\right]
  }{2\sqrt{q} g(\eta_1(x)) g(\eta_2(x))}    \\
&=  \dfrac{  \mathscr{D}f\; \mathscr{M}g -
      \mathscr{D}g\; \mathscr{M}f  }
   { g(\eta_1(x)) g(\eta_2(x)) =
      (\mathscr{M}g)^2 - \mathfrak{q} (\mathscr{D}g)^2  }\, .
\end{split}\end{equation*}


\subsection{$\mathscr{D}-$factorial polynomials, lattices} \  \\

Trying to find a basis of polynomials behaving nicely under $\mathscr{D}$, we look
for a product
\begin{equation} \label{pn}
 \mathsf{p}_n(x)= (x-y_1)(x-y_2)\dotsb (x-y_n)
\end{equation}
  such that
$\mathsf{p}_n(\eta_1(x))$ and $\mathsf{p}_n(\eta_2(x))$ have $n-1$ common factors:
\begin{equation*}\begin{split}
   \mathsf{p}_n(\eta_1(x)) &= (\eta_1(x)-y_1) (\eta_1(x)-y_2) \dotsb (\eta_1(x)-y_{n-1})(\eta_1(x)-y_n) \\
   \mathsf{p}_n(\eta_2(x)) &= (\eta_2(x)-y_1) (\eta_2(x)-y_2) \dotsb (\eta_2(x)-y_{n-1})(\eta_2(x)-y_n)
\end{split}\end{equation*}
Well, well, we can't ask for some $\eta_1(x) -y_j$ to be the same
as some $\eta_2(x) - y_k$ \emph{for all $x$}! This worked only
with $\Delta$, as $\eta_2(x)- \eta_1(x) =$ constant only then.\\
With $D_q$, we found the way out by requesting
$\eta_2(x) -y_{j+1} = $ constant $(\eta_1(x) -y_j)$,
$j= 1,\dotsc,n-1$.\\
And now?

We shall still manage to make  $\mathsf{p}_n(\eta_1(x))$ and
$\mathsf{p}_n(\eta_2(x))$ to vanish \emph{simultaneously} for $n-1$ values of $x$
by having
\begin{equation*}
  y_1= \eta_1(x_1), y_2=\eta_1(x_2) = \eta_2(x_1),\dotsc,
      y_n = \eta_1(x_n) = \eta_2(x_{n-1}).
\end{equation*}

So,
  $\mathsf{p}_n(\eta_1(x))=   0$ at $x=x_1, x_2,\dotsc, x_n$, \\
 $  \mathsf{p}_n(\eta_2(x))= 0$ at $x=x_1,\dotsc, x_{n-1}$, (and $x=x_0$),
 and      \\
 $\mathscr{D}\mathsf{p}_n$

\begin{enumerate}
 \item is a polynomial of degree $n-1$,

 \item vanishes when
$x= x_1, x_2,\dotsc, x_{n-1}$.
\end{enumerate}

So,
\begin{equation}  \label{dhyper1}
(\mathscr{D} \mathsf{p}_n)(x) =
 \mathscr{D} (x-y_1)(x-y_2)\dotsb (x-y_n) = C_n
    (x-x_1)(x-x_2)\dotsb (x-x_{n-1}) = C_n \tilde{\mathsf{p}}_{n-1}(x),
\end{equation}
where $y_j=\eta_1(x_j), y_{j+1}=\eta_2(x_j)$:

\epsfxsize=477bp
\epsfbox[100 150 577 450]{sti1.ps}
\medskip
\begin{center}Figure 1. $x-$ and $y-$ lattices. \end{center}

and where we already encountered the constant $\displaystyle{
 C_n= \dfrac{  (p_0+\sqrt{q_0})^n -
          (p_0-\sqrt{q_0})^n }{2\sqrt{q_0}} }$.


\medskip
Extension to negative indexes:
$$
\text{ if } n<0, \mathsf{p}_n(x) =
             \dfrac{1}{(x-y_{n+1})(x-y_{n+2})\dotsb (x-y_0)}\,.
$$
\textbf{Problem:} we can discuss the effect of $\mathscr{D}$ on
factorial or inverse factorial series [N\"orlund]. What are the
convergence regions of $\mathscr{D}-$factorial series??

\medskip

\emph{\textbf{Exercise.}} Show that
\begin{equation}  \label{mhyper1}   \begin{split}
 \mathscr{M} (x-y_1)(x-y_2)\dotsb (x-y_n) &= [D_n (x-x_n) +E_n]
    (x-x_1)(x-x_2)\dotsb (x-x_{n-1}),\\
   &= D_n \tilde{\mathsf{p}}_n(x) +E_n \tilde{\mathsf{p}}_{n-1}(x)
\end{split}\end{equation}
and look at some values of $D_n$ and $E_n$, $n=0,1,\dotsc$


\bigskip

The points $\{x_i\}$ and $\{y_i\}$ make the \emph{\textbf{lattices}} which are
our main subject of study. The relevant discrete operators act on functions
defined on such lattices.

A set of $x-$ and $y-$ lattices is completely
determined by a starting point $(x_1, y_1)$ on the conic \eqref{conic}.
The, $y_2$ is the second ordinate corresponding to $x_1$, $x_2$ is the
second abscissa corresponding to $y_2$, etc.

{\Small \emph{\textbf{Exercise}} (already stated by people trying to
embarrass the author): what if there is \emph{no} second ordinate,
i.e., if $x_1$ is the abscissa of a vertical tangent to the conic?
}

Note that we should not define $\eta_1$ as the lower part of the
conic as, after some steps, $y_k$ may very well be  the upper
ordinate corresponding to $x_k$. The geometry of the points
$(x_k, y_k)$ supersedes clumsy algebraic definitions. Relevant
items are symmetric functions of $\eta_1$ and $\eta_2$ anyhow.


%\emph{\textbf{Suggestion.}} Some people like to write
%$x_{k+1/2}$ instead of $x_k$, see next section.

\bigskip

\noindent\textbf{Action of $\mathscr{D}$ and $\mathscr{M}$
  on factorial polynomials.}

Let us work \eqref{dhyper1} and \eqref{mhyper1} more
explicitly: as $\mathsf{p}_{n+1}(x) = (x-y_{n+1}) \mathsf{p}_n(x)$, use
\eqref{Dfg}-\eqref{Mfg}:

\begin{equation*} \begin{split}
 \mathscr{D} \mathsf{p}_{n+1}(x) &= \mathscr{M}\mathsf{p}_n(x) + (\mathfrak{p}(x) -
          y_{n+1}) \mathscr{D} \mathsf{p}_n(x) , \\
 \mathscr{M} \mathsf{p}_{n+1}(x) &=
             (\mathfrak{p}(x) -   y_{n+1})   \mathscr{M}\mathsf{p}_n(x)
         +\mathfrak{q}(x) \mathscr{D} \mathsf{p}_n(x) ,
\end{split}\end{equation*}
and put \eqref{dhyper1}-\eqref{mhyper1}:
\begin{equation*} \begin{split}
 C_{n+1} \tilde{\mathsf{p}}_n(x) &= D_n \tilde{\mathsf{p}}_n(x) + E_n \tilde{\mathsf{p}}_{n-1}(x)
    + (\mathfrak{p}(x) -  y_{n+1}) C_n \tilde{\mathsf{p}}_{n-1}(x) , \\
  D_{n+1} \tilde{\mathsf{p}}_{n+1}(x) + E_{n+1} \tilde{\mathsf{p}}_n(x)  &=
     (\mathfrak{p}(x) -  y_{n+1})
      [D_n \tilde{\mathsf{p}}_n(x) + E_n \tilde{\mathsf{p}}_{n-1}(x)]
     +\mathfrak{q}(x)  C_n \tilde{\mathsf{p}}_{n-1}(x) ,
\end{split}\end{equation*}
or, after division by $\tilde{\mathsf{p}}_{n-1}(x)$,
\begin{equation*}\begin{split}
 C_{n+1} (x-x_n)  &= D_n (x-x_n)  + E_n
    + (\mathfrak{p}(x) -  y_{n+1}) C_n, \\
  D_{n+1} (x-x_n)(x-x_{n+1}) + E_{n+1} (x-x_n)   &=
     (\mathfrak{p}(x) -  y_{n+1})
      [D_n (x-x_n)  + E_n ]
     +\mathfrak{q}(x)  C_n,
\end{split}\end{equation*}
we already know $C_n$, check $D_n= C_{n+1} -p_0 C_n$, and find
 $E_n =  - (\mathfrak{p}(x_n) -  y_{n+1}) C_n =
   - C_n \mathfrak{q}(x_n) /(\mathfrak{p}(x_n) -  y_{n+1})$:
\begin{equation} \label{CDE}
 C_n= \dfrac{  (p_0+\sqrt{q_0})^n -
          (p_0-\sqrt{q_0})^n }{2\sqrt{q_0}},
 D_n= \dfrac{  (p_0+\sqrt{q_0})^n +
          (p_0-\sqrt{q_0})^n }{2},
 E_n = C_n \dfrac{y_{n+1} -y_n}{2}.
\end{equation}

\bigskip

So, $\mathscr{D}$ and
    $\mathscr{M}$ transform expansions of $\mathsf{p}$'s as expansions
                                       of $\tilde{\mathsf{p}}$'s.

We will of course need operators making the reverse action.

Let
\begin{equation} \label{tildeD}
 (\widetilde{\mathscr{D}} f)(y) := \dfrac{ f(\xi_2(y)) - f(\xi_1(y)) }
                                         { \xi_2(y) - \xi_1(y) } ,
\end{equation}
where $\xi_1(y)$ and $\xi_2(y)$ are the roots in $x$ of \eqref{conic}
for a given $y$. Then,   \\
$\tilde{\mathsf{p}}_n(\xi_1(y)) = (\xi_1(y)-x_1) \dotsb (\xi_1(y)-x_n)$
vanishes at $y=y_1,\dotsc, y_n$, and \\
$\tilde{\mathsf{p}}_n(\xi_2(y)) = (\xi_2(y)-x_1) \dotsb (\xi_2(y)-x_n)$
vanishes at $y=y_2,\dotsc, y_n$, (and $y=y_{n+1}$):
\begin{equation} \label{tildeDpn}
  \widetilde{\mathscr{D}} (y-x_1)\dotsb (y-x_n) =
 \widetilde{C}_n (y-y_2)\dotsb (y-y_n)
\end{equation}

\subsection{$q$ again} \   \\

Let us try to get the inside secrets of these lattices. Remember
that the two ordinates $y_k$ and $y_{k+1}$ corresponding to
$x_k$ are the two roots of \eqref{conic} when $x=x_k$, so, their
sum is
 $$ y_k +y_{k+1} = 2\mathfrak{p}(x_k) =
              -2(\mathfrak{b} x_k +\mathfrak{d})/\mathfrak{a}.
 $$
Write the same at $k-1$ and sum:
 $$ y_{k-1}+2y_k+y_{k+1} =
          -2(\mathfrak{b} (x_{k-1}+x_k) +2\mathfrak{d})/\mathfrak{a}.
 $$
But $x_{k-1}$ and $x_k$ are the two abscissae corresponding to $y_k$:
look now at $x$ as a root of \eqref{conic} for given $y$:
 $$ x_{k-1}+x_k =
              -2(\mathfrak{b} y_k +\mathfrak{e})/\mathfrak{c}.
 $$
And voila!
\begin{equation} \label{reclat} \begin{split}
    y_{k-1}+(2-4\mathfrak{b}^2/\mathfrak{ac}) y_k+y_{k+1} &=
          4\dfrac{
          \mathfrak{b}\mathfrak{e} -\mathfrak{c}\mathfrak{d} }
              {\mathfrak{ac}},           \\
    x_{k-1}+(2-4\mathfrak{b}^2/\mathfrak{ac}) x_k+x_{k+1} &=
          4\dfrac{
          \mathfrak{b}\mathfrak{d} -\mathfrak{a}\mathfrak{e} }
              {\mathfrak{ac}}.
\end{split}\end{equation}

Remark that $\mathfrak{ac}\neq 0$ if we want the two lattices to exist.

The common exponential part of the solutions of the recurrences
\eqref{reclat} above will of course be called $q^k$ and $q^{-k}$:

\begin{equation} \label{xyq} \begin{split}
 y_k &= \text{ const. }q^k +\text{const. } q^{-k} +
   \dfrac{  \mathfrak{b}\mathfrak{e} -\mathfrak{c}\mathfrak{d} }
              {\mathfrak{ac}- \mathfrak{b}^2},  \\
 x_k &= \text{ const. }q^k +\text{const. } q^{-k} +
   \dfrac{  \mathfrak{b}\mathfrak{d} -\mathfrak{a}\mathfrak{e} }
              {\mathfrak{ac}- \mathfrak{b}^2},
\end{split}\end{equation}
where
\begin{equation} \label{eqq}
  q + q^{-1} =
               4\dfrac{\mathfrak{b}^2}{\mathfrak{ac}}  - 2.
\end{equation}
Remark that the last terms are the coordinates of the \emph{center} of
the conic \eqref{conic}.
\eqref{xyq} is no more valid if $\mathfrak{b}^2 = \mathfrak{ac}$,
i.e., $q=1$.    However, \eqref{reclat} is then easily solved as
\begin{equation} \label{xyq1} \begin{split}
 y_k &= \text{ const. } +\text{const. } k +
  2\dfrac{  \mathfrak{b}\mathfrak{e} -\mathfrak{c}\mathfrak{d} }
              {\mathfrak{ac}} k^2,  \\
 x_k &= \text{ const. } +\text{const. } k +
  2\dfrac{  \mathfrak{b}\mathfrak{d} -\mathfrak{a}\mathfrak{e} }
              {\mathfrak{ac}} k^2.
\end{split}\end{equation}

\medskip

The shape of real conics is linked to $q$:

\begin{center}\begin{tabular}{ccc}
   $q$  &  conic &  Fig.  \\
  $|q| = 1 $     & Ellipse   &     1.           \\
  $q$ real  $\neq \pm 1$      & Hyperbola   &  $2a, 2b$     \\
   $q=1$             &  Parabola    & $2c, 2d$
 \end{tabular}\end{center}

\medskip

Actually, the four constants in \eqref{xyq} or \eqref{xyq1} are not
independent: choose a square root of $\mathfrak{ac}$ and take
$q^{1/2}$ from
\begin{equation*}
   q^{1/2}+q^{-1/2} = 2\dfrac{\mathfrak{b}}{\sqrt{\mathfrak{ac}}}
                            \tag{\ref{eqq}'}
\end{equation*}
then, if $$x_k = A q^k + B q^{-k} +
             \dfrac{\mathfrak{bd-ae}}{\mathfrak{ac-b}^2},$$
$x_{k-1}+x_k = -2(\mathfrak{b} y_k +\mathfrak{e})/\mathfrak{c}$ shows
that $$y_k = -\sqrt{ \dfrac{\mathfrak{c}}{\mathfrak{a}} } \left[
         A q^{k-1/2}   +    B q^{-k+1/2}    \right]
             +\dfrac{\mathfrak{be-cd}}{\mathfrak{ac-b}^2}.$$
Similarly, if $q=q^{1/2}=1$,i.e., $\mathfrak{b} = \sqrt{\mathfrak{ac}}$,
\begin{multline*}
   x_k = A+Bk +2\mathfrak{a}
\dfrac{\mathfrak{d}\sqrt{\mathfrak{c}}
      -\mathfrak{e}\sqrt{\mathfrak{a}}  }
{\mathfrak{ac}} k^2\Rightarrow    \\
 y_k = -\left( \dfrac{ \mathfrak{d}}{2\mathfrak{a}} +
               \dfrac{ \mathfrak{e}}{2\sqrt{\mathfrak{ac}}} \right)
             -\sqrt{ \dfrac{\mathfrak{c}}{\mathfrak{a}} } \left[
         A+B(k-1/2) +2\mathfrak{a}
\dfrac{\mathfrak{d}\sqrt{\mathfrak{c}}
      -\mathfrak{e}\sqrt{\mathfrak{a}}  }
{\mathfrak{ac}} (k-1/2)^2 \right]
\end{multline*}


\noindent
{\Small \emph{\textbf{Funny problem:}} So, $q+q^{-1} \in (-1,1)$ is
for ellipses, $q+q^{-1} \notin (-1,1)$ for hyperbolas, and
$q+q^{-1}=2$, i.e., $q=1$ for parabolas, as expected. But what happens at
$q=-1$??  }

\bigskip

\epsfxsize=477bp
\epsfbox[50  150 527 475]{sti2.ps}
\medskip
\begin{center}Figure 2. Other kinds of lattices. \end{center}

{\Small The figures are picked from \cite{magsnul}.}

A single pair of $(x,y)$ lattices determine the conic \eqref{conic}
if there are at least 5 distinct points. That's why $q=-1$ is so
special: each lattice pair has only 4 points, belonging to ellipses as
well as to hyperbolas (solution to the funny problem).

\subsection{About factorial series as (basic or not) hypergeometric series}
  \    \\

``Old'' plain difference calculus needs factorial series containing products of
the form $\prod_k (x-ak-b)$, i.e., plain hypergeometric series. With ``old''
$q-$difference calculus, $\prod_k (x-aq^k-b)$, so the so-called \emph{\textbf{basic}}
hypergeometric expansions.

\smallskip
(Remember \cite{17} that the ratio of the $(k+1)^{\text{th}}$ term
to the $k^{\text{th}}$ one of ${}_r\phi_s(a_1,\dotsc, a_r;
 b_1,\dotsc,b_s; q,z)$ is
$$
  \dfrac{\text{term}_{k+1}}{\text{term}_k} =
  \dfrac{ (1-a_1 q^k)\dotsb (1-a_r q^k) (-q^k)^{1+s-r} z }
        { (1-b_1 q^k)\dotsb (1-b_s q^k) (1-q^{k+1})  }
$$
[one starts with the $0^{\text{th}}$ term]).

\smallskip

And now, new expansions?

No.

If $q=1$, factorize $x-ak^2 -bk-c$ in $\prod_k (x-ak^2-bk-c)$ as const.
$(k-\xi_1(x)) (k-\xi_2(x))$!  Result is again a plain hypergeometric expansion,
\emph{\textbf{with twice as many parameters}}, involving the algebraic functions
$\xi_1(x)$ and $\xi_2(x)$ in the parameters.

Similarly, if $q\neq 1$, one factorizes $x-A q^k -B q^{-k} -C$ as, say,
 const. $q^{-k} (q^k -\xi_1(x)) (q^k -\xi_2(x))$ to get a basic hypergeometric series.

\smallskip

So, the familiar ${}_2 F$ and ${}_2 \phi$ expansions of the classical orthogonal
polynomials will be succeeded by ${}_4 F$ and ${}_4 \phi$ expansions, very likely.
See also \cite[\S~5, p.~33]{AW}.

\subsection{Parametric representation, the Eastern way} \   \\

Just replace $k$ above by a continuous parameter $s$:

\begin{equation} \label{xys} \begin{split}
  x = x(s)  = A q^s + B q^{-s} +   \text{ const.} , \\
  y = y(s) = -\sqrt{  \dfrac{ \mathfrak{c}}{\mathfrak{a}}} \left(
         A q^{s-1/2}      +B q^{-s+1/2} \right)
             +\text{ const. }
\end{split}\end{equation}
So, the $x-$ and $y-$lattices are basically the same:
$y(s) = \text{ const. } x(s-1/2)+$ const. The $\mathscr{D}$ operator can be written
as a combination of \emph{\textbf{elementary}} difference operators
in terms of this $s$ variable: at $x=x(s)$,
\begin{equation} \label{nikif}
  (\mathscr{D}f)(x(s)) = \dfrac{(\Delta f)(x(s-1/2))  }
                               {\Delta x(s-1/2)  }
                       = \dfrac{(\delta f)(x(s))  }
                               {\delta x(s)  }
\end{equation}
This is the ``Eastern way'' of Nikiforov \emph{et al.}
\cite{24,NS,25}

\smallskip


Considering the limit case $q=1$ and the comfort of real formulas,
6 kinds of lattices are distinguished:

\begin{center}
\begin{tabular}{|c|c|}
 \hline
 ``Eastern''&       \\
  lattice  &  formula       \\
    number &         \\
 \hline\hline
  I    &   $as+b$   \\
  II    &  $as^2 +bs +c$  \\
  III  &  \eqref{xys}, with $AB=0$ \\
  IV  &  \eqref{xys}, with $AB<0$ \\
  V  &  \eqref{xys}, with $AB>0$ \\
  VI  &  \eqref{xys}, with $|q|=1$
 \\
 \hline
\end{tabular}
\end{center}



\subsection{Canonical choices} \   \\

$x-$ and $y-$ lattices are easily transformed under separate
translations and dilations of $x$ and $y$:
\begin{equation} \label{transf} \begin{split}
x=\alpha x'+\beta, y=\gamma y'+\delta &\Rightarrow
    G(x',y') := F(\alpha x'+\beta, \gamma y'+\delta) =0, \\
   & x_k=\alpha x'_k+\beta, y_k=\gamma y'_k+\delta.
\end{split}\end{equation}
The points $(x'_k,y'_k)$ now stay on the conic $G(x',y')=0$ and their
dance in the plane is still homothetic to the $(x_k,y_k)$'s one.

A conic has 5 degrees of freedom, the conics $G$'s associated to $F$
move in a fourdimensional manifold. So, ONE characteristic should be
common to all the $G$'s. Which one? Who did not guess?

Simplest way to get it is to look at large $x$ (and $y$): from
\eqref{conic}, the slopes of the asymptotes of the conic $F=0$
are the two roots of $\mathfrak{a} \sigma^2 +2\mathfrak{b} \sigma +\mathfrak{c}=0$.
Now, from \eqref{transf}, the same for $G$ is given by
 $\mathfrak{a} \beta^2 \sigma'{}^2 +2\mathfrak{b} \alpha\beta \sigma' +
        \mathfrak{c}\alpha^2 =0$. Of course, we find $\displaystyle{
  \sigma' = \dfrac{\alpha}{\beta} \sigma }$ (obvious: $\displaystyle{ \sigma = \lim
  \dfrac{y}{x} = \dfrac{\beta}{\alpha} \lim \dfrac{y'}{x'} =
                 \dfrac{\beta}{\alpha} \sigma' }$).

So, the \emph{ratio of the slopes} of the asymptotes is kept, and this ratio
is (this is your last chance) $q$: let $\rho$ be this ratio, then the
slopes are $\sigma$ and $\rho \sigma$, with sum $\sigma( \rho+1) =
  -2\mathfrak{b}/\mathfrak{a}$ and
 product $\sigma^2 \rho =
  \mathfrak{c}/\mathfrak{a}$, eliminate $\sigma$:\\
$$
  \dfrac{(\rho+1)^2}{\rho} = \rho +2 +\rho^{-1} =
    \dfrac{4 \mathfrak{b}^2/\mathfrak{a}^2}
          {\mathfrak{c}/\mathfrak{a}}   =
   4\dfrac{\mathfrak{b}^2}
          {\mathfrak{ac}} ,
$$
compare with \eqref{eqq}:

\smallskip \emph{\textbf{$q$ is the ratio of the slopes of the
      asymptotes of \eqref{conic}.  }}

(Of course, we may choose beteween a ratio and its inverse. Final formulas
are symmetric in $q$ and $q^{-1}$).

\medskip

As the root of the matter is basically left unchanged by
translations and dilations, one often make special choices for
the conic, such as
\begin{equation} \label{AW1}
  y^2 -(q^{1/2} +q^{-1/2}) xy +x^2  +\dfrac{q+q^{-1}-2}{4} =0,
\end{equation}
which is:
\begin{enumerate}
  \item centered at origin,
  \item symmetric in $x$ and $y$,
  \item inscribed, or circumscribed, to square $[-1,1]\times [-1,1]$:
\end{enumerate}

\psset{unit=1cm}
\begin{pspicture}(17,7)
 \psline{->}(0,3)(7,3)  \uput[90](6.9,3.2){$x$}
 \psline{->}(3,0)(3,6.5) \uput[0](3.2,6.2){$y$}
 \psline(0,0)(6,0)(6,6)(0,6)(0,0)
%\put(2.5,0){\line(1,0){0.5}}
 \psline{->}(8,3)(15,3)  \uput[90](14.9,3.2){$x$}
 \psline{->}(11,0)(11,6)  \uput[0](11.2,5.8){$y$}
 \psline(9.5,1.5)(12.5,1.5)(12.5,4.5)(9.5,4.5)(9.5,1.5)
 \psline[linestyle=dashed,dash=0.2 0.2](8,1.5)(14,1.5)
 \psline[linestyle=dashed](8,4.5)(14,4.5)
 \psline[linestyle=dashed](9.5,0)(9.5,6)
 \psline[linestyle=dashed](12.5,0)(12.5,6)
 \psset{linewidth=0.075}
 \psccurve[linecolor=red](6,5.5)(6,5.75)(5.75,6)(5.5,6)(5,6)(2,4)(0,1)(0,0.5)(0,0.25)(0.25,0)(0.5,0)(1,0)(4,2)(6,5)
 \pscurve[linecolor=red](13,7)(12.75,4.75)(15,5)
 \pscurve[linecolor=red](9,-1)(9.25,1.25)(7,1)
\end{pspicture}

\vspace*{1cm}

\medskip


The related difference operator  is known as the
 \emph{\textbf{Askey-Wilson}} difference operator $\mathcal{D}_q$ in the Western world.

$\mathcal{D}_q$ is often presented as an avatar of the elementary $q-$difference operator
$\delta_q$ through a change of variable:
\begin{equation} \label{AWD}
 (\mathcal{D}_q f)(\cos\theta) =\dfrac{(\delta_q g)(e^{i\theta})}
                                        {\delta_q e^{i\theta}}
  =\dfrac{ g(q^{1/2} e^{i\theta})  - g(q^{-1/2} e^{i\theta}) }
          {  i(q^{1/2} - q^{-1/2}) \sin\theta },
\end{equation}
where $g$ is defined by $g(e^{i\theta} ) = f(\cos\theta)$, i.e.,
$g(z) = f((z+z^{-1}/2)$. So,
$$ (\mathcal{D}_q f)(\cos\theta) =
  \dfrac{ f((q^{1/2} e^{i\theta} + q^{-1/2}e^{-i\theta})/2)
         - f((q^{-1/2} e^{i\theta}+ q^{1/2}e^{-i\theta})/2) }
          {  i(q^{1/2} - q^{-1/2}) \sin\theta },
$$
and one recognizes that $y= (q^{\pm 1/2} e^{i\theta} +
                             q^{\mp 1/2} e^{-i\theta} )/2$ is indeed
related to $x=\cos\theta$ by
$(y-q^{1/2}x)(y-q^{-1/2}x) =
 (y-q^{\pm 1/2}x)(y-q^{\mp 1/2}x) =
 -(q^{\pm 1/2} - q^{\mp 1/2})^2/4$ which is indeed \eqref{AW1}.

\medskip

Being symmetric has some drawbacks, such as forbidding $\eta_1$ and $\eta_2$ to be real
in the hyperbolic case when $x\in (-1,1)$\dots\footnote{However, this strange situation
has triggered interesting investigations by the authors of \cite{BI1}, see next section
on Chebyshev polynomials\dots}. Also, honest real hyperbolas like $y^2-x^2 = 1$ don't
appear real (because $q<0$). And symmetric forms of elementary difference operators, as
$\delta f(x) = (f(x+h/2) - f(x-h/2))/h$ do not supersede the unsymmetric forms
$\Delta$ and $\nabla$. That's why I don't devote my life to symmetric forms.

\medskip


\begin{minipage}[t]{3cm}
\  \par

\begin{pspicture}(3,3)
 \psline{->}(0,0)(3,0)  \uput[90](3,0){$x$}
 \psline{->}(0,0)(0,3) \uput[0](0,3){$y$}
 \psset{linewidth=0.075}
 \pscurve[linecolor=red](2.56,0.36)(1,0)(0.25,0.25)(0,1)(1,3)
\end{pspicture}
\end{minipage}
\begin{minipage}[t]{12.5cm}
\  \par

$q=1$: The \emph{\textbf{Wilson-Racah}} difference operator is
usually written
  $$
     \dfrac{  f( (\sqrt{x}+\alpha/2)^2 )  -
              f( (\sqrt{x}-\alpha/2)^2 )  }{2\alpha\sqrt{x}  }
  $$
(Racah: $\alpha=1$,
 Wilson: $\alpha=i$).

The corresponding parabola is $(y-x)^2 -\alpha^2(x+y)/2
 +\alpha^4/16 =0$. Its axis is $y=x$ and the parabola touches the
$x$ and $y$ axes at $(0,\alpha^2/4)$ and $(\alpha^2/4,0)$.
\end{minipage}

\begin{center}
\begin{tabular}{|c|c||c||c|c|}
 \hline
 ``Western''  & symbol  & conic &``Eastern'' &  canonical \\
              &         &       & lattice    & lattice        \\
     name     &         &       &   number     & \cite[\S 3.4]{25}    \\
 \hline\hline
 AW      &  $\mathcal{D}_q$  &      ell. & VI     &  $\cos(2\omega s)$  \\
         &                   &  hyp.     & IV, V  & $\sinh, \cosh(2\omega s)$   \\
   \hline
 Wilson  &                   &    par. & II  &  $s(s+1)$      \\
          &                  &               &       \\
 \hline
          &   $D_q$         & 2 lines  &  III    &  $e^{2\omega s}$ \\
 \hline
          &   $\Delta$         & 2 parallel lines  &  I    &  $s$ \\
 \hline
\end{tabular}
\end{center}


\subsection{Chebyshev polynomials}  \  \\

Looking for polynomial families with remarkable behaviour under
differencing, one encounters a surprising stability of the
identity
\begin{equation*}\begin{split}
\dfrac{dT_n(x)}{dx} = nU_{n-1}(x), \text{ where }
 T_n(x) &= \dfrac{ (x+\sqrt{x^2-1})^n
                 +(x-\sqrt{x^2-1})^n   }{2}
\\  \text{ and }
 U_{n-1}(x) &= \dfrac{ (x+\sqrt{x^2-1})^n
                 -(x-\sqrt{x^2-1})^n   }{\sqrt{x^2-1} }.
\end{split}\end{equation*}

\smallskip

Apply \eqref{AWD}: $T_n(\cos\theta) = \cos n\theta$, so
$g(z) = (z^n+ z^{-n})/2$,
\begin{equation*} \begin{split}
 (\mathcal{D}_q T_n)(\cos\theta) &=
   \dfrac{ q^{n/2} e^{in\theta}
          +q^{-n/2} e^{-in\theta}
          -q^{-n/2} e^{in\theta}
          -q^{n/2} e^{-in\theta}  } {2i(q^{1/2}-q^{-1/2})\sin\theta} \\
  &= \dfrac{q^{n/2}-q^{-n/2}}{q^{1/2}-q^{-1/2}}
     \dfrac{\sin n\theta}{\sin\theta} \\
  &= \dfrac{q^{n/2}-q^{-n/2}}{q^{1/2}-q^{-1/2}}
        U_{n-1}(\cos\theta).
\end{split}\end{equation*}

The formula is used in \cite{BI1} in order to study $\mathscr{D}$ as
operator on Hilbert spaces.

\medskip

In a noncanonical setting, what is the driving force of such
formulas? Here it is: let the functions $\xi_1$ and $\xi_2$
give the two abscissas for a given ordinate of the conic \eqref{conic}:
$$  F(\xi,y)=0  \Rightarrow  \xi= \xi_1(y) \text{ and }
        \xi=\xi_2(y). $$
$\xi_1$ and $\xi_2$ are the inverse functions to $\eta_1$ and $\eta_2$.

Remark that $\xi_1(x) +\xi_2(x) = -2(\mathfrak{b} x+\mathfrak{e})/\mathfrak{c}$.

Let $f(x) = G(\xi_1(x)) + G(\xi_2(x))$, then
\begin{equation*} \begin{split}
 (\mathscr{D} f)(x) &=
  \dfrac{
            G(\xi_1(\eta_2(x))) + G(\xi_2(\eta_2(x)))
           - G(\xi_1(\eta_1(x))) - G(\xi_2(\eta_1(x)))
       }{\eta_2(x) - \eta_1(x)} \\
  &=\dfrac{
            G(-2(\mathfrak{b} \eta_2(x)+\mathfrak{e})/\mathfrak{c}
              -\eta_2(x))
           -G(-2(\mathfrak{b} \eta_1(x)+\mathfrak{e})/\mathfrak{c}
              -\eta_1(x))
       }{\eta_2(x) - \eta_1(x)} \\
  &\quad \text{ as } \xi_1(\eta_1(x)) =
                   \xi_2(\eta_2(x)) = x, \text{ and from what we
     just saw about } \xi_1+\xi_2 \\
  &= (\mathscr{D}G)(-2(\mathfrak{b} x+\mathfrak{e})/\mathfrak{c} -x).
\end{split}\end{equation*}

\medskip

Hmmm, what if $\mathfrak{b}/\mathfrak{c} =1/2$?

\section{First order difference equations.}
\renewcommand{\mylmark}{$1^{\text{st}}$ order diff. equations.}
\pagestyle{fancy}

\subsection{Summation under $\mathscr{D}$, first difference equations} \  \\

If $\mathscr{D}f=g$, what kind of summation of $g$ shall we
encounter? At $x=x_1$,
$$    f(y_2) - f(y_1) = (y_2-y_1) g(x_1),$$
and we must go further with something $-f(y_2)$:

\begin{equation*} \begin{split}
           f(y_3) -f(y_2) &= (y_3-y_2) g(x_2) \\
           f(y_4) -f(y_3) &= (y_4-y_3)g(x_3) \\
              \cdots \ \    & \ \   \cdots  \\
           f(y_N) -f(y_{N-1}) &= (y_N-y_{N-1}g(x_{N-1}) \\
         \rule[3pt]{106pt}{0.1pt}\  & \  \  \rule[3pt]{86pt}{0.1pt} \\
    f(y_N) - f(y_1) &= \sum_1^{N-1} (y_{k+1}-y_k) g(x_k)
\end{split}\end{equation*}
which we may consider as a plain summation of, say,
$G(x_k) = (y_{k-1}-y_k)g(x_k)$ as, after all,
$          y_{k+1}-y_k       $ is a function of $x_k$.



\subsection{The adjoint of $\mathscr{D}$}  \  \\

Adjoint of $\mathscr{D}$ for the bilinear form
 $\langle f,g\rangle = \sum_k w_k f(y_k) g(y_k)$:

$$
  \langle f, \mathscr{D}g \rangle =
\sum_k  w_k f(y_k) \dfrac{  g(y_{k+1}) -g(y_k)  }{y_{k+1}-y_k}   =
\sum_k  w_k g(y_k)  \left[
         \dfrac{  f(y_{k-1} ) w_{k-1}/w_k  }{y_k-y_{k-1}}   -
         \dfrac{  f(y_k)  }{y_{k+1}-y_k}   \right]
$$

So, $$(\mathscr{D}^\dag  f)(x_k) =
         \dfrac{  f(y_{k-1} ) w_{k-1}/w_k  }{y_k-y_{k-1}}   -
         \dfrac{  f(y_k)  }{y_{k+1}-y_k}
$$

\subsection{First order difference equations}  \   \\


Just as $(\mathscr{D}f)(x) =g(x)$ involves $f(\eta_1(x))$ and
$f(\eta_2(x))$, a first order $\mathscr{D}-$difference equation is
  $$  F( f(\eta_1(x)), f(\eta_2(x)),x)=0,$$
or
  $$   f(\eta_2(x)) = H( f(\eta_1(x)),x),$$
allows the \emph{\textbf{recurrence}}
  $$   f(y_{k+1}) = H( f(y_k),x_k),\qquad k=0,1,\dotsc$$

\smallskip

In $(\mathscr{D}, \mathscr{M})$ form:
$$ G((\mathscr{D}f)(x), (\mathscr{M}f)(x), x)=0. $$

\subsection{Linear first order difference equations}  \   \\


In $(\mathscr{D}, \mathscr{M})$ form:
$$ A(x)(\mathscr{D}f)(x) +B(x) (\mathscr{M}f)(x) +C(x)=0. $$

\smallskip

Becomes

  $$
\left(  \dfrac{A(x)}{\eta_2(x)-\eta_1(x)} +\dfrac{B(x)}{2} \right)
     f(\eta_2(x))
-\left(  \dfrac{A(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{B(x)}{2} \right)
      f(\eta_1(x))
  +C(x) =0.$$

\medskip

\textbf{Solution in factorial series.} For instance, the
$\mathscr{D}-$exponential function, solution of
 $$  \mathscr{D} f  = \omega \mathscr{M}f $$
can be searched as $\displaystyle{\sum_0^\infty  c_n (x-y_1)\dotsb
  (x-y_n)  }$. From \eqref{dhyper1} and \eqref{mhyper1} (p.\pageref{mhyper1}),
$$ c_n C_n = \omega (c_{n-1} D_{n-1} +c_n E_n), \text{ i.e. },
   c_n = \dfrac{ \omega D_{n-1} }{C_n -\omega E_n}\, c_{n-1}, n=1,2,
   \dotsc
$$

\smallskip

For $  A\mathscr{D} f  = \omega \mathscr{M}f $, where $A(x)=A_0x +A_1$,
one finds
$$   c_n = \dfrac{ \omega D_{n-1} -A_0 C_{n-1}}
       {C_n A(x_n) -\omega E_n}\, c_{n-1}, n=1,2,   \dotsc
$$
\medskip

For $A\mathscr{D}f +B\mathscr{M}f +C=0$ with polynomial $A,B,C$:
$$ A \sum c_n \mathscr{D} \mathsf{p}_n +B \sum c_n \mathscr{M} \mathsf{p}_n +C=0,$$
$$A \sum c_n C_n \tilde{\mathsf{p}}_{n-1}
 +B \sum c_n [D_n \tilde{\mathsf{p}}_n +E_n \tilde{\mathsf{p}}_{n-1}] +C=0,$$
expand $A$ and $B$ as Newton interpolants about $x_n$, $x_{n+1}$,\dots
\begin{equation*}  \begin{split}
&\sum_n  c_n C_n \{ A(x_n) +[x_n, x_{n+1}]_A (x-x_n) +\dotsb \} \tilde{\mathsf{p}}_{n-1}   \\
+&\sum_n  c_n D_n \{ B(x_{n+1}) +[x_{n+1}, x_{n+2}]_B (x-x_{n+1}) +\dotsb
           \} \tilde{\mathsf{p}}_n                                      \\
+&\sum_n c_n E_n \{ B(x_n) +[x_n, x_{n+1}]_B (x-x_n) +\dotsb \} \tilde{\mathsf{p}}_{n-1}
+C=0,
\end{split}\end{equation*}
Finally, coefficient of $\tilde{\mathsf{p}}_N$:
\begin{equation}  \label{coeffpn}   \hspace*{-22pt}
\sum_{n=N-\max((\text{degree } A)-1, \text{ degree } B)}^{N+1}
 \hspace*{-43pt}
 c_n\{ C_n [x_n,\dotsc, x_{N+1}]_A
      +D_n [x_{n+1},\dotsc, x_{N+1}]_B
      +E_n [x_n,\dotsc, x_{N+1}]_B \} =0,
\end{equation}
$N\notin[0,\text{ degree } C]$, which is normally a linear
recurrence relation of order $\max(\text{ degree } A, 1+\text{ degree } B)$.

\subsection{Riccati difference equations}  \   \\


\begin{equation} \label{ric}
 A(x)(\mathscr{D}f)(x) =
  B(x)  f(\eta_1(x)) f(\eta_2(x))
+C(x) (\mathscr{M}f)(x) +D(x).
\end{equation}

\smallskip

\begin{equation} \label{ricr}
       f(\eta_2(x)) =
  \dfrac{
        \left( \dfrac{A(x)}{\eta_2(x)-\eta_1(x)} +\dfrac{C(x)}{2} \right)
     f(\eta_1(x))  +D(x)
        }{
            \dfrac{A(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{C(x)}{2}
       -B(x) f(\eta_1(x))
        }
\end{equation}


\medskip

Continued fraction:
\begin{equation*}\begin{split}
     f(y_{k+1}) &=
  \dfrac{
        \left( \dfrac{A(x_k)}{\Delta y_k} +\dfrac{C(x_k)}{2} \right)
     f(y_k)  +D(x_k)
        }{
            \dfrac{A(x_k)}{\Delta y_k} -\dfrac{C(x_k)}{2}
       -B(x_k) f(y_k)
        }        \\
  &= -\dfrac{A(x_k)}{B(x_k)\Delta y_k} -\dfrac{C(x_k)}{2B(x_k)}
 + \dfrac{1}{B(x_k)}\, \dfrac{
        \left( \dfrac{A(x_k)}{\Delta y_k} \right)^2
        -\dfrac{(C(x_k))^2}{4}
       +B(x_k)D(x_k)
        }{
            \dfrac{A(x_k)}{\Delta y_k} -\dfrac{C(x_k)}{2}
       -B(x_k) f(y_k)
        }.
\end{split}\end{equation*}
Remark the importance of $\Omega := A^2 -\mathfrak{q}(C^2 -4BD)$: if $f$ is
finite at a zero of $\Omega$, $f$ can be computed in finite numbers of
steps at each subsequent point of the lattice of this zero.

\medskip

\textbf{Relation with linear second order equation:}

let  $\displaystyle{
            \dfrac{A(x_k)}{y_{k+1}-y_k} -\dfrac{C(x_k)}{2}
       -B(x_k) f(y_k)  = \dfrac{F_k}{F_{k-1}}   }$, then

  $\displaystyle{
 f(y_{k+1}) =
     \dfrac{   -\dfrac{F_{k+1}}{F_k}
     + \dfrac{A(x_{k+1})}{y_{k+2}-y_{k+1}} -\dfrac{C(x_{k+1})}{2}
         }{B(x_{k+1})}\;,
  }$

\smallskip

\begin{equation*} \begin{split}
     B(x_{k+1}) & F_{k+1} =    \\
 & \left\{  B(x_k)
     \left[ \dfrac{A(x_{k+1})}{y_{k+2}-y_{k+1}} -\dfrac{C(x_{k+1})}{2}
     \right]
       - B(x_{k+1})\left[\dfrac{A(x_k)}{y_{k+1}-y_k} -\dfrac{C(x_k)}{2} \right]  \right\}
     F_k \\
  &  -B(x_{k+1}) \left\{
       \left(\dfrac{A(x_k)}{y_{k+1}-y_k}\right)^2
        -\dfrac{(C(x_k))^2}{4} +B(x_k)D(x_k)  \right\}
    F_{k-1}
\end{split}\end{equation*}


\section{Second order difference equations.}
\renewcommand{\mylmark}{$2^{\text{nd}}$ order diff. equations.}
\pagestyle{fancy}


\subsection{}

\begin{equation} \label{op2}
  (\mathscr{L} f)(x_k) =
       \alpha(x_k) f(y_{k+1})
      + \beta(x_k) f(y_k)
      +\gamma(x_k) f(y_{k-1})
\end{equation}


If $\mathscr{L}$ must send polynomials to polynomials,

\begin{equation}  \begin{split}
  \alpha_k  +\beta_k +\gamma_k  &= \text{ polynomial } = p_0(x_k) \\
  \alpha_k y_{k+1} +\beta_k y_k +\gamma_k y_{k-1} &= \text{ polynomial } = p_1(x_k) \\
  \alpha_k y_{k+1}^2 +\beta_k y_k^2 +\gamma_k y_{k-1}^2 &= \text{ polynomial } = p_2(x_k)
\end{split}\end{equation}

Use \eqref{conic}:
$$
\alpha(x_k) F(x,y_{k+1}) +\beta(x_k) F(x,y_k) +\gamma(x_k) F(x,y_{k-1}) =
\mathfrak{a} p_2(x_k) +2(\mathfrak{b} x+\mathfrak{d}) p_1(x_k)
 +(\mathfrak{c} x^2 +2\mathfrak{e}x +\mathfrak{f}) p_0(x_k)
$$

$x=x_k \Rightarrow  \gamma = $

$x=x_{k-1} \Rightarrow  \alpha = $

Hmm, drop that subsection

\subsection{Selfadjointness}

\section{Orthogonal polynomials satisfying $\mathscr{D}-$difference equations.}
\renewcommand{\mylmark}{Orth. polynomials.}
\pagestyle{fancy}

\subsection{Instant Pad\'e-continued fractions- orthogonal polynomials
  theory}  \  \\


Given a real or complex sequence $\displaystyle{\{\mu_m\}_0^\infty  }$,
the three following problems are basically the same ([Chebyshev,
 Frobenius, Pad\'e, Perron, Stieltjes, Szeg\H o, etc., etc.]):

\begin{enumerate}

 \item \emph{\textbf{Continued fraction:}} expand the (formal)
 series $\displaystyle{  \mathcal{S}(z) = \sum_{j=0}^\infty
   \dfrac{\mu_j}{z^{j+1}}   }$ in a Jacobi continued fraction
   $$
      \dfrac{\mu_0}{ z-b_0 -
                           \dfrac{a_1^2}{ z-b_1 -
                                                \dfrac{a_2^2}{ z-b_2 -
                                                      \dotsb }
                                         }
                   }
   $$

 \item \emph{\textbf{(Formal) orthogonal polynomials:}} look for
   polynomials $p_0, p_1,\dotsc$, with $p_n$  of degree $n$:
   $p_n(z) = \gamma_n z^n + \gamma'_n z^{n-1} +\dotsb$, orthonormal
   with respect to the form of moments $\mu_m$:
   $$ \langle p_n, p_m \rangle = \delta_{n,m}, n,m=0,1,\dotsc, \quad
      \text{ where } \langle x^i, x^j \rangle = \mu_{i+j}.
   $$

 \item \emph{\textbf{Pad\'e approximation:}} look for rational functions
  $\tilde{p}_{n-1}/p_n$ (yes, the same $p_n$ than just before) such that
   \begin{equation} \label{Pade}
        \mathcal{S}(z) - \dfrac{\tilde{p}_{n-1}(z) }{p_n(z) } =
         O\left(  \dfrac{1}{z^{2n+1}}  \right)\;, n=0,1,\dotsc
    \end{equation}
\end{enumerate}

\noindent
The problems so stated have each a unique solution\footnote{Unique up to the sign
  of $\gamma_n$.} if all the Hankel determinants     \\
  $\displaystyle{
      H_n = \begin{vmatrix} \mu_0   &   \mu_1  &  \dotsb  &  \mu_n  \\
                            \mu_1   &   \mu_2  &  \dotsb  &  \mu_{n+1} \\
                            \vdots  &  \vdots  &  \dotsb  &  \vdots    \\
                            \mu_n   &  \mu_{n+1}& \dotsb  & \mu_{2n}
            \end{vmatrix}
                }
    $
  do not vanish.

\smallskip

The relations between the solutions of the three problems are seen from the recurrence
relations
\begin{equation}  \label{recur}  \begin{split}
  a_{n+1} p_{n+1}(z) &= (z-b_n) p_n(z) -a_n p_{n-1}(z), \\
  a_{n+1}\tilde{p}_n(z) &= (z-b_n)\tilde{p}_{n-1}(z) -a_n \tilde{p}_{n-2}(z),
\end{split}\end{equation}
$n=1,2,\dotsc$,
\begin{equation} \label{recurinit} \begin{split}
 p_0 = \gamma_0 = 1/\sqrt{\mu_0},   \quad & \quad     \tilde{p}_{-1} =0,  \\
 p_1(z) = \gamma_0 (z-b_0)/a_1,  \quad & \quad
                \tilde{p}_0 = \mu_0 \gamma_1 = \sqrt{\mu_0}/a_1, \\
  \text{ i.e., extended by }  & \eqref{recur} \text{ to $n=0$ by  }   \\
  a_0 p_{-1}(z)\equiv 0,  \quad & \quad   a_0 \tilde{p}_{-2} = -\sqrt{\mu_0}
\end{split}\end{equation}

\medskip

Interesting identities are, among other ones,
 \begin{equation*}
    \gamma_0^2 \gamma_1^2 \dotsb \gamma_n^2 = H_n,
 \end{equation*}
 \begin{equation} \label{agamma}
   a_n = \dfrac{\gamma_{n-1}}{\gamma_n}, \quad n=1,2,\dotsc
 \end{equation}
 \begin{equation} \label{casor}
 p_n \tilde{p}_n -p_{n+1}\tilde{p}_{n-1} = \dfrac{1}{a_{n+1}}, \quad
                     n=0,1,\dotsc
 \end{equation}
\begin{equation}  \label{mainc} \left. \begin{split}
   p_n(z) &= \dfrac{\gamma_0}{a_1 \dotsb a_n} [
         z^n -(b_0  +\dotsb +b_{n-1} )z^{n-1} + \dotsb ] \\
  \tilde{p}_n(z) &=
             \dfrac{1}{\gamma_0 a_1 \dotsb a_{n+1}} [
         z^n -(b_1 +\dotsb +b_n )z^{n-1} + \dotsb ] \\
  \mathcal{S}(z) p_n(z) - \tilde{p}_{n-1}(z) &=
  \dfrac{a_1 \dotsb a_n}{\gamma_0 z^{n+1} }\left[ 1+
     \dfrac{ b_0 + \dotsb + b_n}{z} +\dotsb \right]
\end{split}\right\}\end{equation}


\bigskip

A special family of orthogonal polynomials is considered
essentially known if we have a compact formula for $\mathcal{S}$
(or the moments $\mu_j$) \emph{\textbf{and}} the recurrence
(or continued fraction) coefficients $\{a_n\}$ and $\{b_n\}$.

For instance, $\displaystyle{ \mathcal{S}(z) =
  \int_{-1}^1 \dfrac{dx}{z-x} = \dfrac{2}{z} + \dfrac{2}{3z^3}
  +\dfrac{2}{5z^5} +\dotsb = \log\left( \dfrac{z+1}{z-1}\right)  }$
gives $b_n=0$, $a_n^2 = \displaystyle{ \dfrac{n^2}{4n^2-1} }$
(Legendre).  How on earth is it being done?


\subsection{Semi-classical and Laguerre-Hahn orthogonal polynomials}
  \  \\

Extensions of Gauss and Heine theory of orthogonal polynomials-continued
fractions, emphasizing the role of differential (or difference [in $x$])
equations and relations received the name ``semi-classical''. The choice
of this name \cite{10,22} is perhaps not so fortunate, as a body of
calculations in quantum mechanics bears the same name.

I coined the name ``Laguerre-Hahn'' for what could be the ultimate
big family of orthogonal polynomials satisfying differential or
difference relations. Laguerre \cite{19} started to investigate
continued fractions satisfying differential equations, and showed how
to construct the differential equations satisfied by each
denominator (orthogonal polynomial). Hahn
(\cite{ha1,ha2} and earlier works quoted therein) worked in the
converse direction: what happens when orthogonal polynomials
satisfy linear differential equations of fixed order, with coefficients
of fixed degree? The answer is that they are basically the family that
can be reached by the Laguerre theory. Well, I suppose, because not many
people, and certainly not me, went very far in the study of the most
difficult Hahn's theorems.

\smallskip

Other smart people involved: Bonan-Lubinsky-Nevai \cite{4,5},
the Chudnovsky brothers, Freud,
Gammel-Nuttall, Shohat,\dots

\bigskip
\bigskip
Here is the Laguerre-like theory, in $\mathscr{D}-$language:
we consider Stieltjes functions satisfying a \emph{\textbf{Riccati}}
equation   \eqref{ric}
$$ A(x)(\mathscr{D}\mathcal{S})(x) =
  B(x)  \mathcal{S}(\eta_1(x)) \mathcal{S}(\eta_2(x))
+C(x) (\mathscr{M}\mathcal{S})(x) +D(x)   $$
 with polynomial $A, B, C$, and $D$.
The \emph{\textbf{semi-classical}} special case is when
$B(x)\equiv 0$.

Let $\displaystyle{  \mathcal{S}_m(z) = \dfrac{1}{z-b_m -\dfrac{a_{m+1}^2}
                              {z-b_{m+1} - \dfrac{a_{m+2}^2}{\dotsb}
                                         }}
                }$.

Then, each $\mathcal{S}_m$ satisfies a Riccati equation
\begin{equation} \label{ricm}
 A_m(x)(\mathscr{D}\mathcal{S}_m)(x) =
  B_m(x)  \mathcal{S}_m(\eta_1(x)) \mathcal{S}_m(\eta_2(x))
+C_m(x) (\mathscr{M}\mathcal{S})(x) +D_m(x).
\end{equation}
 similar to
\eqref{ric} with polynomial coefficients with \emph{\textbf{same
degrees}}:   1.) it's true for $m=0$, as $\mathcal{S}=
 \mu_0 \mathcal{S}_0$,\quad 2.) enter
  $\displaystyle{  \mathcal{S}_m(x) =
                           \dfrac{1}{x-b_m -a_{m+1}^2
              \mathcal{S}_{m+1}(x)   }   }$ in \eqref{ricm}.
Things are easier with the recurrent form \eqref{ricr}:
\begin{equation*}   \begin{split}
 \dfrac{1}{\eta_2-b_m -a_{m+1}^2  \mathcal{S}_{m+1}(\eta_2)   }
  &=
  \dfrac{
        \left( \dfrac{A_m(x)}{\eta_2(x)-\eta_1(x)} +\dfrac{C_m(x)}{2} \right)
 \dfrac{1}{\eta_1-b_m -a_{m+1}^2  \mathcal{S}_{m+1}(\eta_1)   }
       +D_m(x)
        }{
            \dfrac{A_m(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{C_m(x)}{2}
 -\dfrac{B_m(x)}{\eta_1-b_m -a_{m+1}^2  \mathcal{S}_{m+1}(\eta_1)   }
        }   \\
  &=
  \dfrac{
         \dfrac{A_m(x)}{\eta_2(x)-\eta_1(x)} +\dfrac{C_m(x)}{2}
 +D_m(x) [ \eta_1-b_m -a_{m+1}^2  \mathcal{S}_{m+1}(\eta_1)   ]  }
 { -B_m(x) + \left(\dfrac{A_m(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{C_m(x)}{2}
                \right)
   [ \eta_1-b_m -a_{m+1}^2  \mathcal{S}_{m+1}(\eta_1)  ] } \;,
\end{split}\end{equation*}
or
$
-a_{m+1}^2 \mathcal{S}_{m+1}(\eta_2) =$ \\
$$
  \dfrac{
     -a_{m+1}^2\left(\dfrac{A_m(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{C_m(x)}{2}
               -(\eta_2 -b_m)D_m(x)\right)
         \mathcal{S}_{m+1}(\eta_1)     +X_m
          }
   {      \dfrac{A_m(x)}{\eta_2(x)-\eta_1(x)} +\dfrac{C_m(x)}{2}
 +( \eta_1-b_m )D_m(x)
            -a_{m+1}^2  \mathcal{S}_{m+1}(\eta_1)     }\;,
$$
where $X_m =$ \\
 $\displaystyle{ \hspace*{-1.75cm}
   \left(\dfrac{A_m(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{C_m(x)}{2}\right)
               (\eta_1 -b_m)
            -B_m(x)
   -\left(\dfrac{A_m(x)}{\eta_2(x)-\eta_1(x)} +\dfrac{C_m(x)}{2}\right)
               (\eta_2 -b_m)
    -(\eta_1-b_m)(\eta_2-b_m)D_m(x)   }$,
so that
\begin{equation*}   \begin{split}
    \dfrac{A_{m+1}(x)}{\eta_2(x)-\eta_1(x)} +\dfrac{C_{m+1}(x)}{2} &=
    \dfrac{A_m(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{C_m(x)}{2}
                  -( \eta_2-b_m )D_m(x)    \\
    \dfrac{A_{m+1}(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{C_{m+1}(x)}{2} &=
    \dfrac{A_m(x)}{\eta_2(x)-\eta_1(x)} +\dfrac{C_m(x)}{2}
                  +( \eta_1-b_m )D_m(x)    \\
   B_{m+1}(x) &= a_{m+1}^2 D_m(x) \\
 a_{m+1}^2 D_{m+1}(x) = -X_m &= A_m(x) +(\eta_1 +\eta_2 -2b_m) C_m(x)/2
   +B_m(x) +(\eta_1-b_m)(\eta_2-b_m) D_m(x)
\end{split}\end{equation*}
The two first equations above simplify as
\begin{equation}   \label{eqA}
  A_{m+1}(x) =  A_m(x) -\dfrac{(\eta_2(x) -\eta_1(x) )^2}{2}\, D_m(x)
             =  A_m(x) -2\mathfrak{q}(x) D_m(x),
\end{equation}
\begin{equation}   \label{eqC}
  C_{m+1}(x) = -C_m(x) -(\eta_1(x) +\eta_2(x) -2b_m) D_m(x)
         = -C_m(x) -2(\mathfrak{p}(x) -b_m) D_m(x)
\end{equation}
showing, together with the already encountered
\begin{equation}   \label{eqB}
   B_{m+1}(x) = a_{m+1}^2 D_m(x),
\end{equation}
and
\begin{equation}   \label{eqD}     \begin{split}
 a_{m+1}^2 D_{m+1}(x) &=  A_m(x) +(\eta_1 +\eta_2 -2b_m) C_m(x)/2
   +B_m(x) +(\eta_1-b_m)(\eta_2-b_m) D_m(x)   \\
 &=  A_m(x) -\dfrac{(\eta_2 -\eta_1)^2}{4}  D_m(x)
          +\dfrac{(\eta_1 +\eta_2 -2b_m) C_m(x)}{2}
   +B_m(x) +\dfrac{(\eta_1+\eta_2-2b_m)^2 D_m(x)}{4} \\
 &=  A_m(x) -\dfrac{(\eta_2 -\eta_1)^2}{4} D_m(x)
          +\dfrac{(\eta_1 +\eta_2 -2b_m) (C_m(x)-C_{m+1}(x))}{2}
   +B_m(x),
\end{split}\end{equation}
that $A_{m+1}, B_{m+1}, C_{m+1}$, and $D_{m+1}$ are
polynomials if $A_m, B_m, C_m$, and $D_m$ are polynomials ($\eta_1+
\eta_2$ and $\eta_1\eta_2$ are polynomials of degrees $\leqslant 1$ and
2 [see p.~\pageref{opdiff}]). What about the degrees? Let
$$ s = \max[ (\text{degree } A_m)-2,
             (\text{degree } B_m)-2,
             (\text{degree } C_m)-1  ]
$$
(the notation was launched by Maroni \cite{22}). Then, degree $D_m \leqslant s$,
from \eqref{ricm}, as $\mathcal{S}_m = O(z^{-1}) \Rightarrow
              \mathscr{D}\mathcal{S}_m = O(z^{-2}) \text{ and }
              \mathscr{M}\mathcal{S}_m = O(z^{-1})$, so the
degrees of $A_{m+1}, B_{m+1}$, and $C_{m+1}$ are bounded by
 $s+2$, $s$, and $s+1$, from the equations above. And the degree of
$D_{m+1}$ \emph{\textbf{must}} be bounded by $s$, as
                      $\mathcal{S}_{m+1} = O(z^{-1})$ still holds.


\bigskip
\eqref{eqA}-\eqref{eqD} can be arranged in a matrix-vector
recurrence with polynomial coefficients. However, painful experiments
lead to prefer (again) the intermediate form in $\eta_1$ and $\eta_2$

\begin{multline} \label{recmatrix}
  \begin{bmatrix}
    -\dfrac{A_{m+1}(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{C_{m+1}(x)}{2} \\
    \dfrac{A_{m+1}(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{C_{m+1}(x)}{2} \\
        B_{m+1}(x) \\  D_{m+1}(x)
  \end{bmatrix}
    =  \\
  \begin{bmatrix}
                &   -1     &        &   \eta_2 -b_m  \\
          -1     &         &        &    \eta_1 -b_m  \\
                &         &        &    a_{m+1}^2   \\
            -\dfrac{\eta_2-b_m}{a_{m+1}^2}  &
            -\dfrac{\eta_1-b_m}{a_{m+1}^2}  &
             \dfrac{1}{a_{m+1}^2}  &
            \dfrac{(\eta_1-b_m)(\eta_2-b_m)}{a_{m+1}^2}
  \end{bmatrix}
  \begin{bmatrix}
    -\dfrac{A_{m}(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{C_{m}(x)}{2} \\
    \dfrac{A_{m}(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{C_{m}(x)}{2} \\
        B_{m}(x) \\  D_{m}(x)
  \end{bmatrix}
\end{multline}

\bigskip
That the degree of $D_m$ stays $\leqslant s$ does not seem obvious from
   \eqref{eqD}. This non-obviousness
is a clue to the calculation of the $a_m$'s and the $b_m$'s!!!


\medskip

Often useful, the sum of \eqref{eqD}:
$$
 a_{m+1}^2 D_{m+1}(x) D_m(x) - a_m^2 D_m(x) D_{m-1}(x)=
  \dfrac{A_m^2(x) -A_{m+1}^2(x)}{(\eta_2 -\eta_1)^2}
   +\dfrac{ C_{m+1}^2(x) -C_m^2(x) }{4},
$$
so
\begin{equation} \label{sumeqD} \begin{split}
 a_m^2 D_m(x) D_{m-1}(x) - B_0(x)  D_0(x)&=
  \dfrac{A_0^2(x) -A_m^2(x)}{(\eta_2 -\eta_1)^2}
   +\dfrac{ C_m^2(x) -C_0^2(x) }{4},  \\
&= A(x)\sum_{j=0}^{m-1} D_j(x) \;
   +\dfrac{ C_m^2(x) -C_0^2(x) }{4}, \text{ if } \eta_1(x)=\eta_2(x),
\end{split}\end{equation}
which shows that the polynomial of degree $\leqslant 2s+4$
\begin{equation} \label{Omega}
  \Omega = A_m^2 -(\eta_2 -\eta_1)^2 \left[ \dfrac{C_m^2}{4} -a_m^2 D_m D_{m-1} \right]
      = A_m^2 -(\eta_2 -\eta_1)^2 \left[ \dfrac{C_m^2}{4} -B_m  D_m  \right]
\end{equation}
does not depend on $m$\dots

\medskip

Summary:

\begin{center}
\begin{tabular}{|c|c|c|}
 \hline
        &    Semi-classical      &   Laguerre-Hahn  \\
        &     $B_0(x)\equiv 0$   &                  \\
 \hline\hline
        &                        &  associate of    \\
  $s=0$ &   classical           & classical if      \\
        &                        & $B_0$ is a constant \\
 \hline
  $s>0$&   General             &  includes transformed \\
       &   semi-classical      &  classical cases \cite{Zhed} \\
 \hline
\end{tabular}
\end{center}


\subsection{Classical orthogonal polynomials from the semi-classical-Laguerre-Hahn
  point of view}
  \  \\

Let us take again the \emph{\textbf{Legendre}} example. Here,
$\eta_1(x)\equiv \eta_2(x) \equiv x$;  $\mathcal{S}' = -2/(z^2-1)$:
$A(z)= A_0(z)= z^2-1$, $B(z)\equiv C(z)\equiv 0$, and $D(z)=-2 \Rightarrow
 D_0(z)= -1$, as $\mu_0 =2$.

\eqref{eqA} $\Rightarrow A_m(x) = x^2-1$ for all $m$, \\
let $C_m(x) = \Gamma_m x+\Gamma'_m$, $D_m = \Delta_m$,
\eqref{eqB} $\Rightarrow B_m = a_m^2 \Delta_m$.
Equations \eqref{eqC} and \eqref{eqD} are \\
$ \Gamma_{m+1} x+\Gamma'_{m+1} = -\Gamma_m x -\Gamma'_m
    -2(x-b_m)\Delta_m, $ \\
$a_{m+1}^2 \Delta_{m+1} = x^2 -1 +(x-b_m)(\Gamma_m x+\Gamma'_m)
 +a_m^2\Delta_{m-1} +(x-b_m)^2 \Delta_m$

\smallskip

\noindent
$x^2:\qquad \Gamma_m + \Gamma_{m+1}= -2\Delta_m \text{ and }
  0=1+\Gamma_m+\Delta_m \Rightarrow \Delta_m = -(2m+1),
  \Gamma_m = 2m,$\\
$x:\ \qquad \Gamma'_m + \Gamma'_{m+1}= -2(2m+1) b_m  \text{ and }
 \Gamma'_m = -2(m+1)b_m \Rightarrow b_m=0$,\\
$x^0:\qquad (2m+3)a_{m+1}^2 -(2m-1)a_m^2 = 1 \Rightarrow
  a_m^2 = m^2/(4m^2-1)$. I told you!

\medskip

\emph{\textbf{Hahn}}: $\eta_1(x)=x$,
      $\eta_2(x)=x+h$. Again, $A_m$ has degree 2, $C_m(x) = \Gamma_m x+ \Gamma'_m$. As
    $D_m$ has degree 0, \eqref{eqA} shows that only the constant coefficient of
    $A_m$, say $\alpha_m$, depends on $m$: let
$A_m(x) = \mathsf{a}x^2 +\mathsf{b} x +\alpha_m$. Then,
     $\alpha_{m+1}-\alpha_m = -h^2 D_m/2$
  from \eqref{eqA};  constant = $2\mathsf{a} \alpha_m -h^2\Gamma_m^2/4$,
i.e., $\alpha_m = h^2 \Gamma_m^2/(8\mathsf{a}) +$ constant, from the
 $x^2$ term of \eqref{Omega}, so,
     $\alpha_{m+1}-\alpha_m = -h^2 D_m/2 = h^2(\Gamma_{m+1}^2
           -\Gamma_m^2)/(8\mathsf{a})$, but, as \eqref{eqC}
  $\Rightarrow \Gamma_m +\Gamma_{m+1} = -2D_m$, we find
 $\Gamma_{m+1} -\Gamma_m = 2\mathsf{a}$, so\\
 $\Gamma_m = 2m\mathsf{a} +\Gamma_0$, and
 $D_m = -(2m+1)\mathsf{a} -\Gamma_0$, $\alpha_m = \alpha_0
                       +\Gamma_0 h^2 m/2 + \mathsf{a}h^2 m^2/2$.

Remark that $\Gamma_m = - D_{m-1/2}$.

Coefficient of $x$ in \eqref{Omega}: $\Omega'(0) = 2\mathsf{b}\alpha_m
 -h^2 \Gamma_m\Gamma'm/2 \Rightarrow
  \Gamma'm = 2(\Omega'(0) -2\mathsf{b}\alpha_m)/(h^2 D_{m-1/2})$,
and \eqref{eqC} at $x=0 \Rightarrow  b_m =
  (\Gamma'_m +\Gamma'_{m+1})/(2D_m) +h/2$.

Finally, $x=0$ in \eqref{Omega} $\Rightarrow
  \displaystyle{ a_m^2 =  \dfrac{ (\alpha_m^2 -\Omega(0) )/h^2
     -(\Gamma'_m)^2/4}{D_m D_{m-1} }   }$.


\bigskip

The full $\mathscr{D}-$classical orthogonal polynomials
(\emph{\textbf{Askey-Wilson orthogonal polynomials}})
will deserve the full next subsection.


\bigskip

\emph{\textbf{A non classical semi-classical example: Freud.}}
$\displaystyle{\mathcal{S}_0(z) = \mu_0^{-1}\int_{-\infty}^\infty (z-x)^{-1}
 e^{-x^4/4}\, dx }  \Rightarrow \mathcal{S}'(z) = -z^3 \mathcal{S}(z)
  +z^2 +a_1^2$. So, $A_0(z)\equiv 1$, $B_0(z)\equiv 0$,
$C_0(z) = -z^3$, and $D_0(z) = z^2+ a_1^2$. Repeated application
of \eqref{eqA}-\eqref{eqD} leads to
$A_m(z)\equiv 1$, $B_m(z) = a_m^2 (z^2 +a_{m-1}^2 +a_m^2)$,
$C_m(z) = -z^3 -2a_m^2z$, and $D_m(z) = z^2 +a_m^2 +a_{m+1}^2$,
\emph{\textbf{and}} $a_m^2( a_{m-1}^2 + a_m^2 + a_{m+1}^2) = m$,
$m=1,2,\dotsc$ ($a_0=0)$.


\subsection{\label{recAW}Recurrence relations  of
the Askey-Wilson polynomials} \   \\

From \eqref{conic}, $\eta_1(x) + \eta_2(x) = 2\mathfrak{p}(x) =
-2(\mathfrak{b} x
       + \mathfrak{d} )/  \mathfrak{a} $,
        $\eta_1(x)  \eta_2(x) = (\mathfrak{c} x^2  +2\mathfrak{e} x
       + \mathfrak{f} )/  \mathfrak{a} $,
$(\eta_2(x) - \eta_1(x))^2 = 4\mathfrak{q}(x) =
4[(\mathfrak{b} x
       + \mathfrak{d} )^2-  \mathfrak{a} (
         \mathfrak{c} x^2  +2\mathfrak{e} x
       + \mathfrak{f} )]/  \mathfrak{a^2} $.

\smallskip

Let us consider first the main coefficients $A_m(x) = \hat{A}_m x^{s+2}+\dotsb$,
 $C_m(x) = \hat{C}_m x^{s+1}+\dotsb$, $D_m(x) = \hat{D}_m x^s+\dotsb$ (this part is
 valid for $s>0$ as well). From \eqref{eqA}-\eqref{eqD},\\
  $\hat{A}_{m+1} = \hat{A}_m -2(\mathfrak{b}^2-\mathfrak{ac}) \hat{D}_m
  /\mathfrak{a}^2$,\\
  $\hat{C}_{m+1} = -\hat{C}_m +2\mathfrak{b} \hat{D}_m /\mathfrak{a}$,\\
  $0 = \hat{A}_m -\mathfrak{b} \hat{C}_m/\mathfrak{a}
   +\mathfrak{c} \hat{D}_m /\mathfrak{a}$.\\

It yields
  $2\mathfrak{b} \hat{C}_{m+1}/\mathfrak{a} =
   \mathfrak{c} (\hat{D}_{m+1} + \hat{D}_m) /\mathfrak{a}$  and
  $2\mathfrak{b} \hat{C}_m/\mathfrak{a} =
   -\mathfrak{c} \hat{D}_{m+1} /\mathfrak{a}
            +(4\mathfrak{b}^2-\mathfrak{ac})\hat{D}_m/\mathfrak{a}^2$,

and so
$$ \mathfrak{ac} \hat{D}_{m+1}
       -(4\mathfrak{b}^2-2\mathfrak{ac})\hat{D}_m
   + \mathfrak{ac} \hat{D}_{m-1}  =0, $$
again solved by a combination of $q^m$ and $q^{-m}$ if
$q\neq 1$, by a combination of 1 and $m$ if $q=1$:

\begin{center}
\begin{tabular}{|c|c|c|}
 \hline
        &    $q\neq 1$           &   $q=1$          \\
 \hline\hline
        &                        &      \\
  $\hat{A}_m $ &   $\dfrac{\mathfrak{c}}{\mathfrak{a}}
                    \dfrac{q^{-1/2} -q^{1/2}}{2}
                   [\alpha q^{m-1/2} -\beta q^{1/2-m} ]$
       &  $\hat{A}_0 = -\dfrac{\mathfrak{c}\alpha}{2\mathfrak{a}}$      \\
        &                        &      \\
 \hline
        &                        &      \\
  $\hat{C}_m$&
             $\sqrt{\dfrac{\mathfrak{c}}{\mathfrak{a}} }
                \hat{D}_{m-1/2}$
     &    $\sqrt{\dfrac{\mathfrak{c}}{\mathfrak{a}} }
                \bigg(\alpha (m-1/2) +\beta \bigg)$ \\
        &                        &      \\
 \hline
        &                        &      \\
  $\hat{D}_m$&
             $\alpha q^m +\beta q^{-m} $
     &   $\alpha m+\beta $           \\
        &                        &      \\
 \hline
\end{tabular}
\end{center}


\smallskip
For classical polynomials ($s=0$), we ``only'' have yet to find
formulas for $C_m(0)=\Gamma'_m$, $a_m^2$ and $b_m$.

From \eqref{eqA}, as $D_m = \hat{D}_m$ is a known constant,
$$ A_m(x) = A_0(x) -2\mathfrak{q}(x)[D_0+ \dotsb +D_{m-1}] =
      A_0(x) -2\mathfrak{q}(x) \left[
       \alpha \dfrac{1-q^m}{1-q} +
       \beta  \dfrac{1-q^{-m}}{1-q^{-1}} \right].
$$
We now need to express $D_{m-1/2}$ as an odd function of $m+\mu$:
 $$ D_{m-1/2} = \mathfrak{D}(m+\mu).$$
This is obvious when $q=1$, $\mathfrak{D}(x) = $ constant $x$, with
$\mu = \beta/\alpha -1/2$. When $q\neq 1$, write
$$ D_{m-1/2} = \alpha q^{-\mu-1/2} q^{m+\mu} +\beta q^{\mu+1/2}
   q^{-m-\mu} = \text{ constant } (q^{m+\mu} -q^{-m-\mu}),$$
provided $q^{2\mu+1} = -\alpha/\beta$.

We then find that $D_0+ \dotsb +D_{m-1}$ is an \emph{even} function
of $m+\mu$:
$$ D_0 +\dotsb + D_{m-1} =
   \dfrac{\alpha}{2}[(m+\mu)^2 -\mu^2] \text{ if } q=1;\;
  \alpha q^{-\mu} \dfrac{ q^\mu +q^{-\mu} -(q^{m+\mu} +q^{-m-\mu}) }
                        {1-q} \text{ if } q\neq 1.
$$

\smallskip

We now look at the
coefficient of $x$ in the expansion of $\Omega/\mathfrak{q}$ about
$x=\infty$, from \eqref{Omega}:
\begin{equation*}
 \left(\dfrac{\Omega}{\mathfrak{q}}\right)_x = \text{ constant } =
 \left(\dfrac{A_m^2}{4\mathfrak{q}}\right)_x -2\hat{C}_m \Gamma'_m,
\end{equation*}
and, as $A_m^2 = A_0^2 -4\mathfrak{q}(D_0+\dotsb +D_{m-1})
                       +4\mathfrak{q}^2(D_0+\dotsb +D_{m-1})^2$,  \\
whence $\displaystyle{ \Gamma'_m = C_m(0) =
 \dfrac{ \text{ constant } - q_1 (D_0+\dotsb +D_{m-1})^2 }
       {2 \hat{C}_m } =
\dfrac{ \text{ a Laurent
  polynomial of degree 2 in } q^m }{\hat{C}_m = \text{ constant }
   D_{m-1/2} }   }$, and, from \eqref{eqC},
$$ b_m = \mathfrak{p}(0) +\dfrac{ C_{m+1}(0) +C_m(0) }{D_m} =
   \dfrac{ \text{ a Laurent
  polynomial of degree 3 in } q^m }{ D_{m-1/2} D_m D_{m+1/2} }.
$$
\textbf{Exercise.} Show that $D_m$ is a factor of the numerator of $b_m$
and can therefore be simplified.

Indeed, $C_m(0) + C_{m+1}(0) = $ odd function of $(m+\mu) +$ odd function
 of $(m+\mu+1) =$ odd function of $(m+\mu+1/2)$.

\medskip
Next,
looking at \eqref{Omega} at a zero $\omega$ of $\Omega$:
$$
 a_m^2 = \dfrac{ C_m^2(\omega) -\dfrac{A_m^2(\omega)}{\mathfrak{q}(\omega)}  }
            {4 D_m D_{m-1} } =
 \dfrac{\text{ Laurent pol. of degree 4 in } q^m}{ D_m D_{m-1/2}^2
   D_{m-1}  }.
$$
Remark that $a_m^2$ is an even function of $m+\mu$.

\smallskip
Let $F_m = [C_m(\omega) +A_m(\omega)/\sqrt{\mathfrak{q}(\omega)}]/(2D_{m-1})$ and
 $G_m = [C_m(\omega) -A_m(\omega)/\sqrt{\mathfrak{q}(\omega)}]/(2D_m)$,
where we manage to have, say, $G_0=0$ (from \eqref{Omega}, as
  $\Omega(\omega)=0$ and $B_0(x)\equiv 0$).

Then, \begin{equation} \label{AWan} a_m^2 = F_m G_m, \end{equation} and
\begin{equation} \label{AWbn}  F_{m+1} +G_m =
   \dfrac{
     C_{m+1}(\omega) +C_m (\omega) + [A_{m+1}(\omega) -A_m(\omega)]
           /\sqrt{\mathfrak{q}(\omega)} } {2D_m} =
 b_m -\mathfrak{p}(\omega) -\sqrt{\mathfrak{q}(\omega)} =
   b_m -\eta(\omega),
\end{equation}
which make four ways (the four zeros of $\Omega$)
 to write $a_m^2$ and $b_m$ as a product $F_m G_m$
and a sum $F_{m+1} +G_m +$ constant.

These four ways are no doubt linked to the four ways to express the
Askey-Wilson polynomials as ${}_4\phi_3$ functions
(Whipple-Sears identities, \cite[p.~6]{AW}).

\medskip

From \cite[etc.]{AW,16,17}, $p_n(x; a,b,c,d |q)$ is symmetric in
 $a, b, c$, and $d$, in
$$
\dfrac{a^n p_n(x; a,b,c,d | q) }{(ab, ac, ad ;q)_n} =       r_n(x) =
 {}_4\phi_3 \left( \left.
    \genfrac{}{}{0pt}{} {q^{-n} , abcdq^{n-1} , a e^{i\theta},
        a e^{-i\theta}  }{ab, ac, ad } \right| q;q \right),\; x=\cos\theta.
$$

$$
2x r_n(x) = \mathsf{A}_n r_{n+1}(x) +[a+a^{-1} -
    (\mathsf{A}_n +\mathsf{C}_n)] r_n(x) +\mathsf{C}_n r_{n-1}(x),
$$
where
$$
  \mathsf{A}_n = \dfrac{ (1-ab q^n) (1-ac q^n) (1-ad q^n) (1-abcd q^n) }
              { a(1-abcd q^{2n-1}) (1-abcd q^{2n}) },
$$
$$
  \mathsf{C}_n = \dfrac{ a(1- q^n) (1-bc q^{n-1}) (1-bd q^{n-1}) (1-cd q^{n-1}) }
              { (1-abcd q^{2n-2}) (1-abcd q^{2n-1}) }.
$$
Here, $b_n =   [a+a^{-1} - (\mathsf{A}_n +\mathsf{C}_n)]/2$, and
 $a_n^2 =  \mathsf{A}_{n-1} \mathsf{C}_n/4$.

These hypergeometric expansions aspects will be considered agaoin in
section \ref{hAW}.


\subsection{Difference relations and equations for
           semi-classical and Laguerre-Hahn
 orthogonal polynomials\label{DLH1}}
  \  \\

Here, we follow Laguerre \cite{19} rather closely:

from \eqref{ric} and \eqref{Pade}, let
 $$ \mathcal{S}(z) =
                         \dfrac{\tilde{p}_{n-1}(z) }{p_n(z) }
             +\varepsilon_n(z),$$
then
\begin{multline*}
   A(x)\left(\mathscr{D}
                   \dfrac{\tilde{p}_{n-1}}{p_n}\right)(x)
          +A(x)(\mathscr{D} \varepsilon_n)(x)
 = B(x)
           \dfrac{\tilde{p}_{n-1}(\eta_1(x)) }{p_n(\eta_1(x)) }
           \dfrac{\tilde{p}_{n-1}(\eta_2(x)) }{p_n(\eta_2(x)) }
 + B(x)
           \dfrac{\tilde{p}_{n-1}(\eta_1(x)) }{p_n(\eta_1(x)) }
            \varepsilon(\eta_2(x))       \\
 + B(x)
           \dfrac{\tilde{p}_{n-1}(\eta_2(x)) }{p_n(\eta_2(x)) }
            \varepsilon(\eta_1(x))
 + B(x)
            \varepsilon(\eta_1(x))
            \varepsilon(\eta_2(x))
+C(x) \left(\mathscr{M}
                   \dfrac{\tilde{p}_{n-1}}{p_n}\right)(x)
          +C(x)(\mathscr{M} \varepsilon_n)(x)
 +D(x).
\end{multline*}

So,
\begin{multline}    \label{Thetan}
 \Theta_n(x) =
 A(x)\dfrac{
      \tilde{p}_{n-1}(\eta_2(x)) p_n(\eta_1(x))
     - \tilde{p}_{n-1}(\eta_1(x)) p_n(\eta_2(x))
     }
     {\eta_2(x) -\eta_1(x)}
 - B(x)
           \tilde{p}_{n-1}(\eta_1(x))
           \tilde{p}_{n-1}(\eta_2(x))         \\
   -C(x)\dfrac{
      \tilde{p}_{n-1}(\eta_2(x)) p_n(\eta_1(x))
     + \tilde{p}_{n-1}(\eta_1(x)) p_n(\eta_2(x))
     }{2}
   -D(x)
       p_n(\eta_1(x))   p_n(\eta_2(x))
\end{multline}
is \  1.) a polynomial (symmetry in $\eta_1$ and $\eta_2$)\qquad 2.)
of degree $\leqslant s$ (as $\varepsilon_n(x) \sim x^{-2n-1}$).

\medskip

This $\Theta_n$ is not a completely new actor. To see things further,
remark that \eqref{Thetan} is a combination of products of solutions of
the recurrence \eqref{recur} at $\eta_1$ and at $\eta_2$.

\smallskip

\emph{We now consider the \textbf{linear} recurrence satisfied by
combinations of such products}, i.e., by combinations of
$$ p_n(\eta_1) p_n(\eta_2),
   p_n(\eta_1) \tilde{p}_{n-1}(\eta_2),
    \tilde{p}_{n-1}(\eta_1) p_n(\eta_2),
   \text{ and }
    \tilde{p}_{n-1}(\eta_1) \tilde{p}_{n-1}(\eta_2) . $$
We just have to consider a product $r_n s_n$, kowing that
$a_{n+1} r_{n+1} = (\eta_1 -b_n) r_n -a_n r_{n-1}$, and
$a_{n+1} s_{n+1} = (\eta_2 -b_n) s_n -a_n s_{n-1}$.
Simplest is to write it as matrix-vector recurrence (\emph{doesn't it
 ring a bell}?)

\def\vrz{\rule{0cm}{1cm}}

\begin{equation*}
  \begin{bmatrix}
 \vrz r_n s_{n+1} \\ \vrz r_{n+1} s_n  \\ \vrz r_n s_n \\ \vrz r_{n+1} s_{n+1}
  \end{bmatrix}
    =
  \begin{bmatrix}
   \vrz     & -\dfrac{a_n}{a_{n+1}}   &   &  \dfrac{\eta_2 -b_n}{a_{n+1}}\\
   \vrz   -\dfrac{a_n}{a_{n+1}}   & &    &  \dfrac{\eta_1 -b_n}{a_{n+1}}  \\
   \vrz             &         &        &    1   \\
   \vrz         -\dfrac{a_n(\eta_2-b_n)}{a_{n+1}^2}  &
            -\dfrac{a_n(\eta_1-b_n)}{a_{n+1}^2}  &
             \dfrac{a_n^2}{a_{n+1}^2}  &
            \dfrac{(\eta_1-b_n)(\eta_2-b_n)}{a_{n+1}^2}
  \end{bmatrix}
  \begin{bmatrix}
 \vrz  r_{n-1} s_n \\ \vrz r_n s_{n-1}  \\ \vrz r_{n-1} s_{n-1} \\ \vrz r_n s_n
  \end{bmatrix}
\end{equation*}





Looks damned close to \eqref{recmatrix}, needs just a little bit of
scaling:


\begin{equation*}
  \begin{bmatrix}
  a_{n+1} r_n s_{n+1} \\  a_{n+1} r_{n+1} s_n  \\
  a_{n+1}^2 r_n s_n \\ \vrz r_{n+1} s_{n+1}
  \end{bmatrix}
    =
  \begin{bmatrix}
        & -1   &   &  \eta_2 -b_n\\
      -1   & &    &  \eta_1 -b_n  \\
                &         &        &    a_{n+1}^2   \\
   \vrz         -\dfrac{\eta_2-b_n}{a_{n+1}^2}  &
            -\dfrac{\eta_1-b_n}{a_{n+1}^2}  &
             \dfrac{1}{a_{n+1}^2}  &
            \dfrac{(\eta_1-b_n)(\eta_2-b_n)}{a_{n+1}^2}
  \end{bmatrix}
  \begin{bmatrix}
   a_n r_{n-1} s_n \\  a_n r_n s_{n-1}  \\
      a_n^2 r_{n-1} s_{n-1} \\ \vrz r_n s_n
  \end{bmatrix}
\end{equation*}

So, each vector
$\displaystyle{  \begin{bmatrix}
    -\dfrac{A_{n}(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{C_{n}(x)}{2} \\
    \dfrac{A_{n}(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{C_{n}(x)}{2} \\
        B_{n}(x) \\  D_{n}(x)
  \end{bmatrix}    }$
will be a fixed (independent of $n$) combination of the four vectors
$$
\begin{bmatrix}
 a_n p_{n-1}(\eta_1) \tilde{p}_{n-1}(\eta_2) \\
 a_n p_{n}(\eta_1) \tilde{p}_{n-2}(\eta_2) \\
 a_n^2 p_{n-1}(\eta_1) \tilde{p}_{n-2}(\eta_2) \\
     p_{n}(\eta_1) \tilde{p}_{n-1}(\eta_2)
\end{bmatrix}   ,
\begin{bmatrix}
 a_n  \tilde{p}_{n-2}(\eta_1) p_n(\eta_2) \\
 a_n  \tilde{p}_{n-1}(\eta_1) p_{n-1}(\eta_2) \\
 a_n^2  \tilde{p}_{n-2}(\eta_1) p_{n-1}(\eta_2) \\
      \tilde{p}_{n-1}(\eta_1) p_n(\eta_2)
\end{bmatrix}   ,
\begin{bmatrix}
 a_n  \tilde{p}_{n-2}(\eta_1) \tilde{p}_{n-1}(\eta_2) \\
 a_n  \tilde{p}_{n-1}(\eta_1) \tilde{p}_{n-2}(\eta_2) \\
 a_n^2  \tilde{p}_{n-2}(\eta_1) \tilde{p}_{n-2}(\eta_2) \\
      \tilde{p}_{n-1}(\eta_1) \tilde{p}_{n-1}(\eta_2)
\end{bmatrix}   , \text{ and }
\begin{bmatrix}
 a_n  p_{n-1}(\eta_1) p_n(\eta_2) \\
 a_n  p_{n}(\eta_1) p_{n-1}(\eta_2) \\
 a_n^2  p_{n-1}(\eta_1) p_{n-1}(\eta_2) \\
      p_{n}(\eta_1) p_n(\eta_2)
\end{bmatrix}   ,
$$
which are $\displaystyle{
\begin{bmatrix}   0 \\  -1 \\ 0  \\ 0 \end{bmatrix}   ,
\begin{bmatrix}   -1 \\  0 \\ 0  \\ 0 \end{bmatrix}   ,
\begin{bmatrix}   0 \\  0 \\ \mu_0  \\ 0 \end{bmatrix}   ,
\begin{bmatrix}   0 \\  0 \\  0      \\  1/\mu_0  \end{bmatrix}
}$ at $n=0$ (from \eqref{recurinit}), so that

{\Small
\begin{equation} \label{ABCD} \hspace*{-1cm} \begin{split}
    -\dfrac{A_{n}(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{C_{n}(x)}{2} &=
  \alpha
         a_n p_{n-1}(\eta_1) \tilde{p}_{n-1}(\eta_2)
 +\beta
         a_n \tilde{p}_{n-2}(\eta_1) p_n(\eta_2)
  +\gamma
         a_n \tilde{p}_{n-2}(\eta_1) \tilde{p}_{n-1}(\eta_2)
  +\delta
         a_n p_{n-1}(\eta_1) p_n(\eta_2)     \\
    \dfrac{A_{n}(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{C_{n}(x)}{2}   &=
  \alpha
         a_n p_n(\eta_1) \tilde{p}_{n-2}(\eta_2)
 +\beta
         a_n \tilde{p}_{n-1}(\eta_1) p_{n-1}(\eta_2)
  +\gamma
         a_n \tilde{p}_{n-1}(\eta_1) \tilde{p}_{n-2}(\eta_2)
  +\delta
         a_n p_{n}(\eta_1) p_{n-1}(\eta_2)
  \\
        B_{n}(x)    &=
  \alpha
         a_n^2 p_{n-1}(\eta_1) \tilde{p}_{n-2}(\eta_2)
 +\beta
         a_n^2 \tilde{p}_{n-2}(\eta_1) p_{n-1}(\eta_2)
  +\gamma
         a_n^2 \tilde{p}_{n-2}(\eta_1) \tilde{p}_{n-2}(\eta_2)
  +\delta
         a_n^2 p_{n-1}(\eta_1) p_{n-1}(\eta_2)
  \\
          D_{n}(x) &=
  \alpha
          p_{n}(\eta_1) \tilde{p}_{n-1}(\eta_2)
 +\beta
          \tilde{p}_{n-1}(\eta_1) p_{n}(\eta_2)
  +\gamma
          \tilde{p}_{n-1}(\eta_1) \tilde{p}_{n-1}(\eta_2)
  +\delta
          p_{n}(\eta_1) p_{n}(\eta_2)
\end{split}\end{equation}
}

with $\displaystyle{
  \alpha=-\dfrac{A_{0}(x)}{\eta_2-\eta_1} +\dfrac{C_{0}(x)}{2},
 \beta = \dfrac{A_{0}(x)}{\eta_2-\eta_1} +\dfrac{C_{0}(x)}{2},
  \gamma= \dfrac{B_0(x)}{\mu_0},
  \delta = \mu_0 D_0(x)  }$.

\smallskip

One may now check that Laguerre's $\Theta_n = -D_n$,
\emph{\textbf{but the truly
useful outcome of \eqref{ABCD} is the production of relations
linking $p_n$ and $p_{n-1}$.}}

Trying to eliminate as many $\tilde{p}$'s as possible, let us take
the first equation of \eqref{ABCD} times $p_n(\eta_1)$ minus the last
equation times $a_n p_{n-1}(\eta_1)$: the contributions of $\alpha$ and
$\delta$ disappear, and
$$
    \left(-\dfrac{A_{n}(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{C_{n}(x)}{2}
    \right) p_n(\eta_1)  -a_n D_n(x) p_{n-1}(\eta_1) =
 \beta a_n X_n p_n(\eta_2)
 +\gamma a_n X_n \tilde{p}_{n-1}(\eta_2)  , $$
with $X_n = \tilde{p}_{n-2}(\eta_1) p_n(\eta_1) -
            \tilde{p}_{n-1}(\eta_1) p_{n-1}(\eta_1)$,
whence the formidable simplification $a_n X_n =-1$, thanks to
\eqref{casor}, and
\begin{equation} \label{pp1}   \hspace*{-8mm}
    \left(-\dfrac{A_{n}(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{C_{n}(x)}{2}
    \right) p_n(\eta_1)  -a_n D_n(x) p_{n-1}(\eta_1) =
   -\left(\dfrac{A_{0}(x)}{\eta_2-\eta_1} +\dfrac{C_{0}(x)}{2}\right)
   p_n(\eta_2)
 -  \dfrac{B_0(x)}{\mu_0}
  \tilde{p}_{n-1}(\eta_2)  .
\end{equation}

Similarly with the second and last equations of \eqref{ABCD}, and with
$\eta_2$ instead of $\eta_1$:
\begin{equation} \label{pp2}  \hspace*{-8mm}
    \left(\dfrac{A_{n}(x)}{\eta_2(x)-\eta_1(x)} -\dfrac{C_{n}(x)}{2}
    \right) p_n(\eta_2)  -a_n D_n(x) p_{n-1}(\eta_2) =
   \left(\dfrac{A_{0}(x)}{\eta_2-\eta_1} -\dfrac{C_{0}(x)}{2}\right)
   p_n(\eta_1)
 -  \dfrac{B_0(x)}{\mu_0}
  \tilde{p}_{n-1}(\eta_1)  .
\end{equation}


\subsection{Difference relations \& equations: the importance of
being semi-classical} \   \\

In the semi-classical case, $B_0=0$, \eqref{pp1}-\eqref{pp2} do not
contain $\tilde{p}_{n-1}$ any more. Summing \eqref{pp1} and
\eqref{pp2} yields a good looking \emph{\textbf{difference
relation}}
\begin{equation} \label{diffrec}
 (A_n(x)+ A_0(x)) (\mathscr{D} p_n)(x)=
 (C_n(x)- C_0(x)) (\mathscr{M} p_n)(x)
 +2a_n D_n(x) (\mathscr{M} p_{n-1})(x),
\end{equation}
but the true useful meaning of \eqref{pp1}-\eqref{pp2} is to show how
$[p_n(\eta_1), p_n(\eta_2)]$ can be determined from
$[p_{n-1}(\eta_1), p_{n-1}(\eta_2)]$, er, oops, it seems that it is
more obvious to get
$[p_{n-1}(\eta_1), p_{n-1}(\eta_2)]$ from
$[p_n(\eta_1), p_n(\eta_2)]$. \\
Ah, this is the subject of a little \emph{\textbf{exercise}}: from
the recurrence relation \eqref{recur},
$-a_n p_{n-1}(\eta) = -(\eta -b_n) p_n(\eta) +a_{n+1} p_{n+1}(\eta)$,
and one can establish, using \eqref{recmatrix}
\begin{equation*} \begin{split}
    \left(-\dfrac{A_{n+1}(x)}{\eta_2(x)-\eta_1(x)} +\dfrac{C_{n+1}(x)}{2}
    \right) p_n(\eta_1)  +a_{n+1} D_n(x) p_{n+1}(\eta_1) &=
   -\left(\dfrac{A_{0}(x)}{\eta_2-\eta_1} +\dfrac{C_{0}(x)}{2}\right)
   p_n(\eta_2) , \\
    \left(\dfrac{A_{n+1}(x)}{\eta_2(x)-\eta_1(x)} +\dfrac{C_{n+1}(x)}{2}
    \right) p_n(\eta_2)  +a_{n+1} D_n(x) p_{n+1}(\eta_2) &=
   \left(\dfrac{A_{0}(x)}{\eta_2-\eta_1} -\dfrac{C_{0}(x)}{2}\right)
   p_n(\eta_1),
\end{split}\end{equation*}
so, $\displaystyle{  a_n D_{n-1}(x)
   \begin{bmatrix}  p_n(\eta_1) \\ p_n(\eta_2) \end{bmatrix} =
   \begin{bmatrix} A_n/\Delta\eta -C_n/2  &  -A_0/\Delta\eta -C_0/2 \\
                   A_0/\Delta\eta-C_0/2 &  -A_n/\Delta\eta -C_n/2
    \end{bmatrix}
   \begin{bmatrix}  p_{n-1}(\eta_1) \\ p_{n-1}(\eta_2) \end{bmatrix}
 }$.

But we also have, from \eqref{pp1}-\eqref{pp2} (remember that $B_0=0$)
$$
   \begin{bmatrix}  p_n(\eta_1) \\ p_n(\eta_2) \end{bmatrix} =
   \begin{bmatrix} -A_n/\Delta\eta -C_n/2  &  A_0/\Delta\eta +C_0/2 \\
                  -A_0/\Delta\eta+C_0/2 &  A_n/\Delta\eta -C_n/2
    \end{bmatrix}^{-1}  a_n D_n(x)
   \begin{bmatrix}  p_{n-1}(\eta_1) \\ p_{n-1}(\eta_2) \end{bmatrix}
$$
the two matrices are found to be the same (use \eqref{Omega}).

\bigskip

Most dramatic consequence of \eqref{pp1}-\eqref{pp2}:

\emph{\textbf{The second order difference equation.}}

Take \eqref{pp1} at $x=x_k$ and \eqref{pp2} at $x=x_{k-1}$, and
eliminate $p_{n-1}(y_k)$ (as $\eta_1(x_k) = \eta_2(x_{k-1}) = x_k$):

\begin{multline} \label{eqdiff}
   D_n(x_{k-1}) \left(  \dfrac{A_0(x_k)}{y_{k+1}-y_k}
                     +\dfrac{C_0(x_k)}{2} \right) p_n(y_{k+1}) \\
 -\left[
   D_n(x_{k-1}) \left(  \dfrac{A_n(x_k)}{y_{k+1}-y_k}
                     +\dfrac{C_n(x_k)}{2} \right)
  + D_n(x_k) \left(  \dfrac{A_n(x_{k-1})}{y_k-y_{k-1}}
                     -\dfrac{C_n(x_{k-1})}{2} \right)
  \right] p_n(y_k) \\
  + D_n(x_k) \left(  \dfrac{A_0(x_{k-1})}{y_k-y_{k-1}}
                     -\dfrac{C_0(x_{k-1})}{2} \right)
   p_n(y_{k-1}) =0.
\end{multline}


\subsection{Difference relations, difference equations, and hypergeometric
expansions of the Askey-Wilson polynomials} \  \label{hAW} \\

We are at last able now to recover the difference equation of the
$\mathscr{D}-$classical orthogonal polynomials (Nikiforov \& \textit{al}.
starting point \cite{24,NS,25}),
 as well as the ${}_4\phi_3$ expansion (Askey \& Wilson
starting point \cite{AW}).

\smallskip

First, the difference equation: when $s=($degree $D_n) = 0$ drop the
same constant $D_n$ from the three terms of \eqref{eqdiff}.
The dependence in $n$ is quite interesting. Remark already that $n$
appears only in the middle term.
 As $C_n$ is
only of degree 1, $C_n(x_k) - C_n(x_{k-1}) = \hat{C}_n (x_k-x_{k-1})$
with the $\hat{C}_n$ of \S~\ref{recAW}. And $A_n(x_p) =
A_0(x_p) +((y_{p+1}- y_p)^2/2) (D_0 + \dotsb + D_{n-1})$ from
\eqref{eqA}, whence the form
 \begin{equation} \label{eqdiffl}
  \mathscr{L} p_n =\lambda_n p_n,
 \end{equation}
with
\begin{multline*}
\mathscr{L}p_n(y_k) =  \dfrac{1}{x_k - x_{k-1} }
   \left(  A_0(x_k)\dfrac{p_n(y_{k+1})-p_n(y_k)}{y_{k+1}-y_k}
         -A_0(x_{k-1}) \dfrac{p_n(y_k) -p_n(y_{k-1})}{y_k-y_{k-1}} \right.   \\
 \left.  +C_0(x_k)\dfrac{p_n(y_{k+1})+p_n(y_k)}{2}
         -C_0(x_{k-1}) \dfrac{p_n(y_k) +p_n(y_{k-1})}{2} \right),
\end{multline*}
i.e.,
\begin{equation}     \label{eqdiffl2}
   \mathscr{L} = \widetilde{\mathscr{D}} A_0  \mathscr{D}
                +\widetilde{\mathscr{D}} C_0  \mathscr{M} ,
 \end{equation}
and
$$
   \lambda_n = -\dfrac{ \mathfrak{b} }{\mathfrak{a} }
      (D_0 +\dotsb + D_{n-1}) +\dfrac{ \hat{C}_n +\hat{C}_0}{2}
$$
(use $y_{k+1} -y_k +(y_k -y_{k-1}) = y_{k+1}+y_k -(y_k + y_{k-1}) =
 (-2\mathfrak{b}/\mathfrak{a}) ( x_k-x_{k-1})$).

\medskip

To explore the \emph{\textbf{hypergeometric expansion}} of $p_n$, look
at the recurrence relation \eqref{recur} (for the monic polynomials
$p_n/\gamma_n$), when \eqref{AWan}-\eqref{AWbn} hold:
$$
p_{n+1}(x)/\gamma_{n+1} = (x-y_1 -F_{n+1} -G_n) p_n(x)/\gamma_n
  -F_n G_n p_{n-1}(x)/\gamma_{n-1},$$
where $x_0=\omega$, $y_1 = \eta(\omega) = \eta_2(\omega)$,
or
$$
p_{n+1}(x)/\gamma_{n+1} +F_{n+1}  p_n(x)/\gamma_n
 +G_n[ p_n(x)/\gamma_n  +F_n  p_{n-1}(x)/\gamma_{n-1}] =
 (x-y_1) p_n(x)/\gamma_n,
$$
showing that $p_n(y_1)/\gamma_n = (-1)^n F_1 \dotsb F_n$, as
$G_n=0$.

Now, \eqref{eqdiffl}-\eqref{eqdiffl2} at $y_1$, i.e., \eqref{eqdiff} at $k=1$:
\begin{equation*} \begin{split}
\left( \dfrac{A_0(x_1)}{y_2-y_1} +\dfrac{C_0(x_1)}{2} \right) p_n(y_2) &=
\left( \dfrac{A_0(x_1)}{y_2-y_1} +
       \dfrac{A_0(x_0)}{y_1-y_0}
       -\dfrac{C_0(x_1) -C_0(x_0)}{2}
          + (x_1-x_0)\lambda_n
 \right) p_n(y_1)  \\
&= \left( \dfrac{A_0(x_1)}{y_2-y_1}
     +C_0(x_0)  -\dfrac{C_0(x_1)}{2} +  (x_1-x_0)\lambda_n
 \right) p_n(y_1)  \\
&= \left( \dfrac{A_0(x_1)}{y_2-y_1}
      + \dfrac{C_0(x_1)}{2} +  (x_1-x_0)(\lambda_n -\hat{C}_0 )
 \right) p_n(y_1),
\end{split} \end{equation*}
having used $2A_0(x_0) = (y_1-y_0)C_0(x_0)$, so,
$$
p_n(y_2) = \left[ 1+
     \dfrac{    (x_1-x_0)(\lambda_n -\hat{C}_0 )  (y_2-y_1) }
           {    A_0(x_1) +(y_2-y_1) C_0(x_1)/2  }
           \right] p_n(y_1)$$
preparing an expansion $p_n(x) = p_n(y_1)[ 1 + c_1 (x-y_1)
 + c_2 (x-y_1)(x-y_2) +\dotsb ]$.

\subsection{Difference equation for Laguerre-Hahn orthogonal
polynomials\label{DLH2}} \   \\

When $B_0(x)\not\equiv 0$, there is no more a second order difference
equation for $p_n$, but one can still do something with the solutions
of \eqref{eqdiff} (Hahn).

Let $\varphi_n$ and $\psi_n$ be two independent solutions of \eqref{eqdiff},
i.e.,

\begin{multline} \label{eqdiffphi}
   D_n(x_{k-1}) \left(  \dfrac{A_0(x_k)}{y_{k+1}-y_k}
                     +\dfrac{C_0(x_k)}{2} \right) \varphi_n(y_{k+1}) \\
 -\left[
   D_n(x_{k-1}) \left(  \dfrac{A_n(x_k)}{y_{k+1}-y_k}
                     +\dfrac{C_n(x_k)}{2} \right)
  + D_n(x_k) \left(  \dfrac{A_n(x_{k-1})}{y_k-y_{k-1}}
                     -\dfrac{C_n(x_{k-1})}{2} \right)
  \right] \varphi_n(y_k) \\
  + D_n(x_k) \left(  \dfrac{A_0(x_{k-1})}{y_k-y_{k-1}}
                     -\dfrac{C_0(x_{k-1})}{2} \right)
   \varphi_n(y_{k-1}) =0,
\end{multline}
and the same for $\psi_n$. Now, nor $\varphi_n$ nor $\psi_n$ need to
be a polynomial. The index $n$ is no more related to a degree of any
kind.

\medskip To be continued\dots

\section{New directions?}
\renewcommand{\mylmark}{New directions?}
\pagestyle{fancy}

\subsection{Chasles-Euler-Poncelet} \   \\

``\textbf{Poncelet Transverse:}   Let a Circle $C_1$ lie inside another
 Circle $C2$. From any point on $C_2$,
 draw a tangent to $C_1$ and extend it to
   $C_2$. From the point, draw another tangent, etc.
      For $n$ tangents, the result is called an $n-$sided Poncelet Transform.

  References      D\"orrie, H.
  \textsl{100 Great Problems of Elementary Mathematics: Their History
   and Solutions}. New York: Dover, p. 192, 1965.''

   \copyright 1996-7 Eric W. Weisstein
    http://www.astro.virginia.edu/\~{}eww6n/math/PonceletTransverse.html

\medskip
For two general conics $\mathscr{C}_1$ and $\mathscr{C}_2$, take a point
$Q_1$ on $\mathscr{C}_2$, call $Q_0$ and $Q_2$ its two Poncelet transforms
defining also $P_1$ and $P_2$ on $\mathscr{C}_1$, etc. We so have two lattices
on the two conics. There is a theorem by Poncelet stating that if the
polygonal line closes on itself after a finite number of steps for some
$Q_1$, the it does so for any starting point on $\mathscr{C}_2$. There must be some
$q-$geometric progression somewhere\dots
\psset{unit=1cm}
\begin{pspicture}(2,-1)(17,7)
%\psline(0,0)(8,0)(8,6)(8,0)(16,0)(16,6)(0,6)(0,0)
\psline(2,0)(2,4.7)    \uput[270](2,0){$Q_1$}   \uput[90](2,4.7){$Q_0$}
    \psline(2,0)(15.7,5.2) \uput[0](15.3,5.2){$Q_2$}
 \uput[180](2,2){$P_1$}   \uput[300](5.8,1.5){$P_2$}
%\psline{->}(0,3)(7,3)  \uput[90](6.9,3.2){$x$}
%\psline{->}(3,0)(3,6.5) \uput[0](3.2,6.2){$y$}
%\psline(0,0)(6,0)(6,6)(0,6)(0,0)
%\put(2.5,0){\line(1,0){0.5}}
%\psline{->}(8,3)(15,3)  \uput[90](14.9,3.2){$x$}
%\psline{->}(11,0)(11,6)  \uput[0](11.2,5.8){$y$}
%\psline(9.5,1.5)(12.5,1.5)(12.5,4.5)(9.5,4.5)(9.5,1.5)
%\psline[linestyle=dashed,dash=0.2 0.2](8,1.5)(14,1.5)
%\psline[linestyle=dashed](8,4.5)(14,4.5)
%\psline[linestyle=dashed](9.5,0)(9.5,6)
%\psline[linestyle=dashed](12.5,0)(12.5,6)
 \uput[180](8,4){$\mathscr{C}_1$}  \uput[150](12,0.4){$\mathscr{C}_2$}
 \psset{linewidth=0.075}
 \psccurve[linecolor=red](2,0)(0,2)(14,6)(16,4)
 \psccurve[linecolor=green](4,1)(2,2)(6,5)(8,4)
\end{pspicture}


As one can always project one of the conics on a circle, we could have a
clue about $\mathscr{D}-$orthogonal polynomials on the circle.
Askey-Wilson polynomials on the circle have been considered in \cite{zhed1}.


\section*{ }
\renewcommand{\mylmark}{References.}
\pagestyle{fancy}

\begin{thebibliography}{99}

\bibitem{1} W.A.\ Al-Salam, Characterization theorems for orthogonal
    polynomials, pp. 1-24 {\it in\/} {\sl Orthogonal Polynomials:
    Theory and Practice\/}, (P.\ Nevai, editor) {\sl NATO ASI Series
   C: Math.\ and Phys.\ Sciences\/} {\bf 294\/}, Kluwer, Dordrecht, 1990.

\bibitem{AW} R.~Askey, J.~Wilson, Some basic hypergeometric
      orthogonal polynomials that generalize Jacobi polynomials,
      {\sl Mem.\ A.\ M.\ S.\/} {\bf 54\/} no.\ 319, 1985.
pp. 1-55.

\bibitem{Ban} G.~Bangerezako,
 Discrete Darboux transformation for discrete polynomials
of hypergeometric type.
     Submitted.

\bibitem{beardon} A.F.~Beardon, Sums of powers of integers,
    \textsl{Amer.\ Math.\ Monthly} \textbf{103} (1996) 201-213.

\bibitem{4}  S.~Bonan, P.~Nevai, Orthogonal polynomials and their
       derivatives,I, {\sl J. Approx. Theory\/} {\bf 40\/} (1984),
       134-147.

\bibitem{5}  S.~Bonan, D.S.~Lubinsky, P.~Nevai, Orthogonal polynomials
 and their
       derivatives,II, {\sl SIAM J. Math. An.\/} {\bf 18\/} (1987),
       1163-1176.


\bibitem{BI1} B.M.~Brown, M.E.H.~Ismail, A right inverse of the
     Askey-Wilson operator, \textsl{Proc.\ AMS} \textbf{123}
    (1995) 2071-2079.

\bibitem{BI2} B.M.~Brown, W.D.~Evans, M.E.H.~Ismail, The
     Askey-Wilson polynomials and $q$-Sturm-Liouville problems,
     \textsl{Math.\ Proc.\ Camb.\ Philos.\ Soc.} \textbf{119}
    (1996) 1-16.

\bibitem{GH} F.~Alberto Gr\"unbaum, L.~Haine, Orthogonal polynomials
  satisfying differential equations: the role of the Darboux
  transformation, pp.~143-154 \emph{in} \textsl{Symmetries and
  Integrability of Difference Equations} (D.~Levi, L.~Vinet, and
  P.~Winternitz, editors), \textsl{CRM (Centre de Recherches
  Math\'ematiques Univ.\ Montr\'eal) Proceedings \& Lecture Notes}
  \textbf{9} (1996), AMS, Providence.

\bibitem{dibuc} A.~Di~Bucchianico, D.~Loeb, A selected survey of
   umbral calculus, \textsl{The Electronic Journal of Combinatorics}
   \textbf{2} (1995), \#DS3,
    \texttt{
    http://www.win.tue.nl/win/math/bs/statistics/bucchianico/hypersyrvey.html}

\bibitem{ha1} W.~Hahn, On differential equations for orthogonal
  polynomials, \textsl{Funk.\ Ekv.} \textbf{21} (1978), 1-9.

\bibitem{ha2} W.~Hahn, \"Uber Orthogonalpolynome, die linearen
  Functionalgleichungen gen\"ugen, pp.~16-35 \emph{in}
  \textsl{Polyn\^omes Orthogonaux et Applications, Proceedings,
  Bar-le-Duc 1984},
(C.BREZINSKI , A.RONVEAUX , A.DRAUX , A.P.MAGNUS , P.MARONI editors),
\textsl{Lecture Notes in Mathematics} {\bf 1171\/}, Springer, Berlin 1985.


\bibitem{10}  E.~Hendriksen, H.~van Rossum, Semi-classical orthogonal
     polynomials,
                                       pp. 354-361 {\it in\/}
      {\sl Polyn\^omes Orthogonaux et Applications, Proceedings,
      Bar-le-Duc 1984\/}, (C.Brezinski \& al., editors),
      {\sl Lecture Notes Math.\/} {\bf 1171\/}, Springer, Berlin 1985.

\bibitem{I} M.E.H.~Ismail, The Askey-Wilson operator and
   summation theorems, pp.~171-178 \emph{in}
  M.E.H.~Ismail \emph{et al.} (ed.), \textsl{Mathematical Analysis,
  Wavelets, and Signal Processing. An International Conference on
  Mathematical Analysis and Signal Processing, Cairo University,
   Cairo, Egypt, January 3-9, 1994}, Contemp.\ Math.\ \textbf{190}, AMS,
   Providence, 1995.

\bibitem{IZ}
     M.E.H.~Ismail, R.~Zhang, Diagonalization
     of certain integral operators, \textsl{Adv.\ Math.}
     \textbf{109} (1994) 1-33.

\bibitem{IRZ}
      M.E.H.~Ismail, M.~Rahman, R.~Zhang,
     Diagonalization of certain integral operators II. Preprint.
To appear in J. Comp. Appl. Math (1996).


\bibitem{16}  R.~Koekoek, R.F.~Swarttouw, The Askey-scheme of hypergeometric
      polynomials and its $q-$analogue,
      Report Fac. of Technical Math. and Informatics 94-05 T.U.\ Delft
      (1994) pp.1--120.
    See also the ``Askey-Wilson
   scheme project'' entry in http://www.cs.vu.nl/\~{}rene/ and also
    \texttt{ftp://unvie6.un.or.at/siam/opsf\_new/koekoek\_swarttouw/koekoek\_swarttouw1.ps}
    and \emph{ibid.}\texttt{
    /koekoek\_swarttouw2.ps}

\bibitem{17} T.H.~Koornwinder,
  {\sl Compact quantum groups and $q$-special functions},
  in {\sl Representations of Lie groups and quantum groups},
  V. Baldoni \& M. A. Picardello (eds.),
  Pitman Research Notes in Mathematics Series 311,
  Longman Scientific \& Technical, 1994,
pp. 46--128.
 \\ See also the part  ``q-Special functions, a tutorial'', in
      \texttt{  ftp://unvie6.un.or.at/siam/opsf\_new/koornwinder/
      koornwinder3/koornwinder3.tex}
   and \emph{ibid.}\texttt{/wilsonqhahn.ps}

\bibitem{19}  E.~Laguerre, Sur la r\'eduction en fractions continues d'une
      fraction qui satisfait \`a une \'equation diff\'erentielle lin\'eaire
      du premier ordre dont les coefficients sont rationnels,
      {\sl J. Math. Pures Appl. (4)\/} {\bf 1\/} (1885), 135-165 =
      pp. 685-711 {\it in\/}
      {\sl Oeuvres\/}, Vol.II, Chelsea, New-York 1972.

\bibitem{lehmer} D.H.~Lehmer, A new approach to Bernoulli
   polynomials, \textsl{Amer.\ Math.\ Monthly} \textbf{95} (1988)
   905-911.

\bibitem{21} A.P.~Magnus, Associated Askey-Wilson polynomials as
   Laguerre-Hahn orthogonal polynomials, pp.~261-278 {\it in\/}
   M.~Alfaro {\it et al.\/}, editors: {\sl Orthogonal Polynomials and
   their Applications, Proceedings, Segovia 1986. Springer Lecture
   Notes Math.\/} {\bf 1329\/}, Springer, Berlin, 1988.

\bibitem{magsnul} A.P.~Magnus,
    Special non uniform lattice ($snul$) orthogonal polynomials
    on discrete dense sets of points,
    \textsl{J.\ Comp.\ Appl.\ Math.} \textbf{65} (1995) 253-265.

\bibitem{22} P.\ Maroni, Une caract\'erisation des  polyn\^omes
  orthogonaux semi-classiques, {\sl C.R.\ Acad.\ Sci.\ Paris ser.\ \/}
   1 {\bf 301\/} (1985) 269-272.


\bibitem{Mi} L.M.~Milne-Thomson, \textsl{The Calculus of Finite
     Differences}, Macmillan, London, 1933.

\bibitem{24} A.F.\ Nikiforov, V.B.\ Uvarov, {\sl Classical orthogonal
  polynomials of a discrete variable on nonuniform lattices\/}
 (in Russian), Keldysh Institute of Applied Mathematics, Preprint
  \#~17, Moscow, 1983.

\bibitem{NS} A.F.~Nikiforov, S.K.~Suslov, Classical orthogonal
   polynomials of a discrete variable on nonuniform lattices,
   \textsl{Letters in Math.\ Phys.} \textbf{11} (1986) 27-34.

\bibitem{25} A.F.~Nikiforov, S.K.~Suslov, V.B.~Uvarov, {\sl Classical
    Orthogonal Polynomials of a Discrete Variable\/}, Springer, Berlin,
    1991.

\bibitem {26} O.Perron, {\sl Die Lehre von den Kettenbr\"uchen\/},
        $2^{\rm nd}$ edition, Teubner, Leipzig, 1929 = Chelsea,

\bibitem{29} M.\ Rahman, S.K.\ Suslov, The Pearson equation and the beta
    integral, {\sl SIAM J.\ Math.\ Anal.\ \/} {\bf 25\/} (1994) 646-693.

\bibitem{30} M.\ Rahman, S.K.\ Suslov, Barnes and Ramanujan-type integrals
    on the $q-$linear lattice,
    {\sl SIAM J.\ Math.\ Anal.\ \/} {\bf 25\/} (1994) 1002-1022.

\bibitem{Smi} D.E.~Smith, \textsl{History of Mathematics, {\Small
   the evolution of arithmetic, geometry, trigonometry, calculating
   devices, algebra, the calculus\dots with a wealth of problems,
   recreations, constructions, applications explained and illustrated}},
   vol.~II, \textsl{Special topics of elemntary mathematics},
   Ginn \& Co., 1925 = Dover, 1958.

\bibitem{svz1} V. Spiridonov, L. Vinet and  A. Zhedanov,
 Periodic reductions of the factorization chain and the Hahn
polynomials, \textsl{J. Phys. A: Math. Gen.} {\bf 27} (1994), L669-L675.

\bibitem{sz} V. Spiridonov and A. Zhedanov, Discrete Darboux transformations,
the discrete-time Toda lattice and the Askey-Wilson polynomials,
\textsl{Methods Appl. Anal.} {\bf 2} (1995), 369-398.

\bibitem{SVZ} Vyacheslav Spiridonov,   Luc Vinet,  Alexei Zhedanov,
 Spectral Transformations, Self-Similar Reductions
and Orthogonal Polynomials, preprint CRM-2464, Montr\'eal, March 1997.

\bibitem{31} H.~Stahl, V.~Totik, {\sl General Orthogonal Polynomials\/},
         ({\sl Encyc.\ Math.\ Appl.\/} {\bf 43\/}), Cambridge U.P.,
         Cambridge, 1992.

\bibitem{zhed1} A.~Zhedanov, On some classes of the polynomials
  orthogonal on arcs of the unit circle connected with symmetric orthogonal
  polynomials on the interval, to appear in \textsl{J.\ Approx.\ Theory}.

\bibitem{Zhed} A.~Zhedanov, Rational spectral transformations and
  orthogonal polynomials, Preprint, 1997.


\end{thebibliography}



\section*{Appendix. About polygonal, or figurate, numbers.}
\renewcommand{\mylmark}{Appendix.}
\pagestyle{fancy}

I just encountered this enthusiastic account and can't resist
to include it:

{\Small

\noindent
sci.math \#225968 (24 + 4158 more)  \\
From: TechMaster Pinyan $<$jefpin@bergen.org$>$ \\
$[1] n-$gonal numbers    \\
Date: Sat Jul 19 09:07:43 EDT 1997

Ok... I am happy.   \\
I have made a great mathematical discovery.

Actually, I'm hyping this up a bit.

I was looking a "Five Hundred Mathematical Challenges" by Edward J.
Barbeau, Murray S. Klamkin, and William O. J. Moser.  In the back, they
have a "tool chest" of mathematical tidbits that can help you solve the
500 questions.  One of them had to do with figurate numbers:
        triangular, square, pentagonal, hexagonal...
And it got me thinking...

        "I know the formula for triangular numbers... $(x/2)(x+1)$...
         square
numbers are obviously $x^2$, and pentagonal numbers are $(x/2)(3x-1)$...
 there
has to be a pattern, so you can find the formula for the $x$th $n-$gonal
number."

There IS, but $x^2$ is a red herring.  Instead, look at $(x/2)(x+1)$ and
$(x/2)(3x-1)$:

       $$ 1x + 1$$
       $$ 3x - 1$$

Would it not follow that for a square number, you would use:

       $$ 2x + 0 $$

I was amazed.  I was thrilled.  I am pathetic :(.  But I did it!  I
figured out the formula.

For the $x$th $n-$gonal number, use this:

     $$   (x/2)[(n-2)x + (4-n)]   $$

Hooray...  this will make my calc teacher happy.  It may not be calculus,
but it's math...

For now, I'm gonna calculate the 453rd 32-gonal number...

----------------     \\
| "I don't contemplate it, I just sit and think about it."    \\
|       - Sonia Balsky             \\
----------------
Jeff "TechMaster" Pinyan | http://www.bergen.org/\~{}jefpin  \\
HTML/CGI Designer and Consultant and JavaScripter     \\
jefpin@bergen.org | TechMasterJeff@juno.com    \\
Got a JavaScript/CGI/Perl question or problem?  Let me know! \\
webXS - the new eZine for WebProgrammers! webXS@juno.com  \\
Visit us @ http://www.bergen.org/~jefpin/webXS

End of article 225968 (of 226703) -- what next? [npq]

\noindent
sci.math \#226021 (23 + 4158 more)   \\
From: ptwahl@aol.com (PTWahl)      \\
$[1] Re: n-$gonal numbers            \\
Date: Sat Jul 19 16:16:59 EDT 1997

Good work.  None of this is new, but there is value in finding things out
for oneself; new perspectives sometimes result in new shoots off an old
root.

For your amusement (I don't need to know the answer!)  I offer these
little tasks:

1.  Note that "$x$" in your formulas "for the $x$th $n-$gonal number" can be a
positive or negative integer.  Investigate what happens when you do that.

2.  We all know that $( x )^2 = ( -x )^2$, for the square numbers ($n=4$).
How about for the triangular numbers ($n=3$)?  Something similar (not quite
the same) is going on.  Discover what it is, and write a formula.

3.  List some hexagonal numbers ($n=6$), say from $x = -5$ to $x=+5$.  Compare
these to the triangular numbers. What is going on here?  Can that happen
for other sets of n-gonal numbers?

You write,  "It may not be calculus, but it's math ..."  But you're just
scratching the surface, so far.  There's calculus around the corner
(formal power series, modular forms, etc.) when you get farther along.
Figurate numbers provide a set of neat examples, all the way along.

There is a lot of fascination, even today, with the figurate numbers.  I
enjoy the topic partly because of the connection between "elementary"
techniques and much deeper results.  For example, there turns out to be a
connection between these two problems:

(A) The number of ways, counting sums in different order as different, to
express a positive integer $N$ as the sum of 24 triangular numbers.  This
turns out to be expressible
in terms of the sum of the 12th powers of the divisors of $N$ and $(N/2)$,
together with Ramanujan's intriguing "Tau" function (essentially the 24th
power of the "eta" function).

(B) The optimal packing of hyper-spheres in 24 dimensions, which among
other applications resulted in an efficient code for correcting errors in
data transmissions such as the Internet.

Good luck!

Patrick T. Wahl       \\
(no institutional affiliation)

End of article 226021 (of 226703) -- what next? [npq]

}

\newpage

Afternote: adjoint in Wolfgang Hahn's \textsl{Lineare Geometrische
Differenzengleichungen}, Bericht Nr.~169 der Mathematisch-Statistischen
Sektion im Forschungszentrum Graz (1981):

the general linear homogeneous system is
 \begin{equation} \label{fPx}
     \boldsymbol{f}(qx) = \boldsymbol{P}(x) \boldsymbol{f}(x),
 \end{equation}
where $\boldsymbol{f}$ is a column vector and $\boldsymbol{P}$ is
a square matrix. If one comes from a scalar equation of order $n$
$$
 p_0(x) f(x)+ p_1(x)f(qx) +\dotsb + p_n(x) f(q^n x) =0,
$$
i.e.,
 \begin{equation} \label{Lf}
   \mathcal{L}f := \sum_{k=0}^n \dfrac{p_k}{p_n} \mathsf{J}^k f =0,
 \end{equation}
(with $\mathsf{J}f(x):= f(qx)$), one then has
$$
 \boldsymbol{f}(x) = \begin{bmatrix} f(x) \\  \mathsf{J}f(x) \\  \vdots \\
                              \mathsf{J}^{n-1}f(x) \\
                              \mathsf{J}^n f(x) \end{bmatrix}
   \;, \qquad \boldsymbol{P} =
         \begin{bmatrix} 0  &  1  &   0  &  \cdots   &  0  \\
                         0  &  0  &   1  &  \cdots   &  0  \\
                            & \cdots&\cdots&\cdots   &    \\
                         0  &  0  &   0  &  \cdots   &  1  \\
   -\dfrac{p_0}{p_n} &
   -\dfrac{p_1}{p_n} &
   -\dfrac{p_2}{p_n} &
         \cdots &
   -\dfrac{p_{n-1}}{p_n}
  \end{bmatrix}
$$


\medskip

One defines the adjoint to \eqref{fPx} to be the equation
 \begin{equation} \label{gQx}
     \boldsymbol{g}(qx) = \boldsymbol{Q}(x) \boldsymbol{g}(x),
 \end{equation}
in such a way that the scalar product $\boldsymbol{g}^T \boldsymbol{f}$
is constant by multiplication by $q$, i.e., $q-$periodic.

As
$   \boldsymbol{g}^T(qx)   \boldsymbol{f}(qx)  =
   \boldsymbol{g}^T(x) \boldsymbol{Q}^T(x)
    \boldsymbol{P}(x) \boldsymbol{f}(x)$,  it follows
   $ \boldsymbol{Q}^T  \boldsymbol{P} = \boldsymbol{E}$,
and the adjoint to \eqref{fPx} is
 \begin{equation} \label{gQx2}
     \boldsymbol{g}(x) = \boldsymbol{P}^T(x) \boldsymbol{g}(qx), \qquad
 \text{ or } \qquad
     \boldsymbol{g}(q^{-1}x) = \boldsymbol{P}^T(q^{-1}x) \boldsymbol{g}(x).
 \end{equation}

\medskip

If the initial equation is scalar, as \eqref{Lf}, the adjoint
\eqref{gQx2} has not the same form, as $\boldsymbol{g}$ is not
compelled to have the form $[g(x), g(qx),\dotsc, g(q^{n-1}x)]^T$, but
we simply have to look for the recurrence satisfied by the
$n^{\text{th}}$ component of $\boldsymbol{g}$. We get this by
eliminating the other components: from \eqref{gQx2},
\begin{align}
g_{i+1}(q^{-1}x) &= g_i(x) -\dfrac{p_i(q^{-1}x)}{p_n(q^{-1}x)} g_n(x),\qquad
          i=0,1,\dotsc,n-1\qquad\qquad (g_0=0) \label{gi}\\
   &  \text{ or} \notag \\
g_i(x) &= g_{i+1}(q^{-1}x) +\dfrac{p_i(q^{-1}x)}{p_n(q^{-1}x)}
                                                      g_n(x), \notag \\
    &= g_{i+2}(q^{-2}x) +\dfrac{p_{i+1}(q^{-2}x)}{p_n(q^{-2}x)} g_n(q^{-1}x)
        +\dfrac{p_i(q^{-1}x)}{p_n(q^{-1}x)} g_n(x), \notag \\
   &  \dotsb \notag \\
  &= g_n(q^{i-n}x) +\dfrac{p_{n-1}(q^{i-n}x)}{p_n(q^{i-n}x)} g_n(q^{i-n+1}x)
      + \dotsb
        +\dfrac{p_i(q^{-1}x)}{p_n(q^{-1}x)} g_n(x)      \notag
\end{align}
and look at $i=0$:
 \begin{equation} \label{Mg}
   (\mathcal{M}g)(x) := g_0(x) =\sum_{k=0}^n
       \dfrac{p_k(q^{-k-1}x)}{p_n(q^{-k-1}x)} (\mathsf{J}^{-k} g)(x) =0,
 \end{equation}
(where $g(x)= g_n(x)$).

\bigskip

Finally, consider now, for any $f$ and $g$,
so that $\mathcal{L}f$ and $g_0=\mathcal{M}g$ need not to vanish,
$$ R(x) = \boldsymbol{g}^T\boldsymbol{f} =
  \sum_{i=0}^{n-1} g_{i+1}(x) f(q^i x)
 =\sum_{i=0}^{n-1} f(q^i x)\sum_{k=i+1}^n \dfrac{p_k(q^{i-k} x)}
                                                {p_n(q^{i-k} x)}
      g(q^{i-k+1} x),
$$
which should be constant ($q-$periodic) when
        $\mathcal{L}f =0$ and  $g_0=\mathcal{M}g=0$. Check:
\begin{equation*}\begin{split}
R(x)-R(qx) &=     \sum_{i=0}^{n-1} g_{i+1}(x) f(q^i x)
                - \sum_{j=0}^{n-1} g_{j+1}(qx) f(q^{j+1} x)       \\
           &=     \sum_{i=0}^{n-1} g_{i+1}(x) f(q^i x)
                - \sum_{i=1}^{n} g_{i}(qx) f(q^{i} x)       \\
           &= g_0(qx)f(x)+ \sum_{i=0}^{n-1}
            [g_{i+1}(x)-g_i(qx)] f(q^i x) \quad -g(qx) f(q^n x) \\
           &= f(x)(\mathcal{M}g)(qx) - \sum_{i=0}^{n-1}
           \dfrac{p_i(x)}{p_n(x)}g(qx) f(q^i x) \quad -g(qx) f(q^n x)
               \qquad \text{ (from }\eqref{gi}) \\
         &= f(x)(\mathcal{M}g)(qx) - g(qx) (\mathcal{L}f)(x).
\end{split}\end{equation*}
Sum this exact difference on a lattice $\{x_0, x_0q,\dotsc, x_0q^{N-1}\}$:
$$
R(x_0) -R(x_0 q^N) =  \sum_{r=0}^{N-1}
 f(x_r)(\mathcal{M}g)(qx_r) - g(qx_r) (\mathcal{L}f)(x_r).
$$
And the sensible duality form is
$$ \langle F,G\rangle = \sum_{r=0}^{N-1} F(x_r) G(q x_r) =
                        \sum_{r=0}^{N-1} F(x_r) G( x_{r+1}).
$$

\bigskip


cf. N.E.~N\"orlund, \textsl{Vorlesungen \"uber Differenzenrechnung},
 Berlin, 1924 = Chelsea, N.Y., 1954 (p.~289);
W.~Hahn, \"Uber die h\"oheren Heineschen Reihen und eine einheitliche
Theorie der sogenannten speziellen Funktionen, \textsl{Math. Nachr.}
\textbf{3} (1950) 257-294.

\smallskip

\noindent
(Of course, N\"orlund's book is about \emph{arithmetic} difference
equations, he finds $\displaystyle{\mu(x) + \sum_{i=1}^n q_{n-i}(x)
 \mu(x+i) =0}$ to be the adjoint of
                    $\displaystyle{u(x+n) + \sum_{i=0}^{n-1} q_i(x)
 u(x+i) =0}$, where $q_i(x) = p_i(x)/p_n(x)$).


\bigskip

\begin{center} \rule{10cm}{0.1mm} \end{center}

\vspace{66pt}

Hilbertian theory: according to Brown, Evans, and Ismail \cite{BI2},
the adjoint of $\mathcal{D}_q$ (see \eqref{AWD}, p.~\pageref{AWD})
with respect to $\langle F,G\rangle = \displaystyle{
  \int_0^\pi F(\cos\theta) \overline{G(\cos\theta)}\,d\theta
    = \int_{-1}^1 F(x)\overline{G(x)}\, (1-x^2)^{-1/2}\,dx  }$
   is \\
$\mathcal{D}^*_q$ defined by
$(\mathcal{D}^*_q)(f(x)) = -(1-x^2)^{1/2}\mathcal{D}_q[ (1-x^2)^{-1/2}
         f(x)]$.

\smallskip

The proof is rather hard.


\newpage

Other afternote: a (slightly) more complete version of the section
\\
\textbf{\ref{DLH2}.
Difference equation for Laguerre-Hahn orthogonal polynomials.}

We work again the general identities of \S~\ref{DLH1}, introducing
first other, hopefully more convenient, notations:

$$
\alpha_n(x) = -\dfrac{A_n(x)n}{\eta_2(x)- \eta_1(x)} +\dfrac{C_n(x)}{2},
\beta_n(x) = \dfrac{A_n(x)n}{\eta_2(x)- \eta_1(x)} +\dfrac{C_n(x)}{2},
\gamma_n(x) = \dfrac{B_n(x)}{\mu_0},
\delta_n(x) = \mu_0  D_n(x),
$$
so that  \eqref{ABCD} becomes

{\Small
\begin{equation*}  \hspace*{-1cm} \begin{split}
    -\beta_n(x)  &=
  \alpha_0
         a_n p_{n-1}(\eta_1) \tilde{p}_{n-1}(\eta_2)
 +\beta_0
         a_n \tilde{p}_{n-2}(\eta_1) p_n(\eta_2)
  +\gamma_0
         a_n \tilde{p}_{n-2}(\eta_1) \tilde{p}_{n-1}(\eta_2)
  +\delta_0
         a_n p_{n-1}(\eta_1) p_n(\eta_2)     \\
    -\alpha_n(x) &=
  \alpha_0
         a_n p_n(\eta_1) \tilde{p}_{n-2}(\eta_2)
 +\beta_0
         a_n \tilde{p}_{n-1}(\eta_1) p_{n-1}(\eta_2)
  +\gamma_0
         a_n \tilde{p}_{n-1}(\eta_1) \tilde{p}_{n-2}(\eta_2)
  +\delta_0
         a_n p_{n}(\eta_1) p_{n-1}(\eta_2)
  \\
        B_{n}(x)    &=
  \alpha_0
         a_n^2 p_{n-1}(\eta_1) \tilde{p}_{n-2}(\eta_2)
 +\beta_0
         a_n^2 \tilde{p}_{n-2}(\eta_1) p_{n-1}(\eta_2)
  +\gamma_0
         a_n^2 \tilde{p}_{n-2}(\eta_1) \tilde{p}_{n-2}(\eta_2)
  +\delta_0
         a_n^2 p_{n-1}(\eta_1) p_{n-1}(\eta_2)
  \\
          D_{n}(x) &=
  \alpha_0
          p_{n}(\eta_1) \tilde{p}_{n-1}(\eta_2)
 +\beta_0
          \tilde{p}_{n-1}(\eta_1) p_{n}(\eta_2)
  +\gamma_0
          \tilde{p}_{n-1}(\eta_1) \tilde{p}_{n-1}(\eta_2)
  +\delta_0
          p_{n}(\eta_1) p_{n}(\eta_2)
\end{split}\end{equation*}
}

We already got \eqref{pp1}-\eqref{pp2} through combinations of the
lines 1 and 4 above, then of lines 2 and 4. We go further, with  now
lines 1 and 3, then lines 2 and 3, always managing to use \eqref{casor}:

\begin{equation*} \begin{split}
 -\beta_n p_n(\eta_1) -a_n D_n p_{n-1}(\eta_1) &=
           -\beta_0 p_n(\eta_2) -\gamma_0  \tilde{p}_{n-1}(\eta_2), \\
 -\alpha_n p_n(\eta_2) -a_n D_n p_{n-1}(\eta_2) &=
           -\alpha_0 p_n(\eta_1) -\gamma_0  \tilde{p}_{n-1}(\eta_1), \\
 -a_n\beta_n \tilde{p}_{n-2}(\eta_2) -B_n \tilde{p}_{n-1}(\eta_2) &=
    -a_n \beta_0 \tilde{p}_{n-2}(\eta_1) -a_n\delta_0 p_{n-1}(\eta_1),\\
 -a_n\alpha_n \tilde{p}_{n-2}(\eta_1) -B_n \tilde{p}_{n-1}(\eta_1) &=
    -a_n \alpha_0 \tilde{p}_{n-2}(\eta_2) -a_n\delta_0 p_{n-1}(\eta_2),
\end{split}\end{equation*}
which shows how the $p$'s and the $\tilde{p}$'s at $\eta_1$ and $\eta_2$
are linked together, in other words: \\
\textbf{the general Laguerre-Hahn orthogonal polynomials enter a first
order difference system of 4 equations}
\begin{multline} \label{DLH4}
\boldsymbol{M}_n \begin{bmatrix} p_n(\eta_2) \\ p_{n-1}(\eta_2) \\
      \tilde{p}_{n-1}(\eta_2) \\  \tilde{p}_{n-2}(\eta_2) \end{bmatrix}
=\boldsymbol{N}_n \begin{bmatrix} p_n(\eta_1) \\ p_{n-1}(\eta_1) \\
      \tilde{p}_{n-1}(\eta_1) \\  \tilde{p}_{n-2}(\eta_1) \end{bmatrix},
\\
\boldsymbol{M}_n = \begin{bmatrix}
      -\beta_0  &           & -\gamma_0  &       \\
     -\alpha_n  & -a_n D_n  &            &       \\
                &           &  -B_n      & -a_n \beta_n
      \end{bmatrix},
\boldsymbol{N}_n = \begin{bmatrix}
   -\beta_n          &     -a_n D_n       &            &         \\
   -\alpha_0       &                    &  -\gamma_0   &         \\
                   &  -a_n \delta_0     &              & -a_n \beta_0 \\
               &                     &  -B_n       &  -a_n \alpha_n
      \end{bmatrix}.
\end{multline}

\bigskip

It should be possible to extract from \eqref{DLH4} a (probably ugly)
\emph{\textbf{scalar difference equation of order 4}} for $p_n$, but it
is much more interesting to investigate the structure of the set of
solutions of \eqref{DLH4}.

What about the matrix in $\boldsymbol{p}(\eta_2) = \boldsymbol{M}_n^{-1}
\boldsymbol{N}_n \boldsymbol{p}(\eta_1)$? Remark first that
$\boldsymbol{M_n}$ and $\boldsymbol{N}_n$ have the same
determinant $a_n^2[\alpha_0 \beta_0 B_n D_n -\gamma_0 \delta_0
\alpha_n \beta_n] =
             a_n^2[\alpha_0 \beta_0 B_n D_n - \alpha_n \beta_n B_0 D_0]$.
\\This asks for a bit more investigation of earlier identities,
the two first equations above \eqref{eqA} are
$$
  \beta_{m+1} = -\alpha_m -(\eta_2 -b_m)D_m, \qquad
  \alpha_{m+1} = -\beta_m -(\eta_1 -b_m)D_m,
$$
and \eqref{eqD} times $\displaystyle{
    \dfrac{B_{m+1}}{a_{m+1}^2} = D_m }$ is rearranged as
\begin{equation*} \begin{split} \hspace*{-50pt}
  D_{m+1} B_{m+1} &= \dfrac{(\eta_2-\eta_1)(\beta_m-\alpha_m)}{2} D_m
   +\dfrac{(\eta_1+\eta_2 -2b_m)(\alpha_m+\beta_m)}{2} D_m
   +B_m D_m +(\eta_1-b_m)(\eta_2-b_m) D_m^2 \\
 &= B_m D_m + [(\eta_1-b_m)\alpha_m +(\eta_2-b_m)\beta_m] D_m
            +(\eta_1-b_m)(\eta_2-b_m) D_m^2 \\
 &= B_m D_m + [(\eta_1-b_m)D_m +\beta_m][(\eta_2-b_m)D_m+\alpha_m]
     -\alpha_m \beta_m \\
  &= B_m D_m +\alpha_{m+1}\beta_{m+1} -\alpha_m\beta_m
\end{split}\end{equation*}
so, we recover  \eqref{Omega} as
$$
  B_m D_m -\alpha_m \beta_m = \dfrac{\Omega}{(\eta_2-\eta_1)^2}.
$$
The common determinant of $\boldsymbol{M}_n$ and $\boldsymbol{N}_n$
is therefore
$$ a_n^2\left[ \alpha_0 \beta_0 \left( \alpha_n \beta_n
           + \dfrac{\Omega}{(\eta_2-\eta_1)^2} \right)
              -\alpha_n \beta_n \left( \alpha_0 \beta_0
           + \dfrac{\Omega}{(\eta_2-\eta_1)^2} \right) \right]
  = a_n^2 (\alpha_0\beta_0 - \alpha_n\beta_n)
            \dfrac{\Omega}{(\eta_2-\eta_1)^2}.
$$
Now, we come to
 $ \boldsymbol{M}_n^{-1}\boldsymbol{N}_n$. I don't try to look for
clever shortcuts, just follow \eqref{DLH4}:
\begin{equation} \label{DLH41} \begin{split}  \hspace*{-44pt}
  \boldsymbol{M}_n^{-1}\boldsymbol{N}_n &=
 \dfrac{(\eta_2-\eta_1)^2}{a_n^2 (\alpha_0\beta_0 -\alpha_n\beta_n)
   \Omega}
 \begin{bmatrix}
  -a_n^2\alpha_0 B_n D_n  & a_n^2 \gamma_0\delta_0\beta_n &
           a_n^2 \alpha_0 \gamma_0 D_n  &  -a_n^2 \gamma_0 \beta_n D_n\\
 a_n \alpha_0 \alpha_n B_n &  -a_n\alpha_0\beta_0    B_n  &
      -a_n \alpha_0 \gamma_0 \alpha_n  & a_n \gamma_0 \alpha_n\beta_n \\
 a_n^2 \delta_0 \alpha_n\beta_n  &  -a_n^2 \beta_0 \delta_0 \beta_n &
   -a_n^2\alpha_0\beta_0 D_n  &  a_n^2 \beta_0 \beta_n D_n \\
 -a_n\delta_0 \alpha_n B_n &  a_n \beta_0\delta_0 B_n  &
   a_n \gamma_0\delta_0 \alpha_n & -a_n \beta_0 B_n D_n
\end{bmatrix}
   \\    &   \hspace*{221pt}
   \times \begin{bmatrix}
   -\beta_n          &     -a_n D_n       &            &         \\
   -\alpha_0       &                    &  -\gamma_0   &         \\
                   &  -a_n \delta_0     &              & -a_n \beta_0 \\
               &                     &  -B_n       &  -a_n \alpha_n
      \end{bmatrix}  \\
&=
 \dfrac{(\eta_2-\eta_1)^2  (\alpha_0\beta_0 -\alpha_n\beta_n) }
 {a_n^2 (\alpha_0\beta_0 -\alpha_n\beta_n)
   \Omega}
 \begin{bmatrix}
  -a_n^2 \alpha_0\beta_n  &  -a_n^3 \alpha_0 D_n &
          -a_n^2 \gamma_0 \beta_n  &  -a_n^3 \gamma_0 D_n \\
  a_n \alpha_0 B_n  &  a_n^2 \alpha_0 \alpha_n &
          a_n \gamma_0 B_n  &  a_n^2 \gamma_0 \alpha_n \\
  a_n^2 \delta_0\beta_n  &  a_n^3 \delta_0 D_n &
          a_n^2 \beta_0 \beta_n  &  a_n^3 \beta_0 D_n \\
  -a_n \delta_0 B_n  &  -a_n^2 \delta_0 \alpha_n &
          -a_n \beta_0 B_n  & -a_n^2 \beta_0 \alpha_n
      \end{bmatrix}.
\end{split}\end{equation}
Hmmm, an array mixing in a very remarkable way (Kronecker product)
things with index $n$ and things with index 0. Let us separate the two
parts of the mixture: consider the  system
\begin{equation} \label{DLH42}
  \begin{bmatrix} \varphi_n(\eta_2) \\  \psi_n(\eta_2) \end{bmatrix}
 = \begin{bmatrix} \beta_n   &   a_n D_n  \\
                    -a_n D_{n-1} & -\alpha_n  \end{bmatrix}
  \begin{bmatrix} \varphi_n(\eta_1) \\  \psi_n(\eta_1) \end{bmatrix}\;.
\end{equation}
Remark that he determinant of the matrix is $a_n^2 D_n D_{n-1}-\alpha_n\beta_n
= B_n D_n -\alpha_n\beta_n = \displaystyle{\dfrac{\Omega}{(\eta_2-\eta_1)^2}}$.

Consider also
\begin{equation} \label{DLH43}
  \begin{bmatrix} \phi_m(\eta_2) \\  \rho_m(\eta_2) \end{bmatrix}
  = \dfrac{(\eta_2-\eta_1)^2}{\Omega}
  \begin{bmatrix} -\alpha_m   &   -\gamma_m  \\
                    \delta_m  & \beta_m  \end{bmatrix}
  \begin{bmatrix} \phi_m(\eta_1) \\  \rho_m(\eta_1) \end{bmatrix}\;.
\end{equation}
Then, one finds (using $a_n D_{n-1} = B_n/a_n$),
$$
 \begin{bmatrix}
    \varphi_n(\eta_2) \phi_0(\eta_2) \\
    \psi_n(\eta_2) \phi_0(\eta_2) \\
    \varphi_n(\eta_2) \rho_0(\eta_2) \\
    \psi_n(\eta_2) \rho_0(\eta_2) \end{bmatrix}
 = \boldsymbol{M}_n^{-1}\boldsymbol{N}_n
 \begin{bmatrix}
    \varphi_n(\eta_1) \phi_0(\eta_1) \\
    \psi_n(\eta_1) \phi_0(\eta_1) \\
    \varphi_n(\eta_1) \rho_0(\eta_1) \\
    \psi_n(\eta_1) \rho_0(\eta_1) \end{bmatrix}
$$
So, if $\displaystyle{
         \begin{bmatrix} \varphi_n^{(1)} \\ \psi_n^{(1)} \end{bmatrix},
         \begin{bmatrix} \varphi_n^{(2)} \\ \psi_n^{(2)} \end{bmatrix};
  \qquad \begin{bmatrix} \phi_n^{(1)} \\ \rho_n^{(1)} \end{bmatrix},
         \begin{bmatrix} \phi_n^{(2)} \\ \rho_n^{(2)} \end{bmatrix}
}$
are fundamental systems of solutions of \eqref{DLH42}-\eqref{DLH43},
any solution of \eqref{DLH4}, and $\displaystyle{
                 \begin{bmatrix} p_n(x) \\ p_{n-1}(x) \\
      \tilde{p}_{n-1}(x) \\  \tilde{p}_{n-2}(x) \end{bmatrix}
   }$ in particular,
 is a combination
$$
\begin{bmatrix}
  \mathsf{A} \varphi_n^{(1)}(x) \phi_0^{(1)}(x)
 +\mathsf{B} \varphi_n^{(1)}(x) \phi_0^{(2)}(x)
 +\mathsf{C} \varphi_n^{(2)}(x) \phi_0^{(1)}(x)
 +\mathsf{D} \varphi_n^{(2)}(x) \phi_0^{(2)}(x)  \\
  \mathsf{A} \psi_n^{(1)}(x) \phi_0^{(1)}(x)
 +\mathsf{B} \psi_n^{(1)}(x) \phi_0^{(2)}(x)
 +\mathsf{C} \psi_n^{(2)}(x) \phi_0^{(1)}(x)
 +\mathsf{D} \psi_n^{(2)}(x) \phi_0^{(2)}(x)  \\
  \mathsf{A} \varphi_n^{(1)}(x) \rho_0^{(1)}(x)
 +\mathsf{B} \varphi_n^{(1)}(x) \rho_0^{(2)}(x)
 +\mathsf{C} \varphi_n^{(2)}(x) \rho_0^{(1)}(x)
 +\mathsf{D} \varphi_n^{(2)}(x) \rho_0^{(2)}(x)  \\
  \mathsf{A} \psi_n^{(1)}(x) \rho_0^{(1)}(x)
 +\mathsf{B} \psi_n^{(1)}(x) \rho_0^{(2)}(x)
 +\mathsf{C} \psi_n^{(2)}(x) \rho_0^{(1)}(x)
 +\mathsf{D} \psi_n^{(2)}(x) \rho_0^{(2)}(x)
\end{bmatrix}
$$

As we are interested in the special combinations leading to
$ \begin{bmatrix} p_n(x) & p_{n-1}(x) &
      \tilde{p}_{n-1}(x) &  \tilde{p}_{n-2}(x) \end{bmatrix}^T$,
let us rename the appropriated special combinations
$ \mathsf{A}  \phi_0^{(1)}
 +\mathsf{B}  \phi_0^{(2)}   $,
$ \mathsf{C}  \phi_0^{(1)}
 +\mathsf{D}  \phi_0^{(2)}  $,
$ \mathsf{A}  \rho_0^{(1)}
 +\mathsf{B}  \rho_0^{(2)}$,
$ \mathsf{C}  \rho_0^{(1)}
 +\mathsf{D}  \rho_0^{(2)}  $:
$   \phi_0^{(1)} , -\phi_0^{(2)}    ,
 \rho_0^{(1)}  , -\rho_0^{(2)}   $ respectively, so that
\begin{equation}  \label{pppp}  \begin{split}
  p_n(x) &= \varphi_n^{(1)}(x) \phi_0^{(1)}(x)
           -\varphi_n^{(2)}(x) \phi_0^{(2)}(x)  ,  \\
  p_{n-1}(x) &= \psi_n^{(1)}(x) \phi_0^{(1)}(x)
           -\psi_n^{(2)}(x) \phi_0^{(2)}(x)  ,  \\
  \tilde{p}_{n-1}(x)  &= \varphi_n^{(1)}(x) \rho_0^{(1)}(x)
           -\varphi_n^{(2)}(x) \rho_0^{(2)}(x)  ,  \\
  \tilde{p}_{n-2}(x)  &= \psi_n^{(1)}(x) \rho_0^{(1)}(x)
           -\psi_n^{(2)}(x) \rho_0^{(2)}(x)  .
\end{split}\end{equation}
The minus sign will make things more comfortable at $n=0$ (and we also
slowly approach important formulas by Hahn, such as
\cite[formula (16), p.4]{ha1}), where the left-hand side is made of
constants. That combinations of products of solutions of \eqref{DLH42}
and \eqref{DLH43} may turn to constants requires a closer look.
Indeed, inverting the matrix of \eqref{DLH43} yields
$$
  \begin{bmatrix} \phi_m(\eta_1) \\  \rho_m(\eta_1) \end{bmatrix}
  =
  \begin{bmatrix} \beta_m   &   \gamma_m  \\
                   -\delta_m  & \alpha_m  \end{bmatrix}
  \begin{bmatrix} \phi_m(\eta_2) \\  \rho_m(\eta_2) \end{bmatrix}\;,
 \text{ \  or\  }
  \begin{bmatrix} \phi_m(\eta_1) \\
   -\dfrac{a_m \rho_m(\eta_1)}{\mu_0} \end{bmatrix}
  =
  \begin{bmatrix} \beta_m   &   -B_m/a_m  \\
                   a_m D_m    & \alpha_m  \end{bmatrix}
  \begin{bmatrix} \phi_m(\eta_2) \\
   -\dfrac{a_m \rho_m(\eta_2)}{\mu_0} \end{bmatrix}    \;,
$$
the adjoint of \eqref{DLH42} (when $m=n$). So, for any solutions of
\eqref{DLH42} and \eqref{DLH43},
$$ \varphi_m \phi_m - a_m \psi_m \rho_m/\mu_0 $$
is constant (on any lattice), or:
\begin{equation} \label{phirho}
\begin{bmatrix}
  \phi_m^{(1)}     &  -a_m \dfrac{\rho_m^{(1)}}{\mu_0} \\
  \phi_m^{(2)}     &  -a_m \dfrac{\rho_m^{(2)}}{\mu_0}
\end{bmatrix} =
 \text{ a constant matrix times }
\begin{bmatrix}
  \psi_m^{(2)} &   -\varphi_m^{(2)} \\
 -\psi_m^{(1)} &    \varphi_m^{(1)}
\end{bmatrix}
 \dfrac{1}{\varphi_m^{(1)} \psi_m^{(2)} - \psi_m^{(1)} \varphi_m^{(2)} }
\end{equation}

\medskip

Let us look more closely at the solutions of \eqref{DLH42}. The scalar
second order equations for $\varphi_n$ and $\psi_n$ are obtained by
performing various eliminations in
\begin{equation*} \begin{split}
\varphi_n(y_{k+1}) &= \beta_n(x_k) \varphi_n(y_k)
                      +a_n D_n(x_k) \psi_n(y_k), \\
\varphi_n(y_k) &= \beta_n(x_{k-1}) \varphi_n(y_{k-1})
                      +a_n D_n(x_{k-1}) \psi_n(y_{k-1}), \\
\psi_n(y_{k+1}) &= -a_n D_{n-1}(x_k) \varphi_n(y_k)
                      -\alpha_n(x_k) \psi_n(y_k), \\
\psi_n(y_k) &= -a_n D_{n-1}(x_{k-1}) \varphi_n(y_{k-1})
                      -\alpha_n(x_{k-1}) \psi_n(y_{k-1}),
\end{split}\end{equation*}
and are
\begin{equation} \label{eqphin2}    \hspace*{-30pt}
D_n(x_{k-1}) \varphi_n(y_{k+1}) +[\alpha_n(x_{k-1}) D_n(x_k)
   -\beta_n(x_k) D_n(x_{k-1}] \varphi_n(y_k)
  +D_n(x_k) \dfrac{\Omega(x_{k-1})}{(y_k-y_{k-1})^2}\varphi_n(y_{k-1}) =0
\end{equation}
and

$$
\hspace*{-55pt}
D_{n-1}(x_{k-1}) \psi_n(y_{k+1}) +[\alpha_{n-1}(x_{k-1}) D_{n-1}(x_k)
   -\beta_{n-1}(x_k) D_{n-1}(x_{k-1}] \psi_n(y_k)
  +D_{n-1}(x_k) \dfrac{\Omega(x_{k-1})}{(y_k-y_{k-1})^2}
  \psi_n(y_{k-1}) =0
$$
showing clearly\footnote{After performing $\alpha_n(x_k) D_{n-1}(x_{k-1})
-\beta_n(x_{k-1})D_{n-1}(x_k) = [\alpha_{n-1}(x_{k-1}) +(y_k-b_{n-1})
 D_{n-1}(x_{k-1})] D_{n-1}(x_k)
                               -[\beta_{n-1}(x_k) +(y_k-b_{n-1})
  D_{n-1}(x_k)]D_{n-1}(x_{k-1}) =
             \alpha_{n-1}(x_{k-1})  D_{n-1}(x_k)
         -\beta_{n-1}(x_k) D_{n-1}(x_{k-1})$.}that $\psi_n$ is a solution
of the equation for $\varphi_{n-1}$.

Shall we conclude $\psi_n^{(1)} = \varphi_{n-1}^{(1)}$ and
 $\psi_n^{(2)} = \varphi_{n-1}^{(2)}$? We can do that for a single value
of $n$, say $n=0$. Then, from the various identities for $\alpha_{n+1}$,
$\beta_{n+1}$, and $D_{n+1}$ above, we see that
\begin{equation*} \begin{split} \hspace*{-22pt}
(\eta_2-b_n)\varphi_n(\eta_2) -a_n\psi_n(\eta_2) &=
  (\eta_2-b_n)\beta_n \varphi_n(\eta_1) +(\eta_2-b_n)a_n D_n \psi_n(\eta_1)
    +a_n^2 D_{n-1}\varphi_n(\eta_1) +a_n\alpha_n \psi_n(\eta_1) \\
           =\beta_{n+1} [(\eta_1-bn)\varphi_n(\eta_1) -&a_n\psi_n(\eta_1)]
  +[(\eta_2-b_n)\beta_n +(\eta_1-b_n)\alpha_n +(\eta_1-b_n)(\eta_2-b_n)
      D_n +a_n^2 D_{n-1}]\varphi_n(\eta_1) \\
 &= \beta_{n+1} [(\eta_1-bn)\varphi_n(\eta_1) -a_n\psi_n(\eta_1)]
       +a_{n+1}^2 D_{n+1}\varphi_n(\eta_1)
\end{split}\end{equation*}
and
\begin{equation*} \begin{split}
 \varphi_n(\eta_2) &= \beta_n \varphi_n(\eta_1) +a_n D_n \psi_n(\eta_1) \\
     &= -D_n [(\eta_1-b_n)\varphi_n(\eta_1) -a_n\psi_n(\eta_1)]
          +[(\eta_1-b_n)D_n +\beta_n] \varphi_n(\eta_1) \\
     &= -D_n [(\eta_1-b_n)\varphi_n(\eta_1) -a_n\psi_n(\eta_1)]
          -\alpha_{n+1}  \varphi_n(\eta_1)
\end{split}\end{equation*}
so, we recover the recurrence
\begin{equation} \label{recphi}
    a_{n+1} \varphi_{n+1}(x) = (x-b_n) \varphi_n(x) -a_n \varphi_{n-1}(x)
\end{equation}
and summarize the whole thing in the \\
\textbf{Theorem.} \textsl{There exist sequences of functions
 $\{\varphi_n^{(1)}\}$ and
 $\{\varphi_n^{(2)}\}$, $n\in\mathbb{Z}$, satisfying simultaneously the difference
equation  \eqref{eqphin2} and the recurrence relation \eqref{recphi}
such that}
\begin{equation*} \begin{split}
 p_n(x) &= \dfrac{1}{\sqrt{\mu_0}}
    \dfrac{\varphi_n^{(1)}(x) \varphi_{-1}^{(2)}(x)
          -\varphi_n^{(2)}(x) \varphi_{-1}^{(1)}(x)}
          {\varphi_0^{(1)}(x) \varphi_{-1}^{(2)}(x)
          -\varphi_0^{(2)}(x) \varphi_{-1}^{(1)}(x)} , \\
 \tilde{p}_{n-1}(x) &= \dfrac{\sqrt{\mu_0}}{a_0}
    \dfrac{\varphi_n^{(1)}(x) \varphi_{-1}^{(2)}(x)
          -\varphi_n^{(2)}(x) \varphi_{-1}^{(1)}(x)}
          {\varphi_0^{(1)}(x) \varphi_{-1}^{(2)}(x)
          -\varphi_0^{(2)}(x) \varphi_{-1}^{(1)}(x)} , \\
      &    \\
    & \hspace*{-33pt}\text{associated polynomials:}  \\
 p_{n-r}^{(r)}(x) &= K_r
    \dfrac{\varphi_n^{(1)}(x) \varphi_{r-1}^{(2)}(x)
          -\varphi_n^{(2)}(x) \varphi_{r-1}^{(1)}(x)}
          {\varphi_0^{(1)}(x) \varphi_{-1}^{(2)}(x)
          -\varphi_0^{(2)}(x) \varphi_{-1}^{(1)}(x)}.
\end{split}\end{equation*}
The factors $\phi$ and $\rho$ in \eqref{pppp} have been replaced by
their expressions in \eqref{phirho}, in such a way that
 $p_{-1} = \tilde{p}_{-1} = p_{-1}^{(r)} =0$.

\bigskip

\textbf{Hahn's theory.}

We just had a little taste of a formidable field of investigation
worked by W.~Hahn:

   Hahn, Wolfgang:

   Ueber lineare Differentialgleichungen, deren Loesungen einer
   Rekursionsformel genuegen.
   \textsl{Math. Nachr.} \textbf{4}, 1-11 (1951).

   Ueber lineare Differentialgleichungen, deren Loesungen einer
   Rekursionsformel genuegen. II.
   \textsl{Math. Nachr.} \textbf{7}, 85-104 (1952).

summarized (?) in

   Hahn, Wolfgang

   Kennzeichnung der linearen Differentialgleichungen, die durch
   Orthogonalpolynome befriedigt werden.
    Bericht No.54. Graz: Mathematisch-statistische Sektion im
   Forschungszentrum Graz. 16 S. (1976).

where we find:

  \emph{\textbf{Theorem.}} \textsl{Let the sequence $\{y_n(x)\}$,
  $n\in\mathbb{Z}$, of functions of the complex variable $x$
 satisfy simultaneously
   \begin{enumerate}
  \item Each $y_n$ is solution of a linear homogeneous differential
   equation of order $k$
     \begin{equation} \label{hahnL}
     L_n y := \sum_{j=0}^k c_{n,j}(x) \dfrac{d^{k-j} y}{dx^{k-j}}(x) =0
    \end{equation}
  in $x$, where the coefficients $c_{n,j}$ are polynomials with
 degrees bounded with respect to $n$; $y_n$ can \emph{not} be written as
  solution of such an equation of order smaller than $k$.
   \item The recurrence relation
    \begin{equation} \label{hahnR}
         y_n(x) = A_n(x) y_{n-1}(x) +B_n(x) y_{n-2}(x)
    \end{equation}
  holds, with coefficients which are polynomials with uniformly
    bounded degrees.
  \item Neither the differential equation \eqref{hahnL} nor any of the two linear
    differential equations whose $A_n y_{n-1}$ and $B_n y_{n-2}$ are
   solutions  have any common solution.
  \end{enumerate}
  Then,
   \begin{enumerate}
    \item $k$ must be even: $k=2r$,
   \item The differential equation \eqref{hahnL} has a fundamental
    system of solutions of the form
     $$
         \phi^{(i)} \varphi^{(1)}_n,  \phi^{(i)} \varphi_n^{(2)}\;,
             \qquad i=1,\dotsc, r,
    $$
      where
              $\varphi^{(1)}_n$ and $ \varphi_n^{(2)}$ are
    solutions of a differential equations of order~2 with polynomial
     coefficients of bounded degree, and where $\{\phi^{(1)},\dotsc,
    \phi^{(r)}\}$ is a fundamental system of solutions of a linear
    differential equation of order~$r$.
   \item If we want $y_n$ to be a polynomial of degree $n$, we can
    only have $k=2$ or $k=4$, and the degrees of $A_n$ and $B_n$ be
    1 and 0. \\
       When $k=4$, one has $y_n= \varphi_n^{(1)} \varphi_{-1}^{(2)}
                               -\varphi_n^{(2)} \varphi_{-1}^{(1)}  $.
   \end{enumerate}
  }

Challenge: extend this to the $\mathscr{D}$ operator!!


\bigskip

\textbf{Stieltjes function:}
 if we succeed to have $\varphi_n^{(2)}(x)/\varphi_n^{(1)}(x) \to 0$
when $n\to\infty$,
 $$ \mathcal{S}(x) = \lim_{n\to\infty} \dfrac{\tilde{p}_{n-1}(x)}{p_n(x)} =
     \dfrac{\varphi_0^{(2)}(x)}{\varphi_{-1}^{(2)}(x)},
 $$
 i.e., a solution of a Riccati equation as ratio of solutions of linear second order
 difference equations, as it should.

\bigskip

See also
   Gruenbaum, F.Alberto; Haine, Luc:
   Some functions that generalize the Askey-Wilson polynomials.
   \textsl{
    Commun. Math. Phys.} \textbf{184}, No.1, 173-202 (1997).

   The author's abstract: We determine all biinfinite tridiagonal
   matrices for which some family of eigenfunctions are also
   eigenfunctions of a second order $q$-difference operator. The solution
   is described in terms of an arbitrary solution of a $q$-analogue of
   Gauss hypergeometric equation depending on five free parameters and
   extends the four dimensional family of solutions given by the
   Askey-Wilson polynomials. There is some evidence that this bispectral
   problem, for an arbitrary order $q$-difference operator, is intimately
   related with some $q$-deformation of the Toda lattice hierarchy and
   its Virasoro symmetries. When tridiagonal matrices are replaced by the
   Schroedinger operator, and $q=1$, this statement holds with Toda
   replaced by KdV. In this context, this paper determines the analogs of
   the Bessel and Airy potentials.

   Gruenbaum, F.Alberto; Haine, Luc:
   On a $q$-analogue of Gauss equation and some $q$-Riccati equations.
    Ismail, Mourad E. H. (ed.) et al., \textsl{Special functions, $q$-series
   and related topics}. Providence, RI: American Mathematical Society,
   \textsl{Fields Inst. Commun}. \textbf{14}, 77-81 (1997).

   Gruenbaum, F.Alberto; Haine, Luc
   A theorem of Bochner, revisited.
    Fokas, A. S. (ed.) et al., \textsl{Algebraic aspects of integrable
   systems: in memory of Irene Dorfman}. Boston, MA: Birkhaeuser, \textsl{Prog.
   Nonlinear Differ. Equ. Appl}. \textbf{26}, 143-172 (1997).

   We pose a bispectral problem associated with the Toda lattice
   hierarchy and describe explicitly all its ``basic solutions'' from
   which all other solutions are expected to be obtained by repeated
   application of the Darboux process. We identify these basic solutions
   as fixed points of (in general) some (sub)algebras of the Virasoro
   algebra formed by the Toda master symmetries. Finally, we illustrate
   by four simple examples our contention that repeated applications of
   the Darboux transformation from some of these fixed points will
   provide integral manifolds of the Toda hierarchy or its master
   symmetries.

\end{document}

