以下の内容はhttps://inamori.hateblo.jp/entry/2024/10/09/210036より取得しました。


MojoでProject Euler 70

https://projecteuler.net/problem=70

 n/\phi(n)が小さいほうから調べればよいです。一番小さくなるのは素数ですが、素数はpermutationになりえないので、素数の平方から調べます。100未満とすると、7×7です。次に小さいのは5×5か7×13です。 25/\phi(25) = 5/4 = 1.25 91/\phi(91) = 91/72 = 1.26388\cdotsなので、次は5×5です。その次は7×13と3×3と3×31ですが、7×13が最も小さくなります。3×3の次は3×3×3もPriorityQueueに入ってきます。
1016で約18秒でした。

import sys


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

fn int_sqrt(n: Int) -> Int:
    var x = 1
    x = (x + n // x) // 2
    while True:
        var x1 = (x + n // x) // 2
        if x1 >= x:
            return x
        x = x1

fn digits(owned n: Int) -> List[Int]:
    var ds = List[Int]()
    while n > 0:
        var d = n % 10
        n //= 10
        ds.append(d)
    return ds


#################### List ####################

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

trait Printable(CollectionElement, Stringable):
    pass

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


#################### prime ####################

fn make_prime_table(N: Int) -> List[Int]:
    var a = initialize_list(N+1, True)
    for p in range(2, N+1):
        if p * p > N:
            break
        elif not a[p]:
            continue
        
        for n in range(p*p, N+1, p):
            a[n] = False
    
    var ps = List[Int]()
    for n in range(2, N+1):
        if a[n]:
            ps.append(n)
    return ps

fn extend_prime_table(first: Int, last: Int, inout ps: List[Int]):
    var a = initialize_list(last - first, True)
    for p in ps:
        if p[] * p[] >= last:
            break
        
        for n in range((first+p[]-1)//p[]*p[], last, p[]):
            a[n-first] = False
    
    for n in range(first, last):
        if a[n-first]:
            ps.append(n)

# N未満の最も大きい素数のインデックス
fn max_prime(N: Int) -> Int:
    var max_p = primes[len(primes)-1]
    if N >= max_p:
        var last = max_p * 3 // 2 + 5
        extend_prime_table(max_p + 1, last, primes)
    var first = 0
    var last = len(primes)
    while first < last - 1:
        var mid = (first + last) // 2
        if primes[mid] >= N:
            last = mid
        else:
            first = mid
    return first


#################### HeapQueue ####################

struct HeapQueue[T: ComparableCollectionElement]:
    var a: List[T]
    
    fn __init__(inout self):
        self.a = List[T]()
    
    fn is_empty(self) -> Bool:
        return len(self.a) == 0
    
    fn pop(inout self) -> T:
        var top = self.a[0]
        self.a[0] = self.a.pop()
        var N = self.a.size
        var i = 0
        while i < N-1:
            var j1 = i*2+1
            var j2 = i*2+2
            var j: Int = 0
            if j1 >= N:
                break
            elif j2 == N:
                j = j1
            else:
                if self.a[j1] <= self.a[j2]:
                    j = j1
                else:
                    j = j2
            if self.a[j] < self.a[i]:
                self.swap(i, j)
                i = j
            else:
                break
        
        return top
    
    fn push(inout self, e: T):
        self.a.append(e)
        var i = self.a.size - 1
        while i > 0:
            var j = (i-1) // 2
            if self.a[j] <= self.a[i]:
                break
            else:
                self.swap(i, j)
                i = j
    
    fn swap(inout self, i: Int, j: Int):
        var tmp = self.a[i]
        self.a[i] = self.a[j]
        self.a[j] = tmp^


#################### State ####################

@value
struct State(Comparable, Stringable):
    var a: List[Int]
    var n: Int
    var phi: Int
    var ratio: Float64
    
    fn __init__(inout self, a: List[Int]):
        self.a = a
        self.n = 1
        self.phi = 1
        self.ratio = 1.0
        for i in range(len(a)):
            var p = primes[a[i]]
            self.n *= p
            if i == 0 or a[i] != a[i-1]:
                self.phi *= p - 1
            else:
                self.phi *= p
        self.ratio = self.n / self.phi
    
    fn __lt__(self, other: State) -> Bool:
        return self.ratio < other.ratio
    
    fn __le__(self, other: State) -> Bool:
        return self.ratio <= other.ratio
    
    fn __eq__(self, other: State) -> Bool:
        return self.ratio == other.ratio
    
    fn __ne__(self, other: State) -> Bool:
        return self.ratio != other.ratio
    
    fn __gt__(self, other: State) -> Bool:
        return self.ratio > other.ratio
    
    fn __ge__(self, other: State) -> Bool:
        return self.ratio >= other.ratio
    
    fn __str__(self) -> String:
        if self.a.size > 0:
            var s = "[" + str(self.a[0])
            for i in range(1, self.a.size):
                s += ", " + str(self.a[i])
            s += "]"
            return s
        else:
            return "[]"
    
    fn is_all_same(self) -> Bool:
        for i in range(1, len(self.a)):
            if self.a[i] != self.a[0]:
                return False
        else:
            return True
    
    fn nexts(self, N: Int) -> List[State]:
        var states = List[State]()
        var L = len(self.a)
        if self.a[L-2] == self.a[L-1]:
            # 7 * 7 -> 7 * 13
            var p = primes[self.a[L-1]]
            var n = self.n // p
            var q = (N - 1) // n
            var i = max_prime(q + 1)
            if i > self.a[L-1]:
                var a = self.a
                a[L-1] = i
                states.append(State(a))
        elif self.a[L-2] < self.a[L-1] - 1:
            # 7 * 13 -> 7 * 11
            var a = self.a
            a[L-1] -= 1
            states.append(State(a))
        
        if self.n <= (N - 1) // primes[self.a[L-1]]:
            # 5 * 7 -> 5 * 7 * 7
            var a = self.a
            a.append(a[L-1])
            var s = State(a)
            states.append(State(a))
        
        if len(self.a) == 2 and self.a[0] == self.a[1] and self.a[0] > 0:
            # 5 * 5 -> 3 * 3
            var a = List[Int](self.a[0]-1, self.a[0]-1)
            states.append(State(a))
        return states
    
    fn is_permutation(self) -> Bool:
        var ds1 = digits(self.n)
        var ds2 = digits(self.phi)
        if len(ds1) != len(ds2):
            return False
        sort(ds1)
        sort(ds2)
        for i in range(len(ds1)):
            if ds1[i] != ds2[i]:
                return False
        else:
            return True


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

var primes: List[Int] = List[Int]()

fn f(N: Int) -> Int:
    var R = int_sqrt(N)
    primes = make_prime_table(N // R)
    var hq = HeapQueue[State]()
    var i = max_prime(R + 1)
    hq.push(State(List[Int](i, i)))
    while not hq.is_empty():
        var state = hq.pop()
        if state.is_permutation():
            return state.n
        var new_states = state.nexts(N)
        for new_state in new_states:
            hq.push(new_state[])
    return 0

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



以上の内容はhttps://inamori.hateblo.jp/entry/2024/10/09/210036より取得しました。
このページはhttp://font.textar.tv/のウェブフォントを使用してます

不具合報告/要望等はこちらへお願いします。
モバイルやる夫Viewer Ver0.14