E - GCD of Subset
Source: atcoder
Difficulty: 475

Problem Statement
You are given a sequence
A = (A_1, A_2, \dots, A_N)
of length
N
and a positive integer
K
(at most
N
).
For each
i = 1, 2, \dots, N
, solve the following problem:
When you choose
K
elements from
A
that include
A_i
, find the maximum possible GCD (greatest common divisor) of those chosen elements.
Constraints
1 \leq K \leq N \leq 1.2 \times 10^6
1 \leq A_i \leq 10^6
All input values are integers.
Input
The input is given from Standard Input in the following format:
N
K
A_1
A_2
\dots
A_N
Output
Print
N
lines. The
j
-th line should contain the answer for
i=j
.
Sample Input 1
5 2
3 4 6 7 12
Sample Output 1
3
4
6
1
6
For
i=1
, choosing
A_1
and
A_3
yields
\gcd(\lbrace 3,6 \rbrace) = 3
, which is the maximum.
For
i=2
, choosing
A_2
and
A_5
yields
\gcd(\lbrace 4,12 \rbrace) = 4
, which is the maximum.
For
i=3
, choosing
A_3
and
A_5
yields
\gcd(\lbrace 6,12 \rbrace) = 6
, which is the maximum.
For
i=4
, choosing
A_4
and
A_2
yields
\gcd(\lbrace 7,4 \rbrace) = 1
, which is the maximum.
For
i=5
, choosing
A_5
and
A_3
yields
\gcd(\lbrace 12,6 \rbrace) = 6
, which is the maximum.
Sample Input 2
3 3
6 10 15
Sample Output 2
1
1
1
Sample Input 3
10 3
414003 854320 485570 52740 833292 625990 909680 885153 435420 221663
Sample Output 3
59
590
590
879
879
590
20
879
590
59

C++ solution:
// competitive-verifier: PROBLEM
#include <array>
#include <cmath>
#include <cstdint>
#include <vector>
/// @brief エラトステネスの篩
/// @see https://qiita.com/peria/items/a4ff4ddb3336f7b81d50
template <int N = (1 << 22)>
struct eratosthenes {
  private:
    static constexpr int SIZE = N / 30 + (N % 30 != 0);
    static constexpr std::array<int, 8> kMod30 = {1, 7, 11, 13, 17, 19, 23, 29};
    static constexpr std::array<int, 8> C1 = {6, 4, 2, 4, 2, 4, 6, 2};
    static constexpr std::array<std::array<int, 8>, 8> C0 = {
        0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 2, 2, 0, 2, 0, 2,
        2, 1, 3, 1, 1, 2, 1, 1, 3, 1, 3, 3, 1, 2, 1, 3, 3, 1, 4, 2, 2, 2,
        2, 2, 4, 1, 5, 3, 1, 4, 1, 3, 5, 1, 6, 4, 2, 4, 2, 4, 6, 1,
    };
    static constexpr std::array<std::array<std::uint8_t, 8>, 8> kMask = {
        0xfe, 0xfd, 0xfb, 0xf7, 0xef, 0xdf, 0xbf, 0x7f, 0xfd, 0xdf, 0xef, 0xfe, 0x7f,
        0xf7, 0xfb, 0xbf, 0xfb, 0xef, 0xfe, 0xbf, 0xfd, 0x7f, 0xf7, 0xdf, 0xf7, 0xfe,
        0xbf, 0xdf, 0xfb, 0xfd, 0x7f, 0xef, 0xef, 0x7f, 0xfd, 0xfb, 0xdf, 0xbf, 0xfe,
        0xf7, 0xdf, 0xf7, 0x7f, 0xfd, 0xbf, 0xfe, 0xef, 0xfb, 0xbf, 0xfb, 0xf7, 0x7f,
        0xfe, 0xef, 0xdf, 0xfd, 0x7f, 0xbf, 0xdf, 0xef, 0xf7, 0xfb, 0xfd, 0xfe,
    };
  public:
    constexpr eratosthenes() {
        prime_number.fill(0xff);
        prime_number[0] = 0xfe;
        if (int r = N % 30) {
            if (r < 7) prime_number[SIZE - 1] = 0x1;
            else if (r < 11) prime_number[SIZE - 1] = 0x3;
            else if (r < 13) prime_number[SIZE - 1] = 0x7;
            else if (r < 17) prime_number[SIZE - 1] = 0xf;
            else if (r < 19) prime_number[SIZE - 1] = 0x1f;
            else if (r < 23) prime_number[SIZE - 1] = 0x3f;
            else if (r < 29) prime_number[SIZE - 1] = 0x7f;
        }
        const std::uint64_t sqrt_x = std::ceil(std::sqrt(N) + 0.1);
        const std::uint64_t sqrt_xi = sqrt_x / 30 + 1;
        for (std::uint64_t i = 0; i < sqrt_xi; ++i) {
            for (std::uint8_t flags = prime_number[i]; flags; flags &= flags - 1) {
                std::uint8_t lsb = flags & (-flags);
                const int ibit = __builtin_ctz(lsb);
                const int m = kMod30[ibit];
                const int pm = 30 * i + 2 * m;
                for (std::uint64_t j = i * pm + (m * m) / 30, k = ibit; j < SIZE;
                     j += i * C1[k] + C0[ibit][k], k = (k + 1) & 7) {
                    prime_number[j] &= kMask[ibit][k];
                }
            }
        }
    }
    /// @brief 素数判定
    bool is_prime(int x) const {
        switch (x % 30) {
            case 1: return prime_number[x / 30] >> 0 & 1;
            case 7: return prime_number[x / 30] >> 1 & 1;
            case 11: return prime_number[x / 30] >> 2 & 1;
            case 13: return prime_number[x / 30] >> 3 & 1;
            case 17: return prime_number[x / 30] >> 4 & 1;
            case 19: return prime_number[x / 30] >> 5 & 1;
            case 23: return prime_number[x / 30] >> 6 & 1;
            case 29: return prime_number[x / 30] >> 7 & 1;
        }
        if (x < 6) {
            if (x == 2) return true;
            if (x == 3) return true;
            if (x == 5) return true;
        }
        return false;
    }
    std::vector<int> prime_numbers(int x) const {
        if (x < 2) return std::vector<int>();
        std::vector<int> res = {2};
        for (int i = 3; i <= x; i += 2) {
            if (is_prime(i)) res.emplace_back(i);
        }
        return res;
    }
  private:
    std::array<std::uint8_t, SIZE> prime_number;
};
#ifdef ATCODER
#pragma GCC target("sse4.2,avx512f,avx512dq,avx512ifma,avx512cd,avx512bw,avx512vl,bmi2")
#endif
#pragma GCC optimize("Ofast,fast-math,unroll-all-loops")
#include <bits/stdc++.h>
#ifndef ATCODER
#pragma GCC target("sse4.2,avx2,bmi2")
#endif
template <class T, class U>
constexpr bool chmax(T &a, const U &b) {
    return a < (T)b ? a = (T)b, true : false;
}
template <class T, class U>
constexpr bool chmin(T &a, const U &b) {
    return (T)b < a ? a = (T)b, true : false;
}
constexpr std::int64_t INF = 1000000000000000003;
constexpr int Inf = 1000000003;
constexpr double EPS = 1e-7;
constexpr double PI = 3.14159265358979323846;
#define FOR(i, m, n) for (int i = (m); i < int(n); ++i)
#define FORR(i, m, n) for (int i = (m)-1; i >= int(n); --i)
#define FORL(i, m, n) for (int64_t i = (m); i < int64_t(n); ++i)
#define rep(i, n) FOR (i, 0, n)
#define repn(i, n) FOR (i, 1, n + 1)
#define repr(i, n) FORR (i, n, 0)
#define repnr(i, n) FORR (i, n + 1, 1)
#define all(s) (s).begin(), (s).end()
struct Sonic {
    Sonic() {
        std::ios::sync_with_stdio(false);
        std::cin.tie(nullptr);
        std::cout << std::fixed << std::setprecision(20);
    }
    constexpr void operator()() const {}
} sonic;
using namespace std;
using ll = std::int64_t;
using ld = long double;
template <class T, class U>
std::istream &operator>>(std::istream &is, std::pair<T, U> &p) {
    return is >> p.first >> p.second;
}
template <class T>
std::istream &operator>>(std::istream &is, std::vector<T> &v) {
    for (T &i : v) is >> i;
    return is;
}
template <class T, class U>
std::ostream &operator<<(std::ostream &os, const std::pair<T, U> &p) {
    return os << '(' << p.first << ',' << p.second << ')';
}
template <class T>
std::ostream &operator<<(std::ostream &os, const std::vector<T> &v) {
    for (auto it = v.begin(); it != v.end(); ++it) os << (it == v.begin() ? "" : " ") << *it;
    return os;
}
template <class Head, class... Tail>
void co(Head &&head, Tail &&...tail) {
    if constexpr (sizeof...(tail) == 0) std::cout << head << '\n';
    else std::cout << head << ' ', co(std::forward<Tail>(tail)...);
}
template <class Head, class... Tail>
void ce(Head &&head, Tail &&...tail) {
    if constexpr (sizeof...(tail) == 0) std::cerr << head << '\n';
    else std::cerr << head << ' ', ce(std::forward<Tail>(tail)...);
}
void Yes(bool is_correct = true) { std::cout << (is_correct ? "Yes\n" : "No\n"); }
void No(bool is_not_correct = true) { Yes(!is_not_correct); }
void YES(bool is_correct = true) { std::cout << (is_correct ? "YES\n" : "NO\n"); }
void NO(bool is_not_correct = true) { YES(!is_not_correct); }
void Takahashi(bool is_correct = true) { std::cout << (is_correct ? "Takahashi" : "Aoki") << '\n'; }
void Aoki(bool is_not_correct = true) { Takahashi(!is_not_correct); }
eratosthenes pn;
int main(void) {
    int n, k;
    cin >> n >> k;
    vector<int> a(n);
    cin >> a;
    int l = 1000000;
    vector<int> b(l + 1);
    for (auto x : a) ++b[x];
    for (int p = 2; p <= l; ++p) {
        if (!pn.is_prime(p))
            continue;
        for (int i = l / p; i >= 1; --i) {
            b[i] += b[i * p];
        }
    }
    vector<int> ans(l + 1, 1);
    for (int i = 2; i <= l; ++i) {
        if (b[i] >= k) ans[i] = i;
    }
    for (int p = 2; p <= l; ++p) {
        if (!pn.is_prime(p))
            continue;
        for (int i = 2; i * p <= l; ++i) {
            chmax(ans[i * p], ans[i]);
        }
    }
    
    rep (i, n) co(ans[a[i]]);
    return 0;
}

Go solution:
package main

import (
	"bufio"
	"fmt"
	"io"
	"os"
)

func main() { solveE(os.Stdin, os.Stdout) }

func solveE(_r io.Reader, _w io.Writer) {
	in, out := bufio.NewReader(_r), bufio.NewWriter(_w)
	defer out.Flush()

	n, k := readTwoNums(in)
	a := readNNums(in, n)
	M := 0
	for _, v := range a {
		M = max(M, v+1)
	}

	cnt := make([]int, M+1)
	cntFact := make([]int, M+1)
	for _, v := range a {
		cnt[v]++
	}
	for i := 1; i < M; i++ {
		for j := i; j < M; j += i {
			cntFact[i] += cnt[j]
		}
	}

	res := make([]int, M+1)
	for i := 1; i < M; i++ {
		if cntFact[i] < k {
			continue
		}
		for j := i; j < M; j += i {
			res[j] = i
		}
	}

	for _, v := range a {
		fmt.Fprintln(out, res[v])
	}
}

func solveD(_r io.Reader, _w io.Writer) {
	in, out := bufio.NewReader(_r), bufio.NewWriter(_w)
	defer out.Flush()

	readNum(in)
	s := readString(in)
	idx1 := []int{}
	for i, ch := range s {
		if ch == '1' {
			idx1 = append(idx1, i)
		}
	}
	m := len(idx1)
	idx2 := make([]int, m)
	idx2[m/2] = idx1[m/2]
	for i := m/2 - 1; i >= 0; i-- {
		idx2[i] = idx2[i+1] - 1
	}
	for i := m/2 + 1; i < m; i++ {
		idx2[i] = idx2[i-1] + 1
	}
	ans := 0
	for i := 0; i < m; i++ {
		ans += abs(idx1[i] - idx2[i])
	}
	fmt.Fprintln(out, ans)
}

func solveC(_r io.Reader, _w io.Writer) {
	in, out := bufio.NewReader(_r), bufio.NewWriter(_w)
	defer out.Flush()

	_, m := readTwoNums(in)
	type pair struct{ x, y int }
	mp := map[pair]struct{}{}
	ans := 0
	for i := 0; i < m; i++ {
		u, v := readTwoNums(in)
		u, v = min(u, v), max(u, v)
		if _, ok := mp[pair{u, v}]; ok || u == v {
			ans++
		} else {
			mp[pair{u, v}] = struct{}{}
		}
	}
	fmt.Fprintln(out, ans)
}

func abs(x int) int {
	if x < 0 {
		return -x
	}
	return x
}

func max(a, b int) int {
	if a >= b {
		return a
	}
	return b
}

func min(a, b int) int {
	if a <= b {
		return a
	}
	return b
}

func readString(reader *bufio.Reader) string {
	s, _ := reader.ReadString('\n')
	for i := 0; i < len(s); i++ {
		if s[i] == '\n' || s[i] == '\r' {
			return s[:i]
		}
	}
	return s
}

func readInt(bytes []byte, from int, val *int) int {
	i := from
	sign := 1
	if bytes[i] == '-' {
		sign = -1
		i++
	}
	tmp := 0
	for i < len(bytes) && bytes[i] >= '0' && bytes[i] <= '9' {
		tmp = tmp*10 + int(bytes[i]-'0')
		i++
	}
	*val = tmp * sign
	return i
}

func readNum(reader *bufio.Reader) (a int) {
	bs, _ := reader.ReadBytes('\n')
	readInt(bs, 0, &a)
	return
}

func readTwoNums(reader *bufio.Reader) (a int, b int) {
	res := readNNums(reader, 2)
	a, b = res[0], res[1]
	return
}

func readThreeNums(reader *bufio.Reader) (a int, b int, c int) {
	res := readNNums(reader, 3)
	a, b, c = res[0], res[1], res[2]
	return
}

func readNNums(reader *bufio.Reader, n int) []int {
	res := make([]int, n)
	x := 0
	bs, _ := reader.ReadBytes('\n')
	for i := 0; i < n; i++ {
		for x < len(bs) && (bs[x] < '0' || bs[x] > '9') && bs[x] != '-' {
			x++
		}
		x = readInt(bs, x, &res[i])
	}
	return res
}

Java solution:
import java.io.*;
import java.util.*;

public class Main {
    static class MyScanner {
        static BufferedReader r;
        static StringTokenizer st;
        public MyScanner() {
            r = new BufferedReader(new InputStreamReader(System.in));
        }
        public String next() {
            try {
                while (st == null || !st.hasMoreTokens()) {
                    st = new StringTokenizer(r.readLine());
                }
                return st.nextToken();
            } catch (Exception e) {
                return null;
            }
        }
        public int nextInt() {
            return Integer.parseInt(next());
        }
        public long nextLong() {
            return Long.parseLong(next());
        }
        public double nextDouble() {
            return Double.parseDouble(next());
        }
    }
    static BufferedWriter bw = new BufferedWriter(new OutputStreamWriter(System.out));
    static MyScanner sc = new MyScanner();
    static int N, K;
    static final int MAXN = 1200010, MAXA = 1000010;
    static int[] a = new int[MAXN], ans = new int[MAXA], cnt = new int[MAXA];
    public static void main(String[] args) throws IOException {
        N = sc.nextInt();
        K = sc.nextInt();
        for (int i = 1; i <= N; i++) {
            a[i] = sc.nextInt();
            cnt[a[i]]++;
        }
        for (int i = 1; i < MAXA; i++) {
            int num = 0;
            for (int j = i; j < MAXA; j += i) {
                num += cnt[j];
            }
            if (num >= K) {
                for (int j = i; j < MAXA; j += i) {
                    ans[j] = i;
                }
            }
        }
        for (int i = 1; i <= N; i++) {
            bw.write(ans[a[i]] + "\n");
        }
        bw.flush();
    }
}

JavaScript solution:
function Main(input) {
  input = input.split("\n").map((line) => line.trim());
  const [N, K] = input[0].split(" ").map(Number);
  const A = input[1].split(" ").map(Number);
  const max = Math.max(...A) + 1;

  const freq = new Array(max).fill(0);
  for (let i = 0; i < N; i++) freq[A[i]]++;
  const f = new Array(max).fill(0);
  for (let i = 1; i < max; i++) {
    for (let j = i; j < max; j += i) f[i] += freq[j];
  }

  const g = Array(max).fill(1);
  for (let i = 2; i < max; i++) {
    if (f[i] < K) continue;
    for (let j = i; j < max; j += i) g[j] = i;
  }

  const ans = [];
  for (let i = 0; i < N; i++) ans.push(g[A[i]]);
  console.log(ans.join("\n"));
}
class Primetools {
  /**
   * @description 高速素数列挙（エラトステネスの篩）, O(NloglogN)
   * @param {number} n　1~nまでの整数の中で素数のリストを返す
   * @returns {number[]} n以下の素数のリスト
   */
  static sievePrimes(n) {
    const isPrime = Array(n + 1).fill(true);
    isPrime[0] = isPrime[1] = false;
    for (let i = 2; i * i <= n; i++) {
      if (!isPrime[i]) continue;
      for (let j = i * i; j <= n; j += i) isPrime[j] = false;
    }
    const primes = [];
    for (let i = 2; i <= n; i++) if (isPrime[i]) primes.push(i);
    return primes;
  }

  /**
   * @description 高速約数列挙（Smallest Prime Factor）, O(NloglogN + NlogN)
   * @param {number} n　1~nまでの整数の中で各整数の約数のリストを返す
   * @returns {number[][]} n以下の各整数の約数のリスト
   * @warning ソートされていないため、必要に応じてソートすること
   */
  static sieveDivisors(n) {
    const spf = new Array(n + 1).fill(0).map((_, i) => i);
    for (let i = 2; i * i <= n; i++) {
      if (spf[i] !== i) continue;
      for (let j = i * i; j <= n; j += i) if (spf[j] === j) spf[j] = i;
    }

    const divisors = new Array(n + 1).fill(0).map(() => []);
    divisors[1] = [1];
    for (let i = 2; i <= n; i++) {
      let num = i;
      const p = spf[num];
      let cnt = 0;
      while (num % p === 0) {
        num /= p;
        cnt++;
      }

      const prevDivs = divisors[num];
      const curDivs = [];
      let p_pow = 1;
      for (let exp = 0; exp <= cnt; exp++) {
        for (const d of prevDivs) curDivs.push(d * p_pow);
        p_pow *= p;
      }
      divisors[i] = curDivs;
    }
    return divisors;
  }

  /**
   * @description 単一数値の素因数分解, O(sqrt(N))
   * @param {number} n　整数N
   * @returns {number[]} 素因数分解の結果
   */
  static primeFactors(n) {
    const res = [];
    for (let i = 2; i * i <= n; i++) {
      while (n % i === 0) {
        res.push(i);
        n /= i;
      }
    }
    if (n !== 1) res.push(n);
    return res;
  }

  /**
   * @description 単一数値の約数列挙, O(sqrt(N))
   * @param {number} n　整数N
   * @returns {number[]} 約数のリスト
   */
  static divisorsOf(n) {
    const res = [];
    for (let i = 1; i * i <= n; i++) {
      if (n % i === 0) {
        res.push(i);
        if (i !== n / i) res.push(n / i);
      }
    }
    return res;
  }

  /**
   * @description 素数判定, O(sqrt(N))
   * @param {number} n　整数N
   * @returns {boolean} 素数かどうか
   */
  static isPrime(n) {
    if (n < 2) return false;
    for (let i = 2; i * i <= n; i++) if (n % i === 0) return false;
    return true;
  }

  /**
   * @description 最大公約数, O(log(max(a,b)))
   * @param {number} a
   * @param {number} b
   * @returns {number} aとbの最大公約数
   */
  static gcd(a, b) {
    if (b === 0) return a;
    return Primetools.gcd(b, a % b);
  }

  /**
   * @description 最小公倍数, O(log(max(a,b)))
   * @param {number} a
   * @param {number} b
   * @returns {number} aとbの最小公倍数
   * @warning a * b がオーバーフローする可能性があるため、適宜注意すること
   */
  static lcm(a, b) {
    return (a * b) / Primetools.gcd(a, b);
  }

  /**
   * @description 拡張ユークリッドの互除法, O(log(max(a,b)))
   * 最大公約数を求めるだけでなく、不定方程式 ax + by = gcd(a,b) の整数解(x,y)を求める
   * @param {number} a
   * @param {number} b
   * @returns {number[]} [gcd(a,b), x, y] を返す
   * @example ax = 1 (mod m) の整数解xを求める場合、[gcd(a,m), x, y] = Primetools.ext_gcd(a, m) として x = x % m とすればよい
   */
  static ext_gcd(a, b) {
    if (b === 0) return [a, 1, 0];
    const [gcd, x, y] = Primetools.ext_gcd(b, a % b);
    return [gcd, y, x - Math.floor(a / b) * y];
  }

  /**
   * @description 約数ゼータ変換, O(NlogN)
   * g = Σ f(d) (dはnの約数) としたとき、fを求める
   * @param {number[]} f 約数ごとの寄与を表す配列
   * @returns {number[]} 約数ゼータ変換後の配列
   */
  static div_zt(f) {
    const n = f.length;
    const g = f.slice();
    for (let i = 1; i < n; i++) {
      for (let j = 2; i * j < n; j++) g[i * j] += f[i];
    }
    return g;
  }

  /**
   * @description メビウス関数 (線形篩), O(N)
   * @param {number} n 0~nまでの整数のメビウス関数のリストを返す
   * @returns {number[]} n以下のメビウス関数のリスト
   */
  static getMobius(n) {
    const mu = Array(n + 1).fill(0);
    const primes = [];
    const isComposite = Array(n + 1).fill(false);
    mu[1] = 1;
    for (let i = 2; i <= n; i++) {
      if (!isComposite[i]) {
        primes.push(i);
        mu[i] = -1;
      }
      for (let j = 0; j < primes.length && i * primes[j] <= n; j++) {
        const x = i * primes[j];
        isComposite[x] = true;
        if (i % primes[j] === 0) {
          mu[x] = 0;
          break;
        } else {
          mu[x] = -mu[i];
        }
      }
    }
    return mu;
  }

  /**
   * @description メビウス逆変換, O(d(n)) ※d(n)はnの約数の個数
   * @param {number[]} g 約数ゼータ変換後の配列
   * @param {number} x 逆変換により求めたい値
   * @param {number[]} mu メビウス関数のリスト
   * @returns {number} f(x)の値
   */
  static mu_inv(g, x, mu) {
    let res = 0;
    for (let i = 1; i * i <= x; i++) {
      if (x % i === 0) res += mu[i] * g[x / i];
      if (i * i !== x) res += mu[x / i] * g[i];
    }
    return res;
  }
}


Main(require("fs").readFileSync(0, "utf8"));

Python solution:
def find_answer(N, K, A):
    M = max(A)

    s = [0] * (M + 1)
    for num in A:
        s[num] += 1

    t = [0] * (M + 1)
    for d in range(1, M + 1):
        for n in range(d, M + 1, d):
            t[d] += s[n]

    u = [0] * (M + 1)
    for d in range(1, M + 1):
        if t[d] < K:
            continue
        for n in range(d, M + 1, d):
            u[n] = max(u[n], d)

    return [u[num] for num in A]


N, K = map(int, input().split())
A = [int(i) for i in input().split()]

answer = find_answer(N, K, A)
print("\n".join(map(str, answer)))

Ruby solution:
#!/usr/local/bin/ruby
require "prime"
class GCDofSubset
  def initialize(*arg)
    @n,@k,@a, = arg
    @m = @a.max
    @f = Array.new(@m+1,0)
  end

  private
  def set_freq
    @a.each{@f[_1]+=1}
  end
  def zeta_gcd(arr,m,primes)
    primes.each do |p|
      if p >= m
        break
      else
        c = (m-1)/p
        while c>=1
          arr[c] += arr[p*c]
          c-=1
        end
      end
    end
  end

  def zeta_max(arr,m,primes)
    primes.each do |p|
      if p>=m
        break
      else
        c = 1
        while c*p < m
          arr[c*p] = [arr[c*p],arr[c]].max
          c+=1
        end
      end
    end
  end

  public

  def ans
    set_freq
    primes = Prime.each(upper_bound=@m).to_a
    zeta_gcd(@f,@m+1,primes)
    h = Array.new(@m+1,1)
    (2..@m).each{h[_1] = _1 if @f[_1]>=@k}
    zeta_max(h,@m+1,primes)
    @a.map{h[_1]}
  end
end
### END: class GCDofSubset
iod = STDIN
n,k = iod.gets.split.map &:to_i
a =  iod.gets.split.map &:to_i
puts GCDofSubset.new(n,k,a).ans
exit 0