chakokuのブログ(rev4)

テック・コミック・DTM・・・ごくまれにチャリ

Pythonでquotanionを扱えるようにクラスを定義してみる

多分Scipyとかにあるんだろうけど、Quotanion(四元数)の演算をプログラムでできるようにPythonでクラスを定義してみた*1
証明とか手で計算していると間違いまくるので、簡単な式変換ならプログラムにやってもらおうという。。
演算子オーバーロード(?)等を活用すると、q1.mul(q2)とせずとも、q1 * q2 とクールに書けるんだろうけど、そのあたり不勉強なので、愚直にメソッドを呼んで計算させる。
注意:強烈なバグが残っている

#!/usr/bin/python3
#
#


class Quotanion:

    def __init__(this,a,b,c,d):
        this.term  = { "1" : a , "i" : b , "j" : c , "k" : d }
        this.keys = ("1", "i", "j", "k")

    def mul(this,arg):
        exp = []
        term  = { "1" : "", "i" : "" , "j" : "" , "k" : "" }
        for key1 in this.keys:
           for key2 in arg.keys:
              new_term = this.term[key1] + arg.term[key2] 
              imagNumber = key1 + key2
              if imagNumber == "11":
                 imagNumber = None
              elif imagNumber == "1i" or imagNumber == "i1" :
                 imagNumber = "i"
              elif imagNumber == "1j" or imagNumber == "j1" :
                 imagNumber = "j"
              elif imagNumber == "1k" or imagNumber == "k1" :
                 imagNumber = "k"
              elif imagNumber == "ij":
                 imagNumber = "k"
              elif imagNumber == "ji":
                 imagNumber = "-k"
              elif imagNumber == "jk":
                 imagNumber = "i"
              elif imagNumber == "kj":
                 imagNumber = "-i"
              elif imagNumber == "ki":
                 imagNumber = "j"
              elif imagNumber == "ik":
                 imagNumber = "-j"

              if not imagNumber:
                 term["1"] += new_term
              elif "i" in imagNumber:
                 if "-" in imagNumber:
                    term["i"] += "-" + new_term
                 else:
                    term["i"] += "+" + new_term
              elif "j" in imagNumber:
                 if "-" in imagNumber:
                    term["j"] += "-" + new_term
                 else:
                    term["j"] += "+" + new_term
              elif "k" in imagNumber:
                 if "-" in imagNumber:
                    term["k"] += "-" + new_term
                 else:
                    term["k"] += "+" + new_term


        q = Quotanion(term["1"], term["i"], term["j"], term["k"])
        return(q)        

    def print(this):
        for key in this.keys[:-1]:
            if key == "1":
                print(f"{this.term[key]} + ",end="")
            else:
                print(f"({this.term[key]}){key} + ",end="")
        key = this.keys[-1]
        print(f"({this.term[key]}){key}")


q1 = Quotanion('a','b','c','d')
q1.print()

q2 = Quotanion('x','y','z','w')
q2.print()

q3 = q1.mul(q2)
q3.print()


# $ python3 qota.py
# a + (b)i + (c)j + (d)k
# x + (y)i + (z)j + (w)k
# ax + (+ay+bx+by+cw-dz)i + (+az-bw+cx+cz+dy)j + (+aw+bz-cy+dx+dw)k


内積はおかしいが、少し改良版。虚数同士の計算はできるが、係数が(a+b)iといった式だとまとに動かない。表面的な実装だとこれぐらいが限界で、式を解釈できるように構文解析が必要だ。lispのS式だったら構文解析なんて面倒な処理いらんけど、pythonベースだからしょうがない。

#!/usr/bin/python3
#
#


# flag for debug print 
DBG = False


class Quotanion:

    def __init__(this,a,b,c,d):
        this.term  = { "1" : a , "i" : b , "j" : c , "k" : d }
        this.type = "Quotanion"

    def set_type(this,type):
        this.type = type

    def size(this):
        exp = "sqrt("
        for key in this.term.keys()[:-1]:
             exp += f"{this.term[key]} + "
        key = this.term.keys()[-1]
        exp += f"{this.term[key]}"
        exp += ")"
        return(exp)

    def conjugate(this):
        term  = { "1" : "", "i" : "" , "j" : "" , "k" : "" }
        term["1"] = this.term["1"]

        for key in ("i","j","k"):
            if "-" in this.term[key]:
                term[key] = this.term[key].replace("-","")
            else:
                term[key] = "-" + this.term[key]
        q = Quotanion(term["1"], term["i"], term["j"], term["k"])
        return(q)        

    def iprod(this,arg):
        term  = { "1" : "", "i" : "" , "j" : "" , "k" : "" }
        term["1"] = this.term["1"] + arg.term["1"]
        term["i"] = this.term["i"] + arg.term["i"]
        term["j"] = this.term["j"] + arg.term["j"]
        term["k"] = this.term["k"] + arg.term["k"]
        q = Quotanion(term["1"], term["i"], term["j"], term["k"])
        q.set_type("innerProduct")
        return(q)        

    def add(this,arg):
        exp = []
        term  = { "1" : "", "i" : "" , "j" : "" , "k" : "" }
        for key in term.keys():
              new_term = this.term[key] 
              if "-" in arg.term[key]:
                   new_term += arg.term[key]
              else:
                   new_term += "+" + arg.term[key]

              term[key] = new_term
        q = Quotanion(term["1"], term["i"], term["j"], term["k"])
        return(q)        

    def mul(this,arg):
        exp = []
        term  = { "1" : "", "i" : "" , "j" : "" , "k" : "" }
        for key1 in this.term.keys():
           for key2 in arg.term.keys():

              if DBG:
                  print(f"{this.term[key1]}{key1} * {arg.term[key2]}{key2} -> ",end="")
              negaFlag = False
              if "-" in this.term[key1]:
                   negaFlag = not negaFlag 
                   new_term = this.term[key1].replace("-","")
              else:
                   new_term = this.term[key1]
              if "-" in arg.term[key2]:
                   negaFlag = not negaFlag
                   new_term += arg.term[key2].replace("-","")
              else:
                   new_term += arg.term[key2]

              image_number = key1 + key2
              if image_number == "11":
                 image_number = None
              elif image_number == "1i" or image_number == "i1" :
                 image_number = "i"
              elif image_number == "1j" or image_number == "j1" :
                 image_number = "j"
              elif image_number == "1k" or image_number == "k1" :
                 image_number = "k"
              elif image_number == "ij":
                 image_number = "k"
              elif image_number == "ji":
                 negaFlag = not negaFlag
                 image_number = "k"
              elif image_number == "jk":
                 image_number = "i"
              elif image_number == "kj":
                 negaFlag = not negaFlag
                 image_number = "i"
              elif image_number == "ki":
                 image_number = "j"
              elif image_number == "ik":
                 negaFlag = not negaFlag
                 image_number = "j"
              elif image_number == "ii" or image_number == "jj" or image_number == "kk":
                 image_number = None
                 negaFlag = not negaFlag

              if DBG:
                  if negaFlag:
                       print(f"-{new_term}{image_number}")
                  else:
                       print(f"{new_term}{image_number}")

              if not image_number:
                 if negaFlag:
                    term["1"] += "-" + new_term
                 else:
                    term["1"] += "+" + new_term
              elif image_number == "i":
                 if negaFlag:
                    term["i"] += "-" + new_term
                 else:
                    term["i"] += "+" + new_term
              elif image_number == "j":
                 if negaFlag:
                    term["j"] += "-" + new_term
                 else:
                    term["j"] += "+" + new_term
              elif image_number == "k":
                 if negaFlag:
                    term["k"] += "-" + new_term
                 else:
                    term["k"] += "+" + new_term
        q = Quotanion(term["1"], term["i"], term["j"], term["k"])
        return(q)        



    def print(this):
        keys = ("1","i","j","k")
        if this.type == "innerProduct":
           print("sqrt(",end="")
           for key in keys[:-1]:
                print(f"{this.term[key]} + ",end="")
           key = keys[-1]
           print(f"{this.term[key]}",end="")
           print(")")

        else:
          for key in keys[:-1]:
            if key == "1":
                print(f"{this.term[key]} + ",end="")
            else:
                print(f"({this.term[key]}){key} + ",end="")
          key = keys[-1]
          print(f"({this.term[key]}){key}")


#
# test
#

q1 = Quotanion('a','b','c','d')
#q1.print()

q2 = Quotanion('x','y','z','w')
#q2.print()


if False:
    print("-----add----")
    q1.print()
    q2.print()
    print("-------------")
    q3 = q1.add(q2)
    print("-----ans----")
    q3.print()


if False:
    print("-----multi----")
    q1.print()
    q2.print()
    print("-------------")
    q3 = q1.mul(q2)
    print("-----ans----")
    q3.print()

#q4 = q1.iprod(q2)
#q4.print()

#print(q1.size())
#print(q2.size())

if False:
    print("----q1----")
    q1.print()
    print("----conjugate of q1----")
    q1c = q1.conjugate()
    q1c.print()


if False:
    q1c = q1.conjugate()
    print("-----multi----")
    q1.print()
    q1c.print()
    print("-------------")
    qq = q1.mul(q1c)
    print("-----ans----")
    qq.print()

a = Quotanion('a','b','c','d')
ac = a.conjugate()

b = Quotanion('x','y','z','w')
bc = b.conjugate()

if False:
    print("--- a(~b) + (b)(~a) -------")
    a.mul(bc).add(b.mul(ac)).print()
#+ax+by+cz+dw++xa+yb+zc+wd + (-ay+bx-cw+dz-xb+ya-zd+wc)i + (-az+bw+cx-dy-xc+yd+za-wb)j + (-aw-bz+cy+dx-xd-yc+zb+wa)k
# 2ax + 2by + 2cz + 2dw  

*1:公式ライブラリとはレベルが違いすぎるので、車輪の再発明とは決して言わない