I input
Solve[35345 x^47 + 1345326 x^2 + 134636246 x + 13451345 == 0, x, Modulus ->2037035976334486086268445688409378161051468393665936250636140449354381299763336706183399101]
to Mathematica and the output is
x->
341367685451396522792304545319361914942948733456457424375640931820820238770112149829924341
I wonder how Mathematica gets this result? Are there special algorithms to solve this problem?
On
As I did a little bit of homework on Cantor-Zassenhaus algorithm I will write up a toy example that can be done with paper and pencil (if you have enough patience and paper). Explaining bits and pieces while I go.
Task: Find the zeros of $f(x)=x^4+13x^2+6x+5$ in the field $F=\Bbb{F}_{17}$.
The first step is to calculate the greatest common divisor $\gcd(f(x),x^{17}-x)$ in the ring $F[x]$. This is fast enough, and it turns out (I set it up that way) that $f(x)\mid x^{17}-x$. This has the following two important consequences.
Consequently $f(x)=(x-a_1)(x-a_2)(x-a_3)(x-a_4)$ for some pairwise distinct $a_1,a_2,a_3,a_4\in F$. We could easily find them by trial and error, but let's not. I want to describe the C-Z idea that also works for primes $q$ much larger than $p=17$. At this point we know (from Chinese remainder theorem) that $$ \begin{aligned} R=F[x]/\langle f(x)\rangle&\cong F[x]/\langle x-a_1\rangle\oplus\cdots F[x]/\langle x-a_4\rangle\\ &\cong F\oplus F\oplus F\oplus F. \end{aligned} $$ The idea is that if we take a random non-zero element $r\in R$, we can easily calculate $r^{(p-1)/2}-1=r^8-1$. If in the above isomorphism the image of $r$ is $(r_1,r_2,r_3,r_4)\in F^4$ is such that some of the (unknown) components $r_i$, $i=1,2,3,4,$ are non-zero quadratic residues, and some are not, then some of the components of $r^8-1,$ i.e.$r_i^8-1, i=1,2,3,4, $ are non-zero and some are not. This happens with probability at least $\approx 1/2$ (for larger primes). When this happens it means that some but not all of the factors $(x-a_i)$ are also factors $r^{(p-1)/2}-1$. Calculating the greatest common divisor of $f(x)$ and the lowest degree polynomial in the coset of $r^8-1$ modulo $f(x)$ then reveals a non-trivial factor. Rinse. Repeat.
In our example case let us try $r=x$ for simplicity. Then it turns out that $$ x^8-1\equiv 12+15x+9x^2+3x^3\pmod{f(x)}. $$ We see (Euclid's algorithm)) that with the remainder $3f_1, f_1(x)=x^3+3x^2+5x+4$, we actually have $f_1(x)\mid f(x)$. A long division shows that $$ f(x)/f_1(x)=x+14. $$ This shows that $a_4=3$ is one of the zeros. This step was already explained in Jack's excellent answer. We can conclude that $3$ is the only zero that is NOT a quadratic residue modulo $17$. After all, if $a$ is a square in $F$, then $x-a\mid x^8-1$.
But we are not done, so let's use another $r$. Let's try $r=x+1+\langle f(x)\rangle$. This time $r^8-1\equiv 6+3x+7x^2+16x^3\pmod{f(x)}$. Another run of Euclid's algorithm gives $$ f_2(x)=\gcd(6+3x+7x^2+16x^3,f(x))=7+6x+x^2. $$ This is a quadratic factor, so our luck is in. Yet another run of Euclid's algorithm shows that $$ \gcd(f_2(x),f_1(x))=x+9, $$ so we see that $a_3=8$ is one of the zeros. By polynomial division (or using Vieta relations) we see that $a_4=3$ is the other zero of $f_2(x)$, so two zeros remain unknown.
We could keep going, but it is worth pointing out that it is better (=computationally simpler) to start working with $$ f_3(x)=f(x)/f_2(x)=8+11x+x^2 $$ instead of $f(x)$.
Let's pick $r=(x+2)$. A calculation shows that $$ (x+2)^8-1\equiv 2+16x=2-x\pmod{f_3(x)}. $$ As $\gcd(x-2,f_3(x))=x-2$, this reveals the zero $a_2=2$. Vieta relations then reveal the final zero $a_1=4$.
Observe that powers in rings such as $R$ can be computed fast with the aid of square-and-multiply. The run time of Euclid's algorithm depends on the degree of the polynomials. Here only in the first step of calculating $\gcd(f(x),x^p-x)$ is there a high degree polynomial. But because it is a binomial, its remainder can be computed (as pointed out by Thomas Andrews) very quickly, again by square-and-multiply.
The non-deterministic element here is in the choice of $r$. If at first you don't succeed, pick another. The factors I used in the above toy example don't look random, but I swear :-) I didn't design the example so that they would work (but in retrospect I could have - I only understood this afterwards).
The factorization of a polynomial over a finite field is a well-understood problem.
The key idea is that $$ x^p-x = \prod_{\xi\in\mathbb{F}_p}(x-\xi),$$ hence the greatest common divisor between $x^p-x$ and $f(x)$ gives all the roots of $f$ in $\mathbb{F}_p$, and there are many efficient methods to compute the $\gcd$ of two polynomials, the classical ones relying on the theory of (sub-)resultants. The complexity essentially depends on a fixed power of the degree of $f$ and on a fixed power of $\log(p)$, hence the problem is very affordable even if the field has a huge characteristic.
Addendum. To recover the roots from $g(x)=\gcd(x^p-x,f(x))$ you can for first get rid of multiple roots by considering $h(x)=g(x)/\gcd(g(x),g'(x))$.
Then you can use a "divide-and-conquer" strategy base on the fact that $x^{\frac{p+1}{2}}-x$ is the product of all the linear factors $(x-\xi)$ with $\xi$ being a quadratic residue, hence you can "split" $h(x)$ into two polynomials, the first one having quadratic residues as roots, the second one having non-quadratic residues as roots. If you know how $p-1$ factors, you can continue this way by separating (now I am assuming that $3\mid(p-1)$) the roots in any factor according to the value of $$\xi^{\frac{p-1}{3}}\in\{1,\omega,\omega^2\}\subset\mathbb{F}_p,$$ till you get the complete list of linear factors.
The probabilistic Cantor-Zassenhaus algorithm do essentially the same without the knowledge of the factorization of $(p-1)$.