对偶四元数——使用python3实现对偶四元数的符号运算 v1.0

v2.0版本：https://blog.csdn.net/yukinoai/article/details/89060211

目录

1 创建一个包含对偶算子的单项式类

2 创建多项式类

3 创建对偶四元数类

4 测试

---

实现对偶四元数简单的符号运算，数值运算

1 创建一个包含对偶算子的单项式类

主要功能：乘法，输出

# coding: utf-8
"""
@ Time: Created on 2019.03.13\n
@ Author: yukino\n
@ Description：Create a monomial class
"""


class Monomial:
    def __init__(self, coef=0, para={}, dual=False):
        self.coef = coef  # 系数
        self.var = para.keys()  # 元
        self.deg = para.values()  # 次数
        self.para = para  # 变量
        self.dual = dual  # 对偶符

    def __mul__(self, monomial):
        """乘法运算符重载"""
        if self.dual and monomial.dual is True:  # 判断是否都是对偶数
            return None  # 如果都是返回空
        else:
            result = Monomial()
            result.coef = self.coef
            result.para = self.para.copy()
            result.coef *= monomial.coef
            for k, v in monomial.para.items():
                if k in result.para.keys():  # 如果元存在则次数相加
                    result.para[k] += v
                else:
                    result.para[k] = v  # 否则在字典中添加新键值对
            result.dual = self.dual or monomial.dual  # 计算对偶符
            return result

    def multiply(self, monomial):
        result = Monomial()
        result.coef = self.coef
        result.para = self.para.copy()
        result.coef *= monomial.coef
        for k, v in monomial.para.items():
            if k in result.para.keys():
                result.para[k] += v
            else:
                result.para[k] = v
        return result

    def __str__(self):
        """输出重载"""
        self.para = dict(sorted(self.para.items(), key=lambda x: x[0]))
        term = ""
        if self.para == {}:  # 如果只有系数
            return str(self.coef)  # 只返回系数
        else:
            for key in self.para:
                if self.para[key] == 1:  # 如果次数为1则不打印幂
                    term += " " + (str(key))
                else:
                    term += " " + (str(key)) + "^" + str(self.para[key])
            if self.coef == 1:  # 如果系数为1则不打印系数
                if self.dual is True:
                    return "@(" + term.strip() + ")"
                else:
                    return term.strip()
            elif self.coef == -1:  # 如果系数为-1则只打印符号
                if self.dual is True:
                    return "@(-" + term.strip() + ")"
                else:
                    return "-" + term.strip()
            else:
                if self.dual is True:
                    return "@(" + str(self.coef) + term.strip() + ")"
                else:
                    return str(self.coef) + term.strip()


'''
# test
a = Monomial(3, {'x': 2, 'y': 3, 'z': 4})
b = Monomial(-3, {'a': 5, 'y': 6, 'z': 7})
c = Monomial(9)
d = b * a
print(d)
'''

2 创建多项式类

主要功能：相加，相减，输出，化简

# coding: utf-8
"""
@ Time: Created on 2019.03.13\n
@ Author: yukino\n
@ Description：Create a polynomial class
"""
from monomial import Monomial


class Polynomial:
    """
    多项式类
    初始化参数：任意数量Monomial类
    """
    def __init__(self, *monomial):
        self.monomials = monomial  # 单项式组成的元组

    def __str__(self):
        """输出重载"""
        term = []
        if self.monomials.__len__() == 1:  # 如果只有一个单项式则不打印括号
            term.append(str(self.monomials[0]))
        else:
            for monomial in self.monomials:
                if monomial not in [None]:  # 如果单项式不为空
                    if monomial.coef < 0:  # 系数小于0则加括号
                        term.append("(" + str(monomial) + ")")
                    elif monomial.coef > 0:
                        term.append(str(monomial))
        term = ' + '.join(term)
        return term

    def __add__(self, polynomial):
        """加法运算符重载"""
        return Polynomial(*(self.monomials + polynomial.monomials))

    def __sub__(self, polynomial):
        """减法运算符重载"""
        minus = []
        for monomial in polynomial.monomials:
            temp = Monomial()
            temp.coef = -monomial.coef
            temp.para = monomial.para.copy()
            temp.dual = monomial.dual
            minus.append(temp)
        return Polynomial(* (self.monomials + tuple(minus)))

    def __mul__(self, polynomial):
        """乘法运算符重载"""
        mult = []
        for monomialP in polynomial.monomials:
            for monomialR in self.monomials:
                mult.append(monomialR * monomialP)
        return Polynomial(* mult)


def simplify(poly):
    """
    多项式化简函数\n
    Input: Polynomial Class\n
    Output: Polynomial Class
    """
    # 分离实部与对偶部
    real = []
    dual = []
    for monomial in poly.monomials:
        if monomial.dual is True:
            dual.append(monomial)
        else:
            real.append(monomial)
    # 简化实部
    resparas = []
    rescoefs = []
    monomials = []
    paras = [monomial.para for monomial in real]
    coefs = [monomial.coef for monomial in real]
    for i in range(len(paras)):
        if not paras[i] in resparas:
            resparas.append(paras[i])
            rescoefs.append(coefs[i])
        else:
            rescoefs[resparas.index(paras[i])] += coefs[i]
    for i in range(len(rescoefs)):
        monomials.append(Monomial(rescoefs[i], resparas[i]))
    # 简化对偶部
    resparas = []
    rescoefs = []
    paras = [monomial.para for monomial in dual]
    coefs = [monomial.coef for monomial in dual]
    for i in range(len(paras)):
        if not paras[i] in resparas:
            resparas.append(paras[i])
            rescoefs.append(coefs[i])
        else:
            rescoefs[resparas.index(paras[i])] += coefs[i]
    for i in range(len(rescoefs)):
        monomials.append(Monomial(rescoefs[i], resparas[i], True))
    simp = Polynomial(*monomials)
    return simp


'''
# test
a = Monomial(3, {'x': 2, 'y': 3, 'z': 4}, True)
b = Monomial(-3, {'a': 5, 'y': 6, 'z': 7})
c = Ploynomial(a, b)
A = Ploynomial(a)
B = Ploynomial(b)
C = A + A
x0 = symbol('x0', True)
D = simplify(C)
print(A)
print(B)
print(C)
print(D)
'''

3 创建对偶四元数类

主要功能：相加，相减，输出，化简

# coding: utf-8
"""
@ Time: Created on 2019.03.13\n
@ Author: yukino\n
@ Description：Create a dual-quaternion class
"""
from polynomial import Polynomial, simplify
from monomial import Monomial


class DualQuaternion:
    """
    对偶四元数类\n
    初始化参数：List[8 个 Ploynomial 类]
    """
    def __init__(self, quaterlist):
        """构造函数, 1个由8个多项式类构成的列表"""
        self.real = quaterlist[:4]
        self.dual = quaterlist[4:]
        self.all = quaterlist
        
    def __str__(self):
        """输出重载"""
        q = ['', 'i', 'j', 'k']
        result = []
        for i in range(4):
            if self.real[i].monomials[0].coef != 0:  # 如果只有系数0则不输出
                result.append('(' + str(self.real[i]) + ')' + q[i])
            if self.dual[i].monomials[0].coef != 0:
                result.append('(' + str(self.dual[i]) + ')' + q[i])
        term = ' + '.join(result)
        return term

    def __add__(self, dualquater):
        """对偶四元数加法重载"""
        result = []
        self = self.all
        dualquater = dualquater.all
        for i in range(8):
            result.append(self[i] + dualquater[i])
        return DualQuaternion(result)

    def __sub__(self, dualquater):
        """对偶四元数减法重载"""
        result = []
        self = self.all
        dualquater = dualquater.all
        for i in range(8):
            result.append(self[i] - dualquater[i])
        return DualQuaternion(result)

    def __mul__(self, dualquater):
        """对偶四元数乘法重载"""
        q = self.all
        p = dualquater.all
        x0 = q[0]*p[0] - q[1]*p[1] - q[2]*p[2] - q[3]*p[3]
        x1 = q[0]*p[1] + q[1]*p[0] + q[2]*p[3] - q[3]*p[2]
        x2 = q[0]*p[2] - q[1]*p[3] + q[2]*p[0] + q[3]*p[1]
        x3 = q[0]*p[3] + q[1]*p[2] - q[2]*p[1] + q[3]*p[0]
        y0 = q[0]*p[4] - q[1]*p[5] - q[2]*p[6] - q[3]*p[7] + q[4]*p[0] - q[5]*p[1] - q[6]*p[2] - q[7]*p[3]
        y1 = q[0]*p[5] - q[1]*p[4] - q[2]*p[7] - q[3]*p[6] + q[4]*p[1] - q[5]*p[0] - q[6]*p[3] - q[7]*p[2]
        y2 = q[0]*p[6] - q[1]*p[7] - q[2]*p[4] - q[3]*p[5] + q[4]*p[2] - q[5]*p[3] - q[6]*p[0] - q[7]*p[1]
        y3 = q[0]*p[7] - q[1]*p[6] - q[2]*p[5] - q[3]*p[4] + q[4]*p[3] - q[5]*p[2] - q[6]*p[1] - q[7]*p[0]
        return DualQuaternion([x0, x1, x2, x3, y0, y1, y2, y3])
    
    def toArray(self):
        """
        转换为列表[字符串]格式\n
        Output: List[str]
        """
        return [str(i) for i in self.all]

    def dualqsimp(self):
        """
        对偶四元数类简化函数\n
        Output: Polynomial Class
        """
        result = []
        for i in self.all:
            result.append(simplify(i))
        return DualQuaternion(result)


def symbol(sym, dual=False):
    """
    符号定义函数，定义为多项式类\n
    Input: str\n
    Output: Polynomial Class
    """
    return Polynomial(Monomial(1, {sym: 1}, dual))


def number(num, dual=False):
    """
    数字定义函数，定义为多项式类\n
    Input: int or float\n
    Output: Polynomial Class
    """
    return Polynomial(Monomial(num, dual=dual))


def symbols(*syms):
    """一次定义多个符号的函数，方便定义符号，由于使用全局变量，应谨慎使用"""
    names = globals()
    for sym in syms:
        if sym[0] == '@':
            names[sym[1:]] = Polynomial(Monomial(1, {sym[1:]: 1}, True))
        else:
            names[sym] = Polynomial(Monomial(1, {sym: 1}))


def quaterSym(*syms):
    """
    生成对偶四元数类所需的列表参数\n
    Input: 8 str, int or float\n
    Output: List of Polynomial Class
    """
    result = []
    for i in range(4):
        if isinstance(syms[i], str):
            result.append(symbol(syms[i]))
        else:
            result.append(number(syms[i]))
    for i in range(4, 8):
        if isinstance(syms[i], str):
            result.append(symbol(syms[i], True))
        else:
            result.append(number(syms[i], True))
    return result

4 测试

# coding: utf-8
"""
@ Time: Created on 2019.03.13\n
@ Author: yukino\n
@ Description：For debugging programs
"""
import daulquaternion as dq

# 创建对偶四元数参数，可以是符号也可以是数字
e = dq.quaterSym('x0', 'x1', 'x2', 'x3', 'y0', 'y1', 'y2', 'y3')
g = dq.quaterSym(1, 2, 3, 4, 0, 0, 0, 0)
p = dq.quaterSym(1, -2, -3, -4, 0, 0, 0, 0)

q1 = dq.DualQuaternion(p)
q2 = dq.DualQuaternion(g)
q3 = dq.DualQuaternion(e)
q4 = q1 * q2  # 数字运算
q5 = (q2 * q2).dualqsimp()
q6 = q1 * q3  # 符号运算
print(q1)
print(q3)
print(q4.dualqsimp())
print(q5)
print(q6.dualqsimp().toArray())

结果如下：