Very nice. In terms of personal preferences, I would like it better if the user weren't responsible for evaluating d[] in order and up to a sufficiently large limit. It's hard to do that without marring the elegance. And it's not that inconvenient to pre-evaluate d[].

One way:

Clear[a, b, d];
b[1] = a[1] = 5;
d[1] = 0; (* a recursion stopper; value is irrelevant *)
d[n_] := (d[n - 1];
   With[{m = a[n - 1]}, d[n] = a[m + 1] = a[m] + 1]); 
a[n_] := a[n] = a[n - 1] + 6;
b[n_] := (d[n - 1]; a[n]);

Then b[100] yields 504.

Similar, but without recursion-limit issues (could use a variable to keep track of the max n on which d[] has been evaluated instead of the number of down-values):

Clear[a, b, d];
b[1] = a[1] = 5;
d[n_] := With[{m = a[n - 1]}, d[n] = a[m + 1] = a[m] + 1]; 
a[n_] := a[n] = a[n - 1] + 6;
b[n_] := (d /@ Range[Length@DownValues[d] + 1, n - 1]; a[n]);

Even less elegant but more efficient:

Clear[a, b, d];
b[1] = a[1] = 5;
d[1] = 0;
d[n_] := With[{m = a[n - 1]}, d[n] = a[m + 1] = a[m] + 1]; 
a[n_] := a[n] = a[n - 1] + 6;
bdef = Function[b[n_] := (Replace[
      Hold[d[#]] //. 
       Hold[d[k_]] /; d[k] < Min[6 n, 5 n + 50] :> Hold[d[k + 1]],
      Hold[d[l_]] /; l != # :> bdef[l]]; a[n])];
bdef[1];

b[100]
(*  504  *)

Table[b[k] - 5 k, {k, 1000000}] // Max // AbsoluteTiming
(*  {6.60977, 33}  *)

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]

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].