快速计算 (a^b) mod n

需求：计算(a^d) mod n，其中a、d、n都是超大数，比如 :

a=29033793948611024383985226589380327268410365915345045221422018356137431660932535186202642526599446262022568978622521153097299581568679434816825712461058592712562397012627020575723923589246658426796073196472141192591117528362339550638668003488192557371657081780684417377855541134663485614245972808435461007434
d=222432460867810830547268623836384189159087209512350317868716836088090819194176604003161569409739044308971343425911343340527149421824735849336260097881116010135275955029797307117840761229961305637595400241553817104545776652334554538872529742802718662643385866831969462915848856054189567689190451958864625233
n=97017536198236097887027810122515410619659485697283600160144570596909548609965615008038739116646997559005845248037584107688077248878768279884800118900545837204807618610350247395446183630604490529251059570366856350131995095330465544789339594760729309518819101385700485986350835403338285335501732735378685171873

直接计算程序卡死。

由于(a^d) mod n可以展开为 [(a%n)(a%n)...(a%n)] % n ， (a+b) mod n可以展开为(a%n+b%n) % n。同时，由于在本地测试，发现 a^1000秒回，所以把 a^d 拆成 [(a^1000)^x]*(a^(d%1000))，其中x=d//1000。该过程的Python实现代码如下：

def powM(a, d, n):
    # step1, set result = (a^d) mod n = ((a^5)^x * (a^y)) % n = ((((a^5) %n )^m) * ((a^y)%n)) % n
    # step2, set b = (a^5) % n, c = (a^y)%n
    # step3, get result = ((b^x)*c) % n
    # if m is bigger than 5, then, we can repeate step1~step3
    # finally we get the result
    step = 5
    if d > step:
        x = d // step
        y = d % step
        # d % step 余数处理
        c = (a**y ) %n
        b = (a**step) % n
        if x > step:
            return (powM(b, x, n)*c)%n
        else:
            return (((b**x)%n)*c) % n
    else:
        return (a**d) % n

可以使用上述a, d, n进行测试，效率大幅提高。
经在服务器上测试，发现step=5时，效率最高