finding a^b^c^... mod m

Question

a^bc mod m = a^{bc mod n} mod m, where n = φ(m) Euler's totient function.

If m is prime, then n = m-1.

Edit: as Nabb pointed out, this only holds if a is coprime to m. So you would have to check this first.

The answer does not contain full formal mathematical proof of correctness. I assumed that it is unnecessary here. Besides, it would be very illegible on SO, (no MathJax for example).
I will use (just a little bit) specific prime factorization algorithm. It's not best option, but enough.

tl;dr

We want calculate a^x mod m. We will use function modpow(a,x,m). Described below.

1. If x is small enough (not exponential form or exists p^x | m) just calculate it and return
2. Split into primes and calculate p^x mod m separately for each prime, using modpow function
   1. Calculate c' = gcd(p^x,m) and t' = totient(m/c')
   2. Calculate w = modpow(x.base, x.exponent, t') + t'
   3. Save pow(p, w - log_p c', m) * c' in A table
3. Multiple all elements from A and return modulo m

Here pow should look like python's pow.

Main problem:

Because current best answer is about only special case gcd(a,m) = 1, and OP did not consider this assumption in question, I decided to write this answer. I will also use Euler's totient theorem. Quoting wikipedia:

Euler's totient theorem:
If n and a are coprime positive integers, then where φ(n) is Euler's totient function.

The assumption numbers are co-primeis very important, as Nabb shows in comment. So, firstly we need to ensure that the numbers are co-prime. (For greater clarity assume x = b^(c^...).) Because , where we can factorize a, and separately calculate q1 = (p1^alpha)^x mod m,q2 = (p2^beta)^x mod m... and then calculate answer in easy way (q1 * q2 * q3 * ... mod m). Number has at most o(log a) prime factors, so we will be force to perform at most o(log a) calculations.

In fact we doesn't have to split to every prime factor of a (if not all occur in m with other exponents) and we can combine with same exponent, but it is not noteworthy by now.

Now take a look at (p^z)^x mod m problem, where p is prime. Notice some important observation:

If a,b are positive integers smaller than m and c is some positive integer and , then true is sentence .

Using the above observation, we can receive solution for actual problem. We can easily calculate gcd((p^z)^x, m). If x*z are big, it is number how many times we can divide m by p. Let m' = m /gcd((p^z)^x, m). (Notice (p^z)^x = p^(z*x).) Let c = gcd(p^(zx),m). Now we can easily (look below) calculate w = p^(zx - c) mod m' using Euler's theorem, because this numbers are co-prime! And after, using above observation, we can receive p^(zx) mod m. From above assumption wc mod m'c = p^(zx) mod m, so the answer for now is p^(zx) mod m = wc and w,c are easy to calculate.

Therefore we can easily calculate a^x mod m.

Calculate a^x mod m using Euler's theorem

Now assume a,m are co-prime. If we want calculate a^x mod m, we can calculate t = totient(m) and notice a^x mod m = a^(x mod t) mod m. It can be helpful, if x is big and we know only specific expression of x, like for example x = 7^200.

Look at example x = b^c. we can calculate t = totient(m) and x' = b^c mod t using exponentiation by squaring algorithm in Θ(log c) time. And after (using same algorithm) a^x' mod m, which is equal to solution.

If x = b^(c^(d^...) we will solve it recursively. Firstly calculate t1 = totient(m), after t2 = totient(t1) and so on. For example take x=b^(c^d). If t1=totient(m), a^x mod m = a^(b^(c^d) mod t1), and we are able to say b^(c^d) mod t1 = b^(c^d mod t2) mod t1, where t2 = totient(t1). everything we are calculating using exponentiation by squaring algorithm. Note: If some totient isn't co-prime to exponent, it is necessary to use same trick, as in main problem (in fact, we should forget that it's exponent and recursively solve problem, like in main problem). In above example, if t2 isn't relatively prime with c, we have to use this trick.

Calculate φ(n)

Notice simple facts:

1. if gcd(a,b)=1, then φ(ab) = φ(a)*φ(b)
2. if p is prime φ(p^k)=(p-1)*p^(k-1)

Therefore we can factorize n (ak. n = p1^k1 * p2^k2 * ...) and separately calculate φ(p1^k1),φ(p2^k2),... using fact 2. Then combine this using fact 1. φ(n)=φ(p1^k1)*φ(p2^k2)*...

It is worth remembering that, if we will calculate totient repeatedly, we may want to use Sieve of Eratosthenes and save prime numbers in table. It will reduce the constant.

python example: (it is correct, for the same reason as this factorization algorithm)

def totient(n) :          # n - unsigned int
    result = 1
    p = 2                 #prime numbers - 'iterator'
    while p**2 <= n :
        if(n%p == 0) :    # * (p-1)
            result *= (p-1)
            n /= p
        while(n%p == 0) : # * p^(k-1)
            result *=  p
            n /= p
        p += 1
    if n != 1 :
        result *= (n-1)
    return result         # in O(sqrt(n))

Case: abc mod m

Cause it's in fact doing the same thing many times, I believe this case will show you how to solve this generally.
Firstly, we have to split a into prime powers. Best representation will be pair <number, exponent>.
c++11 example:

std::vector<std::tuple<unsigned, unsigned>> split(unsigned n) {
  std::vector<std::tuple<unsigned, unsigned>> result;
  for(unsigned p = 2; p*p <= n; ++p) {
    unsigned current = 0;
    while(n % p == 0) {
      current += 1;
      n /= p;
     }
    if(current != 0)
     result.emplace_back(p, current);
   }
  if(n != 1)
   result.emplace_back(n, 1);
  return result;
 }

After split, we have to calculate (p^z)^(b^c) mod m=p^(z*(b^c)) mod m for every pair. Firstly we should check, if p^(z*(b^c)) | m. If, yes the answer is just (p^z)^(b^c), but it's possible only in case in which z,b,c are very small. I believe I don't have to show code example to it.
And finally if p^(z*b^c) > m we have to calculate the answer. Firstly, we have to calculate c' = gcd(m, p^(z*b^c)). After we are able to calculate t = totient(m'). and (z*b^c - c' mod t). It's easy way to get an answer.

function modpow(p, z, b, c, m : integers) # (p^z)^(b^c) mod m
    c' = 0
    m' = m
    while m' % p == 0 :
        c' += 1
        m' /= p
    # now m' = m / gcd((p^z)^(b^c), m)
    t = totient(m')
    exponent = z*(b^c)-c' mod t
    return p^c' * (p^exponent mod m')

And below Python working example:

def modpow(p, z, b, c, m) : # (p^z)^(b^c) mod m
    cp = 0
    while m % p == 0 :
        cp += 1
        m /= p              # m = m' now
    t = totient(m)
    exponent = ((pow(b,c,t)*z)%t + t - (cp%t))%t
                            # exponent = z*(b^c)-cp mod t
    return pow(p, cp)*pow(p, exponent, m)

Using this function, we can easily calculate (p^z)^(b^c) mod m, after we just have to multiple all results (mod m), we can also calculate everything on an ongoing basis. Example below. (I hope I didn't make mistake, writing.) Only assumption, b,c are big enough (b^c > log(m) ak. each p^(z*b^k) doesn't divide m), it's simple check and I don't see point to make clutter by it.

def solve(a,b,c,m) : # split and solve
    result = 1
    p = 2            # primes
    while p**2 <= a :
        z = 0
        while a % p == 0 :
                     # calculate z
            a /= p
            z += 1
        if z != 0 :
            result *=  modpow(p,z,b,c,m)
            result %= m
        p += 1
    if a != 1 :      # Possible last prime
        result *= modpow(a, 1, b, c, m)
    return result % m

Looks, like it works.
DEMO and it's correct!

1. Since for any relationship a=x^y, the relationship is invariant with respect to the numeric base you are using (base 2, base 6, base 16, etc).

2. Since the mod N operation is equivalent to extracting the least significant digit (LSD) in base N

3. Since the LSD of the result A in base N can only be affected by the LSD of X in base N, and not digits in higher places. (e.g. 34*56 = 30*50+30*6+50*4+4*5 = 10*(3+50+3*6+5*4)+4*6)

Therefore, from LSD(A)=LSD(X^Y) we can deduce

LSD(A)=LSD(LSD(X)^Y)

Therefore

A mod N = ((X mod N) ^ Y) mod N

and

(X ^ Y) mod N = ((X mod N) ^ Y) mod N)

Therefore you can do the mod before each power step, which keeps your result in the range of integers.

---

This assumes a is not negative, and for any x^y, a^y < MAXINT

---

This answer answers the wrong question. (alex)