polynomial.py

# 多項式のクラス
class Polynomial():
    def __init__(self,array):        
        self.__zero = (array[0]-array[0])   # 任意の係数の "0" を作るためのアクロバティックな処理
        # 多項式の次数を確認し、次数以上の係数を取り除く前処理
        d = 0
        for i in range(len(array)):
            if array[i] != self.__zero: 
                d = i
        size = d+1
        self.__array = []
        for i in range(size):
            self.__array.append(array[i])
    
    # 多項式の次数
    def deg(self):
        #print(self.__array)
        return len(self.__array)-1
    
    # 最高次係数
    def leading_coefficient():
        return self.__array[self.deg()]

    # -f を計算
    def __neg__(self):
        new_arr = []
        for i in range(self.deg()+1):
            new_arr.append(-self.__array[i])
        return self.__class__(new_arr)

    # f == g ? を計算    
    def __eq__(self, other):
        if self.deg() != other.deg():
            return False
        for i in range(self.deg()+1):
            if self.__array != other.__array:
                return False
        return True

    # f + g を計算
    def __add__(self, other):
        new_arr = []
        a_size = self.deg()+1
        b_size = other.deg()+1
        size = max(a_size, b_size)
        for i in range(size):
            a = self.__zero
            b = self.__zero
            if i < len(self.__array):
                a = self.__array[i]
            if i < len(other.__array):
                b = other.__array[i]            
            new_arr.append(a + b)
        return self.__class__(new_arr)

    # f - g を計算
    def __sub__(self, other):
        new_arr = []
        a_size = self.deg()+1
        b_size = other.deg()+1
        size = max(a_size, b_size)
        for i in range(size):
            a = self.__zero
            b = self.__zero
            if i < len(self.__array):
                a = self.__array[i]
            if i < len(other.__array):
                b = other.__array[i]            
            new_arr.append(a - b)
        return self.__class__(new_arr)

    # f * g を計算
    def __mul__(self, other):
        a_size = self.deg()+1
        b_size = other.deg()+1
        
        mul_size = self.deg() + other.deg() + 1
        mul_arr = [self.__zero for _ in range(mul_size)]

        for i in range(a_size):
            for j in range(b_size):
                mul_arr[i+j] += self.__array[i] * other.__array[j]
        return self.__class__(mul_arr)
    
    # f / g と f % g を同時に計算
    def div_mod(self, other):
        f = self.__array
        g = other.__array
        
        # 商を保存するための配列
        q_size = len(f) - len(g) + 1
        q_arr = [self.__zero for _ in range(q_size)]
        if len(q_arr) == 0:
            q_arr.append(self.__zero)
        
        # あまりを保存しておくための配列
        r_arr = [f[i] for i in range(len(f))]

        # 割り算の実行
        for i in range(q_size):
            j = len(r_arr) - 1 - i
            deg_g = other.deg()
            
            q = r_arr[j] / g[deg_g]     # 体係数だと r = g*q なる q in K^* がとれる
            q_arr[q_size - i - 1] = q
            for k in range(len(g)):
                r_arr[j-k] -= q*g[deg_g - k]
        
        return self.__class__(q_arr), self.__class__(r_arr)
    

    # f / g を計算
    def __truediv__(self, other):
        q,r = self.div_mod(other)
        return q

    # f // g を計算
    def __floordiv__(self, other):
        q,r = self.div_mod(other)
        return q

    # f % g を計算
    def __mod__(self, other):
        q,r = self.div_mod(other)
        return r

    # 多項式を文字列に変換するメソッド
    def __str__(self):
        res = []
        arr = self.__array
        for i in range(len(arr)):
            j = len(arr) - 1 - i
            res.append("{}*X^{}".format(arr[j],j))
        return " + ".join(res)
    
    # 多項式 f(X) に X = a を代入して得られる値を返す関数
    def apply(self, a):
        res = self.__zero
        for i in range(self.deg()+1):
            res += self.__array[i] * (a ** i)
        return res
    
    # 最大公約元 gcd(f,g) を返す関数
    def gcd(self, other):
        zero = self - self
        
        f = self
        g = other
        while(True):
            q = f/g
            r = f - g*q
            #print("{} = ({})*({}) + {}".format(f, q, g, r))
            if r == zero:   # 多項式としての 0 と比較
                return g
            f = g
            g = r
    
    # 拡張ユークリッドの互除法
    def euclidean(self, other):
        zero = self - self     # "0" を作る
        one = zero
        if self != zero:
            one  = self / self     # "1" を作る
        else:
            one  = other / other   # "1" を作る            
        
        f = self
        g = other    
        q_array = []

        # ユークリッドの互除法
        while(True):
            q = f/g
            r = f - g*q
            #print("euc: {} = ({}) * ({}) + {}".format(f, q, g, r))
            if r == zero:
                break
            q_array.append(q)
            f = g
            g = r

        # 拡張ユークリッドの互除法
        n = len(q_array)
        x = [one, zero]    # ここで "0, 1" を使う 
        y = [zero, one]    # ここで "0, 1" を使う
        for i in range(len(q_array)):
            x.append(q_array[i]*x[i+1] + x[i]) 
            y.append(q_array[i]*y[i+1] + y[i]) 
        
        res_x = x[len(x)-1]
        res_y = -y[len(y)-1]
        
        if len(q_array) % 2 == 0:
            res_x = -res_x
            res_y = -res_y
        
        return res_x,res_y,g
    
    # 中国剰余定理: 
    # m と n が互いに素なとき
    #   x == a (mod m)
    #   x == b (mod n)
    # なる x を出力する
    @staticmethod
    def crt4(a, b, m, n):
        s,t,g = m.euclidean(n)
        #print("test2: ({}) * ({}) + ({}) * ({}) = {}".format(m,s,n,t,g))
        d = (b - a)/g
        s = d*s
        #t = d*t
        x = a + m*s   # == b - n*t
        #print("test",a+m*s, b-n*t)
        return x % (m*n)


    # 中国剰余定理:
    # left  := [a_1, ..., a_n]
    # right := [m_1, ..., m_n]
    # m_1, ..., m_n が互いに素であるとき
    #
    #   x == a_1 (mod m_1)
    #     ...
    #   x == a_n (mod m_n)
    #
    # なる x を出力する
    @staticmethod
    def crt(left, right):
        n = len(left)
        a = left[0]
        m = right[0]
        for i in range(n-1):
            a = Polynomial.crt4(a, left[i+1], m, right[i+1]) #% (m * right[i+1])
            m = m * right[i+1]
        return a #% m



# 以下、テスト
from fractions import Fraction as Q   # Q係数多項式を考える

f = Polynomial([Q(1),Q(4),Q(6),Q(4),Q(1)])
g = Polynomial([Q(2),Q(2),Q(3)])
x = Polynomial([Q(0)])
print("f :=", f)
print("deg(f) =", f.deg())
print("g :=", g)
print("deg(g) =", g.deg())
print("x :=", x)
print("deg(x) =", x.deg())

print("")

print("# 加・減・乗算のテスト")
print("f+g =", f+g)
print("f-g =", f-g)
print("f*g =", f*g)

print("")

print("# 商・あまりの計算のテスト")
q = f/g
r = f%g
print("q := f/g =", q)
print("r := f%g =", r)

print("h := q*g+r =", q*g + r)
print("f == h:", f == q*g + r)

print("")

print("g/f =", g/f)
print("g%f =", g%f)

print("")

print("# 代入計算のテスト")
print("f(-1) =", f.apply(Q(-1)))

print("")

print("# ユークリッドの互除法のテスト：")
print("gcd(f, g) =", f.gcd(g))
x,y,d = f.euclidean(g)
print("({}) * ({}) + ({}) * ({}) = {}".format(f,x,g,y,d))


print("")
print("### 中国剰余定理（多項式）のテスト：")
print("# f(1) = 4, f(2) = 5, f(3) = 7 となるような多項式 f(X) を求めよ：")
###
# f(1) = 4
# f(2) = 5
# f(3) = 7
# となるような多項式関数 f を求めよ

a1 = Polynomial([Q(4)])
a2 = Polynomial([Q(5)])
a3 = Polynomial([Q(7)])

m1 = Polynomial([Q(-1),Q(1)])
m2 = Polynomial([Q(-2),Q(1)])
m3 = Polynomial([Q(-3),Q(1)])

func = Polynomial.crt([a1,a2,a3], [m1,m2,m3])

print("f(X) =",func)
print("  f(1) =",func.apply(Q(1)))
print("  f(2) =",func.apply(Q(2)))
print("  f(3) =",func.apply(Q(3)))

print("")
# 参考：「月を入力すると日を返す多項式」と中国剰余定理
# https://tsujimotter.hatenablog.com/entry/chinese-reminder-theorem
print("### 中国剰余定理（多項式）のテスト：")
print("# f(1) = 31, f(2) = 28, f(3) = 31, ..., f(12) = 31 となるような多項式 f(X) を求めよ：")
a_list = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
m_list = [1,  2,  3,  4,  5,  6,  7,  8,  9,  10, 11, 12]

func = Polynomial.crt([Polynomial([Q(i)]) for i in a_list], [Polynomial([Q(-i), Q(1)]) for i in m_list])
print("f(X) =",func)
for i in m_list:
    print("  f({}) = {}".format(i, func.apply(Q(i))))