Recursive function with if-statement

I am trying to represent the following function definition in Mathematica:

\begin{align*} f(1)&=1 \\ f(2n)&= \begin{cases}f(n) & \text{if}\space n\equiv0\pmod{2} \\ 2f(n) & \text{if}\space n\equiv1\pmod{2}\end{cases} \\ f(2n+1)&=\begin{cases}2f(n)+1 & \text{if}\space n\equiv0\pmod{2} \\ f(n) & \text{if}\space n\equiv1\pmod{2}\end{cases} \end{align*}

I have tried both the following definitions, but both of them return f[2011] as their output when I input f[2011] into the kernel.

Using if-statements:

f[1] := 1
f[2 n_] := If[Mod[n, 2] == 0, f[n], 2 f[n]]
f[2 n_ + 1] := If[Mod[n, 2] == 0, 2 f[n] + 1, f[n]]


And using switch-statements:

f[1] := 1
f[2 n_] := Switch[Mod[n, 2], 0, Evaluate[f[n]], 1, 2 Evaluate[f[n]]]
f[2 n_ + 1] := Switch[Mod[n, 2], 0, 2 Evaluate[f[n]] + 1, 1, Evaluate[f[n]]]


Neither of these definitions seem to give me what I'm after, so I'd appreciate any help.

-

I'd have done something like this:

f[1] = 1;
f[n_Integer?EvenQ] := f[n] =
Block[{m = n/2}, (1 + Boole[Mod[m, 2] == 1]) f[m]];
f[n_Integer?OddQ] := f[n] =
Block[{m = (n - 1)/2, t}, t = Boole[Mod[m, 2] == 0]; (1 + t) f[m] + t]


The use of both := and = in the odd and even cases is sometimes referred to as memoization; there are a number of threads on this site (e.g. the one just linked to) that discuss that in more detail.

In any event,

AbsoluteTiming[f[2011]]
{0.000785, 21}


If one wanted to compress things as with the other answers given,

f[1] = 1;
f[n_Integer] := f[n] =
Block[{m = Quotient[n, 2], t}, t = Mod[m, 2];
If[EvenQ[n], (1 + t) f[m], (2 - t) f[m] - t + 1]]

-
Thank you very much! I'll look up memoization when I get time. Cheers. –  Shaktal Jun 13 '12 at 14:29
@Shaktal: for further fun, execute ??f after that last evaluation of f[2011]. –  Guess who it is. Jun 13 '12 at 15:07

While you've already been told how you can do it, you haven't been told yet why your way doesn't work.

When you write f[something] := something_else what you define is not actually a function, but a pattern replacement rule. Let's look for example at your first try:

f[1] := 1
f[2 n_] := If[Mod[n, 2] == 0, f[n], 2 f[n]]
f[2 n_ + 1] := If[Mod[n, 2] == 0, 2 f[n] + 1, f[n]]


Here, the first definition tells Mathematica that whenever it finds the expression f[1], it shall replace it with 1. This part of your definition works as expected.

The second definition tells Mathematica that if it finds any expression which matches the pattern f[2 n_] then it shall replace it with the expression If[Mod[n, 2] == 0, f[n], 2 f[n]] where n has been replaced by whatever expression matched the n_ part. OK, so what does it mean to match that expression? Well, it means that your expression (after being evaluated) has exactly this form. So for example, assuming x was not assigned a value, then if you type f[2 x] then Mathematica will note that if in the pattern f[2 n_] n_ is replaced by x, then it will get f[2 x], therefore it has a match, and therefore replaces the right hand side and continues evaluation.

So what happens with e.g. f[4]? You might be inclined to say that if you replace n with 2, you'd get f[4] from f[2 n_], but that is not true. What you get by this replacement is f[2*2]. Further evaluation of that expression of course would eventually give f[4], but that's not how Mathematica works (and BTW, further evaluation would also include evaluation of f[4] so you'd get an endless loop otherwise, except if making an extra exception for that). Now f[4] is clearly not the same expression as f[2 2], and therefore you don't get a match. Since you don't get a match for the other two definitions, no further evaluation is possible, and you get f[4] as result.

Note that f[2*2] would not invoke that definition either because the argument is evaluated before trying to match the function, and since that evaluation results in 4, there's again no match. You could, however, explicitly tell Mathematica to not evaluate that argument before pattern matching, by using f[Unevaluated[2*2]]. In that case, you will get a match for the second definition, with n being 2. Thus the right hand side would evaluate to f[2] which would be the final result because again there's no matching rule.

Note however that even if Mathematica worked differently, your definition would still not be meaningful because there's no built-in assumption that a variable called n has to be an integer (and any such assumption would be a bad idea anyway, since n could e.g. describe the refractive index of a material).

-
Ahh, a very interesting read, thanks for enlightening me! I'll bear this in mind for future reference, thanks! –  Shaktal Jun 13 '12 at 17:59

A more streamlined version (after re-writing the functions involved in equivalent forms):

h[1] = 1;
h[n_?IntegerQ] :=  h[n] =
If[OddQ[n],
2 h[(n - 1)/2] + 1 - Mod[(n - 1)/2, 2] (h[(n - 1)/2] + 1),
h[n/2] + Mod[n/2, 2] h[n/2]]

-
Even more streamlined: h[n_Integer] := h[n] = Module[{nh = n ~Quotient~ 2}, Module[{m = n ~Mod~ 2, mh = nh ~Mod~ 2}, m (1 - mh) + h[nh] ((1 + m)(1 + mh) - 3 m mh)]] (indeed, given that now there's only a single recursive call, I think memoization doesn't help much for this version). –  celtschk Jun 13 '12 at 16:13
I omitted the h[1] = 1, but with that added, it should be identical. –  celtschk Jun 13 '12 at 16:29
@celtschk I made a ninja edit of my comment; it originally said that I was about to post a similar solution (to your comment version). Then I realized they were not the same. –  Mr.Wizard Jun 13 '12 at 16:31
Ah, OK, I understood your comment to say that my code does not give the same result. –  celtschk Jun 13 '12 at 16:33
@celtschk, really nice! –  kguler Jun 13 '12 at 17:00

My version:

f[1] = 1
f[n_Integer?Positive] := f[n] =
With[{x = n ~Mod~ 2},
If[OddQ[# + x], 2 f@# + x, f@#] &[(n - x)/2]
]


maximally terse (code golf):

f@1=1;f@n_:=If[OddQ@+##,2f@#+#2,f@#]&[(n-#)/2,#]&@Mod[n,2]

-
It's not maximally terse unless you write it as f@n_ := f@n = ... ;) –  The Toad Jun 13 '12 at 16:25
@R.M I almost did but I thought it was too much. Now you've thrown down the gauntlet. :-> –  Mr.Wizard Jun 13 '12 at 16:26
@R.M there, how's that? :D –  Mr.Wizard Jun 13 '12 at 16:29
@R.M no comment on my Golf? I'm hurt. –  Mr.Wizard Jun 14 '12 at 14:06
Heh, it's no secret that I'm not the biggest fan of golf ;) I'm in it purely for the infix jokes :D nevertheless, that's pretty neat! –  The Toad Jun 14 '12 at 15:00

After taking a step back from the given definition, it seems to me that the Mod[...,2]' aren't quite capturing the real periodicity interval of the four case distinctions that are required here.

A more natural approach would be to use Mod[...,4] and make the case distinction based on the four possible values of the remainder.

If $m$ is an integer, the four cases are that the argument is equal to $2(2m)$ (even), $2(2m)+1$ (odd), $2(2m+1)$ (even) or $2(2m+1)+1$ (odd). The remainders of these under division by $4$ are $0$, $1$, $2$ and $3$. And here, $2m$ is even while $(2m+1)$ is odd. Then you can write down the cases in a (to me) more readable way:

Clear[f, m, rule]

rule = {
{m_, 0} :> f[2 m],
{m_, 1} :> 2 f[2 m] + 1,
{m_, 2} :> 2 f[2 m + 1],
{m_, 3} :> f[2 m + 1]
};

f[0] = 0;
f[x_] := QuotientRemainder[x, 4] /. rule


The recursion is now both readable and fast, even without memoization. With recursion relations, readability is more important to me than terseness.

A little more explanation (so that the answer looks less terse...): I implemented the "If statement" in the form of (delayed) rules that depend on the output of the function QuotientRemainder[x, 4]. That output is a list whose first argument is the integer $m$ I introduced earlier, and the second argument is the remainder on which the case distinction is based.

Therefore, the left-hand sides of the rules (e.g., {m_, 0} :> f[2 m]) are tuples, where now m_ appears instead of m because it is a placeholder for the actual value Quotient[x, 4] when f[x] is called. On the right-hand side of the rules, 2 m or 2 m +1 appears in place of the even or odd integer n in the original formulation.

P.S.

One could also implement conditions using the Switch statement, as in

Clear[h];
h[0] = 0;
h[x_] := Switch[#2,
0, h[2 #1],
1, 2 h[2 #1] + 1,
2, 2 h[2 #1 + 1],
3, h[2 #1 + 1]
] & @@
QuotientRemainder[x, 4]


This follows the same idea as in my first approach, but it doesn't rely on rules. Which of these is clearer may be a matter of taste.

-

Here's one using conditions on the SetDelayed:

f[1] = 1;
f[n_?EvenQ] := (f[n] = f[n/2])           /; Mod[n/2, 2] == 0;
f[n_?EvenQ] := (f[n] = 2 f[n/2])         /; Mod[n/2, 2] == 1;
f[n_?OddQ]  := (f[n] = 2 f[(n-1)/2] + 1) /; Mod[(n-1)/2, 2] == 0;
f[n_?OddQ]  := (f[n] = f[(n-1)/2])       /; Mod[(n-1)/2, 2] == 1;

-
Probably! But I was aiming for something similar to how to original equations were written. –  wxffles Jun 13 '12 at 20:05
Absolutely. Which brings to mind a question which I might ask. –  wxffles Jun 13 '12 at 20:17
Would that question relate to this? –  Mr.Wizard Jun 13 '12 at 20:19
Yes, but not nearly that complicated. All I want is a space width the same as a character width. –  wxffles Jun 13 '12 at 20:23
I presume you mean in a standard Input cell rather than a Code cell? –  Mr.Wizard Jun 13 '12 at 20:24

I am late but just for fun another way:

rcv[1] = 1
rcv[n_?IntegerQ] := (rcv[n] = rcv[Floor[n/2]]) /;
Or[Mod[n, 4] == 0, Mod[n, 4] == 3]
rcv[n_?IntegerQ] := rcv[n] = 2 rcv[Floor[n/2]] + 2 - Mod[n, 4]


Comparing with Mr. Wizard f (I could have used any of the answers):

f[1] = 1
f[n_Integer?Positive] :=
f[n] = With[{x = n~Mod~2},
If[OddQ[# + x], 2 f@# + x, f@#] &[(n - x)/2]]


testing:

Table[f[j], {j, 1000}] == Table[rcv[j], {j, 1000}]


yields True`

-