euler/967.cpp

/* Generated by GPT-5 on 2025-11-05. Tested on Ubuntu 25.04. Runtime on Lenovo ThinkPad P14s: 0m0.559s */

#include <bits/stdc++.h>
using i128 = __int128_t;
using u128 = __uint128_t;
using ll = long long;

struct Node {
  unsigned long long prod;
  int t;   // (-1)^{|U|}
  long long s0, s1, s2; // s-vector
};

void extend_vec_mod3(int cls, long long &s0, long long &s1, long long &s2){
  long long a=s0,b=s1,c=s2;
  if(cls==1){ s0 = a - c; s1 = b - a; s2 = c - b; }
  else      { s0 = a - b; s1 = b - c; s2 = c - a; }
}

void gen_half(const std::vector<int>& primes, unsigned long long N, std::vector<Node>& out){
  struct Frame { int i; unsigned long long prod; int t; long long s0,s1,s2; };
  out.clear(); out.reserve(40000);
  std::vector<Frame> st;
  st.push_back({0,1, +1, 1,0,0});
  while(!st.empty()){
    auto [i,prod,t,s0,s1,s2] = st.back(); st.pop_back();
    out.push_back({prod,t,s0,s1,s2});
    for(int j=i;j<(int)primes.size();++j){
      unsigned long long p=primes[j];
      __uint128_t np = (u128)prod * p;
      if(np > N) break;
      long long ns0=s0,ns1=s1,ns2=s2;
      extend_vec_mod3(primes[j]%3, ns0, ns1, ns2);
      st.push_back({j+1,(unsigned long long)np,-t,ns0,ns1,ns2});
    }
  }
}

int main(){
  auto sieve = [](int B){
    std::vector<int> ps; std::vector<char> s(B+1,1); s[0]=s[1]=0;
    for(int i=2;i*i<=B;++i) if(s[i]) for(long long j=1LL*i*i;j<=B;j+=i) s[(int)j]=0;
    for(int i=2;i<=B;++i) if(s[i] && i!=3) ps.push_back(i);
    return ps;
  };

  auto solve = [&](unsigned long long N, int B){
    auto allp = sieve(B);
    // split
    int mid = allp.size()/2;
    std::vector<int> L(allp.begin(), allp.begin()+mid);
    std::vector<int> R(allp.begin()+mid, allp.end());

    std::vector<Node> A,Bv;
    gen_half(L,N,A);
    gen_half(R,N,Bv);

    std::sort(Bv.begin(), Bv.end(), [](const Node& x, const Node& y){ return x.prod < y.prod; });

    int m = (int)Bv.size();
    std::vector<unsigned long long> bp(m);
    std::vector<__int128_t> pref0(m+1,0), pref1(m+1,0), pref2(m+1,0);
    for(int i=0;i<m;++i){
      bp[i]=Bv[i].prod;
      pref0[i+1]=pref0[i] + (__int128_t)Bv[i].t * Bv[i].s0;
      pref1[i+1]=pref1[i] + (__int128_t)Bv[i].t * Bv[i].s1;
      pref2[i+1]=pref2[i] + (__int128_t)Bv[i].t * Bv[i].s2;
    }

    auto sum_range = [&](int l, int r){ // [l,r)
      return std::array<__int128_t,3>{ pref0[r]-pref0[l], pref1[r]-pref1[l], pref2[r]-pref2[l] };
    };

    __int128_t ans = 0;
    for(const auto& a : A){
      unsigned long long limit = N / a.prod;
      int tL = a.t;
      long long x0=a.s0, x1=a.s1, x2=a.s2;
      // process right items with b <= limit in blocks where floor(N/(a.prod*b)) is constant
      int r = std::upper_bound(bp.begin(), bp.end(), limit) - bp.begin();
      int i = 0;
      while(i < r){
        unsigned long long v = (unsigned long long)( (unsigned long long)(N / a.prod) / bp[i] );
        if(v==0){ break; }
        unsigned long long maxb = (unsigned long long)((N / a.prod) / v);
        int j = std::upper_bound(bp.begin()+i, bp.begin()+r, maxb) - bp.begin();
        auto sums = sum_range(i,j);
        __int128_t comb = (__int128_t)x0*sums[0] + (__int128_t)x1*sums[2] + (__int128_t)x2*sums[1];
        ans += (__int128_t)tL * (__int128_t)v * comb;
        i = j;
      }
    }
    // add the empty/empty pair is included (A contains U=∅ and B contains ∅), good.
    unsigned long long out = (unsigned long long)ans; // result is nonnegative and fits 64-bit
    return out;
  };

  // Examples:
  std::cout << solve(10ULL, 4) << "\n";   // 5
  std::cout << solve(10ULL, 10) << "\n";  // 3
  std::cout << solve(100ULL, 10) << "\n"; // 41

  // The asked value:
  std::cout << solve(1000000000000000000ULL, 120) << "\n";
  return 0;
}