JavaScript (Babel Node), 59 bytes

n=>g=([c,a,b])=>(a=a?g(a):c,b=b?g(b):c,[a+b,a**b,a%b,n][c])

Try it online!

This is using @xnor's approach, but without eval. I'm using Javascript because Python is too careful with its arithmetic (there are exceptions when doing mod-by-zero).

The simple datatype D is a BigInt or an array of D. Given a program p (a D) and an input n (a BigInt), the expression e(n)(p) interprets the program with input n. Programs are interpreted as follows:

• [0, a, b] evaluates a and b, returning their sum
• [1, a, b] evaluates a and b, returning a to the power of b
• [2, a, b] evaluates a and b, returning a modulo b
• [3] returns n

For example, the following program outputs the constant value of two (using @xnor's observations):

[0, [1, [2, [1, [3], [3]],
            [1, [3], [3]]],
        [2, [1, [3], [3]],
            [1, [3], [3]]]],
    [1, [2, [1, [3], [3]],
            [1, [3], [3]]],
        [2, [1, [3], [3]],
            [1, [3], [3]]]]]

It is such a surprising result that all recursively enumerable functions can be represented by such a program! This is Corollary 3 of the following paper, which @xnor cited.

Marchenkov, S. S. (2007). Superpositions of elementary arithmetic functions. Journal of Applied and Industrial Mathematics, 1(3), 351–360. doi:10.1134/s1990478907030106