Group Abstract Group Abstract

Message Boards Message Boards

1
|
88 Views
|
2 Replies
|
2 Total Likes
View groups...
Share
Share this post:

Self-referential recurrence

Brute force computation of the sequence A000966

5,11,17,23,29,30,36,42,48,54,60,61,67,73,79,85,91,92,98,104,110,116,122,123,129,135, …

is straightforward. Define e, f , and g:

e = n |-> (5^n-1)/4;
f = n |-> (1-x^(e[n]-1))/(1-x^e[n-1]);
t = n |-> x^(e[n]-6);  

Then the series generating recurrence is

a[2] = 1;
a[n_/;n>2]:= a[n ]= f[n]a[n-1] + t[n] 

And we can compute the sequence—up to some order—as follows:

Position[Rest@CoefficientList[x^5 a[4] + O[x]^200, x],1] // Flatten

It is clear that the differences of the terms in this sequence is either 6 or 1. What is interesting—and does not appear to have been noted before—is that the position of the 1’s is self-referentially indexed by the sequence itself: the first 1 is in the position 5 of the difference sequence, the second 1 is in position 11, and so on. So there should be a “simple” way of “self-referentially” computing this series, starting with 5 and adding 6 or 1 (in positions depending upon the coefficients already determined). However, I’ve not been able to write elegant code to do this yet—so looking for solutions!

POSTED BY: Paul Abbott
2 Replies

I'm not sure any two people have exactly the same notions of elegance, but it is self-referential without infinite recursion:

A000966 // ClearAll;
Begin["A000966`"];
ClearAll[iA000966, plusX];
(* make a plus-x def & queue a future plus-one *)
plusX[n_, x_] :=
  iA000966[n] :=  (* SetDelayed[] delays future plus-one *)
   iA000966[n] =  (* memoize def for n *)
    With[{res = iA000966[n - 1] + x},
     plusX[res, 1]; res  (* queue def for n=res; return res *)
     ];(* init *)
A000966[0] = iA000966[0] = (plusX[5, 1]; 5);
(* make a plus-six def & queue a future plus-one *)
iA000966[n_] := (plusX[n, 6]; iA000966[n]);
(* must define sequence from bottom up *)
A000966[n_Integer?Positive] /; IntegerQ[A000966- 1]] := iA000966[n]; 
End[];

test = Position[Rest@CoefficientList[x^5 a[6] + O[x]^1, x], 1] // Flatten;

Table[A000966[k], {k, 0, Length[test] - 1}] == test
(*True  *)

You have to define the sequence sequentially from the beginning or a plus-one definition may be skipped. The IntegerQ[] test makes sure the previous term has been defined; if not, it triggers a definition of the previous term, first checking the one before that; and so on back down to the beginning if necessary. The double initialization A000966[0] = iA000966[0] =... makes sure there is a beginning at which to stop.

One might think a more straightforward code like the following would work:

A000966[n_Integer?Positive] /; IntegerQ[A000966[n - 1]] :=
 (plusX[n, 6]; A000966[n])

But if A000966[30] is the first evaluation performed, its value turns out to be 180. If instead, one evaluates Do[A000966[k], {k, 30}], then its value turns out to be the correct 155. The problem is this: The definition of A000966[30] has started when A000966[5] is first called. The definition of A000966[5] queues a plus-one definition for A000966[30]. But this definition is overwritten by a plus-six definition when the recursive checking returns to finish the definition for A000966[30] that had started. So in the actual code used, the recursive checking is done on the symbol A000966 and the appropriate plus-one and plus-six definitions are made for the symbol iA000966.

A minor weakness is that the recursive checking might cause $RecursionLimit error. For instance if the first call were A000966[2025].

POSTED BY: Michael Rogers
recur = RSolveValue[{f[n] == f[n - 1] + f[n - 5] - f[n - 6], 
    f[1] == 6, f[2] == 11, f[3] == 17, f[4] == 23, f[5] == 29, 
    f[6] == 30}, f[n], n];

It's long and ugly. Shorter but still ugly:

srecur = FullSimplify[recur]

(* (4^-n (6965 I 2^(1/2 + 2 n) Sqrt[5 + Sqrt[5]] + 
     3115 I 2^(1/2 + 2 n) Sqrt[5 (5 + Sqrt[5])] + 
     15 4^(1 + n)
       n (246 + 110 Sqrt[5] + 
        Root[1280 + 2292560 #1^2 + #1^4 &, 2]) + 
     2 (-1 - Sqrt[5] + Root[80 + 20 #1^2 + #1^4 &, 1])^
      n (1040 + 465 Sqrt[5] + 
        Root[15125 + 15278650 #1^2 + #1^4 &, 2]) + 
     2 (35 4^n (123 + 55 Sqrt[5]) + (-1 + Sqrt[5] + 
           Root[80 + 20 #1^2 + #1^4 &, 3])^
         n (4025 + 1800 Sqrt[5] + 
           Root[120125 + 18562450 #1^2 + #1^4 &, 2])) - 
     2 (-1 + Sqrt[5] + I Sqrt[2 (5 + Sqrt[5])])^
      n (1990 + 890 Sqrt[5] + 
        Root[162000 + 67522500 #1^2 + #1^4 &, 1]) + (-1 - Sqrt[5] + 
        Root[80 + 20 #1^2 + #1^4 &, 2])^
      n (615 + 275 Sqrt[5] + 
        Root[15842000 + 76903700 #1^2 + #1^4 &, 4])))/(25 (123 + 
     55 Sqrt[5] + Root[80 + 573140 #1^2 + #1^4 &, 2])) *)

Table[N@srecur, {n, 30}] // Chop

(* Out[10]= {6., 11., 17., 23., 29., 30., 35., 41., 47., 53., 54., 59., \
65., 71., 77., 78., 83., 89., 95., 101., 102., 107., 113., 119., \
125., 126., 131., 137., 143., 149.} *)
POSTED BY: Daniel Lichtblau
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard