Group Abstract Group Abstract

Message Boards Message Boards

3
|
2.3K Views
|
8 Replies
|
11 Total Likes
View groups...
Share
Share this post:

Self-referential recurrence

POSTED BY: Paul Abbott
8 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[n - 1]] := iA000966[n]; 
End[];

test = Position[Rest@CoefficientList[x^5 a[6] + O[x]^1000, 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
POSTED BY: Daniel Lichtblau

The recurrence is not correct—and it cannot be that simple: the positions of 1 in the differences has to equal the sequence itself, but your answer fails to do that for n=31. Note that there is no simple recurrence given for the sequence A000966—so if you could find one, that would be great!

POSTED BY: Paul Abbott
POSTED BY: Michael Rogers

Thanks, Michael:

I ended up coding it this way:

a[1]=5; 
d[n_]:= a[n+1] = a[n] + 1; (* unit differences *)
a[n_]:= a[n] = a[n-1] + 6; 

Then, computing the unit differences, starting from n=1 and working upwards, causes dynamic programming to compute all the required undetermined a[n] via the second recurrence:

(n |-> d[a[n]]) /@ Range[40]; 
a /@ Range[a[40]+ 1]
POSTED BY: Paul Abbott

Oops. I simply misunderstood. Agreed, a closed recurrence seems unlikely. (A recurrence of misunderstanding, that’s a different matter.)

POSTED BY: Daniel Lichtblau

Sorry about the 6/5 typo. Below is a better method anyway.

ff = RSolveValue[{f[1] == 5, f[n] == f[n - 1] + Piecewise[{{1, Mod[n, 6] == 0}}, 6]}, f[n], n]

(* Out[53]= 30 + Floor[1/6 (-6 + n)] + 6 Floor[1/6 (-5 + n)] + 
 6 Floor[1/6 (-4 + n)] + 6 Floor[1/6 (-3 + n)] + 
 6 Floor[1/6 (-2 + n)] + 6 Floor[1/6 (-1 + n)] *)

Check it:

Table[ff, {n, 30}]

(* Out[54]= {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, 141, 147, 153, 154} *)
POSTED BY: Daniel Lichtblau

Thanks, Daniel. However, f[1] = 5 and one should get f[7] = 36.

I have noted that the (direct) OEIS code by Noe for A000966 can be simplified:

u = Union @ Accumulate @ IntegerExponent[Range[1000], 5]; 
Complement[Range[Last @ u], u]
POSTED BY: Paul Abbott
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard