Now the task is to find and prove a proper upper bound (inequality) for the score. The score might be dependent on the initial number, thus could be used for "branch cutting" when searching for the highest possible score. Or a simple constant is also acceptable.
We formulate the process as following: Let $n_0$ be the initial number, $n_i$ is the number obtained after $i$th operation. Each operation corresponds to a prime number $p_i$ (possibly negative) where $|p_i| \neq |p_j|, \forall i \neq j$. Easy to see that $|p_i p_{i+1}| \leq \max(n_i, n_{i+1})$. Then (notice the index, $p_i := n_i - n_{i-1}$) $$ \begin{align} n_{i+1} &= n_i + p_{i+1} \leq n_i + |p_{i+1}| \leq \max(n_i + n_{i+1}/|p_i|, n_i + n_i/|p_i|) \ &\text{Because } (|p_i|+1)/|p_i| \leq |p_i|/(|p_i|-1), \text{ we always have:} \ n_{i+1} &\leq \frac{|p_i|}{|p_i|-1}n_i \leq \left(\prod_{k=0}^i \frac{|p_k|}{|p_k|-1}\right)n_0 \end{align} $$
Notice $x/(x-1)$ is a decreasing function for $x \geq 1$, and most importantly, $|p_i| \leq n_{i+1}$ (This one needs proof, but loose. If we disable subtraction, then this is plainly true.).
Now, instead of bound the final number using initial number, we bound the initial number (lower bound) using the final number:
$$ \frac{y}{\prod_{q_i < y, q_i \text{ is prime}} \frac{q_i}{q_i-1}} = \left(\prod_{q_i < y, q_i \text{ is prime}} \left(1-\frac{1}{q_i}\right)\right)y \leq x $$
Where $y$ is the final number and $x$ is the initial. Notice we are not using the $p_i$ (operations, which is unknown) here. We are using all the prime numbers below $y$ in the product. We could optimize this result by optimizing $\alpha$ in $|p_i| \leq \alpha n_{i+1}$ which is set as 1 above.
Using Merten's third theorem, for large $y$, $$ \prod_{p \leq y} \left(1 - \frac{1}{p}\right) \sim \frac{e^{-\gamma}}{\ln y} $$ where $\gamma$ is the Euler-Mascheroni constant.
Thus a lower bound (very loosely) for initial number $x$ is approximately $\displaystyle \frac{e^{-\gamma}y}{\ln y}$.

import matplotlib.pyplot as plt
import numpy as np
from math import isqrt
from sympy import sieve
import math
def product_factor(y):
"""
Compute the product of p/(p-1) for all primes p <= y.
Utilizes sympy's sieve for efficient prime enumeration.
"""
if y < 2:
return 1
primes = list(sieve.primerange(2, y + 1))
product = 1.0
for p in primes:
product *= p / (p - 1)
return product
def smallest_initial(y):
"""
Compute the smallest initial number x such that
x >= y / (product of p/(p-1) for all primes p <= y)
"""
if y == 0:
return 0
prod = product_factor(y)
return y / prod
def visualize():
x = np.arange(10, 1001) # Final numbers from 10 to 1000
# Calculate y values: smallest initial x needed for each y
y = [smallest_initial(i) for i in x]
# Create the plot
plt.figure(figsize=(12, 6))
plt.plot(x, y, label=r'$x \geq \frac{y}{\prod_{p \leq y} \frac{p}{p-1}}$')
plt.title('Lower Bound of Initial Number vs Final Number')
plt.xlabel('Final Number (y)')
plt.ylabel('Minimum Initial Number (x)')
plt.grid(True)
plt.legend()
plt.tight_layout()
# Show the plot
plt.show()
if __name__ == "__main__":
visualize()