多項式の解き方をpythonコードで実装してみる。

多項式分解のコードをDKA法で書いてみました。行列の固有多項式の計算用のコードと合わせて載せてあります。最初は稚拙なもので、このコードでは15桁くらいまで桁が上がると時間がかかったり、うまく解を求められなくなったりします。

DKA法による多項式の因数分解の他に、できた多項式の解の重複の確認のために、行列の性質を使っています。ある行列Aと単位行列の整数乗で構成される行列が変数の方程式を考えた際に、対角成分に方程式の解となる数値を置いた行列を代入すると、０行列になります。また、対角成分の右隣に1を置いて同じく代入すると、同じ行の対角成分が両方方程式の解を為す場合にのみ、0行列になる性質があります。この性質を利用して、求めた解の重複の確認をおこなっています。

また、重複の確認を行った後は、計算した解から方程式を作り、元の式から割り算をすることで、元の方程式の簡略化をおこなっています。

出力はこんな感じ。流石に何回かトライしないと失敗する。

import numpy as np
import scipy.linalg as linalg
import random

def solv(a0,b0,c0):
    D=b0**2-4*a0*c0
    if D<0:
        D=complex(0,(-D)**0.5)
    x_1=(-b0+D)/2/a0
    x_2=(-b0-D)/2/a0
    return x_1,x_2

import random
def warizan(keisuu,bunbo,krank,brank):
    p=keisuu
    b=bunbo
    jisuuk=krank
    jisuub=brank
    ci=np.array([0.0]*(jisuuk-jisuub+1),dtype=complex)
    ama=np.array([0.0]*(jisuub),dtype=complex)
    #print(jisuuk-jisuub)
    for i in range(jisuuk-jisuub+1):
        ci[i]=p[i]/b[0]
    #print(ci,ama,b,jisuuk,jisuub)
    for i in range(jisuub):
        ama[i]=p[i+jisuuk-jisuub+1]
    for i in range(1,jisuuk+1):
        
        for j in range(1,jisuub+1):
            if i>=j :
        
                if i<=jisuuk-jisuub:
                    ci[i]-=ci[i-j]*b[j]/b[0]
                elif i>jisuuk-jisuub and i-j<jisuuk-jisuub+1:
                    ama[i-jisuuk+jisuub-1]-=ci[i-j]*b[j]
    return ci,ama,jisuuk-jisuub
    
def dainyu(xval,keisuu):
    x=xval
    pkei=keisuu
    rui=1
    kekka=0
    for i in range(len(pkei)):
        kekka+=pkei[len(pkei)-i-1]*rui
        rui*=xval
    return kekka

def bibun(xval,keisuu):
    x=xval
    pkei=keisuu
    rui=1
    kekka=0
    for i in range(1,len(pkei)):
        kekka+=pkei[len(pkei)-i-1]*rui*i
        rui*=xval
    return kekka
def dkahou(keisuu,keisanjougen=13):
    jougen=keisanjougen
    seido=1.0e-10
    klist=keisuu
    jisuu=len(klist)-1
    bibklist=[0]*(jisuu-1)
    #print(jisuu,klist)
    #if jisuu<jougen:
    jougen=jisuu
    for i in range(jisuu-1):
        bibklist[i]=(jisuu-i)*klist[i]
    i0=0
    waribase=[[]]*jisuu
    amaribase=[[]]*jisuu
    tugir=jisuu
    tugik=klist
    bk=[1,klist[1]/jisuu]
    while tugir>=1:
        tugik,amari0,tugir=warizan(tugik,bk,tugir,1)
        waribase[i0]=tugik
        amaribase[i0]=amari0
        i0+=1
    #print(bk)
    #print(waribase)
    #print(amaribase)
    kumitatek=[0]*(jisuu+1)
    kumitatek[0]=1
    for i in range(jisuu):
        kumitatek[i+1]=-abs(amaribase[jisuu-i-1][0])
    #print(kumitatek)
    r=1
    for i in range(30):
        rsabun=dainyu(r,kumitatek)
        r-=rsabun/bibun(r,kumitatek)
        if rsabun<seido:
            break
        #print(r)
    shokiti=[0]*jougen
    if jougen<0:
        for i in range(jougen):
        #shokiti[i]=-klist[1]/jisuu/klist[0]+r*complex(np.sin(2*np.pi/jougen*i+np.pi/2/jougen),np.cos(2*np.pi/jougen*i+np.pi/2/jougen))
            randoms=random.random()
            shokiti[i]=-klist[1]/jisuu/klist[0]+r*complex(np.sin(i/jougen*2*np.pi+np.pi/2/jougen),np.cos(i/jougen*2*np.pi+np.pi/2/jougen))
    else:
        for i in range(jougen):
            randoms=random.random()
            shokiti[i]=-klist[1]/jisuu/klist[0]+r*complex(np.cos(2*randoms*np.pi),np.sin(2*randoms*np.pi))*random.random()
    #dka
    kouho=shokiti
    normlist=[0]*jougen
    norm=1
    #print(kouho)
    for tr in range(80):
        sabun=np.array([[0]*jougen]*jougen,dtype=complex)
        norm=0
        for i in range(jougen):
            for j in range(jougen):
                sabun[i,j]=kouho[i]-kouho[j]
        #print(sabun)
        for i in range(jougen):
            kansuu=dainyu(kouho[i],klist)
            newkei=1
            for j in range(jougen):
                if j!=i:
                    newkei*=sabun[i,j]
            
            #print(newkei,kansuu)
            kouho[j]-=kansuu/newkei/klist[0]/pow(newkei,jisuu/jougen-1)
            #print(kansuu/newkei)
            norm+=abs(kansuu)
            #print(kouho)
        
        if norm<seido:
            
            break
        #print(norm)
        #print(kouho)
    if np.isnan(kouho).any()==True:
        for i in range(len(kouho)):
            if np.isnan(kouho[i])==True:
                kouho[i]=0
    if norm>seido:
        for tr in range(80):
            norm=0
            for i in range(jougen):
                dif=dainyu(kouho[i],klist)
                kouho[i]-=dif/bibun(kouho[i],klist)
                norm+=abs(dif)
            if norm<seido:
                break
    if np.isnan(kouho).any()==True:
        for i in range(len(kouho)):
            if np.isnan(kouho[i])==True:
                kouho[i]=0        

    #print (kouho)
    for i in range(jougen):
        normlist[i]=abs(dainyu(kouho[i],klist))
        if norm>=seido and normlist[i]<1:
            for tr in range(60):
                #kouho[i]+=0
                kouho[i]-=dainyu(kouho[i],klist)/bibun(kouho[i],klist)/20
            normlist[i]=abs(dainyu(kouho[i],klist))
#    print(normlist)
    return kouho,norm,normlist
    

def kaikeisan(keisuu):
    motometai=keisuu
    seido=1.0e-5
    mrank=len(motometai)
    kailist,norm,norml=dkahou(motometai)
    #print(kailist,norm)
    kailen=len(kailist)
    d=np.zeros((kailen,kailen),dtype=complex)
    d[0,0]=kailist[0]
    for i in range(kailen-1):
        d[i,i]=kailist[i+1]
        d[i,i+1]=1
    b=np.eye(kailen,dtype=complex)*motometai[mrank-1]
    a=np.eye(kailen,dtype=complex)
    for i in range(1,mrank-1):
        a=a@d
        b+=a*motometai[mrank-1-i]
    
    normchouhuku=linalg.norm(b)
    #print(normchouhuku)
    kaikakutei=[]
    if normchouhuku>seido:
        d=np.eye(kailen,dtype=complex)*kailist[0]
        b=np.eye(kailen,dtype=complex)*motometai[mrank-1]
        a=np.eye(kailen,dtype=complex)
        for i in range(1,mrank):
            a=a@d
            b+=a*motometai[mrank-1-i]
#            print(b)
        normchouhuku=linalg.norm(b)
        row=0
        #print(normchouhuku)
        rowmax=kailen
        if normchouhuku>seido:
            for j in range(1,kailen):
                d=np.eye(kailen,dtype=complex)*kailist[j]
                b=np.eye(kailen,dtype=complex)*motometai[mrank-1]
                a=np.eye(kailen,dtype=complex)
                for i in range(1,mrank):
                    a=a@d
                    b=b+a*motometai[mrank-1-i]
                normchouhuku=linalg.norm(b)
                
                if normchouhuku>seido:
                    rowmax-=1
                else:
                    rowmax-=1
                    break
        
        #print(normchouhuku)
        for j in range(1,rowmax):
            if normchouhuku<seido:
                d[row,row+1]=1
                row+=1
            d[row,row]=kailist[j]
            
            #print(d)
            b=np.eye(kailen,dtype=complex)*motometai[mrank-1]
            a=np.eye(kailen,dtype=complex)
            for i in range(1,mrank):
                a=a@d
                b=b+a*motometai[mrank-1-i]
            normchouhuku=linalg.norm(b)
            #print(normchouhuku)
        for j in range(row):
            kaikakutei.append(d[j,j])
    else:
        kaikakutei=kailist
    #print(kaikakutei)
    
    return kaikakutei
    
def kakezan(keisuu1,keisuu2):
    kakutei=[0]*(len(keisuu1)+len(keisuu2)-1)
    for i in range(len(keisuu1)):
        for j in range(len(keisuu2)):
            kakutei[i+j]+=keisuu1[i]*keisuu2[j]
    return kakutei


def kontansaku(keisuu):
    kakuteikeisuu=[]
    kyuukai=keisuu
    subkyuukai=kyuukai
    flag=0
    tr0=0
    while flag==0 and tr0<100:
        tr0+=1
        print(tr0)
        if flag==1:
            return kakuteikeisuu
        kaikeisuu=kaikeisan(subkyuukai)
        #print(subkyuukai,'subkyu')
        warilist=[1]

        if kaikeisuu==[]:
            continue
        if len(kaikeisuu)>=1:
            for i in range(len(kaikeisuu)):
                for tr in range(100):
                    kaikeisuu[i]-=dainyu(kaikeisuu[i],subkyuukai)/bibun(kaikeisuu[i],subkyuukai)
                    if abs(dainyu(kaikeisuu[i],subkyuukai))<1.0e-14:
                        kakekeisuu=[1,-kaikeisuu[i]]
                        warilist=kakezan(warilist,kakekeisuu)
                        kakuteikeisuu.append(kaikeisuu[i])
                        break

        #print(kaikeisuu,warilist,'warimoto')
        if len(warilist)>1:
            subkyuukai,dummy1,dummy2=warizan(subkyuukai,warilist,len(subkyuukai)-1,len(warilist)-1)
        ##print(subkyuukai,'kai')
        print(len(subkyuukai))
        if len(subkyuukai)==2:
            tmpkeisuu=kyuukai
            for i in range(len(kakuteikeisuu)):
                warilist=[1,-kakuteikeisuu[i]]

                tmpkeisuu,dummy1,dummy2=warizan(tmpkeisuu,warilist,len(tmpkeisuu)-1,1)
            kakuteikeisuu.append(-tmpkeisuu[1]/tmpkeisuu[0])
            flag=1
            break
        if len(subkyuukai)==1:
            kakuteikeisuu.append(0)
            flag=1
            break
        if subkyuukai==[] :
            flag=1
            break   
    return kakuteikeisuu



def gyouretukeisuu(gyouretu):
    gyousuu=gyouretu.shape[0]
    if gyousuu!=gyouretu.shape[1]:
        print('not seihougyouretu error')
            
    a=np.eye(gyousuu,dtype=complex)
    b=np.zeros((gyousuu,gyousuu),dtype=complex)
    for i in range(gyousuu):
        m=np.ones((gyousuu,1))
    
        a=a@gyouretu
    
        b[0:gyousuu,gyousuu-1-i:gyousuu-i]=a@m
    #print(b,'b')
    #print(n)
    p=-b**-1@np.ones((gyousuu,1),dtype=complex)#a,b,1
    print(p)
    keisuu=[0]*(gyousuu+1)
    for i in range(gyousuu):
        keisuu[i]=p[i,0]/p[0,0]
    keisuu[gyousuu]=1
    return keisuu

n=[[3,4,1,0],[1,4,-3,2],[-2,-8,1,5],[2,3,1,9]]
print(np.matrix(n),'matrix')
keisuu=gyouretukeisuu(np.matrix(n).reshape(4,4))
keisuu=[1,4,2,8,5,7,1,4,2,8,5,7,1,5,9+1.0j,3,9,2,-4,5,4,6]
print(keisuu,'keisuu')
konlist1=kontansaku(keisuu)
konlist2=[0]*len(konlist1)
for i in range(len(konlist1)):
    konlist2[i]=dainyu(konlist1[i],keisuu)
print(konlist1,konlist2)