
(* Asymptotic solutions of P-finite recurrences *) 
(* by M. Kauers, started 08/09/25 *)

(* ------------------------------------------------------------------------- *)

(* LICENCE:                                                                  *)
(*                                                                           *)
(*  This software package is free; you can redistribute it and/or modify it  *)
(*  under the terms of the GNU General Public License as published by the    *)
(*  Free Software Foundation; either version 2 of the License, or (at your   *)
(*  option) any later version. The package is distributed in the hope that   *)
(*  it will be useful, but without any warranty; without even the implied    *)
(*  warranty of merchantability or fitness for a particular purpose. See the *)
(*  GNU General Public License for more details.                             *)

(* ------------------------------------------------------------------------- *)

(* bugs:

   Hypergeometric2F1[1/2,-1/2,1,x]^2 : higher order terms lead to exception

   problem with higher order terms 
4*n^4*(2 + n)^4*f[n] - 2*(1 + 2*n)*(3 + 2*n)*(81 + 423*n + 888*n^2 + 930*n^3 + 515*n^4 + 144*n^5 + 16*n^6)*
  f[1 + n] + (1 + n)*(3 + n)*(1 + 2*n)*(3 + 2*n)*(5 + 4*n)*(7 + 4*n)*(9 + 4*n)*(11 + 4*n)*f[2 + n]

 *)

Asymptotics`Private`Version = "0.11 (2012-07-19)";

BeginPackage["RISCComb`"]

{`Infolevel}

EndPackage[]

BeginPackage["Asymptotics`", {"RISCComb`"}]

`Asymptotics::usage = "Asymptotics[rec, f[n]] computes the dominant terms of the series solutions " <>
  "of the recurrence rec.  The recurrence rec must be an expression of the form poly*f[n] + poly*f[n+1] " <>
  "+ poly*f[n+2] + ... + poly*f[n+r] where r is an explicit positive integer and 'poly' stands for some " <>
  "expressions which depend polynomially on n. (In other words, rec must be a P-finite recurrence.)\n\n" <>
  "Options:\n" <>
  "Order -> [Integer] ... computes also higher order terms up to the specified order\n" <>
  "Return -> [String] ... determines the output format: \"internal\" returns a machine friendly output " <>
  "and any other value returns a human friendly output\n\n" <>
  "Example:\n" <> 
  "Asymptotics[(n+3)f[n] - (3n+1)f[n+1]+(2n+2)f[n+2], f[n], Order -> 1]";

`FurtherTerms::usage = "FurtherTerms[rec, f[n], init, ord] computes higher order terms of a partially " <>
  "given series solution of a recurrence equation.\n\n" <>
  "Arguments:\n" <>
  "rec... homogeneous linear recurrence equation in f[n] with polynomial coefficients\n" <> 
  "f[n]... function symbol in terms of which the recurrence is specified\n" <>
  "init... partial series solution of 'rec' in internal formart as produced by the Asymptotics command\n" <>
  "ord... total number of terms desired in the output (i.e. ord is the sum of the number of terms already " <>
  "available in 'init' and the number of terms to be computed\n\n" <>
  "Example:\n" <>
  "FurtherTerms[(n+3)f[n] - (3n+1)f[n+1]+(2n+2)f[n+2], f[n], {0,1,{},-4,{1,38,1045,26110}}, 10]";

`GuessAsymptotics::usage = "GuessAsymptotics[init, rec, f[n]] determines empirically the dominant term for " <>
  "the asymptotic behaviour of the sequence determined by the given recurrence and initial values. This is " <>
  "often faster than using the Asymptotics command, but it fails if the specified sequence does not have " <> 
  "a \"pure\" dominant term, e.g., (-2)^n + 2^n is not good.";

Begin["`Private`"]

Off[Solve::svars]

HomogEquation2CoefficientList[eq_, u_[x_]] :=
  Module[ {r, i, j, rMin, rMax, myEq},
    r = Cases[{eq}, u[x+i_.] -> i, Infinity];
    rMin = Min[r]; rMax = Max[r];
    myEq = eq /. x -> x - rMin /. i_ == j_ -> i - j // Together // Numerator;
    If[ (myEq /. u[_] -> 0) =!= 0, Throw["Inhomogeneous equation encountered."] ];
    Return[Table[Coefficient[myEq, u[x+i]], {i, 0, rMax - rMin}]]
  ];

EExpansion[k_Integer, j_] := 
  Which[
    k < 0, 0,
    k == 0, 1,
    k <= Length[ECoefficients], ECoefficients[[k]] /. X -> j,
    True, 
      Module[ {n, i},
        ECoefficients = 
          Rest[CoefficientList[Sum[Series[(1/n Log[1+n X]-X)^i, {n,0,k}]/i!, {i,0,k}], n]];
        Last[ECoefficients] /. X -> j
      ]
   ];
ECoefficients = {-X^2/2};

(* a number d such that expr = O(n^d), works for expressions involving 
   fractional exponents; may overshoot *)
DegreeBound[expr_, n_] := 
  Which[
    FreeQ[expr, n], 0,
    PolynomialQ[expr, n], Exponent[expr, n],
    MatchQ[expr, _^_?NumberQ], DegreeBound[expr[[1]], n]*expr[[2]],
    Head[expr] === Plus, Max[DegreeBound[#, n]& /@ List@@expr],
    Head[expr] === Times, Plus@@(DegreeBound[#, n]& /@ List@@expr),
    True, Throw["illegal expression encountered"]
  ];

(* a leading coefficient extraction which (hopefully) also works for 
   expressions involving radicals *)
LeadingTerm[expr_, n_] :=
  Module[ {s, d},
    If[ PuiseuxRationalFunctionQ[expr, n], Return[MySeries[expr, n, 1]] ];
    (* old code ahead *)
    d = DegreeBound[Together[expr, Extension -> Automatic], n];
    If[ d === -Infinity, Return[ Together[expr, Extension -> Automatic]] ];
    s = Series[expr, {n, Infinity, 1 - Min[0, Ceiling[d]]}]; (* trustworthy? *)
    Which[
      FreeQ[s, n], s,
      Head[s] === SeriesData && Length[s[[3]]] >= 1, s[[3, 1]]n^(-s[[4]]/s[[6]]),
      True, Throw["extraction of leading coefficient failed. "]
    ]
  ];
LeadingCoefficient[expr_, n_] := (LeadingTerm[expr, n] /. n -> 1);

(* checks whether expr is a quotient of two puiseux polynomials *)
PuiseuxRationalFunctionQ[expr_, n_] :=
  Switch[ Head[expr],
    Times|Plus, And@@(PuiseuxRationalFunctionQ[#, n]& /@ List@@expr),
    Power, IntegerQ[Last[expr]] || First[expr]==n,
    _, FreeQ[expr, n]
  ];

(* expansion at infinity. assumption: rat is a quotient of two puiseux polynomials.
   accurracy gives the number of terms to be computed. 
   The result will be a puiseux polynomial. *)
MySeries[rat_, n_, accurracy_Integer] :=
  Module[{num, den, dnum, dden, term},
    If[accurracy === 0, Return[0]];
    {num, den} = {Numerator[#], Denominator[#]}&[Together[rat, Extension -> Automatic]];
    If[num == 0, Return[0]];
    {dnum, dden} = Exponent[#, n]& /@ {num, den};
    term = Coefficient[Expand[num], n, dnum]/Coefficient[Expand[den], n, dden] n^(dnum-dden);
    Return[term + MySeries[rat - term, n, accurracy - 1]];
  ];

(* given a recurrence for f[n] as output by HomogEquation2CoefficientList, 
   this returns a recurrence for f[n] / (n/e)^(gamma n) rho^n Exp[mu.n^(Range[Length[mu]]/(Length[mu] + 1))] n^alpha, for the given values of gamma, rho, mu, and alpha *)
Substitute[coeffs_, n_, asympt_, accurracy_Integer] := 
  Module[ {k, i, j, r, newcoeffs, q, gamma, rho, mu, alpha},
    r = Length[coeffs] - 1; 
    {gamma, rho, mu, alpha} = asympt; q = Length[mu];
    newcoeffs = coeffs*Table[
      Sum[EExpansion[k, j]/n^k, {k, 0, Denominator[gamma]accurracy}]^gamma *
      (n^(j gamma) + Sum[Binomial[j gamma, k]j^k n^(j gamma - k), {k, 1, Denominator[gamma]accurracy}]) *
      (1 + Sum[(mu.(Sum[Binomial[#, i] j^i n^(# - i), {i, 1, (q+1)accurracy - k}]& 
                  /@ (Range[q]/(q+1))))^k/k!, {k, 1, (q+1)accurracy}]) *
      rho^j *
      (1 + Sum[Pochhammer[alpha - k + 1, k]/k!(j/n)^k, {k, 1, accurracy}]) 
    , {j, 0, r}];
    ((*Normal[Series[#, {n, Infinity, accurracy}]]&*)
     MySeries[#, n, accurracy]& /@ Together[newcoeffs/LeadingTerm[First[coeffs], n]]) //. (n^a_)^b_ -> n^(a b)
  ];

(* asymptotic solutions of the form (n/E)^(gamma n) Exp[mu[1] n^(1/q) + mu[2] n^(2/q) + mu[q-1] n^(1-1/q)] rho^n n^alpha 
   returned in the form {gamma, {mu[1], mu[2], ..., mu[q-1]}, rho, alpha} *)
Options[Asymptotics] = {Order -> 0, Return -> "expression", Infolevel -> 0, Rational -> False, 
  ToNumberField -> False, Skip -> {}};
Asymptotics[rec_, f_[n_], opts___] :=
  Module[ {sols, sol, coeffs, mycoeffs, degs, i, j, k, kk, l, r, c, d, max, L, info, rat, alg, skip,
           gamma, rho, mu, alpha, beta, lambda, poly, exp, head, terms, ret},

    {terms, ret, info, rat, alg, skip} = {Order, Return, Infolevel, Rational, ToNumberField, Skip}
                                        /. {opts} /. Options[Asymptotics];
    If[ Head[skip] == String, skip = {skip} ];
    coeffs = HomogEquation2CoefficientList[rec, f[n]];
    r = Length[coeffs] - 1;
    sols = {};

    (* 1. Determine all possible gammas *)
    If[ MemberQ[skip, "gamma"], 
      If[info > 0, Print["Taking default gammas: 0"] ];
      sols = {0}
    , If[info > 0, Print["Determine all possible gammas"] ];
      degs = Select[Exponent[coeffs, n] + gamma Range[0, Length[coeffs] - 1], FreeQ[#, -Infinity]&];
      sols = {gamma, degs} 
          /. Flatten[Table[Solve[degs[[i]] == degs[[j]], gamma], 
                           {i, 1, Length[degs]}, {j, i + 1, Length[degs]}], 2];
      sols = First /@ Select[sols, Length[Position[#[[2]], Max[#[[2]]]]]>1&];
      sols = Union[sols];
    ];

    If[rat, sols = Select[sols, Element[#, Rationals]&]];
    
    If[info > 1, Print["  ", Length[sols], " gammas isolated."]];

    (* 2. For each gamma, determine all possible rhos *)
    If[ MemberQ[skip, "rho"],
      If[info > 0, Print["Taking default rho: 1"] ];
      sols = Append[#, 1]& /@ sols
    , If[info > 0, Print["For each gamma, determine all possible rhos"] ];
      Do[
        gamma = sols[[l]];
        If[info > 2, Print["    considering gamma=", InputForm[gamma]]];
        poly = LeadingCoefficient[Plus@@Substitute[coeffs, n, {gamma, rho, {}, 0}, 1], n];
        If[info > 3, Print["      lc=", InputForm[Numerator[Factor[poly]] /. rho -> Symbol["rho"]]]];
        poly = Table[i[j[#], k], {k, 1, Exponent[#, rho]}]& /@ First /@ FactorList[poly];
        poly = Flatten[poly] /. {i -> Root, j -> Function, rho -> #};
        If[rat, poly = Select[poly, Element[#, Rationals]&]];
        sols[[l]] = {gamma, #}& /@ Union[Flatten[DeleteCases[poly, 0]]];
      , {l, 1, Length[sols]}];
      sols = Union[Flatten[sols, 1]];
    ];

    If[info > 1, Print["  ", Length[sols], " pairs (gamma,rho) isolated."]];

    If[ MemberQ[skip, "subdominant"], 
       sols = Extract[sols, Position[#, Max[#], 1]&[Abs[Last /@ sols]]];
       If[info > 1, Print["keeping ", Length[sols]] ];
    ];

    If[alg, sols = ToNumberField /@ sols]; (* assuming no algebraic numbers are in the input *)

    (* 3. Determine all possible exponential parts *)
    If[ MemberQ[skip, "mu"],
      If[info > 0, Print["Taking default mu: {}"] ];
      sols = Append[#, {}]& /@ sols
    , If[info > 0, Print["For each (gamma,rho), determine all possible exponential parts"]];
      Do[
        {gamma, rho} = sols[[l]]; mu = {}; k = r; max = 1;
        If[info > 2, Print["    considering {gamma,rho}=", InputForm[{gamma,rho}]]];
        mycoeffs = Substitute[coeffs, n, {gamma, rho, mu, alpha}, 2k];
        While[ FreeQ[head = LeadingTerm[Plus@@mycoeffs, n], alpha],
          If[info > 3, Print["      compensating exponential term introduced"]];
          While[ head === 0, 
            head = LeadingTerm[Plus@@(mycoeffs = Substitute[coeffs, n, {gamma, rho, mu, alpha}, ++k]), n]];
          If[ !FreeQ[ head, alpha ], Break[] ];
          While[ True, 
            poly = Expand[Plus@@(mycoeffs.Table[1 + Sum[(lambda beta j n^(beta-1))^i/i!, {i, 1, 2k}], {j, 0, r}])];
            exp = Select[Union[Cases[poly, n^u_ -> u, Infinity]], !FreeQ[#, beta]&];
            sol = First[Solve[# == Exponent[head, n]]]& /@ exp;
            sol = Extract[sol, Position[Max /@ (exp /. sol), Exponent[head, n]]];
            If[ Length[sol] > 0, Break[] ];
            mycoeffs = Substitute[coeffs, n, {gamma, rho, mu, alpha}, ++k];
          ];
          sol = Min[beta /. sol];
          poly = lambda /. Solve[LeadingCoefficient[poly /. beta -> sol, n]==0, lambda];
          kk = LCM[Denominator[sol], Length[mu]+1]/(Length[mu]+1) - 1;
          mu = Flatten[{{Table[0, {kk}], #}& /@ mu, Table[0, {kk}]}];
          mu[[sol*(Length[mu]+1)]] = First[poly];
          mycoeffs = Substitute[coeffs, n, {gamma, rho, mu, alpha}, ++k];
        ];
        Clear[k];
        sols[[l]] = Table[{gamma, rho, FullSimplify[Table[Exp[2*Pi*I*i*k/(Length[mu]+1)], {i, 1, Length[mu]}]]mu}, {k, 0, Length[mu]}];
      , {l, 1, Length[sols]}];
      sols = Union[Flatten[sols, 1]];
    ];

    If[info > 1, Print["  ", Length[sols], " partial solutions (gamma,rho,mu) isolated."]];

    (* 4. For each gamma and rho, determine all possible alphas *)    
    If[ MemberQ[skip, "alpha"], 
      If[info > 0, Print["Taking default alpha: {}"] ];
      sols = Append[#, 0]& /@ sols
    , If[info > 0, Print["For each (gamma,rho,mu), determine all possible alphas"]];
      Do[
        {gamma, rho, mu} = sols[[l]]; poly = 0; k = r;
        If[info > 2, Print["    Considering (gamma,rho,mu)=", InputForm[sols[[l]]]]];
        While[ poly === 0,
          poly = LeadingCoefficient[Plus@@Substitute[coeffs, n, {gamma, rho, mu, alpha}, ++k], n]; 
        ];
        If[ FreeQ[poly, alpha], sols[[l]] = {}; Continue[] ];
        poly = First /@ FactorList[poly];
        poly = Table[i[j[#], k], {k, 1, Exponent[#, alpha]}]& /@ poly;
        poly = Flatten[poly] /. {i -> Root, j -> Function, alpha -> #};
        If[ Length[poly] <= 2, poly = RootReduce[poly] ];
        sols[[l]] = {gamma, rho, mu, #}& /@ poly;
      , {l, 1, Length[sols]}];
      sols = Union[Flatten[sols, 1]];
    ];

    If[info > 1, Print["  ", Length[sols], " partial solutions (gamma,rho,mu,alpha) isolated."]];

    (* 5. determine log terms and the requested number of terms of the expansion *)
    If[info > 0, Print["For each (gamma,rho,mu,alpha), determine all possible log terms"]];

    If[ Length[sols] < r || terms > 0, 
      Do[
        sols[[l]] = {FurtherTerms[rec, f[n], Append[sols[[l]], {1}], terms]};
        alpha = Most[First[sols[[l]]]];
        Catch[
          Do[AppendTo[sols[[l]], FurtherTerms[rec, f[n], Join[alpha, Table[{0},{k}], {{1}}], Max[1, terms]]]
          , {k, 1, r}]];
      , {l, 1, Length[sols]}];
      sols = Flatten[sols, 1];
    , 
      sols = Append[#, {1}]& /@ sols;
    ];

    If[info > 1, Print["  ", Length[sols], " solutions identified."]];

    If[ ret == "internal", Return[sols] ];

    (* 6. convert to expressions and return *)
    (n/E)^(#[[1]]n) *
      Exp[#[[3]].n^(Range[Length[#[[3]]]]/(Length[#[[3]]]+1))] #[[2]]^n n^#[[4]](
        Thread[Drop[#, 4] . 
        (n^(-Range[0, Length[Last[#]] - 1]/LCM[Denominator[#[[1]]], 1+Length[#[[3]]]]))]
        . Log[n]^Range[0, Length[#] - 5])& /@ sols
  ];

(* Ex: FurtherTerms[red, f[n], {gamma,phi,{},alpha,{1}}, 10] *)
FurtherTerms[rec_, f_[n_], asympt_List, desiredorder_Integer] :=
  Module[ {subs, terms, ord, myrec, poly, q, r, m, i, j, k, l, c, lc, L, sys, sol, add, o},
    subs = Take[asympt, 4]; terms = Drop[asympt, 4];
    ord = desiredorder * (q = LCM[Denominator[asympt[[1]]], Length[asympt[[3]]] + 1]);
    myrec = HomogEquation2CoefficientList[rec, f[n]]; r = Length[myrec] - 1; 
    myrec = Substitute[myrec, n, subs, ord + (add = r + 1)];
    Do[
      poly = myrec.Table[Sum[c[i, j] *
                 n^(-i/q)(1+Sum[Binomial[-i/q, k] (l/n)^k, {k, 1, q (o + add)}]) * (* = (n+l)^(-i/q) *)
                 (L - Sum[(-l/n)^k/k, {k, 1, o + add}])^(j-1) (* = Log[n+l]^(j-1) *)
              , {i, 0, o}, {j, 1, Length[terms]}], {l, 0, r}]
              /. (c[i_, j_] /; i < Length[terms[[j]]]) :> terms[[j, i + 1]];
      sys = {};
      (*
      Print["o=", o];
      Print[Collect[Expand[poly],n] /. {c -> Symbol["c"], u -> Symbol["u"], v -> Symbol["v"], w -> Symbol["w"]}];
      poly = poly /. {u->1, v->1, w->1};
      Print["lc=", LeadingTerm[poly, n]];
      *)
      Do[
        If[ o + 1 - Length[terms[[j]]] >= 0, 
          Nest[(AppendTo[sys, LeadingCoefficient[#, n]]; # - LeadingTerm[#, n])&, 
               Coefficient[poly, L, j - 1], o + 1 - Length[terms[[j]]]]]
      , {j, 1, Length[terms]}];
      sol = Solve[Thread[sys==0], Reverse[Flatten[Table[c[i,j], {j, 1, Length[terms]}, {i, 0, o}]]]];
      If[ Length[sol]===0, Throw["failed."] ];
      terms = Table[c[i,j], {j, 1, Length[terms]}, {i, 0, o}]
         /. First[sol] /. (c[i_, j_] /; i < Length[terms[[j]]]) :> terms[[j, i + 1]] /. c[__] -> 0;
    , {o, 1, ord}];
    Return[Join[subs, terms]];
  ];

(* given initial values and a recurrence, this defines a function c that 
   gives the n-th term of the sequence in linear time. *)
Rec2Def[c_Symbol, data_, rec_, f_[n_]] :=
  Module[ {u0,u1,u2,u3,kkk,myrec,ord},

     ClearAll[c];
     Do[c[kkk]=data[[kkk+1]], {kkk,0,Length[data]-1}];
     ord = Max[Cases[rec, f[n+kkk_.] -> kkk, Infinity]];
     If[ord > Length[data], Throw["insufficient initial values."]];
     myrec = Together[f[n] /. First[Solve[0==rec /. n -> n - ord, f[n]]]];
     myrec = (Factor/@Collect[Numerator[myrec], f[_]])/Denominator[myrec] 
        /. {f -> c, n -> kkk};
     u0[u2[c[kkk_Integer], kkk >= Length[data]], u1[c[kkk], u3[myrec]]] 
       /. {u0 -> SetDelayed, u1 -> Set, u2 -> Condition, u3 -> Together};
  ]; 

GuessAsymptotics[init_List, rec_, f_[n_]] :=
  Module[ {c, d, r, coeffs, poly, rho, i, j, k, roughrho, alpha, roughalpha, print},
     (* +++ under construction +++ *)
     print[__] = Null; (* print = Print *)
     print = Print;

     (* Assuming gamma=1 and mu={} and no logs, so the asymptotics is of the form c*rho^n*n^alpha *)

     (* Idea: compute the possible rho's symbolically (easy) and match them to the numerics *)
     (* so as to find the right one. Then determine symbolically the corresponding alphas *)
     (* and match them against the numerics to find the right one. *)
     (* Finally, compute a numerical estimate for the multiplicative constant c. *)

     Rec2Def[c, init, rec, f[n]];
     coeffs = HomogEquation2CoefficientList[rec, f[n]];
     r = Length[coeffs] - 1; d = Max[Exponent[coeffs, n]];
     If[ d > r + 5,
        coeffs = PolynomialQuotient[#, n^(d - r), n]& /@ coeffs;
     ];
     Clear[d];

     (* pasted from Asymptotics: exact rho's *)
     poly = LeadingCoefficient[Plus@@Substitute[coeffs, n, {0, rho, {}, 0}, 1], n];
     poly = Table[i[j[#], k], {k, 1, Exponent[#, rho]}]& /@ First /@ FactorList[poly];
     poly = Flatten[poly] /. {i -> Root, j -> Function, rho -> #};
     rho = Union[Flatten[DeleteCases[poly, 0]]];
     print["candidates for rho: ", N[rho]];

     (* match against numerics *)
     Block[ {$MaxExtraPrecision = 1000}, 
       rho = Extract[rho, Position[#, Min[#]]&[ Sum[ Abs[2c[2n+1]/c[2n] - c[n+1]/c[n] - rho]^2, {n, 20, 25}] ] ];
     ];
     If[Length[rho] =!= 1, Print[rho]; Return[$Failed]];
     rho = ToNumberField[First[rho]]; 
     print["picking ", InputForm[rho]];

     (* roots of the indicial equation for this particular rho *)
     poly = LeadingCoefficient[Plus@@Substitute[coeffs, n, {0, rho, {}, alpha}, 4], n];
     If[ FreeQ[poly, alpha], Return[$Failed] ];
     poly = First /@ FactorList[poly];
     poly = Table[i[j[#], k], {k, 1, Exponent[#, alpha]}]& /@ poly;
     poly = Flatten[poly] /. {i -> Root, j -> Function, alpha -> #};
     If[ Length[poly] <= 2, poly = RootReduce[poly] ];
     alpha = poly;
     print["candidates for alpha: ", InputForm[alpha]];

     (* match against numerics *)
     alpha = Extract[alpha, Position[#, Min[#]]&[ Sum[(Abs[n(c[n+1] - rho c[n]) / c[n]] - alpha)^2, {n, 1000, 1015}] ]];
     If[Length[alpha] =!= 1, Return[$Failed]];
     alpha = First[alpha];
     print["picking ", InputForm[alpha]];

     (* approximate multiplicative constant *)
     Clear[k];
     d[k_Integer] := c[k]/(rho^k k^alpha);
     const = N[(-d[n] + 14*d[2*n] - 56*d[4*n] + 64*d[8*n])/21 /. n -> 20, 20];

     Return[ const RootReduce[rho]^n n^alpha ]
  ];

End[]

EndPackage[]


If[$Notebooks,
   CellPrint[Cell[#, "Print", FontColor -> RGBColor[0, 0, 0], 
                              CellFrame -> 0.5, 
                              Background -> RGBColor[0.796887, 0.789075, 0.871107]]]&,
   Print
  ]["Asymptotics  Package by Manuel Kauers " <> 
    "\[LongDash] \[Copyright] RISC Linz \[LongDash] V " <> Asymptotics`Private`Version];
