solve_mod() takes long time and no result return when computing big numbers

edit

Solution via linear algebra modulo 14282531:

M = Matrix( Zmod(14282531), [[9677035, 7162589], [5302728, 8472081]])
print( M.solve_right( vector([4477085, 11061226]) ) )

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,