Tonelli-Shanks algorithm in Python

I am trying to implement quadratic sieve factorization algorithm, but got stuck at finding modular squares with Tonelli-Shanks algorithm.

For example, if n = 87463, output of Tonelli-Shanks should be:

p = 2 -> x = 1  
p = 3 -> x = 1, 2  
p = 13 -> x = 5, 8  
p = 17 -> x = 7, 10  
p = 19 -> x = 5, 14  
p = 29 -> x = 12, 17  

But I get sometimes correct answers, and sometimes absolute rubbish:

tonelli_shanks(87463, 2)  
(87463, -87461)  
tonelli_shanks(87463, 3)  
(87463, -87460)  
tonelli_shanks(87463, 13)  
(8L, 5L)  
tonelli_shanks(87463, 17)  
(11, 6)  
tonelli_shanks(87463, 19)  
(5L, 14L)  
tonelli_shanks(87463, 29)  
(12L, 17L)   

I am following Tonelli-Algorithm description on Wikipedia but still cannot figure out the spot in the tonelli_shanks function where I messed up. Could someone help find the problem?

from math import *

def bound(num):
    o1 = 0.22
    return int(exp((0.5 + o1) * (log(num) * log(log(num))) ** 0.5))

def atkin_sieve(bound):
    primes = [2, 3, 5]
    sieve = [False] * (bound + 1)
    factor = int(sqrt(bound)) + 1
    for i in range(1, factor):
        for j in range(1, factor):
            n = 4 * i ** 2 + j ** 2
            if n <= bound and (n % 12 == 1 or n % 12 == 5):
                sieve[n] = not sieve[n]
            n = 3 * i ** 2 + j ** 2
            if n <= bound and n % 12 == 7:
                sieve[n] = not sieve[n]
            if i > j:
                n = 3 * i ** 2 - j ** 2
                if n <= bound and n % 12 == 11:
                    sieve[n] = not sieve[n]
    for i in range(5, factor):
        if sieve[i]:
            for j in range(i ** 2, bound, i ** 2):
                sieve[j] = False
    for i in range(7, bound):
        if sieve[i]: primes.append(i)
    return primes

def binary_jacobi(a, n):
    t = 1
    while a != 0:
        while a % 2 == 0:
            a >>= 1
            if n % 8 == 3 or n % 8 == 5: t = -t
        if a < n:
            a, n = n, a
            if a % 4 == 3 and n % 4 == 3: t = -t
        a = (a - n) >> 1
        if n % 8 == 3 or n % 8 == 5: t = -t
    if n == 1: return t
    else: return 0


def factor_base(primes, n):
    return [p for p in primes if binary_jacobi(n, p) == 1]

def stein_gcd(u, v):
    g = 1
    while u % 2 == 0 and v % 2 == 0:
        u >>= 1
        v >>= 1
        g <<= 1
    while u != 0:
        if u % 2 == 0: u >>= 1
        elif v % 2 == 0: v >>= 1
        else:
            t = abs((u - v) >> 2)
            if u >= v: u = t

            else: v = t
    return g * v

def modpower(a, k, n):
    """compute a ** k % n"""
    if k == 0: return 1
    if k == 1: return a % n
    q, r = divmod(k ,2)
    if r == 1:
        return modpower(a * a % n, q, n) * a % n
    else:
        return modpower(a * a % n, q, n)

def tonelli_shanks(n, p):
    h = p >> 1
    Q = p - 1
    S = 0
    while not Q & 1:
        Q >>= 1
        S += 1
    if S == 1:
        R = modpower(n, p + 1 >> 2, p)
        return R, p - R
    z = next(a for a in atkin_sieve(p) if binary_jacobi(p, a) == -1)
    c = modpower(z, Q, p)
    R = modpower(n, Q + 1 >> 1, p)
    if p - R > h: R -= p    
    t = modpower(n, Q, p)
    if p - t > h: t -= p
    M = S
    while t % p != 1:
        i = 0
        while t % p != 1:
            t **= 2
            i += 1
        b = modpower(c, 2 ** (M - i - 1), p)    
        bb = b * b % p
        if p - bb > h: bb -= p                    
        R = (R * b) % p        
        t = (t * bb) % p
        if p - t > h: t -= p        
        c = (bb) % p
        if p - c > h: c -= p        
        M = i        
    return R, p - R

Update 0: Apparently there an error in modpower function, instead of

if k == 0: return 1
if k == 1: return a

There should be:

if k == 0: return 1
if k == 1: return a % n

There is also an error in the stein_gcd method but it is not used yet; don't pay attention to it.

After modpower fix, tonelli_shanks almost works, but for

tonelli_shanks(87463, 17)

it still produces incorrect result - (11, 6) instead of (7, 10)

Update 1:

I found and error, which was in the place where I couldn't even think I might make a mistake:

        i = 0
        while t % p != 1:
            t **= 2
            i += 1

Here I search for the lowest i, s.t. t²ⁱ mod p = 1. But while searching for i, I modify t, which I should not do.

The correct version is with

        i = 0
        temp = t
        while temp % p != 1:
            temp **= 2
            i += 1