description: Let's call a set of positive integers a_1, a_2, ..., a_k quadratic if the product of the factorials of its elements is a square of an integer, i. e. ∏_{i=1}^{k} a_i! = m^2, for some integer m.
  
  You are given a positive integer n.
  
  Your task is to find a quadratic subset of a set 1, 2, ..., n of maximum size. If there are multiple answers, print any of them.
  
  Input
  
  A single line contains a single integer n (1 ≤ n ≤ 10^6).
  
  Output
  
  In the first line, print a single integer — the size of the maximum subset. In the second line, print the subset itself in an arbitrary order.
  
  Examples
  
  Input
  
  
  1
  
  
  Output
  
  
  1
  1 
  
  
  Input
  
  
  4
  
  
  Output
  
  
  3
  1 3 4 
  
  
  Input
  
  
  7
  
  
  Output
  
  
  4
  1 4 5 6 
  
  
  Input
  
  
  9
  
  
  Output
  
  
  7
  1 2 4 5 6 7 9 
solution: #include <bits/stdc++.h>
  using ULL = unsigned long long;
  std::vector<int> spfS(int N) {
    std::vector<int> sp(N);
    std::iota(sp.begin(), sp.end(), 0);
    for (int i = 2; i * i < N; ++i)
      if (sp[i] == i) {
        for (int j = i * i; j < N; j += i)
          if (sp[j] == j) {
            sp[j] = i;
          }
      }
    return sp;
  }
  std::vector<int> spf(int N) {
    std::vector<int> sp(N), p{0, 2};
    p.reserve(N);
    for (int i = 2; i < N; i += 2) sp[i] = 2;
    for (int i = 1; i < N; i += 2) sp[i] = i;
    for (int i = 3; i < N; i += 2) {
      if (sp[i] == i) p.emplace_back(i);
      for (int j = 2, np = (int)p.size(); j < np && p[j] <= sp[i] && i * p[j] < N;
           ++j) {
        sp[i * p[j]] = p[j];
      }
    }
    return sp;
  }
  std::mt19937_64 rnd64(
      std::chrono::steady_clock::now().time_since_epoch().count());
  void solve() {
    int n;
    std::cin >> n;
    auto sp = spfS(n + 1);
    std::vector<ULL> hsh(n + 1);
    for (int i = 2; i <= n; ++i) {
      if (sp[i] == i) {
        hsh[i] = rnd64();
      } else {
        hsh[i] = hsh[i / sp[i]] ^ hsh[sp[i]];
      }
    }
    ULL ans = 0;
    for (int i = n; i > 1; i -= 2) ans ^= hsh[i];
    if (ans == 0) {
      std::cout << n << '\n';
      for (int i = 1; i <= n; ++i) std::cout << i << " \n"[i == n];
      return;
    }
    std::map<ULL, int> mp;
    for (int i = 2; i <= n; ++i) hsh[i] ^= hsh[i - 1];
    for (int i = 1; i <= n; ++i) mp[hsh[i]] = i;
    if (auto it = mp.find(ans); it != mp.end()) {
      std::cout << n - 1 << '\n';
      for (int i = 1; i <= n; ++i) {
        if (i != it->second) std::cout << i << ' ';
      }
      std::cout << '\n';
      return;
    }
    for (int i = 2; i <= n; ++i) {
      ULL need = ans ^ hsh[i];
      auto it = mp.find(need);
      if (it != mp.end()) {
        std::cout << n - 2 << '\n';
        for (int j = 1; j <= n; ++j) {
          if (j != it->second && j != i) {
            std::cout << j << ' ';
          }
        }
        std::cout << '\n';
        return;
      }
    }
    assert(n % 4 == 3);
    std::cout << n - 3 << '\n';
    std::cout << 1;
    for (int i = 3; i < n; ++i) {
      if (i != n / 2) std::cout << ' ' << i;
    }
    std::cout << '\n';
  }
  int main() {
    std::cin.tie(nullptr)->sync_with_stdio(false);
    int cas = 1;
    while (cas--) {
      solve();
    }
    return 0;
  }
  