Algebra-with-Python/各章代码/utility.py

import numpy as np
from fractions import Fraction as Q
def exp(a):
    n=len(a)
    s=''
    start=0
    while start<n and a[start]==0:#查找第一个非零项下标start
        start+=1
    if start>=n:#零多项式
        s+='0'
    else:#非零多项式
        if start==0:#非零常数项
            s+='%s'%a[start]
        if start==1:#非零1次项
            s+='%s∙x'%a[start]
        if start>1:#首项次数大于1
            s+='%s∙x^%d'%(a[start],start)
        for i in range(start+1,n):#首项后各项
            if i==1:#1次项
                if a[i]>0:#正1次项
                    s=s+'+%s∙x'%a[i]
                if a[i]<0:#负1次项
                    s=s+'%s∙x'%a[i]
            else:#非1次项
                if a[i]>0:#正项
                    s=s+'+%s∙x**%d'%(a[i],i)
                if a[i]<0:#负项
                    s=s+'%s∙x**%d'%(a[i],i)
    return s
class myPoly:#多项式类
    def __init__(self, coef):#初始化
        c=np.trim_zeros(coef,trim='b')#删除高次零系数
        self.coef=np.array(c)#设置多项式系数
        self.degree=(self.coef).size-1#设置多项式次数
    def __eq__(self, other):#相等
        if self.degree!=other.degree:#次数不等
            return False
        return (abs(self.coef-other.coef)<1e-10).all()
    def __str__(self):#输出表达式
        return exp(self.coef)        
    def __add__(self, other):#运算符“+”
        n=max(self.degree,other.degree)+1#系数个数
        a=np.append(self.coef,[0]*(n-self.coef.size))#补长
        b=np.append(other.coef,[0]*(n-other.coef.size))#补长
        return myPoly(a+b)#创建并返回多项式和
    def __rmul__(self, k):#右乘数k
        c=self.coef*k#各系数与k的积
        return myPoly(c)
    def __neg__(self):#负多项式
        return (-1)*self
    def __sub__(self, other):#运算符“-”
        return self+(-other)
    def __mul__(self,other):#运算符“*”
        m=self.degree
        n=other.degree
        a=np.append(self.coef,[0]*n)
        b=np.append(other.coef,[0]*m)
        c=[]
        for s in range(1, m+n+2):
            cs=(a[:s]*np.flip(b[:s])).sum()
            c.append(cs)
        return myPoly(c)
    def val(self,x):#多项式在x处的值
        n=self.degree#次数
        power=np.array([x**k for k in range(n+1)])#x各次幂
        v=((self.coef)*power).sum()#多项式的值
def P1(A, i, j, row=True):
    if row:
        A[[i,j],:]=A[[j,i],:]
    else:
        A[:,[i,j]]=A[:,[j,i]]
def P2(A,i,k, row=True):
    if row:
        A[i]=k*A[i]
    else:
        A[:,i]=k*A[:,i]
def P3(A,i,j,k, row=True):
    if row:
        A[j]+=k*A[i]
    else:
        A[:,j]+=k*A[:,i]
def Aij(A,i,j):
    up=np.hstack((A[:i,:j],A[:i,j+1:]))
    lo=np.hstack((A[i+1:,:j],A[i+1:,j+1:]))
    M=np.vstack((up,lo))
    return ((-1)**(i+j))*np.linalg.det(M)
def adjointMatrix(A):
    n,_=A.shape
    Am=np.zeros((n,n))
    for i in range(n):
        for j in range(n):
            Am[j,i]=Aij(A,i,j)
    return Am
def cramer(A, b):
    n=len(A)#未知数个数
    A=np.hstack((A, b.reshape(n,1)))#增广矩阵
    x=[]#解初始化为空
    detA=np.linalg.det(A[:,:n])#系数行列式
    if detA!=0:#有解
        for i in range(n):#计算每一个未知数的值
            A[:,[i,n]]=A[:,[n,i]]#用b替换A的第i列
            detAi=np.linalg.det(A[:,:n])#行列式
            x.append(detAi/detA)#第i个未知数的值
            A[:,[i,n]]=A[:,[n,i]]#还原系数矩阵
    return np.array(x)
def rowLadder(A, m, n):
    rank=0#系数矩阵秩初始化为0
    zero=m#全零行首行下标
    i=0#当前行号
    order=np.array(range(n))#未知数顺序
    while i<min(m,n) and i<zero:#自上向下处理每一行
        flag=False#B[i,i]非零标志初始化为F
        index=np.where(abs(A[i:,i])>1e-10)#当前列A[i,i]及其以下的非零元素下标
        if len(index[0])>0:#存在非零元素
            rank+=1#累加秩
            flag=True#B[i,i]非零标记
            k=index[0][0]#非零元素最小下标
            if k>0:#若非第i行，交换第i，k+i行
                P1(A,i,i+k)
        else:#A[i:,i]内全为0
            index=np.where(abs(A[i,i:n])>1e-10)#当前行A[i,i:n]的非零元素下标           
            if len(index[0])>0:#存在非零元素，交换第i，k+i列
                rank+=1
                flag=True
                k=index[0][0]
                P1(A,i,i+k,row=False)#列交换
                order[[i, k+i]]=order[[k+i, i]]#跟踪未知数顺序
        if flag:#A[i,i]不为零，A[i+1:m,i]消零
            P2(A,i,1/A[i,i])
            for t in range(i+1, zero):
                P3(A,i,t,-A[t,i])
            i+=1#下一行
        else:#将全零元素行交换到矩阵底部
            P1(A,i,zero-1)
            zero-=1#更新全零行首行下标
    return rank, order
def simplestLadder(A,rank):
    for i in range(rank-1,0,-1):#主对角线上方元素消零
        for k in range(i-1, -1,-1):
            P3(A,i,k,-A[k,i])
def mySolve(A,b):
    m,n=A.shape                                     #系数矩阵结构
    b=b.reshape(b.size, 1)                          #常量列向量
    B=np.hstack((A, b))                             #构造增广矩阵
    r, order=rowLadder(B, m, n)                     #消元
    X=np.array([])                                     #解集初始化为空
    index=np.where(abs(B[:,n])>1e-10)               #常数列非零元下标
    nonhomo=index[0].size>0                         #判断是否非齐次
    r1=r                                            #初始化增广矩阵秩
    if nonhomo:                                     #非齐次
        r1=np.max(index)+1                          #修改增广阵秩
    solvable=(r>=r1)                                #判断是否可解
    if solvable:                                    #若可解
        simplestLadder(B, r)#回代
        X=np.vstack((B[:r,n].reshape(r,1),
                            np.zeros((n-r,1))))
        if r<n:
            x1=np.vstack((-B[:r,r:n],np.eye(n-r)))
            X=np.hstack((X,x1))
        X=X[order]
    return X
def matrixSolve(A,B):
    m,n=A.shape
    _,n1=B.shape
    C=np.hstack((A, B))
    r, order=rowLadder(C, m, n)
    simplestLadder(C, r)
    solution=[]
    X=np.array([])
    index=np.where(abs(C[:,n:])>1e-10)
    r1=max(index[0])+1
    if r==r1:
        X=np.vstack((C[:r,n:],
                     np.zeros((n-r,n1))))
        X=X[order,:]
    return X
def maxIndepGrp(A):
    B=A.copy()
    m,n=B.shape
    r,order=rowLadder(B,m,n)
    simplestLadder(B,r)
    return r,order,B[:r,r:]
def orthogonalize(A):
    m,n=A.shape
    B=A.copy()
    for i in range(1,n):
        for j in range(i):
            B[:,i]-=(np.dot(B[:,j],A[:,i])/
                     np.dot(B[:,j],B[:,j]))*B[:,j]
    return B
def unitization(A):
    _,n=A.shape
    for i in range(n):
        A[:,i]/=np.linalg.norm(A[:,i])
def symmetrization(A):
    n,_=A.shape
    for i in range(n):
        for j in range(i+1,n):
            A[i,j]=(A[i,j]+A[j,i])/2
            A[j,i]=A[i,j]