#include <bits/stdc++.h>
using namespace std;
using vi = vector<int>;
using pi = pair<int, int>;
using vs = vector<string>;
using vpi = vector<pi>;
using ll = long long int;
template <typename T>
string tostr(T x) {
  ostringstream oss;
  oss << x;
  return oss.str();
}
template <typename T1, typename T2, typename T3>
T1 modpow(T1 a, T2 p, T3 mod) {
  T1 ret = 1;
  a %= mod;
  for (; p > 0; p /= 2) {
    if (p & 1) ret = 1ll * ret * a % mod;
    a = 1ll * a * a % mod;
  }
  return ret;
}
inline int scan() {
  int x;
  scanf("%d", &x);
  return x;
}
int dx[]{-1, 0, 1, 0, 1, 1, -1, -1};
int dy[]{0, -1, 0, 1, 1, -1, 1, -1};
constexpr auto PI = 3.14159265358979323846L;
constexpr auto oo = numeric_limits<int>::max() / 2 - 10;
constexpr auto eps = 1e-6;
constexpr auto mod = 1000000007;
const int ND = 2e5 + 10;
int par[ND], con[ND], res[ND];
bool loop[ND];
int f[ND];
void init(int n) {
  for (int i = 1; i <= n; ++i) par[i] = i;
}
int get(int x) {
  if (x != par[x]) par[x] = get(par[x]);
  return par[x];
}
bool merge(int x, int y) {
  x = get(x), y = get(y);
  if (x == y) return false;
  par[x] = y;
  return true;
}
void solve() {
  int n = scan();
  init(n);
  for (int i = 1; i <= n; ++i) {
    f[i] = scan();
    res[i] = f[i];
  }
  for (int i = 1; i <= n; ++i) {
    if (merge(i, f[i])) continue;
    res[i] = i;
  }
  int root = -1;
  for (int i = 1; i <= n and root == -1; ++i) {
    if (res[i] == i and f[i] == i) root = i;
  }
  for (int i = 1; i <= n and root == -1; ++i) {
    if (res[i] == i) root = i;
  }
  for (int i = 1; i <= n; ++i) {
    if (res[i] == i) res[i] = root;
  }
  int changes = 0;
  for (int i = 1; i <= n; ++i) {
    if (f[i] != res[i]) ++changes;
  }
  cout << changes << "\n";
  for (int i = 1; i <= n; ++i) cout << res[i] << " \n"[i == n];
}
int main() {
  solve();
  return 0;
}
