Problem 718: Unreachable Numbers

Problem 718 source code

C++

#include <cassert>
#include <cstdint>
#include <functional>
#include <iostream>
#include <limits>
#include <queue>
#include <utility>
#include <vector>

namespace {

using i128 = __int128_t;
using u64 = std::uint64_t;

constexpr u64 kMod = 1'000'000'007ULL;

u64 pow_u64(u64 base, int exp) {
    u64 result = 1;
    while (exp > 0) {
        if (exp & 1) {
            result *= base;
        }
        base *= base;
        exp >>= 1;
    }
    return result;
}

u64 mod_pow(u64 base, u64 exp) {
    u64 result = 1;
    base %= kMod;
    while (exp > 0) {
        if (exp & 1ULL) {
            result = static_cast<u64>((static_cast<i128>(result) * base) % kMod);
        }
        base = static_cast<u64>((static_cast<i128>(base) * base) % kMod);
        exp >>= 1ULL;
    }
    return result;
}

u64 G_mod(const int p) {
    const u64 A = pow_u64(17, p);
    const u64 B = pow_u64(19, p);
    const u64 C = pow_u64(23, p);

    const u64 inf = std::numeric_limits<u64>::max() / 4;
    std::vector<u64> dist(static_cast<std::size_t>(A), inf);

    using Node = std::pair<u64, int>;
    std::priority_queue<Node, std::vector<Node>, std::greater<Node>> pq;
    dist[0] = 0;
    pq.push({0, 0});

    while (!pq.empty()) {
        const auto [d, u] = pq.top();
        pq.pop();
        if (d != dist[static_cast<std::size_t>(u)]) {
            continue;
        }

        int v = static_cast<int>((static_cast<u64>(u) + B) % A);
        u64 nd = d + B;
        if (nd < dist[static_cast<std::size_t>(v)]) {
            dist[static_cast<std::size_t>(v)] = nd;
            pq.push({nd, v});
        }

        v = static_cast<int>((static_cast<u64>(u) + C) % A);
        nd = d + C;
        if (nd < dist[static_cast<std::size_t>(v)]) {
            dist[static_cast<std::size_t>(v)] = nd;
            pq.push({nd, v});
        }
    }

    u64 sum_w_mod = 0;
    u64 sum_w2_mod = 0;
    for (const u64 w : dist) {
        const u64 wm = w % kMod;
        sum_w_mod += wm;
        if (sum_w_mod >= kMod) {
            sum_w_mod -= kMod;
        }
        sum_w2_mod = (sum_w2_mod + static_cast<u64>((static_cast<i128>(wm) * wm) % kMod)) % kMod;
    }

    const u64 inv2 = (kMod + 1) / 2;
    const u64 inv12 = mod_pow(12, kMod - 2);
    const u64 Am = A % kMod;
    const u64 invA = mod_pow(Am, kMod - 2);

    u64 genus = static_cast<u64>((static_cast<i128>(sum_w_mod) * invA) % kMod);
    genus = (genus + kMod - static_cast<u64>((static_cast<i128>((Am + kMod - 1) % kMod) * inv2) % kMod)) % kMod;

    u64 gaps = static_cast<u64>((static_cast<i128>(sum_w2_mod) * inv2) % kMod);
    gaps = static_cast<u64>((static_cast<i128>(gaps) * invA) % kMod);
    gaps = (gaps + kMod - static_cast<u64>((static_cast<i128>(sum_w_mod) * inv2) % kMod)) % kMod;
    const u64 A2m = static_cast<u64>((static_cast<i128>(Am) * Am) % kMod);
    gaps = (gaps + static_cast<u64>((static_cast<i128>((A2m + kMod - 1) % kMod) * inv12) % kMod)) % kMod;

    const u64 shift = (A % kMod + B % kMod + C % kMod) % kMod;
    const u64 tri = static_cast<u64>((static_cast<i128>(shift) * ((shift + kMod - 1) % kMod) % kMod) * inv2 % kMod);

    u64 ans = tri;
    ans = (ans + static_cast<u64>((static_cast<i128>(genus) * shift) % kMod)) % kMod;
    ans = (ans + gaps) % kMod;
    return ans;
}

}  // namespace

int main() {
    assert(G_mod(1) == 8'253);
    assert(G_mod(2) == 60'258'000);

    std::cout << G_mod(6) << '\n';
    return 0;
}

Python

import heapq

def solve():
    MOD = 1000000007
    p = 6
    A = 17**p; B = 19**p; C = 23**p

    def mod_pow(base, exp):
        r = 1; base %= MOD
        while exp > 0:
            if exp & 1: r = r * base % MOD
            base = base * base % MOD; exp >>= 1
        return r

    INF = float('inf')
    dist = [INF] * A; dist[0] = 0
    pq = [(0, 0)]
    while pq:
        d, u = heapq.heappop(pq)
        if d != dist[u]: continue
        for step in [B, C]:
            v = (u + step) % A; nd = d + step
            if nd < dist[v]:
                dist[v] = nd; heapq.heappush(pq, (nd, v))

    sw = sw2 = 0
    for w in dist:
        wm = w % MOD
        sw = (sw + wm) % MOD
        sw2 = (sw2 + wm * wm) % MOD

    inv2 = (MOD + 1) // 2
    inv12 = mod_pow(12, MOD - 2)
    Am = A % MOD; invA = mod_pow(Am, MOD - 2)

    genus = sw * invA % MOD
    genus = (genus - (Am - 1) * inv2) % MOD

    gaps = sw2 * inv2 % MOD * invA % MOD
    gaps = (gaps - sw * inv2) % MOD
    A2m = Am * Am % MOD
    gaps = (gaps + (A2m - 1) * inv12) % MOD

    shift = (A % MOD + B % MOD + C % MOD) % MOD
    tri = shift * ((shift - 1) % MOD) % MOD * inv2 % MOD

    ans = (tri + genus * shift + gaps) % MOD
    return str(ans)

if __name__ == '__main__':
    print(solve())

Java

import java.util.Arrays;
import java.util.PriorityQueue;

public class Euler718 {
    static final long kMod = 1000000007L;

    static long powU64(long base, int exp) {
        long result = 1;
        while (exp > 0) {
            if ((exp & 1) != 0) {
                result *= base;
            }
            base *= base;
            exp >>= 1;
        }
        return result;
    }

    static long modPow(long base, long exp) {
        long result = 1;
        base %= kMod;
        while (exp > 0) {
            if ((exp & 1) != 0) {
                result = (result * base) % kMod;
            }
            base = (base * base) % kMod;
            exp >>= 1;
        }
        return result;
    }

    static class Node implements Comparable<Node> {
        long d;
        int u;

        Node(long d, int u) {
            this.d = d;
            this.u = u;
        }

        @Override
        public int compareTo(Node other) {
            return Long.compare(this.d, other.d);
        }
    }

    public static String solve() {
        int p = 6;
        long A = powU64(17, p);
        long B = powU64(19, p);
        long C = powU64(23, p);

        long[] dist = new long[(int) A];
        Arrays.fill(dist, Long.MAX_VALUE / 4);

        PriorityQueue<Node> pq = new PriorityQueue<>();
        dist[0] = 0;
        pq.add(new Node(0, 0));

        while (!pq.isEmpty()) {
            Node node = pq.poll();
            long d = node.d;
            int u = node.u;

            if (d != dist[u]) {
                continue;
            }

            int v = (int) ((u + B) % A);
            long nd = d + B;
            if (nd < dist[v]) {
                dist[v] = nd;
                pq.add(new Node(nd, v));
            }

            v = (int) ((u + C) % A);
            nd = d + C;
            if (nd < dist[v]) {
                dist[v] = nd;
                pq.add(new Node(nd, v));
            }
        }

        long sumWMod = 0;
        long sumW2Mod = 0;
        for (long w : dist) {
            long wm = w % kMod;
            sumWMod = (sumWMod + wm) % kMod;
            sumW2Mod = (sumW2Mod + (wm * wm) % kMod) % kMod;
        }

        long inv2 = (kMod + 1) / 2;
        long inv12 = modPow(12, kMod - 2);
        long Am = A % kMod;
        long invA = modPow(Am, kMod - 2);

        long genus = (sumWMod * invA) % kMod;
        genus = (genus + kMod - (((Am + kMod - 1) % kMod * inv2) % kMod)) % kMod;

        long gaps = (sumW2Mod * inv2) % kMod;
        gaps = (gaps * invA) % kMod;
        gaps = (gaps + kMod - ((sumWMod * inv2) % kMod)) % kMod;
        long A2m = (Am * Am) % kMod;
        gaps = (gaps + (((A2m + kMod - 1) % kMod * inv12) % kMod)) % kMod;

        long shift = (A % kMod + B % kMod + C % kMod) % kMod;
        long tri = (shift * ((shift + kMod - 1) % kMod) % kMod * inv2) % kMod;

        long ans = tri;
        ans = (ans + (genus * shift) % kMod) % kMod;
        ans = (ans + gaps) % kMod;

        return Long.toString(ans);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}