This reduces to the so called modular square root problem. That is, solve for x, such that
mod(x^2,p) == r
where p and r are given values. The solution can be gained from the Shanks-Tonelli algorithm, here:
A problem is this is more difficult if p is composite.
But here, p is indeed prime, so a solution may exist. That is not always true, as you seem to understand. Your problem is to solve for b, such that:
b^2 = (9361)^3 + 23698* 9361 + 9684 (mod 36683)
First, expand the right hand side, as
rhs = mod((9361)^3 + 23698* 9361 + 9684,p)
Note that p is a prime of the form 4*k+3.
In that case, a simple, direct solution exists as:
x = powermod(sym(rhs),(p+1)/4,p)
Is this a valid solution?
mod(x^2,p) == rhs
ans = 
Yes. Only in the case where p is a prime of the form 4*k+1 do you need to dive into Shanks-Tonelli. As I recall, it is not that difficult to write. I've got the code somewhere, but it is implemented using my VPI toolbox. I might have written a version that works under sym. Not sure about that. Anyway, this is irrelevant, as long as p is a prime of the form 4*k+3.
Again, if p is composite, things get messier yet.
