The radical function of integers: finding a reference and proof of a relation We can write each positive integer greater than 1 as the product of its prime factors: $n=\prod{p_i^{e_i}}$. The radical function is defined as $\operatorname{rad}(n)=\prod_i{p_i}$. For instance $\operatorname{rad}(20)=2 \cdot 5=10$. It is convenient to define $\operatorname{rad}(0)=0$ and $\operatorname{rad}(1)=1$. 
For non-zero $n$ this implies that $\operatorname{rad}(n)$ divides $n$, and $\frac{n}{\operatorname{rad}(n)}=\prod_{p \mid n}{p_i^{e_i-1}}$.
We know several properties and relations of the radical function, such as:

*

*For squarefree integers $n$ we have $\operatorname{rad}(n)=n$. Therefore $\operatorname{rad}(n)=0,1$ implies $n=0,1$.

*$\operatorname{rad}(mn) \le \operatorname{rad}(m)\operatorname{rad}(n)$, and only if $m,n$ are coprime then $\operatorname{rad}(mn)=\operatorname{rad}(m)\operatorname{rad}(n)$.

*$\operatorname{rad}(m) \mid n$ implies $\operatorname{rad}(mn)=\operatorname{rad}(n)$.

*$n$ divides $\operatorname{rad}(n)^{|n|}$.

*$\operatorname{rad}(n) = \sum_{d|n}{\phi(d) \cdot \mu(d)^2}$ with $\phi$ the Euler totient function.

*$\operatorname{rad}(n)$ is the period of the sequence $(k^n \mod{n})_{k=1}^{k=\infty}$.

*Using the floor function we have $\frac{n}{\operatorname{rad}(n)} = \sum_{k=1}^{k=n}{\lfloor \frac{k^n}{n} \rfloor-\lfloor \frac{k^n-1}{n} \rfloor}$.

*Using the Dirichlet series we have for $0<t$ and $t+1<s$ that $\sum_{n=1}^{\infty}\frac{\operatorname{rad}(n)^t}{n^s}=\prod_{p,prime}{(1+\frac{p^t}{p^s-1})}$.

*Let $a,b$ be coprime integers with $1<a,b$ and $r_b^2 \nmid b$. Let $s,t$ such that $sr_a^2+tr_b^2=1$. Let $(r)_x=r \pmod{x}$. Then
$b+(r_ar_b)^2=(bs)_{(r_b)^2} \cdot r_a^2+(bt)_{(r_a)^2} \cdot r_b^2$ where $sr_a^2+tr_b^2=1$. (Fermat)

In particular, I am looking for a reference and proof of 9.
Thanks for any help.
 A: Proposition 9 is incorrect. Let $a = 5$ and $b = 3 \times 2^9 = 1536$. Then $\mathrm{gcd}(a, b) = 1$, $r_b^2 \nmid b$ and one can compute $(r_a, r_b, s, t, (bs)_{r_b^2}, (bt)_{r_a^2}) = (5, 6, -23, 16, 24, 1)$. But we have:
$$b + (r_ar_b)^2 = 2436 \not= 636 = (bs)_{r_b^2}.r_a^2 + (bt)_{r_a^2}.r_b^2$$
Is this the only counterexample that exists? How could I reach it?
I explain my unsuccessful approach to prove the statement. Using the division algorithm, we can write
\begin{align}
bs = kr_b^2 + (bs)_{r_b^2} && (1)\\
bt = qr_a^2 + (bt)_{r_a^2} && (2)
\end{align}
By multiplying the equation $(1)$ by $r_a^2$ and the other by $r_b^2$ and then adding them up, we get
\begin{align}
b(sr_a^2 + tr_b^2) &= (k + q)(r_ar_b)^2 + (bs)_{r_b^2}.r_a^2 + (bt)_{r_a^2}.r_b^2\\
\stackrel{sr_a^2 + tr_b^2 = 1}{\implies} b &= (k + q)(r_ar_b)^2 + (bs)_{r_b^2}.r_a^2 + (bt)_{r_a^2}.r_b^2
\end{align}
So if the statement is true, we must have $k + q = -1$. Notice that $k = \lfloor \frac{bs}{r_b^2} \rfloor$ and $q = \lfloor \frac{bt}{r_a^2} \rfloor$. If we multiply the equation $sr_a^2 + tr_b^2 = 1$ by $\frac{b}{(r_ar_b)^2}$, we get
\begin{align}
\frac{bs}{r_b^2} + \frac{bt}{r_a^2} = \frac{b}{(r_ar_b)^2}.
\end{align}
The integer part of the LHS approximately equals $k + q$ (the |difference| is at most $1$) but we can choose $(a, b)$ such that the RHS go to infinity. The statement seems true (in the first examples) because for small $\max\{|a|, |b|\}$ we often have $0 < $ RHS $< 1$ and $k + q = -1$ but the issue appears when we select bigger numbers.
For calculations and finding other counterexamples, I wrote some dirty code in online-python that you can use.
def isPrime(n):
    if n == 1:
        return 0
    elif n == 2:
        return 1
    else:
        for i in range(2, n//2+1):
            if n % i == 0:
                return 0
    return 1

def rad(n):
    r = 1
    if n == 1:
        return 1
    else:
        for i in range(2, n+1):
            if n % i == 0 and isPrime(i) == 1:
                r = r*i
    return r

def sandt(a, b):
    ra2 = rad(a)**2
    rb2 = rad(b)**2
    t = 1
    s = 0
    while s == 0:
        if (1 - t*rb2)%(ra2) == 0:
            s = (1 - t*rb2)//(ra2)
        else:
            t = t + 1
    return (s, t)
    
def isCorrect(a, b):
    ra2 = rad(a)**2
    rb2 = rad(b)**2
    s = sandt(a, b)[0]
    t = sandt(a, b)[1]
    print(str(b) + "*" + str(s) + " = " + str((b*s)//rb2) + "*" + str(rb2) + " + " + str((b*s)%rb2))
    print(str(b) + "*" + str(t) + " = " + str((b*t)//ra2) + "*" + str(ra2) + " + " + str((b*t)%ra2))
    return b + ra2*rb2 == ((b*s)%rb2)*ra2 + ((b*t)%ra2)*rb2
    
print(isCorrect(5, 1536))

