#553667 (C++14) No.1276 3枚のカード

• 問題: No.1276 3枚のカード
• ユーザー: logx
• 提出日時: 2020-09-16 19:28:21
• 言語: C++14 (gcc 13.3.0 + boost 1.87.0)
• 結果: AC
• 実行時間: 297 ms / 2,000 ms
• コード長: 9,561 bytes
• コンパイル時間: 2,046 ms
• コンパイル使用メモリ: 191,908 KB
• 実行使用メモリ: 6,948 KB
• 最終ジャッジ日時: 2024-06-22 03:28:31
• 合計ジャッジ時間: 12,716 ms
• ジャッジサーバーID （参考情報）: judge3 / judge2

ファイルパターン	結果
other	AC * 61

ソースコード

#include <bits/stdc++.h>
#ifndef ONLINE_JUDGE
#define _GLIBCXX_DEBUG
#endif
using namespace std;
#define rep(i, n) for (int i = 0; i < (int)(n); i++)
#define rep2(i, s, n) for (int i = s; i < (int)(n); i++)
#define Clear(a) a = decltype(a)()
#define all(a) a.begin(),a.end()
#define rall(a) a.rbegin(),a.rend()
#define vec vector
typedef long long ll;
typedef pair<ll,ll> P;
//const ll big=998244353;
const ll big=1000000007LL;
const ll INF=1e18;
const int dx[4]={1,0,-1,0};
const int dy[4]={0,1,0,-1};
ll max(ll x,ll y){
    if(x>y)return x;
    else return y;
}
ll min(ll x,ll y){
    if(x<y)return x;
    else return y;
}
ll expm(ll x,ll y){
    if(y==0)return 1;//0^0=1
    if(x==1||x==0)return x;
    if(y%2==1)return (expm(x,y-1)*x)%big;
    ll t=expm(x,y/2);
    return (t*t)%big;
}
ll exp(ll x,ll y){
    if(y==0)return 1;//0^0=1
    if(x==1||y==0)return x;
    if(y%2==1)return exp(x,y-1)*x;
    ll t=exp(x,y/2);
    return t*t;
}

//約数和のためのやつ
#define u64 unsigned long long
#define u32 unsigned
 
struct u128 {
  u64 l;
  u64 h;
  u128() : l(0), h(0) {}
  u128(u64 l) : l(l), h(0) {}
  u128(u64 l, u64 h) : l(l), h(h) {}
 
  static u128 mul64(u64 a, u64 b) {
    u32 a_hi = a >> 32, a_lo = u32(a);
    u32 b_hi = b >> 32, b_lo = u32(b);
    u64 p0 = (u64)(a_lo) * b_lo;
    u64 p10 = (u64)(a_hi) * b_lo + (p0 >> 32);
    u64 p11 = (u64)(a_lo) * b_hi + p10;
    u64 p2 = ((u64)(p11 < p10) << 32) + (u64)(a_hi) * b_hi + (p11 >> 32);
    return u128(u32(p0) | (p11 << 32), p2);
  }
 
  static void divmod(u128 a, u32 d, u128& q, u32& r) {
    u32 a3 = a.h >> 32, a2 = a.h;
    u32 a1 = a.l >> 32, a0 = a.l;
 
    u64 t = a3;
    u32 q3 = t / d; t = ((t % d) << 32) | a2; 
    u32 q2 = t / d; t = ((t % d) << 32) | a1;
    u32 q1 = t / d; t = ((t % d) << 32) | a0;
    u32 q0 = t / d; r = t % d;
    q = u128(q0 | ((u64)(q1) << 32), q2 | ((u64)(q3) << 32));
  }
  bool is_zero() const {
    return l == 0 && h == 0;
  }
  bool operator >= (const u128& rhs) const {
    return h > rhs.h || (h == rhs.h && l >= rhs.l);
  }
  u128& operator += (const u64 rhs) {
    u64 old_l = l;
    l += rhs;
    h += (l < old_l);
    return *this;
  }
  u128& operator += (const u128& rhs) {
    u64 old_l = l;
    l += rhs.l;
    h += (l < old_l) + rhs.h;
    return *this;
  }
  u128& operator -= (const u128& rhs) {
    u64 old_l = l;
    l -= rhs.l;
    h -= (l > old_l) + rhs.h;
    return *this;
  }
  u128 operator + (const u64 rhs) const {
    return u128(*this) += rhs;
  }
  u128 operator + (const u128& rhs) const {
    return u128(*this) += rhs;
  }
  u128 operator - (const u128& rhs) const {
    return u128(*this) -= rhs;
  }
  std::string to_string() const {
    static const u32 ten9 = u32(1e9);
    if (is_zero()) {
      return "0";
    }
    char str[41];
    int i = 0;
 
    u128 n = *this;
    u32 r;
    while (1) {
      divmod(n, ten9, n, r);
      if (n.is_zero()) {
        while (r > 0) str[i++] = r % 10, r /= 10;
        break;
      } else {
        for (int j = 0; j < 9; ++j) str[i++] = r % 10, r /= 10;
      }
    }
    std::string ret;
    while (i) ret += '0' + str[--i];
    return ret;
  }
};
 
u32 isqrt(u64 n) {
  if (n >= (u64)(UINT32_MAX) * UINT32_MAX) {
    return UINT32_MAX;
  } else {
    u32 ret = sqrtl(n);
    assert((u64)(ret) * ret <= n);
    assert((u64)(ret + 1) * (ret + 1) > n);
    return ret;
  }
}
 
inline u64 div64_32(u64 a, u32 b) {
  u32 r, ql, qh;
  u32 al = a, ah = a >> 32;
  qh = ah / b; ah %= b;
  __asm__ (
    "divl %4\n\t"
    : "=a"(ql), "=d"(r)
    : "0"(al), "1"(ah), "rm"(b)
  );
  return ((u64)(qh) << 32) | ql;
}
 
ll sigma0_sum(ll n) {
  // 1 から N までの約数の個数の総和を返す。
  // 計算量: O(n^(1/3+))
  u64 N=(u64)(n);
  const u32 v = isqrt(N);
  const u32 w = pow(N, 0.35);
 
  u64 x = N / v;
  u32 y = N / x + 1;
 
  std::stack< std::pair<u32, u32> > stac;
  stac.emplace(1, 0);
  stac.emplace(1, 1);
 
  auto outside = [N] (u64 x, u32 y) {
    return x * y > N;
  };
 
  auto cut_off = [N] (u64 x, u32 ldx, u32 ldy) {
    return u128::mul64(x, x * ldy) >= u128::mul64(N, ldx);
  };
 
  u128 ret = 0;
  while (1) {
    u32 ldx, ldy;
    std::tie(ldx, ldy) = stac.top(); stac.pop();
 
    while (outside(x + ldx, y - ldy)) {
      ret += x * ldy + (u64)(ldy + 1) * (ldx - 1) / 2;
      x += ldx; y -= ldy;
    }
 
    if (y <= w) {
      break;
    }
 
    u32 udx = ldx, udy = ldy;
    while (1) {
      std::tie(ldx, ldy) = stac.top();
      if (outside(x + ldx, y - ldy)) break;
      udx = ldx, udy = ldy;
      stac.pop();
    }
 
    while (1) {
      u32 mdx = ldx + udx;
      u32 mdy = ldy + udy;
      if (outside(x + mdx, y - mdy)) {
        ldx = mdx, ldy = mdy;
        stac.emplace(ldx, ldy);
      } else {
        if (cut_off(x + mdx, ldx, ldy)) {
          break;
        }
        udx = mdx, udy = mdy;
      }
    }
  }
  for (u32 i = y - 1; i > 0; --i) {
    ret += div64_32(N, i);
  }
  ret = ret + ret - (u64)(v) * v;
  return stoll(ret.to_string().c_str());
}



template <ll mod> struct Fp{
    ll x;
    constexpr Fp(ll x=0) noexcept : x ((x%mod+mod)%mod){}
    /*-------ここから演算子-------*/
    constexpr Fp operator-() const noexcept{
        return Fp(-x);
    }
    constexpr Fp& operator+=(const Fp &a) noexcept{
        if((x+=a.x)>=mod)x-=mod;
        return *this;
    }
    constexpr Fp& operator-=(const Fp &a) noexcept{
        if((x+=mod-a.x)>=mod)x-=mod;
        return *this;
    }
    constexpr Fp& operator*=(const Fp &a) noexcept{
        (x*=a.x)%=mod;
        return *this;
    }
    constexpr Fp& operator/=(const Fp &a) noexcept{
        return (*this)*=a.inv();
    }
    constexpr Fp operator+(const Fp &a)const noexcept{
        Fp res(*this);
        return res+=a;
    }
    constexpr Fp operator-(const Fp &a)const noexcept{
        Fp res(*this);
        return res-=a;
    }
    constexpr Fp operator*(const Fp &a)const noexcept{
        Fp res(*this);
        return res*=a;
    }
    constexpr Fp operator/(const Fp &a)const noexcept{
        Fp res(*this);
        return res/=a;
    }

    constexpr bool operator==(const Fp &a)const noexcept{
        if((*this).x==a.x)return true;
        else return false;
    }
    constexpr bool operator!=(const Fp &a)const noexcept{
        return !((*this)==a);
    }
    constexpr bool operator<(const Fp &a)const noexcept{
        if((*this).x<a.x)return true;
        else return false;
    }
    constexpr bool operator>(const Fp &a)const noexcept{
        if((*this).x>a.x)return true;
        else return false;
    }
    constexpr bool operator<=(const Fp &a)const noexcept{
        return !((*this)>a);
    }
    constexpr bool operator>=(const Fp &a)const noexcept{
        return !((*this)<a);
    }
    /*-------ここまで演算子-------*/

    constexpr Fp pow(const ll &t)const noexcept{
        if(t==0)return (Fp)1;
        if(t%2==0){
            Fp k=pow(t/2);
            return k*k;
        }else{
            if(t>0)return (*this)*pow(t-1);
            else return pow(t+1)/(*this);
        }
    }
    constexpr Fp inv()const noexcept{
        return (Fp)forinv((*this).x,mod).first;
    }


    friend ostream& operator<<(ostream &os,const Fp &m)noexcept{
        os << m.x;
        return os;
   }
    friend istream& operator>>(istream &is,Fp &m)noexcept{
        is >> m.x;
        m.x%=mod;
        if(m.x<0)m.x=(m.x+mod)%mod;
        return is;
    }

private:
    //Euclidの互除法
    constexpr pair<ll,ll> forinv(ll a,ll b)const noexcept{
        if(b==0)return pair<ll,ll>(1,0);

        pair<ll,ll> ans = forinv(b,a%b);
        pair<ll,ll> t=pair<ll,ll>(ans.second,ans.first-a/b*ans.second);
        return t;
    }
};

using mint=Fp<big>;

mint fi(ll x){
    //約数の個数関数
    map<ll,ll> primes;
    for(ll p=2;p*p<=x;p++){
        if(x%p==0){
            while(x%p==0){
                x/=p;
                primes[p]++;
            }
        }
    }
    if(x!=1)primes[x]++;
    mint res=1;
    for(auto p : primes){
        res*=(p.second+1);
    }
    return res;
}
mint f(ll x){//非 約数の個数関数
    return (mint)x-fi(x);
}

int main(){
    ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    cout << fixed << setprecision(10);
    /*--------------------------------*/
    
    ll n;cin >> n;
/*  //愚直解
    mint ans=-(n+1)*n;
    rep2(i,1,n+1){
        ans+= -fi(i)*(n/i) - (n/i)*(n/i) + (n/i)*(n+3);
    }
    cout << ans << endl;*/

    mint cnt=-n*(n+1);

    //vector<pair<ll,ll>> で√nくらいの個数,[何回現れるか , 値]を持っておく→sum([n/i]) 1 to k(k<=√nくらい？)も前計算しておきたいかも…？
    vector<pair<ll,ll>> noveri;
    int s=1;
    while(s<=n){
        pair<ll,ll> ne={0,n/s};
        ll k=n/s;
        if(k==1){ne.first=n-n/2;noveri.push_back(ne);break;}
        ne.first=n/k - n/(k+1);
        noveri.push_back(ne);
        s=(n+k)/k;//(n+1)/kの切り上げ
    }
    //次に個数の累積和みたいなやつを作る
    vector<ll> sum;
    sum.push_back(0);
    rep(i,noveri.size()){
        sum.push_back(sum[i] + noveri[i].first );//sum[i]=Σnoveri[j] ,j=0 to i-1
    }

    rep(k,noveri.size()){
        ll num=noveri[k].first;//その値を取る個数
        ll val=noveri[k].second;//その時見ている[n/i]の値
        cnt+=(-(mint)val*val+(mint)(n+3)*val)*num;
        //[n/i]の値を取る間のσ(i)の和

        if(k==0)cnt-=(mint)sigma0_sum(sum[k+1]) * (mint)val;
        else cnt-=(mint)(sigma0_sum(sum[k+1])- sigma0_sum(sum[k]))*(mint)val;
    }
    cout << cnt << endl;
}