AGARTHA IS REAL (Editoryal)

Solution

(We'll be referring to the current value of $n$ after all the operations as $\text{cur}$)

$\text{cur}$ can't become negative at any time, because then both M-operations and S-operations make it go even further in the negatives so it won't ever reach $0$. That means in order to make $\text{cur}$ equal to $0$, we first have to make it a multiple of $z$ (\text{cur}\mod z = 0) and then keep performing S-operations until it's $0$.

Let's first deal with the $-1$ cases. As we said earlier, $\text{cur}$ can't be negative, so if $n$ is initially negative then the answer is $-1$.

The other $-1$ case is if we can't make make $\text{cur}$ a multiple of $z$. It's enough to check the value (n \cdot x^z \cdot y^z)\mod z, if it's not $0$ then we can't make $n$ a multiple of $z$.

Proof: Since there are only $z$ different cur\mod z values, multiplying by more than $z$ times will end us up in a cycle.

Now that we're done with the $-1$ cases, we will always have an answer.

We can clearly see that in the optimal solution we won't ever do an M-operation unless \text{cur} < z is true, because if \text{cur} \ge z we can keep doing S-operations until \text{cur} < z, that's because \text{cur}\mod z = (\text{cur}-z)\mod z, and since not doing the S-operations before the M-operation will require us to do more S-operations in the future which is clearly not optimal.

Let's solve the case where x=y:

Claim: We won't do more than $\log (z)$ M-operations.

Proof: For $\text{cur}$ to be a multiple of $z$, $\text{cur}$ must have all the prime factors of $z$. When we do an M-operation we're adding prime factors to $\text{cur}$, and the maximum difference between the number of occurrences of a prime factor in $z$ and the number of occurrences of the same prime factor in $\text{cur}$ is $\log (z)$, so we won't need more than $\log (z)$ M-operations to make $\text{cur}$ have all the prime factors of $z$.

That means for x=y we can just keep doing M-operations while \text{cur} < z and if $\text{cur}$ ever becomes larger than $z$ we can make $\text{cur}$ equal to \text{cur}\mod z with a cost of \left \lfloor{\frac{\text{cur}}{z}} \right \rfloor (by repeatedly doing S-operations), We stop when \text{cur}=0.

Now to solve the full problem, Let's first define the state ($i,j$) as the state where we've multiplied $i$ times by $x$ and $j$ times by $y$ and $\text{cur}$ must equal to (n \cdot x^i \cdot y^j)\mod z (we'll refer to the value of $\text{cur}$ in this state as \text{cur}[i][j]).

We can use a dynamic programming approach:

Base case:

Transition:

The answer is the minimum \text{dp}[i][j] value where \text{cur}[i][j] = 0.

Now since we won't make more than $z$ M-operations, the complexity for the dp is \mathcal{O}(z^2), which is way too large.

As we said in the x=y solution, We don't need more than $\log (z)$ M-operations to make $\text{cur}$ a multiple of $z$. That means we can always get an answer less than 10^{11}. Since the maximum answer we could get if we do the M-operations arbitrary is \log(z) \cdot \left \lfloor{\frac{\max(x,y) \cdot (z-1)}{z}}\right \rfloor + \log(z) \cdot (\log(z)+1) which is strictly less than 10^{11}.

Since the answer is always less than 10^{11}, We won't do more than 5\cdot10^5 M-operations because the cost for them alone would be larger than 10^{11}.

That optimizes the dp solution to \mathcal{O}(10^{11}) by making 5 \cdot 10^5 the limit for $i,j$, but that's still too slow.

Let's observe the state ($i,j$) where both $i$ and $j$ are strictly larger than $\log (z)$, it's guaranteed that \text{cur}[i-1][j]=\text{cur}[i][j-1]=0, and since \min(\text{dp}[i-1][j],\text{dp}[i][j-1]) \le \text{dp}[i][j], we don't care about the state ($i,j$) because we already have a better answer.

Now we can run the previous dp two times, the first time with i \le \log(z) , j \le 5 \cdot 10^5, and the second time with i \le 5 \cdot 10^5 , j \le \log(z).

So the total complexity of the code is \mathcal{O}(5 \cdot 10^5 \cdot \log(z)) which is fast enough to pass.

#include<bits/stdc++.h>
#define ll long long
#define d3 ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
#define cans cout << ans << "\n";
 
using namespace std;

const ll N = 500000 , lg = 35;

ll n , x , y , z , ans = LLONG_MAX;
ll dp[lg][N];
ll md[2][N];

ll fp(ll a , ll b)
{
    if(b==0)
        return 1;
    ll c = fp(a,b/2);
    if(b%2)
        return c*c%z*a%z;
    else
        return c*c%z;
}

ll Mcost(ll a)
{
    return a*(a+1)/2;
}

void rec(ll a , ll b)
{
    dp[0][0]=n/z;
    md[0][0]=md[1][0]=1;
    for(int i = 1 ; N>i ; i++)
    {
        md[0][i]=md[0][i-1]*a%z;
        md[1][i]=md[1][i-1]*b%z;
    }
    for(int i = 1 ; lg>i ; i++)
    {
        ll p = n*md[0][i-1]%z;
        p*=a;
        dp[i][0]=dp[i-1][0]+p/z;
    }
    for(int i = 1 ; N>i ; i++)
    {
        ll p = n*md[1][i-1]%z;
        p*=b;
        dp[0][i]=dp[0][i-1]+p/z;
    }
    for(int i = 1 ; lg>i ; i++)
    {
        for(int j = 1 ; N>j ; j++)
        {
            ll p = n*md[0][i-1]%z*md[1][j]%z;
            p*=a;
            dp[i][j]=dp[i-1][j]+p/z;
            p = n*md[0][i]%z*md[1][j-1]%z;
            p*=b;
            dp[i][j]=min(dp[i][j],dp[i][j-1]+p/z);
        }
    }
    for(int i = 0 ; lg>i ; i++)
    {
        for(int j = 0 ; N>j ; j++)
        {
            if(n*md[0][i]%z*md[1][j]%z==0)
                ans=min(ans,dp[i][j]+Mcost(i+j));
        }
    }
}

int main()
{
    d3
    cin >> n >> x >> y >> z;
    ll e = n*fp(x,100)%z*fp(y,100)%z;
    if(n<0||e!=0)
    {
        cout << -1;
        return 0;
    }
    rec(x,y);
    rec(y,x);
    cans
}