MojoでProject Euler 35

https://projecteuler.net/problem=35

エラトステネスの篩を使えば簡単ですが、これだとせいぜい10^9までです。
2桁以上だと1, 3, 7, 9しか使えないので、これらからなる整数を作ってMiller-Rabinで素数判定すればいいです。
さらに攻めると、その整数の中で回転しても最小のものだけ作ればよいです。1を含むなら最上位は1が続いて次が3, 7, 9で、次が1で、という具合になります。これを各パートの長さで表します。2, 1, 1, 1なら、11317などで、確実に最小になります。逆に1, 1, 2, 1なら13117などで、必ず最小にはなりません。1, 1, 1, 2なら13177や17133があって最小になるか分からないので回して最小になるものだけ残します。パートは偶数個あります。なぜなら、奇数だと最後が1になって、最小でなくなるからです。例外はRepunitと呼ばれる1が並んだ整数です。
10^15までなら約30秒でした。
いちおう10^18まで計算したのですが、新たなCircular Primeは出てきませんでした。19桁のRepunitが素数のはずなので、それがでてくるはずですが、19桁だとオーバーフローするので計算は難しいです。

from collections import Set
import sys


#################### library ####################

fn initialize_list(N: Int, init: Bool) -> List[Bool]:
    var a = List[Bool](capacity=N)
    for n in range(N):
        a.append(init)
    return a

fn copy_list(a: List[Int]) -> List[Int]:
    var b = List[Int]()
    for e in a:
        b.append(e[])
    return b

fn sub_list(a: List[Int], first: Int, last: Int) -> List[Int]:
    var b = List[Int]()
    for i in range(first, last):
        b.append(a[i])
    return b

fn connect_list(a: List[Int], b: List[Int]) -> List[Int]:
    var a1 = copy_list(a)
    a1.extend(b)
    return a1

fn print_list(a: List[Int]):
    if a.size > 0:
        var s = str(a[0])
        for i in range(1, a.size):
            s += " " + str(a[i])
        print(s)
    else:
        print()

fn sum_list(a: List[Int]) -> Int:
    var s = 0
    for e in a:
        s += e[]
    return s

fn add(a: Int, b: Int, d: Int) -> Int:
    var c = a + b
    if c < 0:   # overflow
        c -= d
    elif c >= d:
        c -= d
    return c

fn mul2(a: Int, b: Int, d: Int, e1: Int, e2: Int) -> Int:
    var n = 0
    var a1 = a
    for i in range(0, e2, e1):
        var b1 = (b >> i) & ((1 << e1) - 1)
        n += a1 * b1 % d
        a1 = (a1 << e1) % d
    return n % d

fn mul3(a: Int, b: Int, d: Int, e1: Int, e2: Int) -> Int:
    var n = 0
    var a1 = a
    for i in range(0, e2, e1):
        var b1 = (b >> i) & ((1 << e1) - 1)
        n = add(n, a1 * b1 % d, d)
        a1 = (a1 << e1) % d
    return n

fn mul4(a: Int, b: Int, d: Int) -> Int:
    var n = 0
    var a1 = a
    for i in range(0, 63):
        var b1 = (b >> i) & 1
        if b1 == 1:
            n = add(n, a1, d)
        a1 = add(a1, a1, d)
    return n

fn mul(a: Int, b: Int, d: Int) -> Int:
    if d < 1 << 31:
        return a * b % d
    elif d < 1 << 40:
        var b1 = b & ((1 << 20) - 1)
        var b2 = b >> 20
        var n1 = a * b1
        var n2 = (a * b2 % d) << 20
        return (n1 + n2) % d
    elif d < 1 << 45:
        return mul2(a, b, d, 15, 45)
    elif d < 1 << 48:
        return mul2(a, b, d, 12, 48)
    elif d < 1 << 50:
        return mul2(a, b, d, 10, 50)
    elif d < 1 << 54:
        return mul3(a, b, d, 6, 54)
    elif d < 1 << 55:
        return mul3(a, b, d, 5, 55)
    elif d < 1 << 56:
        return mul3(a, b, d, 4, 56)
    elif d < 1 << 57:
        return mul3(a, b, d, 3, 57)
    elif d < 1 << 60:
        return mul3(a, b, d, 2, 60)
    else:
        return mul4(a, b, d)

fn pow(n: Int, e: Int, d: Int) -> Int:
    if e == 0:
        return 1
    elif e == 1:
        return n
    elif (e & 1) == 1:
        return mul(n, pow(n, e-1, d), d)
    else:
        var m = pow(n, e>>1, d)
        return mul(m, m, d)

fn Miller_test(n: Int, owned d: Int, a: Int) -> Bool:
    var y = pow(a, d, n)
    if y == 1:
        return True
    d <<= 1
    while d < n:
        if y == n - 1:
            return True
        y = mul(y, y, n)
        d <<= 1
    return False

fn Miller_Rabin(n: Int) -> Bool:
        var d = (n - 1) >> 1
        while (d & 1) == 0:
            d >>= 1
        
        if not Miller_test(n, d, 2):
            return False
        if not Miller_test(n, d, 3):
            return False
        if not Miller_test(n, d, 5):
            return False
        if not Miller_test(n, d, 7):
            return False
        if not Miller_test(n, d, 11):
            return False
        if not Miller_test(n, d, 13):
            return False
        if not Miller_test(n, d, 17):
            return False
        return True

fn is_prime(n: Int) -> Bool:
    if n == 2:
        return True
    elif (n & 1) == 0 or n == 1:
        return False
    else:
        for p in range(3, 20, 2):
            if p * p > n:
                return True
            elif n % p == 0:
                return False
        
        return Miller_Rabin(n)


#################### process ####################

fn divide(L: Int, d: Int) -> List[List[Int]]:
    var vs = List[List[Int]]()
    if d == 1:
        for L1 in range(1, L):
            for L2 in range(1, L-L1+1):
                var v = List[Int]()
                v.append(L1)
                v.append(L2)
                vs.append(v)
    elif d % 2 == 1:
        var vs1 = divide(L-(d-1)*2, 1)
        var vs2 = divide(L-2, d-1)
        for v1 in vs1:
            var L1 = sum_list(v1[])
            for v2 in vs2:
                var L2 = sum_list(v2[])
                if L1 + L2 <= L:
                    var v = connect_list(v1[], v2[])
                    vs.append(v)
    else:
        var half = d // 2
        var vs1 = divide(L-half, half)
        for v1 in vs1:
            var L1 = sum_list(v1[])
            for v2 in vs1:
                var L2 = sum_list(v2[])
                if L1 + L2 <= L:
                    var v = connect_list(v1[], v2[])
                    vs.append(v)
    return vs

# 1234, 4 -> 2341
fn rotate(n: Int, e: Int) -> Int:
    var d = n // 10**(e-1)
    return (n - d * 10**(e-1)) * 10 + d

fn fill_digits357(first: Int, e: Int) -> List[Int]:
    if e == 1:
        var ns = List[Int]()
        for n in range(first, 10, 2):
            if n != 5:
                ns.append(n)
        return ns
    elif e % 2 == 1:
        var ns1 = fill_digits357(first, 1)
        var ns2 = fill_digits357(first, e-1)
        var ns = List[Int]()
        for n1 in ns1:
            for n2 in ns2:
                ns.append(n1[]*10**(e-1)+n2[])
        return ns
    else:
        var ns1 = fill_digits357(first, e//2)
        var ns = List[Int]()
        for n1 in ns1:
            for n2 in ns1:
                ns.append(n1[]*10**(e//2)+n2[])
        return ns

# rotateしたとき確実に最小になるかどうか
fn div_type(ds: List[Int]) -> Int:
    var head = ds[0]
    var type = 2
    for i in range(2, ds.size, 2):
        if ds[i] > head:
            return 0    # 最小にならない
        elif ds[i] == head:
            type = 1
    return type

fn fill_digits_core(ds: List[Int], first: Int) -> List[Int]:
    var ns = List[Int]()
    if ds.size == 2:
        var head = first * (10**ds[0] - 1) // 9 * 10**ds[1]
        var first1 = 7 if first == 3 else first + 2
        for n in fill_digits357(first1, ds[1]):
            ns.append(head + n[])
    else:
        var mid = 2 if ds.size % 4 == 2 else ds.size // 2
        var ds1 = sub_list(ds, 0, mid)
        var ds2 = sub_list(ds, mid, ds.size)
        var ns1 = fill_digits_core(ds1, first)
        var ns2 = fill_digits_core(ds2, first)
        for n1 in ns1:
            var e = sum_list(ds2)
            var head = n1[] * 10**e
            for n2 in ns2:
                ns.append(head + n2[])
    return ns

# [2, 1] -> [113, 117, 119, 337, 339, 779]
fn fill_digits(ds: List[Int]) -> List[Int]:
    var ns = List[Int]()
    for first in range(1, 10, 2):
        if first == 5:
            continue
        var ns1 = fill_digits_core(ds, first)
        ns.extend(ns1)
    return ns

fn is_smallest(n: Int, e: Int) -> Bool:
    var n1 = n
    for _ in range(e-1):
        n1 = rotate(n1, e)
        if n1 < n:
            return False
    return True

fn is_circular_prime(n: Int, e: Int) -> Bool:
    if n % 3 == 0:
        return False
    
    var n1 = n
    var ns = List[Int]()
    var bs = List[Bool]()
    for e1 in range(e):
        ns.append(n1)
        for p in range(7, 20, 2):
            if p * p > n1:
                bs.append(True)
                break
            elif n1 % p == 0:
                return False
        bs.append(False)
        if e1 != e - 1:
            n1 = rotate(n1, e)
    
    for i in range(e):
        if not bs[i] and not Miller_Rabin(ns[i]):
            return False
    
    return True

fn collect_rotated_numbers(n: Int, e: Int) -> Set[Int]:
    var s = Set[Int]()
    var n1 = n
    s.add(n)
    for _ in range(e-1):
        n1 = rotate(n1, e)
        s.add(n1)
    return s^

fn f(E: Int) -> Int:
    var N = 10**E
    var counter = 4     # 2, 3, 5, 7
    
    var repunit = 1
    for e in range(2, E+1):
        repunit = repunit * 10 + 1
        if not is_prime(e):
            continue
        if is_prime(repunit):
            print(repunit)
            counter += 1
    
    for d in range(1, E//2+1):
        for ds in divide(E, d):
            var e = sum_list(ds[])
            var type = div_type(ds[])
            if type == 0:
                continue
            var ns = fill_digits(ds[])
            for n in ns:
                if type == 1 and not is_smallest(n[], e):
                    continue
                elif is_circular_prime(n[], e):
                    var ns1 = collect_rotated_numbers(n[], e)
                    for n1 in ns1:
                        print(n1[])
                    counter += len(ns1)
    return counter

fn main() raises:
    var args = sys.argv()
    var E = atol(args[1])
    print(f(E))