Let $$$k$$$ be the maximum value of $$$a_i$$$, $$$b_i$$$ on input ($$$10^8$$$). Let $$$D(k)$$$ be the maximum number of divisors a number in range $$$[1, \ldots, k]$$$ can have and $$$P(k)$$$ the maximum number of prime divisors such number can have.
Let's think about how to solve one query (with $$$M$$$ coins) for now. Assume that in the optimal solution, we never swap $$$a_1$$$ and $$$b_1$$$. In the end, we will just run the same solution on input where $$$a_1$$$ and $$$b_1$$$ is swapped and $$$M$$$ is decreased by $$$c_1$$$, and then we will take the maximum of the two values these two runs will find.
So now, in the optimal solution, the gcd of $$$a$$$ is always a divisor of $$$a_1$$$, and the $$$\gcd$$$ of $$$b$$$ is a divisor of $$$b_1$$$. Let's start by factorizing the two numbers in $$$O(\sqrt k)$$$.
Now we will precalculate something in order to be efficiently able to determine for every $$$X$$$, a divisor of $$$a$$$, and $$$Y$$$, a divisor of $$$b_1$$$, whether we can have $$$\gcd(a_1, \ldots, a_n) = X$$$ and $$$gcd(b_1, ..., b_n) = Y$$$ simultaneously (and also the minimum cost of making it so). Then we can obviously answer the query by finding the best pair with sum at most $$$M$$$.
Let's create new two dimensional arrays $$$P$$$ and $$$S$$$ of size $$$D(a_1) \times D(b_1)$$$. We will use $$$P$$$ to be able to tell the number of indexes $$$i$$$ such that we have either $$$X | a_i$$$ and $$$Y | b_i$$$ or $$$X | b_i$$$ and $$$Y | a_i$$$. If this count won't be $$$n$$$, then obviously we can't have $$$X$$$ and $$$Y$$$ as $$$\gcd$$$'s of $$$a$$$ and $$$b$$$. Also, we will use $$$S$$$ to tell us the sum of costs of all swaps we need to perform to have $$$X | \gcd(a_1, \ldots, a_n)$$$ and $$$Y | \gcd(b_1, ..., b_n)$$$.
Now how to make two such arrays efficiently? It is obvious that if the pair of $$$\gcd$$$ s $$$(X, Y)$$$ is consistent with some indexes in the original array, then for every pair $$$(x, y)$$$ such that $$$x|X$$$ and $$$y|Y$$$, this pair of $$$\gcd$$$s is also consistent with those indexes (and maybe even more, also maybe some swaps just became unnecessary, but the point is, it doesn't get worse). So if we want to add some value to a pair, we also want to get it added to all its divisors. That's why, in order to calculate the arrays efficiently, we will first add some values on some positions and then do something like $$$2D$$$ prefix sums — for every cell $$$(x, y)$$$ we will sum the values for all pairs $$$(X, Y)$$$ such that $$$x | X$$$ and $$$y | Y$$$ and update its current value with it. Assuming this is going to happen in the end, let's look at every $$$i$$$ and consider what pairs are "good" for this index — with or without the swap:
a) If $$$X$$$ divides $$$a_i$$$ and $$$Y$$$ divides $$$b_i$$$: For this type of pairs, we don't need to make any swaps on this index. Let's add $$$1$$$ to $$$P[\gcd(a_1, a_i)][\gcd(b_1, b_i)]$$$ to indicate that for all $$$(X, Y)$$$ such that $$$X|a_i$$$ and $$$Y|b_i$$$, we don't have to perform any swaps at the position $$$i$$$, the index is good as it is.
b) $$$X$$$ divides $$$b_i$$$ and $$$Y$$$ divides $$$a_i$$$: In this case we will add $$$1$$$ to $$$P[\gcd(a_1, b_i)][\gcd(b_1, a_i])]$$$ and $$$c_i$$$ to $$$S[\gcd(a_1, b_i)][\gcd(b_1, a_i)]$$$ to indicate that if we pick $$$X$$$, $$$Y$$$ such that $$$X|b_i$$$ and $$$Y|a_i$$$, we can make index $$$i$$$ good if we swap $$$a_i$$$ and $$$b_i$$$.
c) $$$X$$$ and $$$Y$$$ both divide both $$$a_i$$$ and $$$b_i$$$: To avoid counting both of the previous cases and therefore overcounting, we will add $$$-1$$$ to $$$P[\gcd(a_1, a_i, b_i)][\gcd(b_1, a_i, b_i)]$$$ and $$$-c_i$$$ to $$$S[\gcd(a_1, a_i, b_i)][\gcd(b_1, a_i, b_i)]$$$ (we have to undo paying for the swap since in this case we actually don't have to pay for it, but it falls under the case b) too).
This step can be done in $$$O(n \log k)$$$.
Now let's fix the arrays $$$P$$$ and $$$S$$$ so they store the actual values, not just the values we need to add.
We will go through all primes $$$p$$$ dividing $$$a_1$$$ and update $$$P[X][Y]$$$ with $$$P[p \cdot X][Y]$$$ and $$$S[X][Y]$$$ with $$$S[p \cdot X][Y]$$$, similarly for all primes dividing $$$b_1$$$. If we make those updates in the right order, we achieve that $$$P[x][y]$$$ is the sum of all original values $$$P[X][Y]$$$ for all the pairs $$$(X, Y)$$$ such that $$$x|X$$$ and $$$y|Y$$$ like we wanted (and we can do the same for $$$S$$$).
By careful precalculation of divisors and their order while factorizing, we can do this step in $$$O(D(k)^2 P(k))$$$. Some efficient implementations with extra log might pass also, but you will have to be more careful.
For multiple queries, after precalculating the possible sums of $$$\gcd$$$ s and their costs, you can sort them and use binary search to answer the queries.
Time complexity: $$$O(n \log k + \sqrt k + D(k)^2 P(k) + D(k)^2 \log D(k) + q \log D(k))$$$ Memory complexity: $$$O(n + D(k)^2 + P(k))$$$
Implementation in C++: 261999743
Fun fact: While making the original proposal, I invented a version of this problem with lower constraints and thought it could be nice Div2C. Few days later, while I was on a walk, I realized this solution exists, so we decided to try to propose it as the hardest problem.