# How to compute GCD with a MILP?

How would one formulate a linear optimisation program that computes the Greatest Common Divisor (GCD) of a set of $$n$$ numbers $$\{a_1,...,a_n \}$$?

I can only think of a non linear formulation :

$$\max\; \{ p \}$$ subject to \begin{align*} a_i &= p \cdot q_i \quad &\forall i=1,...,n\\ q_i &\in \mathbb{N} \quad &\forall i=1,...,n\\ p &\in \mathbb{N} \end{align*}

One could linearize the product, but I am wondering if there is a better way to this (with an optimization program).

Let $$m = \min_{i \in \{1,\dots,n\}} a_i$$. For $$p\in\{1,\dots,m\}$$, let binary variable $$z_p$$ indicate whether $$p$$ is the GCD. The problem is to maximize $$\sum_{p=1}^m p\ z_p$$ subject to \begin{align} \sum_{p=1}^m z_p &= 1 \\ (a_i-p\ a_i)(1 - z_p) \le a_i - p\ q_i &\le a_i (1 - z_p) &&\text{for i \in \{1,\dots,n\}} \\ q_i &\in [0,a_i] \cap \mathbb{N} &&\text{for i \in \{1,\dots,n\}}\\ z_p &\in \{0,1\} &&\text{for p \in \{1,\dots,m\}}\\ \end{align}