fork download
copy 
  1. #include <cassert>
  2. #include <cstdio>
  3. #include <cstring>
  4. #include <cstdlib>
  5. #include <cmath>
  6. #include <ctime>
  7. #include <cctype>
  8.  
  9. #include <algorithm>
  10. #include <random>
  11. #include <bitset>
  12. #include <queue>
  13. #include <functional>
  14. #include <set>
  15. #include <map>
  16. #include <vector>
  17. #include <chrono>
  18. #include <iostream>
  19. #include <iomanip>
  20. #include <limits>
  21. #include <numeric>
  22.  
  23. typedef long long ll;
  24. typedef long double ld;
  25.  
  26. // quaternion from yhx-12243: https://y...content-available-to-author-only...b.io/OI-transit/records/cf1375I.html
  27. struct quaternion {
  28. 	ld s, x, y, z;
  29. 	quaternion (ld s0 = 0, ld x0 = 0, ld y0 = 0, ld z0 = 0) : s(s0), x(x0), y(y0), z(z0) {}
  30. 	inline quaternion operator - () const {return quaternion(-s, -x, -y, -z);}
  31. 	inline quaternion conj() const {return quaternion(s, -x, -y, -z);}
  32. 	inline quaternion operator + (const quaternion &B) const {return quaternion(s + B.s, x + B.x, y + B.y, z + B.z);}
  33. 	inline quaternion operator - (const quaternion &B) const {return quaternion(s - B.s, x - B.x, y - B.y, z - B.z);}
  34. 	inline quaternion operator * (const ld k) const {return quaternion(s * k, x * k, y * k, z * k);}
  35. 	inline quaternion operator * (const quaternion &B) const {
  36. 		return quaternion(
  37. 			s * B.s - x * B.x - y * B.y - z * B.z,
  38. 			s * B.x + x * B.s + y * B.z - z * B.y,
  39. 			s * B.y - x * B.z + y * B.s + z * B.x,
  40. 			s * B.z + x * B.y - y * B.x + z * B.s
  41. 		);
  42. 	}
  43. 	inline operator void * () const {return s || x || y || z ? (void *)this : 0;}
  44. 	inline ld norm2() const {return s * s + x * x + y * y + z * z;}
  45. 	inline ll intnorm2() const {return llroundl(norm2());}
  46. 	inline quaternion inv() const {return conj() * (1.l / norm2());}
  47. 	inline quaternion round() const {
  48. 		quaternion A(roundl(s), roundl(x), roundl(y), roundl(z)), B(floorl(s) + .5l, floorl(x) + .5l, floorl(y) + .5l, floorl(z) + .5l );
  49. 		return (*this - A).norm2() <= (*this - B).norm2() ? A : B;
  50. 	}
  51. 	inline friend quaternion Ldiv(const quaternion &A, const quaternion &B) {return (B.inv() * A).round();}
  52. 	inline friend quaternion Rdiv(const quaternion &A, const quaternion &B) {return (A * B.inv()).round();}
  53. 	inline friend quaternion Lmod(const quaternion &A, const quaternion &B) {return A - B * Ldiv(A, B);}
  54. 	inline friend quaternion Rmod(const quaternion &A, const quaternion &B) {return A - Rdiv(A, B) * B;}
  55. 	// find (Hurwitz) quaternion d where A = d u, B = d v.
  56. 	inline friend quaternion Lgcd(quaternion A, quaternion B) {
  57. 		static quaternion tmp;
  58. 		for (; B; tmp = Lmod(A, B), A = B, B = tmp);
  59. 		return A;
  60. 	}
  61. 	inline bool is_int() const {return fabsl(s - roundl(s)) < 1e-10l && fabsl(x - roundl(x)) < 1e-10l && fabsl(y - roundl(y)) < 1e-10l && fabsl(z - roundl(z)) < 1e-10l;}
  62. };
  63.  
  64. const int N = 3000;
  65.  
  66. int R;
  67. int sqfree[N];
  68.  
  69. ll ans;
  70.  
  71. void test(quaternion q) {
  72. 	if (q.intnorm2() > R) return;
  73. 	if (q == quaternion()) return;
  74. 	quaternion unit, res;
  75.  
  76. 	unit = quaternion(0, 1, 0, 0); res = q * unit * q.conj();
  77. 	int a = std::abs(res.x) + std::abs(res.y) + std::abs(res.z);
  78. 	unit = quaternion(0, 0, 1, 0); res = q * unit * q.conj();
  79. 	int b = std::abs(res.x) + std::abs(res.y) + std::abs(res.z);
  80. 	unit = quaternion(0, 0, 0, 1); res = q * unit * q.conj();
  81. 	int c = std::abs(res.x) + std::abs(res.y) + std::abs(res.z);
  82.  
  83. 	for (int i = 1; i <= (R - 1) / std::max({a, b, c}); ++i)
  84. 		if (sqfree[i])
  85. 			ans += ll(R - i * a) * (R - i * b) * (R - i * c);
  86. }
  87.  
  88. int main() {
  89. 	std::cin >> R;
  90.  
  91. 	std::fill(sqfree + 1, sqfree + R + 1, 1);
  92. 	for (int i = 2; i <= R; i += 2) sqfree[i] = 0;
  93. 	for (int i = 2; i * i <= R; ++i)
  94. 		for (int j = i * i; j <= R; j += i * i)
  95. 			sqfree[j] = 0;
  96.  
  97. 	int rt = ceil(sqrt(R)) + 1;
  98. 	for (int s = -rt; s <= rt; ++s) for (int x = -rt; x <= rt; ++x) for (int y = -rt; y <= rt; ++y) for (int z = -rt; z <= rt; ++z) {
  99. 		test(quaternion(s, x, y, z));
  100. 		test(quaternion(s + .5, x + .5, y + .5, z + .5));
  101. 	}
  102.  
  103. 	std::cout << ans / 24 << '\n';
  104.  
  105. 	return 0;
  106. }
  107.  
Success #stdin #stdout 7.67s 5548KB
comments (?)
stdin
copy 
2021
stdout
copy 
759516000281546
https://ideone.com/reAYcy
language:
C++14 (gcc 8.3)
created:
4 years ago
visibility:
public

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!