import math
import functools

def factors(n):
    while n > 1:
        for i in range(2, int(n) + 1):
            if n % i == 0:
                n /= i
                yield i
                break


#https://oeis.org/search?q=3%2C+5%2C+9%2C+13%2C+21%2C+29%2C+45&language=english&go=Search
A060482 = [1, 2, 3, 5, 9, 13, 21, 29, 45, 61, 93, 125, 189, 253, 381, 509]
A136252 = [ 1, 3, 5, 9, 13, 21, 29, 45, 61, 93, 125, 189, 253, 381, 509, 765]

#https://oeis.org/search?q=5%2C+11%2C+17%2C+35%2C+53%2C+107&language=english&go=Search
A060647 = [1, 3, 5, 11, 17, 35, 53, 107, 161, 323, 485, 971, 1457, 2915, 4373, 8747, 13121, 26243, 39365]

#https://oeis.org/search?q=3%2C+9%2C+21%2C+45%2C+93&sort=&language=english&go=Search
A068156 = [1, 3, 9, 21, 45, 93, 189, 381, 765, 1533, 3069, 6141, 12285, 24573, 49149, 98301, 196605]

#https://oeis.org/A000225
A000225 = [0, 1, 3, 7, 15, 31, 63, 127, 255, 511, 1023, 2047, 4095, 8191, 16383, 32767, 65535, 131071]

#base
b = 3
#power of base
p = 5
#power of step
q = 9
#number of steps
n = b**p

#sum_{x=1..n}(x^q mod n - n/2)
s = 0
#sum_{x=1..n}(x^q mod n)/n
s2 = 0

for x in range(1,n+1):
    s += ((x**q) % n - n/2)
    s2 += (x**q) % n

# sum s
print("s=sum_{x=1..n}(x^q mod n - n/2):")
print(s)

#coplement
coplement = s*(-2)
print("s*(-2):")
print(coplement)

print("sum_{x=1..n}(x^q mod n)/n:")
print(s2/n)
check = 1/2*(n-math.sqrt(n))
print("1/2*(n-sqrt(n)):")
print(check)

#check 2^p, q=2
'''
print("1/2*(n-A060482[p])")
check = 1/2*(n-A060482[p])
print(check)
print("1/2*(n-A136252[p-1])")
check = 1/2*(n-A136252[p-1])
print(check)
'''

#check 3^p, q=2
'''
print("1/2*(n-A060647[p+1]/3)")
check = 1/2*(n-A060647[p+1]/3)
print(check)
'''

#check 4^p, q=2
'''
print("1/2*(n-A068156[p])")
check = 1/2*(n-A068156[p])
print(check)
'''
#factorize complement
compl_factors_list = list(factors(coplement))
print("factors of s*(-2):")
print(compl_factors_list)
print(len(compl_factors_list))

#product of non base remainder
non_base_factors = [x for x in compl_factors_list if x != b]
if len(non_base_factors) > 0:
    remainder = functools.reduce((lambda x, y: x * y), non_base_factors)
    print("non base remainder:")
    print(remainder)