Problem 699: Triffle Numbers

Problem 699 source code

C++

#include <algorithm>
#include <chrono>
#include <cstdint>
#include <functional>
#include <future>
#include <iostream>
#include <string>
#include <thread>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>

using std::uint64_t;
using u128 = __uint128_t;

namespace {

constexpr uint64_t kDefaultN = 100000000000000ULL;

uint64_t gcd_u64(uint64_t a, uint64_t b) {
    while (b) {
        uint64_t t = a % b;
        a = b;
        b = t;
    }
    return a;
}

uint64_t mul_mod_u64(uint64_t a, uint64_t b, uint64_t mod) {
    return static_cast<uint64_t>((static_cast<u128>(a) * b) % mod);
}

uint64_t pow_mod_u64(uint64_t a, uint64_t d, uint64_t mod) {
    uint64_t r = 1;
    while (d) {
        if (d & 1) r = mul_mod_u64(r, a, mod);
        a = mul_mod_u64(a, a, mod);
        d >>= 1;
    }
    return r;
}

bool is_prime_u64(uint64_t n) {
    if (n < 2) return false;
    static uint64_t small_primes[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
    for (uint64_t p : small_primes) {
        if (n == p) return true;
        if (n % p == 0) return false;
    }

    uint64_t d = n - 1, s = 0;
    while ((d & 1) == 0) {
        d >>= 1;
        ++s;
    }

    auto witness = [&](uint64_t a) -> bool {
        if (a % n == 0) return false;
        uint64_t x = pow_mod_u64(a, d, n);
        if (x == 1 || x == n - 1) return false;
        for (uint64_t i = 1; i < s; ++i) {
            x = mul_mod_u64(x, x, n);
            if (x == n - 1) return false;
        }
        return true;
    };

    static uint64_t bases[] = {2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL};
    for (uint64_t a : bases) {
        if (witness(a)) return false;
    }
    return true;
}

uint64_t splitmix64(uint64_t& x) {
    uint64_t z = (x += 0x9e3779b97f4a7c15ULL);
    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL;
    z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL;
    return z ^ (z >> 31);
}

thread_local uint64_t rng_state =
    static_cast<uint64_t>(std::chrono::high_resolution_clock::now().time_since_epoch().count()) ^
    0x9e3779b97f4a7c15ULL;

uint64_t rand_u64(uint64_t lo, uint64_t hi) {
    uint64_t r = splitmix64(rng_state);
    return lo + (hi > lo ? (r % (hi - lo + 1)) : 0);
}

uint64_t pollard_rho(uint64_t n) {
    if ((n & 1ULL) == 0) return 2;
    if (n % 3ULL == 0) return 3;

    while (true) {
        uint64_t c = rand_u64(1, n - 1);
        uint64_t x = rand_u64(0, n - 1);
        uint64_t y = x;
        uint64_t d = 1;

        auto f = [&](uint64_t v) { return (mul_mod_u64(v, v, n) + c) % n; };

        while (d == 1) {
            x = f(x);
            y = f(f(y));
            uint64_t diff = (x > y) ? (x - y) : (y - x);
            d = gcd_u64(diff, n);
        }
        if (d != n) return d;
    }
}

void factor_u64(uint64_t n, std::vector<uint64_t>& fac) {
    if (n == 1) return;
    if (is_prime_u64(n)) {
        fac.push_back(n);
        return;
    }
    uint64_t d = pollard_rho(n);
    factor_u64(d, fac);
    factor_u64(n / d, fac);
}

struct FactorCache {
    std::unordered_map<uint64_t, std::vector<uint64_t>> mp;

    std::vector<uint64_t> get(uint64_t n) {
        auto it = mp.find(n);
        if (it != mp.end()) return it->second;

        std::vector<uint64_t> f;
        if (n > 1) {
            std::vector<uint64_t> tmp;
            factor_u64(n, tmp);
            std::sort(tmp.begin(), tmp.end());
            tmp.erase(std::unique(tmp.begin(), tmp.end()), tmp.end());
            f = std::move(tmp);
        }
        mp.emplace(n, f);
        return f;
    }
};

int v3_u64(uint64_t x) {
    int c = 0;
    while (x % 3ULL == 0) {
        x /= 3ULL;
        ++c;
    }
    return c;
}

uint64_t sigma_prime_power(uint64_t p, int e, uint64_t p_pow) {
    u128 next_pow = static_cast<u128>(p_pow) * p;
    u128 sig = (next_pow - 1) / (p - 1);
    return static_cast<uint64_t>(sig);
}

std::string to_string_u128(u128 x) {
    if (x == 0) return "0";
    std::string s;
    while (x > 0) {
        int digit = static_cast<int>(x % 10);
        s.push_back(static_cast<char>('0' + digit));
        x /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

struct SolverA {
    uint64_t A;
    uint64_t limit;
    int maxv3;
    bool validate;

    uint64_t sum_m = 0;

    std::unordered_set<uint64_t> visited_m;
    FactorCache cache;
    std::vector<uint64_t> used;

    SolverA(uint64_t A_, uint64_t limit_, int maxv3_, bool validate_)
        : A(A_), limit(limit_), maxv3(maxv3_), validate(validate_) {
        visited_m.reserve(1 << 16);
        cache.mp.reserve(1 << 16);
        used.reserve(32);
    }

    bool is_used(uint64_t p) const {
        for (uint64_t x : used) {
            if (x == p) return true;
        }
        return false;
    }

    void dfs(uint64_t m, uint64_t num, uint64_t den, int v3sig, uint64_t sigma_m) {
        if (m > limit) return;
        if (!visited_m.insert(m).second) return;

        if (den == 1 && v3sig <= maxv3) {
            sum_m += m;

            if (validate) {
                uint64_t g = gcd_u64(m, sigma_m);
                uint64_t r = m / g;
                if (A % r != 0) {
                    std::cerr << "Validation failed: A%r!=0, A=" << A << " m=" << m << "\n";
                    std::abort();
                }
                if (v3_u64(sigma_m) != v3sig) {
                    std::cerr << "Validation failed: v3 mismatch, m=" << m << "\n";
                    std::abort();
                }
                u128 prod = static_cast<u128>(A) * sigma_m;
                if (static_cast<uint64_t>(prod % m) != 0ULL) {
                    std::cerr << "Validation failed: (A*sigma_m)%m!=0, m=" << m << "\n";
                    std::abort();
                }
            }
        }

        std::vector<uint64_t> cand = cache.get(num);
        cand.push_back(2);
        std::sort(cand.begin(), cand.end());
        cand.erase(std::unique(cand.begin(), cand.end()), cand.end());

        for (uint64_t p : cand) {
            if (p == 3) continue;
            if (is_used(p)) continue;

            uint64_t p_pow = p;
            for (int e = 1;; ++e) {
                if (static_cast<u128>(m) * p_pow > limit) break;

                uint64_t sig = sigma_prime_power(p, e, p_pow);
                int v3f = v3_u64(sig);
                if (v3sig + v3f <= maxv3) {
                    uint64_t denom_factor = p_pow;
                    uint64_t sig_factor = sig;

                    uint64_t num1 = num;
                    uint64_t den1 = den;

                    uint64_t g1 = gcd_u64(num1, denom_factor);
                    num1 /= g1;
                    denom_factor /= g1;

                    uint64_t g2 = gcd_u64(sig_factor, den1);
                    sig_factor /= g2;
                    den1 /= g2;

                    u128 new_num = static_cast<u128>(num1) * sig_factor;
                    u128 new_den = static_cast<u128>(den1) * denom_factor;

                    uint64_t new_m = static_cast<uint64_t>(static_cast<u128>(m) * p_pow);
                    uint64_t new_sigma_m = static_cast<uint64_t>(static_cast<u128>(sigma_m) * sig);

                    used.push_back(p);
                    dfs(new_m, static_cast<uint64_t>(new_num), static_cast<uint64_t>(new_den),
                        v3sig + v3f, new_sigma_m);
                    used.pop_back();
                }

                if (p_pow > limit / p) break;
                p_pow *= p;
            }
        }
    }

    void run() {
        dfs(1, A, 1, 0, 1);
    }
};

u128 computeT(uint64_t N, bool validate, bool multithread) {
    std::vector<std::pair<int, uint64_t>> tasks;
    u128 pow3 = 1;
    for (int a = 1;; ++a) {
        pow3 *= 3;
        if (pow3 > N) break;
        tasks.push_back({a, static_cast<uint64_t>(pow3)});
    }

    auto solve_one = [&](int a, uint64_t p3) -> u128 {
        uint64_t limit = N / p3;
        u128 t = static_cast<u128>(p3) * 3 - 1;
        uint64_t A = static_cast<uint64_t>(t / 2);
        SolverA solver(A, limit, a - 1, validate);
        solver.run();
        return static_cast<u128>(p3) * solver.sum_m;
    };

    u128 total = 0;
    if (!multithread || tasks.size() <= 1) {
        for (const auto& task : tasks) total += solve_one(task.first, task.second);
        return total;
    }

    std::vector<std::future<u128>> fut;
    fut.reserve(tasks.size());
    for (const auto& task : tasks) {
        fut.push_back(std::async(std::launch::async, solve_one, task.first, task.second));
    }
    for (auto& f : fut) total += f.get();
    return total;
}

}  // namespace

int main(int argc, char** argv) {
    uint64_t N = kDefaultN;
    bool validate = true;
    bool multithread = true;

    for (int i = 1; i < argc; ++i) {
        std::string s = argv[i];
        if (s.rfind("--N=", 0) == 0) {
            N = std::stoull(s.substr(4));
        } else if (s.rfind("--validate=", 0) == 0) {
            validate = (std::stoi(s.substr(11)) != 0);
        } else if (s.rfind("--mt=", 0) == 0) {
            multithread = (std::stoi(s.substr(5)) != 0);
        }
    }

    if (validate) {
        u128 t100 = computeT(100ULL, false, false);
        if (static_cast<uint64_t>(t100) != 270ULL) {
            std::cerr << "Checkpoint failed: T(100) expected 270, got "
                      << static_cast<uint64_t>(t100) << "\n";
            return 1;
        }
        u128 t1e6 = computeT(1000000ULL, false, false);
        if (static_cast<uint64_t>(t1e6) != 26089287ULL) {
            std::cerr << "Checkpoint failed: T(1e6) expected 26089287, got "
                      << static_cast<uint64_t>(t1e6) << "\n";
            return 1;
        }
    }

    u128 ans = computeT(N, validate, multithread);
    std::cout << to_string_u128(ans) << "\n";
    return 0;
}

Python

import sys
import multiprocessing
import multiprocessing.pool
import math

def gcd(a, b):
    while b:
        a, b = b, a % b
    return a

def is_prime(n):
    if n < 2: return False
    small_primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
    for p in small_primes:
        if n == p: return True
        if n % p == 0: return False
        
    d = n - 1
    s = 0
    while (d & 1) == 0:
        d >>= 1
        s += 1
        
    def witness(a):
        if a % n == 0: return False
        x = pow(a, d, n)
        if x == 1 or x == n - 1: return False
        for _ in range(1, s):
            x = (x * x) % n
            if x == n - 1: return False
        return True
        
    bases = [2, 325, 9375, 28178, 450775, 9780504, 1795265022]
    for a in bases:
        if witness(a): return False
    return True

import random

def pollard_rho(n):
    if n & 1 == 0: return 2
    if n % 3 == 0: return 3
    
    while True:
        c = random.randint(1, n - 1)
        x = random.randint(0, n - 1)
        y = x
        d = 1
        
        def f(v):
            return (v * v + c) % n
            
        while d == 1:
            x = f(x)
            y = f(f(y))
            d = gcd(abs(x - y), n)
            
        if d != n: return d

def factor(n, fac):
    if n == 1: return
    if is_prime(n):
        fac.append(n)
        return
    d = pollard_rho(n)
    factor(d, fac)
    factor(n // d, fac)

class FactorCache:
    def __init__(self):
        self.mp = {}
        
    def get(self, n):
        if n in self.mp: return self.mp[n]
        f = []
        if n > 1:
            tmp = []
            factor(n, tmp)
            f = sorted(list(set(tmp)))
        self.mp[n] = f
        return f

def v3(x):
    c = 0
    while x % 3 == 0:
        x //= 3
        c += 1
    return c

def sigma_prime_power(p, e, p_pow):
    next_pow = p_pow * p
    return (next_pow - 1) // (p - 1)

class SolverA:
    def __init__(self, A, limit, maxv3, validate=False):
        self.A = A
        self.limit = limit
        self.maxv3 = maxv3
        self.validate = validate
        self.sum_m = 0
        self.visited_m = set()
        self.cache = FactorCache()
        self.used = []
        
    def dfs(self, m, num, den, v3sig, sigma_m):
        if m > self.limit: return
        if m in self.visited_m: return
        self.visited_m.add(m)
        
        if den == 1 and v3sig <= self.maxv3:
            self.sum_m += m
            
        cand = self.cache.get(num)[:]
        cand.append(2)
        cand = sorted(list(set(cand)))
        
        for p in cand:
            if p == 3: continue
            if p in self.used: continue
            
            p_pow = p
            e = 1
            while True:
                if m * p_pow > self.limit: break
                
                sig = sigma_prime_power(p, e, p_pow)
                v3f = v3(sig)
                if v3sig + v3f <= self.maxv3:
                    denom_factor = p_pow
                    sig_factor = sig
                    
                    num1 = num
                    den1 = den
                    
                    g1 = gcd(num1, denom_factor)
                    num1 //= g1
                    denom_factor //= g1
                    
                    g2 = gcd(sig_factor, den1)
                    sig_factor //= g2
                    den1 //= g2
                    
                    new_num = num1 * sig_factor
                    new_den = den1 * denom_factor
                    
                    new_m = m * p_pow
                    new_sigma_m = sigma_m * sig
                    
                    self.used.append(p)
                    self.dfs(new_m, new_num, new_den, v3sig + v3f, new_sigma_m)
                    self.used.pop()
                    
                if p_pow > self.limit // p: break
                p_pow *= p
                e += 1

def solve_one(args):
    a, p3, N = args
    limit = N // p3
    t = p3 * 3 - 1
    A = t // 2
    solver = SolverA(A, limit, a - 1, False)
    solver.dfs(1, A, 1, 0, 1)
    return p3 * solver.sum_m

def compute_T(N):
    tasks = []
    pow3 = 1
    a = 1
    while True:
        pow3 *= 3
        if pow3 > N: break
        tasks.append((a, pow3, N))
        a += 1
        
    threads = multiprocessing.cpu_count() or 1
    total = 0
    if threads <= 1 or len(tasks) <= 1:
        for t in tasks:
            total += solve_one(t)
    else:
        with multiprocessing.Pool(threads) as pool:
            results = pool.map(solve_one, tasks)
        total = sum(results)
    return total

def solve():
    N = 100000000000000
    ans = compute_T(N)
    return str(ans)

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

Java

import java.math.BigInteger;
import java.util.*;
import java.util.concurrent.*;

public class Euler699 {
    static long gcd(long a, long b) {
        while (b != 0) {
            long t = a % b;
            a = b;
            b = t;
        }
        return a;
    }

    static long mulMod(long a, long b, long mod) {
        return BigInteger.valueOf(a).multiply(BigInteger.valueOf(b)).mod(BigInteger.valueOf(mod)).longValue();
    }

    static long powMod(long a, long d, long mod) {
        long r = 1;
        while (d > 0) {
            if ((d & 1) != 0)
                r = mulMod(r, a, mod);
            a = mulMod(a, a, mod);
            d >>= 1;
        }
        return r;
    }

    static boolean isPrime(long n) {
        if (n < 2)
            return false;
        long[] smallPrimes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 };
        for (long p : smallPrimes) {
            if (n == p)
                return true;
            if (n % p == 0)
                return false;
        }

        long d = n - 1;
        int s = 0;
        while ((d & 1) == 0) {
            d >>= 1;
            s++;
        }

        long[] bases = { 2L, 325L, 9375L, 28178L, 450775L, 9780504L, 1795265022L };
        for (long a : bases) {
            if (a % n == 0)
                continue;
            long x = powMod(a, d, n);
            if (x == 1 || x == n - 1)
                continue;
            boolean composite = true;
            for (int i = 1; i < s; i++) {
                x = mulMod(x, x, n);
                if (x == n - 1) {
                    composite = false;
                    break;
                }
            }
            if (composite)
                return false;
        }
        return true;
    }

    static long pollardRho(long n) {
        if ((n & 1) == 0)
            return 2;
        if (n % 3 == 0)
            return 3;

        ThreadLocalRandom rand = ThreadLocalRandom.current();
        while (true) {
            long c = rand.nextLong(1, n);
            long x = rand.nextLong(0, n);
            long y = x;
            long d = 1;

            while (d == 1) {
                x = (mulMod(x, x, n) + c) % n;
                long ty = (mulMod(y, y, n) + c) % n;
                y = (mulMod(ty, ty, n) + c) % n;

                long diff = Math.abs(x - y);
                d = gcd(diff, n);
            }
            if (d != n)
                return d;
        }
    }

    static void factor(long n, List<Long> fac) {
        if (n == 1)
            return;
        if (isPrime(n)) {
            fac.add(n);
            return;
        }
        long d = pollardRho(n);
        factor(d, fac);
        factor(n / d, fac);
    }

    static class FactorCache {
        Map<Long, List<Long>> mp = new HashMap<>();

        List<Long> get(long n) {
            if (mp.containsKey(n))
                return mp.get(n);
            List<Long> f = new ArrayList<>();
            if (n > 1) {
                List<Long> tmp = new ArrayList<>();
                factor(n, tmp);
                Collections.sort(tmp);
                for (long p : tmp) {
                    if (f.isEmpty() || f.get(f.size() - 1) != p) {
                        f.add(p);
                    }
                }
            }
            mp.put(n, f);
            return f;
        }
    }

    static int v3(long x) {
        int c = 0;
        while (x % 3 == 0) {
            x /= 3;
            c++;
        }
        return c;
    }

    static long sigmaPrimePower(long p, int e, long pPow) {
        BigInteger nextPow = BigInteger.valueOf(pPow).multiply(BigInteger.valueOf(p));
        return nextPow.subtract(BigInteger.ONE).divide(BigInteger.valueOf(p - 1)).longValue();
    }

    static class SolverA {
        long A;
        long limit;
        int maxv3;
        long sumM = 0;
        HashSet<Long> visitedM = new HashSet<>();
        FactorCache cache = new FactorCache();
        List<Long> used = new ArrayList<>();

        SolverA(long a, long limit, int maxv3) {
            this.A = a;
            this.limit = limit;
            this.maxv3 = maxv3;
        }

        void dfs(long m, long num, long den, int v3sig, long sigmaM) {
            if (m > limit)
                return;
            if (!visitedM.add(m))
                return;

            if (den == 1 && v3sig <= maxv3) {
                sumM += m;
            }

            List<Long> cand = new ArrayList<>(cache.get(num));
            cand.add(2L);
            Collections.sort(cand);
            List<Long> uniqCand = new ArrayList<>();
            for (long p : cand) {
                if (uniqCand.isEmpty() || uniqCand.get(uniqCand.size() - 1) != p) {
                    uniqCand.add(p);
                }
            }

            for (long p : uniqCand) {
                if (p == 3)
                    continue;
                if (used.contains(p))
                    continue;

                long pPow = p;
                for (int e = 1;; e++) {
                    if (BigInteger.valueOf(m).multiply(BigInteger.valueOf(pPow))
                            .compareTo(BigInteger.valueOf(limit)) > 0)
                        break;

                    long sig = sigmaPrimePower(p, e, pPow);
                    int v3f = v3(sig);

                    if (v3sig + v3f <= maxv3) {
                        long denomFactor = pPow;
                        long sigFactor = sig;

                        long num1 = num;
                        long den1 = den;

                        long g1 = gcd(num1, denomFactor);
                        num1 /= g1;
                        denomFactor /= g1;

                        long g2 = gcd(sigFactor, den1);
                        sigFactor /= g2;
                        den1 /= g2;

                        long newNum = num1 * sigFactor;
                        long newDen = den1 * denomFactor;

                        long newM = m * pPow;
                        long newSigmaM = sigmaM * sig;

                        used.add(p);
                        dfs(newM, newNum, newDen, v3sig + v3f, newSigmaM);
                        used.remove(used.size() - 1);
                    }

                    if (pPow > limit / p)
                        break;
                    pPow *= p;
                }
            }
        }
    }

    public static String solve() {
        long N = 100000000000000L;
        List<long[]> tasks = new ArrayList<>();
        long pow3 = 1;
        int a = 1;
        while (true) {
            pow3 *= 3;
            if (pow3 > N)
                break;
            tasks.add(new long[] { a, pow3 });
            a++;
        }

        int threads = Runtime.getRuntime().availableProcessors();
        if (threads < 1)
            threads = 1;
        ExecutorService pool = Executors.newFixedThreadPool(threads);
        List<Future<Long>> futures = new ArrayList<>();

        for (long[] task : tasks) {
            futures.add(pool.submit(() -> {
                int t_a = (int) task[0];
                long t_p3 = task[1];
                long limit = N / t_p3;
                long t_val = t_p3 * 3 - 1;
                long A = t_val / 2;
                SolverA solver = new SolverA(A, limit, t_a - 1);
                solver.dfs(1L, A, 1L, 0, 1L);
                return t_p3 * solver.sumM;
            }));
        }

        long total = 0;
        for (Future<Long> f : futures) {
            try {
                total += f.get();
            } catch (Exception e) {
            }
        }
        pool.shutdown();
        return Long.toString(total);
    }

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