quadratic_integer.py

# coding: utf-8

# 2020/01/19
# 平方因子を持たない0,1を除く有理整数 m について、
# 2次体 Q(√m) の整数環 Z[ω] のイデアルの積を計算できるプログラム
#
# ただし、ωはmに応じて次のルールで計算される:
#   ω = √m           (m ≡ 2, 3 (mod 4))
#   ω = (1+√m)/2     (m ≡ 1 (mod 4))
#
# 整数環Z[ω]の任意の整数 a+bω（a,bは有理整数）は、Pythonのリスト [a,b] の形式で表現される
# 標準基底で表されたZ[ω]のイデアル I = [a, b+cω]（a,b,cは有理整数）は、Pythonのリスト [a, b, c] の形式で表現される

import math

m = -14

# 整数環Z[ω]の整数 a+bω（a,bは有理整数）（作りかけ）
class QuadInteger:
    def __init__(self, a, b, m):
        self.a = a
        self.b = b
        self.m = m

    def __str__(self):
        str_a = "{}".format(self.a)
        if self.a < 0:
            str_a = "({})".format(self.a)

        str_b = "{}".format(self.b)
        if self.a < 0:
            str_b = "({})".format(self.b)

        if m % 4 == 1:
            return "{}+{}*ω".format(str_a, str_b)
        else: 
            return "{}+{}*√({})".format(str_a, str_b, self.m)



# Z[ω]の整数を文字列にする
def disp(a):
    str_a = "{}".format(a[0])
    if a[0] < 0:
        str_a = "({})".format(a[0])

    str_b = "{}".format(a[1])
    if a[1] < 0:
        str_b = "({})".format(a[1])

    if m % 4 == 1:
        return "{}+{}*ω".format(str_a, str_b)
    else: 
        return "{}+{}*√({})".format(str_a, str_b, m)


# Z[ω]の整数 a,b の差を計算する
def sub(a, b):
    return [a[0]-b[0], a[1]-b[1]]

# Z[ω]の整数 a,b の積を計算する
def mul(a, b):
    if m % 4 == 1:
        dd = (1 - m)//4
        return [a[0]*b[0]-dd*a[1]*b[1], a[0]*b[1]+a[1]*b[0]-a[1]*b[1]]
    else:
        return [a[0]*b[0]+m*a[1]*b[1], a[0]*b[1]+a[1]*b[0]]

# Z[ω]の整数 a の共役を計算する
def conjugate(a):
    if m % 4 == 1:
        return [a[0]+a[1], -a[1]]
    else:
        return [a[0], -a[1]]

# Z[ω]の整数 a のノルムを計算する
def norm(a):
    return mul(a, conjugate(a))[0]



# Z[ω]のイデアル ideal を表示する
def disp_ideal(ideal):
    a = ideal[0]
    b = ideal[1]
    c = ideal[2]
    return "[{}, {}]".format(a, disp([b, c]) )

# 標準基底で表された2つのイデアル 
# I1 = [a1, b1+c1ω], I2 = [a2, b2+c2ω]
# の積を計算する
def mul_ideals(ideal1, ideal2):
    x = [ideal1[0], 0]
    y = [ideal1[1], ideal1[2]]
    z = [ideal2[0], 0]
    w = [ideal2[1], ideal2[2]]

    print("*** CALL MUL_IDEALS(I1, I2) ***")
    print("*** I1 = {}".format(disp_ideal(ideal1)))
    print("*** I1 = {}".format(disp_ideal(ideal2)))
    # イデアル I1 * I2 の生成元4つを計算する
    a = mul(x,z)
    b = mul(y,z)
    c = mul(x,w)
    d = mul(y,w)

    # 生成元は、ωの係数が非負になるように正規化する（(-1)倍する）
    if a[1] < 0:
        a = mul(a, [-1, 0]) 
    if b[1] < 0:
        b = mul(b, [-1, 0]) 
    if c[1] < 0:
        c = mul(c, [-1, 0]) 
    if d[1] < 0:
        d = mul(d, [-1, 0]) 

    res = [a, b, c, d]

    # イデアル I1 * I2 の生成元が2個になるまで削減する
    while len(res) > 2:
        n = len(res)
        print("***   {}".format(res))

        # ピボットを選択: ωの係数を比較して、一番小さいものを探す
        for i in range(2,n):
            if res[1][1] > res[i][1]:
                tmp = res[1]
                res[1] = res[i]
                res[i] = tmp
        
        # ピボットの設定: ωの0でない最小の係数をピボットとする
        pivot_index = 1
        for i in range(1,n):
            if res[i][1] > 0:
                pivot_index = i
                break
        pivot = res[pivot_index]

        # ピボット(pivot)を定数倍して引く
        for i in range(pivot_index+1,n):
            k = res[i][1]//pivot[1]
            res[i] = sub( res[i], mul(pivot, [k, 0]) )
            res[i][0] = res[i][0] % res[0][0]
        
        d = res[0][0]   # 計算結果が整数になった生成元の最大公約数を計算する
        for i in range(1, n):
            if res[i][1] == 0 and res[i][0] != 0:
                d = math.gcd(d, res[i][0])
        new_res = [ [d, 0] ]

        for i in range(pivot_index, n):
            if res[i][1] != 0:
                new_res.append(res[i])
        res = new_res

    print("***   {}".format(res))
    
    if len(res) == 2:
        # 生成元が2個以上
        result = [res[0][0], res[1][0], res[1][1]]
    else:
        # 生成元が1個
        result = [res[0][0], 0, res[0][0]]   # a[1,ω] = [a, aω]
    
    print("*** I1 * I2 = {} ***".format(disp_ideal(result)))
    return result


a = [1, 1]  # 1 + 1*sq(-14)
b = [16, 1] # 16 + 1*sq(-14)
print(" MUL_NUMBERS:")
print( "{} * {} = {}".format(disp(a), disp(b), disp(mul(a, b))) )     # 2 + 17*sq(-14)
print("")

#'''
ideal1 = [3, 1, 1]    # (a, b+cω) = (3, 1+ω)
ideal2 = [9, 7, 1]    # (a, b+cω) = (9, 7+ω)
#'''

#ideal1 = [2, 1, 1]    # (a, b+cω) = (3, 1+ω)
#ideal2 = [3, 1, 1]    # (a, b+cω) = (9, 7+ω)

print(" MUL_IDEALS:")
print( "{} * {} = {}".format( disp_ideal(ideal1), disp_ideal(ideal2), disp_ideal(mul_ideals(ideal1, ideal2)) ) )
print("")


ideal1 = [2, 0, 1]    # (a, b+cω) = (2, ω)
ideal2 = [2, 0, 1]    # (a, b+cω) = (2, ω)
print(" MUL_IDEALS:")
print( "{} * {} = {}".format( disp_ideal(ideal1), disp_ideal(ideal2), disp_ideal(mul_ideals(ideal1, ideal2)) ) )
print("")


q = [3, 1, 1]   # [3, 1+√-14]
print("Q   = {}".format( disp_ideal(q) ))
print("")

q2 = mul_ideals(q, q)
print("Q^2 = {}".format( disp_ideal(q2) ))
print("")
q3 = mul_ideals(q2, q)
print("Q^3 = {}".format( disp_ideal(q3) ))
print("")
q4 = mul_ideals(q3, q)
print("Q^4 = {}".format( disp_ideal(q4) ))
print("")