numbers theory

#!/usr/bin/python
# -*- coding: utf-8 -*-
#
#    This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>.

from collections import Counter
from sys import stdin, stdout, stderr
from math import sqrt

def is_prime(n):
    return not any(map(lambda i: not n % i, range(2, n)))

def gcd(a, b, inform = True):
    if inform:
        stdout.write('gcd(%d, %d) = ' % (a, b, ))
    while b:
        if inform:
            stdout.write('gcd(%d, %d) = ' % (b, a % b, ))
        a, b = b, a % b
    if inform:
        print a
    return a

def bit_form(n):
    if not n:
        return []
    if n % 2:
        return [1] + map(lambda x: x*2, bit_form(n / 2))
    else:
        return map(lambda x: x*2, bit_form(n / 2))

def factorization(n):
    _n = n
    factors = Counter()
    i = 2
    while n > 1 and i*i <= n:
        if n % i:
            i += 1
        else:
            n /= i
            factors.update([i])
    if n > 1:
        factors.update([n])
    print '%d = %s' % (_n, ' * '.join(map(lambda factor: '%d^%d' % (factor, factors[factor], ), factors)))
    return factors

def fast_pow_mod(a, b, p):
    print 'Считаем %d^%d по модулю %d' % (a, b, p, )
    deg = 2
    prev = a
    while deg <= b:
        print '%d^%d ≡ %d*%d ≡ %d ≡ %d' % (a, deg, prev, prev, prev*prev, prev*prev % p, )
        prev = prev*prev % p
        deg *= 2
    bits = bit_form(b)
    print '%d = %s' % (b, '+'.join(map(str, bits)))
    stdout.write('%d^%d = %s' % (a, b, ' * '.join(map(lambda b: '%d^%d' % (a, b), bits))))
    cur = 1
    while bits:
        cur *= a**bits[0] % p
        del bits[0]
        stdout.write(' = %s' % (' * '.join([str(cur)] + map(lambda b: '%d^%d' % (a, b), bits))))
        if cur != cur % p:
            cur %= p
            stdout.write(' = %s' % (' * '.join([str(cur)] + map(lambda b: '%d^%d' % (a, b), bits))))
    print
    return cur

def phi(n):
    print 'Ищем φ(%d)' % (n, )
    factors = factorization(n)
    res = 1
    for f in factors:
        c = f**(factors[f]-1) * (f - 1)
        res *= c
        print 'φ(%d^%d) = %d' % (f, factors[f], c, )
    print 'φ(%d) = %d' % (n, res, )
    return res

def primitive_root(p, k):
    print 'Ищем примитивный корень в Z_{%d^%d}' % (p, k, )
    def primitive_root_p(p):
        print 'Ищем примитивный корень в Z_%d' % (p, )
        print 'φ(%d) = %d' % (p, p-1, )
        factors = factorization(p-1)
        print 'Простые делители %d: %s' % (p-1, ', '.join(map(str, factors)), )
        check_degs = map(lambda n: (p-1) / n, factors)
        print 'Надо проверять степени %s' % (', '.join(map(str, check_degs)), )
        for candidate in range(2, p):
            if any(map(lambda deg: candidate**deg % p == 1, check_degs)):
                continue
            print 'Проверяем %d' % (candidate, )
            for deg in check_degs:
                if fast_pow_mod(candidate, deg, p) == 1:
                    continue
            print '%d - примитивный корень по модулю %d' % (candidate, p, )
            return candidate

    p_pr = primitive_root_p(p)
    if k == 1:
        return p_pr
    print 'Проверим, является ли %d примитивным корнем по модулю %d^2' % (p_pr, p, )
    if fast_pow_mod(p_pr, p-1, p*p) != 1:
        print '%d - примитивный корень mod %d^2 (а значит, и %d^%d)' % (p_pr, p, p, k, )
        return p_pr
    else:
        print '%d не подошел; значит, %d -  примитивный корень mod %d^2 (а значит, и %d^%d)' % (p_pr, p_pr+p, p, p, k, )
        return p_pr + p

def solve_congruence_modulo(k1, m1, k2, m2):
    print 'Решаем сравнение: x ≡ %d (%d), x ≡ %d (%d)' % (k1, m1, k2, m2, )
    cur = k1
    while cur % m2 != k2:
        #print '%d mod %d == %d; не подходит' % (cur, m2, cur % m2, )
        cur += m1
    print '%d + %d*%d = %d ≡ %d (%d); подходит' % (k1, m1, (cur - k1)/m1, cur, cur % m2, m2, )
    return cur

def group_decomposition(m):
    print 'Разложение группы Z_{%d}^*' % (m, )
    factors = factorization(m)
    print
    g = {}
    r = {}
    if 2 in factors:
        if factors[2] == 1:
            r[-1] = 1
            r[0] = 1
        else:
            r[-1] = 2
            r[0] = 2**(factors[2] - 2)
        print 'Ищем примитивные корни для 2^%d' % (factors[2], )
        g[-1] = solve_congruence_modulo(-1, 2**factors[2], 1, m / 2**factors[2])
        g[0] = solve_congruence_modulo(5, 2**factors[2], 1, m / 2**factors[2])
        print
    i = 0
    for factor in factors:
        if factor != 2:
            i += 1
            f = factor**factors[factor]
            r[i] = phi(f)
            g[i] = solve_congruence_modulo(primitive_root(factor, factors[factor]), f, 1, m/f)
            print
    print ('Z_{%d}^* = ' % (m, )) + ' * '.join(map(lambda i: '<%d>_{%d}' % (g[i], r[i], ), sorted(g.keys())))

def sqr_solves_count(a, b, c, m):
    print 'Проверяем, сколько решений имеет сравнение %d*x^2 + %d*x + %d ≡ 0 (mod %d)' % (a, b, c, m, )
    if not m % 2 and b % 2:
        print 'b делится на 2, m - нет. Не знаю, что делать'
        return None
    factors = factorization(m)
    d = (b*b - 4*a*c)# % m
    print 'Дискриминант равен %d' % (d, )
    for factor in factors:
        if factors[factor] >= 2 and not d % factor:
            print 'Дискриминант делится на %d, которое входит в %d в степени выше первой. Не знаю, что делать' % (factor, m, )
            return None
    for factor in factors:
        if factor != 2:
            if (d**((factor-1)/2) + 1) % factor == 0:
                p = fast_pow_mod(d, (factor - 1) / 2, factor)
                print '%d - не квадратичный вычет по модулю %d, решений нет' % (d, factor, )
                return 0
    res = 1
    if 2 in factors:
        if factors[2] > 2:
            if d % 8 == 1:
                print 'По модулю %d есть 4 решения' % (2**factors[2], )
                res *= 4
            else:
                print 'А вот нету решений по модулю %d, пичалька:(' % (2**factors[2], )
                return 0
        else:
            print 'Хм, модуль делится на 2, но не делится на 8. ХЗ что делать'
            return None
    for factor in factors:
        if factor != 2:
            if not d % factor:
                print 'Дискриминант делится на %d, которое входит в m в первой степени. Число решений не меняется' % (factor, )
            else:
                print 'Проверяем количество решений по модулю %d' % (factor**factors[factor], )
                fast_pow_mod(d, phi(factor - 1) / 2, factor)
                print 'Астрологи объявили %d квадратичным вычетом по модулю %d. Количество решений удваивается.' % (d, factor, )
                res *= 2
    print 'Итого решений: %d' % (res, )
    return res

def str_polynomial(coef, letter):
    coef = coef[::-1]
    letter = str(letter)
    return ' + '.join(reversed(map(lambda deg: ('%d*%s^%d' % (coef[deg], letter, deg)), xrange(len(coef)))))

def polynomial_solves_number(coef, m):
    def f(x):
        res = 0
        for c in coef:
            res = res*x + c
        return res
    def f_p(x):
        res = 0
        _coef = coef[::-1]
        return sum(map(lambda n: n * x**(n-1) * _coef[n], xrange(1, len(coef))))
    print 'Ищем число решений %s = 0 (mod %d)' % (str_polynomial(coef, 'x'), m, )
    factors = factorization(m)
    solutions_diff_modules = {}
    for factor in factors:
        solutions = []
        print 'Считаем число решений по модулю %d' % (factor**factors[factor], )
        print 'Ищем решения по модулю %d (переборчик) и пытаемся их поднять' % (factor, )
        print
        for i in xrange(factor):
            if f(i) % factor:
                print 'f(%d) = %d, не решение' % (i, f(i), )
            else:
                print 'f(%d) = %d, решение' % (i, f(i), )
                c_f_p = f_p(i)
                print 'f\'(%d) = %d' % (i, c_f_p, )
                if c_f_p % factor != 0:
                    print 'Это решение поднимается единственным образом'
                    cur_sol = i
                    for cur_deg in xrange(2, factors[factor]+1):
                        for t in xrange(0, p):
                            if (f(cur_sol)/(factor**cur_deg) + t*c_f_p) % factor == 0:
                                print 'f(%d) / %d^%d + f\'(%d)*%d = 0 (p), f(%d) = 0 (p^%d+1)' % (cur_sol, factor, cur_deg, cur_sol, t, cur_sol + t*factor**cur_deg, )
                                solutions += 1
                else:
                    print 'Производная равна нулю... придется париться с подъемом:('
                    cur_candidats = [i]
                    cur_deg = 1
                    while cur_deg < factors[factor]:
                        print 'Пытаемся поднять решения со степени %d до %d' % (cur_deg, cur_deg + 1, )
                        next_candidats = []
                        for candidat in cur_candidats:
                            print 'Пробуем поднять %d: f(%d) = %d' % (candidat, candidat, f(candidat), )
                            if not f(candidat) % (factor ** (cur_deg + 1)):
                                print 'Подходит, значит, добавляем кучу решений'
                                for t in xrange(0, factor):
                                    next_candidats.append(candidat + t * factor**cur_deg)
                            else:
                                print 'Не подходит. Дальше не поднимаем.'
                        cur_deg += 1
                        cur_candidats = next_candidats
                        print 'Решения mod %d^%d: %s', (factor, cur_deg, ', '.join(map(str, next_candidats)), )
                print
        print

class RationalNumber:
    def __reduce(self):
        if self.nominator == 0:
            return self
        g = gcd(self.nominator, self.denominator, False)
        if g < 0:
            g *= -1
        self.nominator /= g
        self.denominator /= g
        return self

    def __init__(self, nominator, denominator):
        self.nominator = nominator
        self.denominator = denominator
        self.__reduce()

    def __neg__(self):
        return self.__class__(-self.nominator, self.denominator)
    def __pos__(self):
        return self
    def __abs__(self):
        return self.__class__(abs(self.nominator), self.denominator)
    def __float__(self):
        return float(self.nominator) / self.denominator
    def __str__(self):
        return '(%d/%d)' % (self.nominator, self.denominator, )
    def __add__(self, other):
        return self.__class__(self.nominator * other.denominator + self.denominator * other.nominator, self.denominator * other.denominator).__reduce()
    def __sub__(self, other):
        return self.__class__(self.nominator * other.denominator - self.denominator * other.nominator, self.denominator * other.denominator).__reduce()
    def __mul__(self, other):
        return self.__class__(self.nominator * other.nominator, self.denominator * other.denominator).__reduce()
    def __div__(self, other):
        return self.__class__(self.nominator * other.denominator, self.denominator * other.nominator).__reduce()

def purely_periodic_continued_fraction(coefs):
    print 'Раскрываем понемногу'
    nomenator_free = 1
    nomenator_x = coefs[-1]
    denomenator_free = 0
    denomenator_x = 1
    for coef in coefs[-2::-1]:
        print '(%d*x + %d)/(%d*x + %d)' % (nomenator_x, nomenator_free, denomenator_x, denomenator_free)
        nomenator_x, nomenator_free, denomenator_x, denomenator_free = coef*nomenator_x+denomenator_x, coef*nomenator_free+denomenator_free, nomenator_x, nomenator_free
        print '(%d*x + %d)/(%d*x + %d)' % (nomenator_x, nomenator_free, denomenator_x, denomenator_free)
    a = denomenator_x
    b = denomenator_free - nomenator_x
    c = -nomenator_free
    print '%d*x^2 + %d*x + %d = 0' % (a, b, c, )
    d = b*b - 4*a*c
    print 'd = %d' % (d, )
    print 'x = (%d ± sqrt(%d))/%d' % (-b, d, 2*a)

def number_to_periodic_continued_fraction(normal_part, root, root_coef, denominator, prec = 0.001):
    approx_root = max(filter(lambda n: n*n <= root, xrange(root+1)))
    real_value = (float(normal_part) + root_coef * sqrt(root)) / denominator
    print 'sqrt(%d) ≈ %d' % (root, approx_root, )
    calculated = {}
    a_coefs = []
    i = 0
    while (normal_part, root_coef, denominator) not in calculated:
        #calculated.append((normal_part, root_coef, denominator))
        a = int((normal_part + sqrt(root)*root_coef) / denominator)
        a_coefs.append(a)
        calculated[(normal_part, root_coef, denominator)] = i
        print 'a_%d = %d' % (i, a, )
        normal_part, root_coef, denominator = denominator*(normal_part - a*denominator), -root_coef*denominator, (normal_part - a*denominator)**2 - root_coef**2 * root
        g = gcd(root_coef, gcd(normal_part, denominator))
        normal_part /= g
        root_coef /= g
        denominator /= g
        print 'α_%d = (%d + %d*sqrt(%d))/%d' % (i, normal_part, root_coef, root, denominator, )
        i += 1
    print 'Предпериод: %s' % (', '.join(map(str, a_coefs[:calculated[(normal_part, root_coef, denominator)]])), )
    print 'Период: %s' % (', '.join(map(str, a_coefs[calculated[(normal_part, root_coef, denominator)]:])), )
    all_as = a_coefs[:calculated[(normal_part, root_coef, denominator)]] + a_coefs[calculated[(normal_part, root_coef, denominator)]:] * 100
    ps = {-2: 0, -1: 1}
    qs = {-2: 1, -1: 0}
    a = all_as
    a_coefs = all_as
    for i in xrange(0, 50):
        ps[i] = a[i]*ps[i-1] + ps[i-2]
        qs[i] = a[i]*qs[i-1] + qs[i-2]
    real_values = map(lambda n: float(ps[n])/qs[n], xrange(30))
    first_ok = min(filter(lambda n: abs(real_value - real_values[n]) <= prec, xrange(0, 30)))
    print 'i\ta\tp\tq\tq_n(q_{n+1}+q_n)\tq_n*q_{n+1}'
    for i in xrange(first_ok + 2):
        print '%d\t%d\t%d\t%d\t%d\t\t\t%d' % (i, a[i], ps[i], qs[i], qs[i]*(qs[i+1]+qs[i]), qs[i]*qs[i+1])