What I'm actually trying to do is efficiently calculate Faulhaber(N, E) modulo M (or at least Faulhaber(N, E) - the modulo operation could be worked in later), where N, E, and M can be arbitrarily large. From what I have read, it seems that this is generally expressed in terms of either a polynomial to the degree of N+1 (memory-expensive) or some recursive relationship (time-expensive). Is there any way that this can be done that requires neither very much memory nor time?
