fork(5) download
copy 
  1. #include <cstdio>
  2. #include <cassert>
  3. #include <cmath>
  4.  
  5. #include <functional>
  6. #include <stack>
  7.  
  8. using u64 = unsigned long long;
  9. using u32 = unsigned;
  10.  
  11. struct u128 {
  12.   u128() : l(0), h(0) {}
  13.   u128(u64 l) : l(l), h(0) {}
  14.   u128(u64 l, u64 h) : l(l), h(h) {}
  15.  
  16.   static u128 mul64(u64 a, u64 b) {
  17.     u32 a_hi = a >> 32, a_lo = u32(a);
  18.     u32 b_hi = b >> 32, b_lo = u32(b);
  19.     u64 p0 = u64(a_lo) * b_lo;
  20.     u64 p10 = u64(a_hi) * b_lo + (p0 >> 32);
  21.     u64 p11 = u64(a_lo) * b_hi + p10;
  22.     u64 p2 = (u64(p11 < p10) << 32) + u64(a_hi) * b_hi + (p11 >> 32);
  23.     return u128(u32(p0) | (p11 << 32), p2);
  24.   }
  25.  
  26.   static void divmod(u128 a, u32 d, u128& q, u32& r) {
  27.     u32 a3 = a.h >> 32, a2 = a.h;
  28.     u32 a1 = a.l >> 32, a0 = a.l;
  29.  
  30.     u64 t = a3;
  31.     u32 q3 = t / d; t = ((t % d) << 32) | a2; 
  32.     u32 q2 = t / d; t = ((t % d) << 32) | a1;
  33.     u32 q1 = t / d; t = ((t % d) << 32) | a0;
  34.     u32 q0 = t / d; r = t % d;
  35.     q = u128(q0 | (u64(q1) << 32), q2 | (u64(q3) << 32));
  36.   }
  37.   bool is_zero() const {
  38.     return l == 0 && h == 0;
  39.   }
  40.   bool operator >= (const u128& rhs) const {
  41.     return h > rhs.h || (h == rhs.h && l >= rhs.l);
  42.   }
  43.   u128& operator += (const u64 rhs) {
  44.     u64 old_l = l;
  45.     l += rhs;
  46.     h += (l < old_l);
  47.     return *this;
  48.   }
  49.   u128& operator += (const u128& rhs) {
  50.     u64 old_l = l;
  51.     l += rhs.l;
  52.     h += (l < old_l) + rhs.h;
  53.     return *this;
  54.   }
  55.   u128& operator -= (const u128& rhs) {
  56.     u64 old_l = l;
  57.     l -= rhs.l;
  58.     h -= (l > old_l) + rhs.h;
  59.     return *this;
  60.   }
  61.   u128 operator + (const u64 rhs) const {
  62.     return u128(*this) += rhs;
  63.   }
  64.   u128 operator + (const u128& rhs) const {
  65.     return u128(*this) += rhs;
  66.   }
  67.   u128 operator - (const u128& rhs) const {
  68.     return u128(*this) -= rhs;
  69.   }
  70.   std::string to_string() const {
  71.     static const u32 ten9 = u32(1e9);
  72.     if (is_zero()) {
  73.       return "0";
  74.     }
  75.     char str[41];
  76.     int i = 0;
  77.  
  78.     u128 n = *this;
  79.     u32 r;
  80.     while (1) {
  81.       divmod(n, ten9, n, r);
  82.       if (n.is_zero()) {
  83.         while (r > 0) str[i++] = r % 10, r /= 10;
  84.         break;
  85.       } else {
  86.         for (int j = 0; j < 9; ++j) str[i++] = r % 10, r /= 10;
  87.       }
  88.     }
  89.     std::string ret;
  90.     while (i) ret += '0' + str[--i];
  91.     return ret;
  92.   }
  93.   u64 l, h;
  94. };
  95.  
  96. u32 isqrt(u64 n) {
  97.   if (n >= u64(UINT32_MAX) * UINT32_MAX) {
  98.     return UINT32_MAX;
  99.   } else {
  100.     u32 ret = sqrtl(n);
  101.     assert(u64(ret) * ret <= n);
  102.     assert(u64(ret + 1) * (ret + 1) > n);
  103.     return ret;
  104.   }
  105. }
  106.  
  107. inline u64 div64_32(u64 a, u32 b) {
  108.   u32 r, ql, qh;
  109.   u32 al = a, ah = a >> 32;
  110.   qh = ah / b; ah %= b;
  111.   __asm__ (
  112.     "divl %4\n\t"
  113.     : "=a"(ql), "=d"(r)
  114.     : "0"(al), "1"(ah), "rm"(b)
  115.   );
  116.   return (u64(qh) << 32) | ql;
  117. }
  118.  
  119. u128 sigma0_sum(u64 N) {
  120.   // 1 から N までの約数の個数の総和を返す。
  121.   // 計算量: O(n^(1/3+))
  122.  
  123.   const u32 v = isqrt(N);
  124.   const u32 w = pow(N, 0.35);
  125.  
  126.   u64 x = N / v;
  127.   u32 y = N / x + 1;
  128.  
  129.   std::stack< std::pair<u32, u32> > stac;
  130.   stac.emplace(1, 0);
  131.   stac.emplace(1, 1);
  132.  
  133.   auto outside = [N] (u64 x, u32 y) {
  134.     return x * y > N;
  135.   };
  136.  
  137.   auto cut_off = [N] (u64 x, u32 ldx, u32 ldy) {
  138.     return u128::mul64(x, x * ldy) >= u128::mul64(N, ldx);
  139.   };
  140.  
  141.   u128 ret = 0;
  142.   while (1) {
  143.     u32 ldx, ldy;
  144.     std::tie(ldx, ldy) = stac.top(); stac.pop();
  145.  
  146.     while (outside(x + ldx, y - ldy)) {
  147.       ret += x * ldy + u64(ldy + 1) * (ldx - 1) / 2;
  148.       x += ldx; y -= ldy;
  149.     }
  150.  
  151.     if (y <= w) {
  152.       break;
  153.     }
  154.  
  155.     u32 udx = ldx, udy = ldy;
  156.     while (1) {
  157.       std::tie(ldx, ldy) = stac.top();
  158.       if (outside(x + ldx, y - ldy)) break;
  159.       udx = ldx, udy = ldy;
  160.       stac.pop();
  161.     }
  162.  
  163.     while (1) {
  164.       u32 mdx = ldx + udx;
  165.       u32 mdy = ldy + udy;
  166.       if (outside(x + mdx, y - mdy)) {
  167.         ldx = mdx, ldy = mdy;
  168.         stac.emplace(ldx, ldy);
  169.       } else {
  170.         if (cut_off(x + mdx, ldx, ldy)) {
  171.           break;
  172.         }
  173.         udx = mdx, udy = mdy;
  174.       }
  175.     }
  176.   }
  177.   for (u32 i = y - 1; i > 0; --i) {
  178.     ret += div64_32(N, i);
  179.   }
  180.   ret = ret + ret - u64(v) * v;
  181.   return ret;
  182. }
  183.  
  184. u128 spiral_walk(u64 N) {
  185.   /*
  186.   横 w, 縦 h の格子状の道における距離は (w + 1) * (h + 1) - 1　となる。
  187.   したがって、正の整数 n に対して F(n) = (n + 1 の約数の個数) - 2 となる。
  188.   これを 1 から n まで足し合わせたものが G(n) であるから、
  189.   以下のような式で書くことができる。
  190.   */
  191.   return sigma0_sum(N + 1) - 1 - N - N;
  192. }
  193.  
  194. int main() {
  195.   u64 N;
  196.   while (~scanf("%llu", &N)) {
  197.     // G(1e19) = 419035480966899467511 ; 0.59 秒
  198.     // G(18446744065119617022) = 784279019989592462219 ; 0.71 秒
  199.     assert(1 <= N && N <= 0xFFFFFFFDFFFFFFFFull - 1);
  200.     const u128 ans = spiral_walk(N);
  201.     puts(ans.to_string().c_str());
  202.   }
  203. }
Success #stdin #stdout 0.64s 3480KB
comments (?)
stdin
copy 
1
1000
1000000
1000000000
1000000000000
1000000000000000
10000000000000000000
stdout
copy 
0
5076
11970037
18877697665
25785452449093
32693207724724373
419035480966899467511
https://ideone.com/lZnZTU
language:
C++14 (gcc 8.3)
created:
7 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!