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 namespace std;
  const int N = (int)1e6 + 50;
  mt19937_64 rng(23333);
  int n;
  unsigned long long h[N], v[N];
  int bad[N];
  int p[N];
  void output(vector<int> res) {
    for (int x : res) bad[x] = 1;
    cout << n - (int)(res).size() << '\n';
    for (int i = 1; i < (n + 1); ++i) {
      if (!bad[i]) cout << i << " ";
    }
    exit(0);
  }
  int main() {
    ios::sync_with_stdio(false);
    cin.tie(NULL);
    cin >> n;
    h[0] = 1;
    unsigned long long fp = 0;
    iota(p, p + n + 1, 0);
    for (int i = 2; i < (n + 1); ++i) {
      if (p[i] == i) {
        for (int j = i * 2; j <= n; j += i) p[j] = i;
      }
    }
    h[0] = 1;
    for (int i = 1; i < (n + 1); ++i) {
      if (i == 1)
        v[i] = 1;
      else if (p[i] == i)
        v[i] = rng();
      else
        v[i] = v[i / p[i]] ^ v[p[i]];
      h[i] = h[i - 1] ^ v[i];
      fp ^= h[i];
    }
    if (fp == 0) output({});
    map<unsigned long long, int> mp;
    for (int i = 1; i < (n + 1); ++i) mp[h[i] ^ fp] = i;
    if (mp.find(0) != mp.end()) {
      output({mp[0]});
    }
    for (int i = 1; i < (n + 1); ++i) {
      if (mp.find(h[i]) != mp.end()) {
        output({i, mp[h[i]]});
      }
    }
    output({2, (n - 1) / 2, n});
  }
  