The algorithm *is* horrible ! From the documentation :

Warning

The current implementation splits the modulus into prime powers, then naively enumerates all possible solutions (starting modulo primes and then working up through prime powers), and finally combines the solution using the Chinese Remainder Theorem. The interface is good, but the algorithm is very inefficient if the modulus has some larger prime factors! Sage does have the ability to do something much faster in certain cases at least by using Groebner basis, linear algebra techniques, etc. But for a lot of toy problems this function as is might be useful. At least it establishes an interface.

And, in the ticket documenting its initial implementation, William Stein confesses its youthful sin :

```
WARNING:
Currently this naively enumerates all possible solutions.
The interface is good, but the algorithm is horrible if the
modulus is at all large! Sage *does* have the ability to do
something much faster in certain cases at least by using
the Chinese Remainder Theorem, Groebner basis, linear algebra
techniques, etc. But for a lot of toy problems this function
as is might be useful. At least it establishes an interface.
```

Indeed, trying to reproduce your problem on a moderatly fast laptop, I haven't got an answer in 20 minutes...

**EDIT :** in fact, this ended up in a crash :

```
sage: %time Sol = solve_mod([9677035*u + 7162589*v == 4477085 ,5302728*u + 8472081*v == 11061226], 14282531)
Killed
Process Sage exited abnormally with code 137
```

Since

```
sage: p = 14282531 ; p.is_prime()
True
```

a better solution is to work directly in the ring of polynomials of the finite field `GF(p)`

:

```
sage: R.<x, y> = GF(p)[]
sage: J = R.ideal(9677035*x + 7162589*y - 4477085 ,5302728*x + 8472081*y - 11061226)
sage: J.dimension()
0
```

Your problem has a finite number of solutions ; they are

```
sage: %time J.variety()
CPU times: user 38.7 ms, sys: 23.7 ms, total: 62.4 ms
Wall time: 87.9 ms
[{y: 3263554, x: 5754431}]
```

in fact a unique solution (obtained in a reasonable time...).

HTH,

**EDIT :**

Following @Max Alekseyev's hint, one can indeed use Pari/GP. The problem is to get Pari/GP's *column* vectors. Sage doesn't really have *row* and *column* vectors, just vectors (which, IMHO, is mathematically sounder that "row" and "column" vectors, which are a figment of programmer's imaginations...). In Sage, a `n`

-dimensional "column" vector is but a (`n`

, 1) matrix .. which is *not* accepted by Pari nor GP as a "column" vector.

The best I can come up with is the following horror :

```
sage: %time PSol= pari("matsolvemod(%s,%s~, %s~)"%(M._pari_init_(), vector([p]*2)._pari_init_(), B._pari_init_())) ; PSol
CPU times: user 244 µs, sys: 25 µs, total: 269 µs
Wall time: 273 µs
[5754431, 3263554]~
sage: PSol.sage()
[5754431, 3263554]
```

which is indeed faster (and more general) than the `GF(p)`

solution (nor all moduli ahave to be equal), and the unicity can be checked with :

```
sage: %time pari("matsolvemod(%s,%s~, %s~, flag=1)"%(M._pari_init_(), vector([p]*2)._pari_init_(), B._pari_init_()))
CPU times: user 271 µs, sys: 0 ns, total: 271 µs
Wall time: 276 µs
[[5754431, 3263554]~, [14282531, 0; 0, 14282531]]
```

Don't ask me what is the last part of this answer... ;-).

HTH,