もう少し修正を施したプログラムで因数分解の時間を測ってみる。
前回のコードから、もう少し修正をして、500項までの多項式であれば何分か単位の時間でほぼ全ての解に因数分解することができるようになりました。
前回のプログラムのポイントは、
xの何乗の計算を、できるだけ少ない回数で計算するために、xの(2のn乗)乗を先に計算してから全ての指数分のxの項を計算すること。
ニュートン法を一度に複数の項について適用できるようにするために、30個くらいの数値について、複素数平面でいう角度を制限したニュートン法を適用したこと。
解を求める毎に、元の式の係数を組み立て除法して、求まった解を含む項を除く作業を追加したこと。
この3点になります。ちなみに、修正の結果、速度を稼ぐために解の正確性を10項程切り捨てた方が、他の解の正確性も高まるとわかったため、今回は最後の方の解はあまり正しくない結果を出すようになっています。
100項のランダム多項式の解計算の結果です。12秒ほどで計算が可能です。
400項ほどの多項式の解の計算結果です。100秒ほどです。
400項のランダム多項式の解の計算結果です。400秒ほど。項数が増えると計算量が増えて、当たり前に時間がかかるようになるのですが、案外計算量は増えていないようですね。
コードを置いておきます。
import numpy as np
import scipy.linalg as linalg
import random
import math
import threading
import copy
import concurrent.futures
import multiprocessing
from multiprocessing import Pool
class Keisan():
lock=threading.Lock()
def __init___(self,name):
self.name='keisan'
def warizan(self,keisuu,bunbo,krank,brank):
with self.lock:
p=keisuu
b=bunbo
while b[0]==0 and len(b)>1:
del b[0]
brank-=1
jisuuk=krank
jisuub=brank
if krank<brank:
ci=[0]
ama=keisuu
else:
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 kakezan(self,keisuu1,keisuu2):
with self.lock:
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 bibunkeisuu(self,keisuu):
pkei=keisuu
rank=len(pkei)
kekka=[0]*(rank-1)
for i in range(rank-1):
kekka[i]=pkei[i]*(rank-i-1)
return kekka
def tasizan(self,keisuu1,keisuu2):
with self.lock:
if len(keisuu1)>=len(keisuu2):
rank=len(keisuu1)
kekka=[0]*rank
for i in range(rank):
kekka[i]=keisuu1[i]
for i in range(len(keisuu2)):
kekka[rank-i-1]+=keisuu2[len(keisuu2)-1-i]
else:
rank=len(keisuu2)
kekka=[0]*rank
for i in range(rank):
kekka[i]=keisuu2[i]
for i in range(len(keisuu1)):
kekka[rank-i-1]+=keisuu1[len(keisuu1)-1-i]
return kekka
class ButterEnzan():
def __init__(self,length):
self.lock=threading.Lock()
self.value=length
self.keta=length
self.log2bord=int(np.log2(length))+1
#print(self.log2bord)
self.yousosu=2**(self.log2bord-1)+1
self.bangou=[0]*(self.yousosu*2-1)
self.ruix=[0]*(self.log2bord+1)
self.nii=[2**i-1 for i in range(self.log2bord+1)]
#print(self.log2bord)
for i in range(self.log2bord):
self.bangou[self.nii[i]]=2**(self.log2bord-i-1)
for j in range(self.nii[i]):
self.bangou[j+self.nii[i]+1]=self.bangou[self.nii[i]]+self.bangou[j]
self.bangou=np.array([int(self.bangou[i])for i in range(len(self.bangou))])
#print(self.basho)
#print(self.bangou)
def dainyu2(self,xval,keisuu0):#maxlength<=self.ke
with self.lock:
x=xval
pkei=copy.copy(keisuu0)
maxketa=len(keisuu0)-1
if maxketa>self.keta:
print('error')
elif maxketa==0:
return keisuu0[0]
elif maxketa==0:
return keisuu0[0]*x+keisuu0[1]
maxlog2=int(np.log2(maxketa))+1
ruilist=[0]*self.yousosu
kekka=pkei[maxketa]
#print(pkei)
#print(maxlog2,self.log2bord,self.ruix,pkei,maxketa)
for i in range(self.log2bord):
k=self.nii[i]
if self.bangou[k]>maxketa:
continue
ruilist[k]=x**(2**(self.log2bord-i-1))
kekka+=pkei[maxketa-self.bangou[k]]*ruilist[k]
for j in range(k+1,2*k+1):
if i<self.log2bord-1 and self.bangou[j]<=maxketa:
ruilist[j]=ruilist[k]*ruilist[j-k-1]
kekka+=pkei[maxketa-self.bangou[j]]*ruilist[j]
elif self.bangou[j]<=maxketa and i==self.log2bord-1:
kekka+=pkei[maxketa-self.bangou[j]]*ruilist[k]*ruilist[j-k-1]
# for i in range(maxketa+1):
# kekka+=pkei[maxketa-self.bangou[self.basho[i]]]*ruilist[self.basho[i]]
#print(ruilist,pkei,self.ruix)
return kekka
def bibun2(self,xval,keisuu0):#maxlength<=self.keta
with self.lock:
x=xval
pkei=copy.copy(keisuu0)
maxketa=len(keisuu0)-1
if maxketa>self.keta:
print('error')
elif maxketa==0:
return 0
elif maxketa==1:
return keisuu0[0]
maxlog2=int(np.log2(maxketa))+1
ruilist=[0]*self.yousosu
kekka=pkei[maxketa-1]
#print(self.log2bord,maxlog2)
for i in range(self.log2bord):
k=self.nii[i]
if self.bangou[k]>maxketa-1:
continue
ruilist[k]=x**(2**(self.log2bord-i-1))
kekka+=pkei[maxketa-self.bangou[k]-1]*ruilist[k]*(self.bangou[k]+1)
for j in range(k+1,2*k+1):
if i<self.log2bord-1 and self.bangou[j]<=maxketa-1:
ruilist[j]=ruilist[k]*ruilist[j-1-k]
kekka+=pkei[maxketa-self.bangou[j]-1]*ruilist[j]*(self.bangou[j]+1)
elif self.bangou[j]<=maxketa-1 and i==self.log2bord-1:
kekka+=pkei[maxketa-self.bangou[j]-1]*ruilist[k]*ruilist[j-1-k] *(self.bangou[j]+1)
#print(ruilist,ruix,maxketa,pkei)
return kekka
def dainyuandbibun(self,xval,keisuu0):
with self.lock:
x=xval
pkei=copy.copy(keisuu0)
maxketa=len(keisuu0)-1
if maxketa>self.keta:
print('error')
elif maxketa==0:
return keisuu0[0],0
elif maxketa==1:
return keisuu0[0]*xval+keisuu0[1],keisuu0[0]
maxlog2=int(np.log2(maxketa))+1
ruilist=[0]*self.yousosu
kekkab=pkei[maxketa-1]
kekka=pkei[maxketa]
#print(self.log2bord,maxlog2)
for i in range(self.log2bord):
k=self.nii[i]
if self.bangou[k]>maxketa:
continue
ruilist[k]=x**(2**(self.log2bord-i-1))
if self.bangou[k]<maxketa:
kekkab+=pkei[maxketa-self.bangou[k]-1]*ruilist[k]*(self.bangou[k]+1)
kekka+=pkei[maxketa-self.bangou[k]]*ruilist[k]
for j in range(k+1,2*k+1):
if i<self.log2bord-1 and self.bangou[j]<=maxketa:
ruilist[j]=ruilist[k]*ruilist[j-k-1]
if self.bangou[j]<maxketa:
kekkab+=pkei[maxketa-self.bangou[j]-1]*ruilist[j]*(self.bangou[j]+1)
kekka+=pkei[maxketa-self.bangou[j]]*ruilist[j]
elif self.bangou[j]<=maxketa and i==self.log2bord-1:
if self.bangou[j]<maxketa:
kekkab+=pkei[maxketa-self.bangou[j]-1]*ruilist[k]*ruilist[j-k-1] *(self.bangou[j]+1)
kekka+=pkei[maxketa-self.bangou[j]]*ruilist[k]*ruilist[j-k-1]
#print(ruilist,ruix,maxketa,pkei)
return kekka,kekkab
def patrol_tansaku(keisuu0,rate,bunkatusuu,enzan):
#with threading.Lock():
kotaelist=[]
print('tansaku start')
keisuu=copy.copy(keisuu0)
enz=enzan
if len(keisuu)<=1:
print('keisuu ended')
if len(keisuu)==0:
return []
return []
for j in range(bunkatusuu+1):
rad=j/bunkatusuu*2*np.pi+np.random.uniform(-np.pi/bunkatusuu,np.pi/bunkatusuu)+rate*8/np.pi
flag=0
kakudop=complex(np.cos(rad),np.sin(rad))
kakudon=complex(np.cos(rad),-np.sin(rad))
x=0
for i in range(40):
fx,bx=enz.dainyuandbibun(x,keisuu)
sabun=(fx/bx)*kakudon/2
if sabun.real>-abs((x*kakudon).real):
x-=sabun.real*kakudop
else:
flag=1
#print(x,fx,sabun,bx)
if abs(sabun.real)<0.1:
#print(x,fx,bx,sabun,'sabun0')
if abs(fx)<1.0e-2:
kotaelist.append(x)
#print('appended')
flag=1
break
else:
for k in range(60):
sabun=(fx/bx)/2
x-=sabun
fx,bx=enz.dainyuandbibun(x,keisuu)
#sabun=fx/bx
#print(x,fx,bx,sabun,'sabun')
if abs(fx)<1.0e-2:
#print('appended')
flag=1
break
kotaelist.append(x)
if flag==1:
flag=0
break
if kotaelist==[]:
kotaelist==[complex(random.random(),random.random())]
return kotaelist
def kenshou_kotae(kotaelist0,keisuu0,xkeisuu0,kakuteil,enzan):
keisan=Keisan()
kakuteilist=copy.copy(kakuteil)
kotaelist=copy.copy(kotaelist0)
enz=enzan
keisuu=copy.copy(keisuu0)
xkeisuu=copy.copy(xkeisuu0)
for i in range(len(kotaelist)):
for j in range(80):
fx0,fx0b=enz.dainyuandbibun(kotaelist[i],keisuu)
fx0abs=abs(fx0)
#if len(kotaelist)>0:
#print(kotaelist[i],fx0,'hantei',(fx0/fx0b),len(kakuteilist))
if fx0abs<1.0e-7:# and fx1abs<1.0e-12:
kakuteilist.append(kotaelist[i])
keisuu,amari1,dummy2=keisan.warizan(keisuu,[1,-kotaelist[i]],len(keisuu)-1,1)
#amari=[0]*len(keisuu)
# for k in range(len(keisuu)):
# amari[len(keisuu)-k-1]=-amari1[0]/kotaelist[i]**(k+1).
# keisuu=keisan.tasizan(keisuu,amari)
#wari=np.linalg.norm(keisuu)
#keisuu,dummy1,dummy2=keisan.warizan(keisuu,[wari],len(keisuu)-1,0)
print(kotaelist[i],fx0,len(kakuteilist),'. succeed. ')
break
if np.isnan(kotaelist[i]) or np.isinf(kotaelist[i]):
break
if fx0abs>1.0e-5:
kotaelist[i]-=(fx0/fx0b)/2
continue
if (fx0abs<=1.0e-5):
kotaelist[i]-=(fx0/fx0b)/2
continue
#if (fx0abs>100 or fx1abs>100) and j>20:
#break
#if fx0abs>1 and fx1abs>1 and j>50:
#break
return kakuteilist,keisuu
def seiri_kakutei(kakuteilist0,keisuu0,xkeisuu0,enzan,prelength0):
#with threading.Lock():
keisan=Keisan()
enz=enzan
xkeisuu=copy.copy(xkeisuu0)
kakuteilist=copy.copy(kakuteilist0)
keisuu=copy.copy(keisuu0)
prelength=prelength0
printamari=0
kakutei2=[]
kakutei3=[]
print('organizing...')
if np.isnan(any(keisuu)):
keisuu=copy.copy(xkeisuu)
return [],keisuu,prelength,[],xkeisuu
if kakuteilist==[]:
print('blank list')
return [],keisuu,0,[],xkeisuu
if True:
keisuu=copy.copy(xkeisuu)
bunbolist=[1]
for i in range(len(kakuteilist)):
fx0,fx0b=enz.dainyuandbibun(kakuteilist[i],keisuu)
fx1,fx1b=enz.dainyuandbibun(kakuteilist[i],xkeisuu)
for j in range(60):
if abs(fx1)<1.0e-11 and abs(fx0)<1.0e-4 and np.isnan(kakuteilist[i])==False:
kakutei3.append(kakuteilist[i])
xkeisuu,amarilist,dummy2=keisan.warizan(xkeisuu,[1,-kakuteilist[i]],len(xkeisuu)-1,1)
#wari=xkeisuu[len(xkeisuu)-1]
#xkeisuu,dummy1,dummy2=keisan.warizan(xkeisuu,[wari],len(xkeisuu)-1,0)
keisuu,amarilist,dummy2=keisan.warizan(keisuu,[1,-kakuteilist[i]],len(keisuu)-1,1)
#print(amarilist[0],'amari')
#wari=np.linalg.norm(keisuu)
#keisuu,dummy1,dummy2=keisan.warizan(keisuu,[wari],len(keisuu)-1,0)
break
if abs(fx0)>2.0e-15:
kakuteilist[i]-=(fx0/fx0b)/2
fx0,fx0b=enz.dainyuandbibun(kakuteilist[i],keisuu)
fx1,fx1b=enz.dainyuandbibun(kakuteilist[i],xkeisuu)
else:
kakutei2.append(kakuteilist[i])
keisuu,amarilist,dummy2=keisan.warizan(keisuu,[1,-kakuteilist[i]],len(keisuu)-1,1)
printamari=amarilist[0]
#print(amarilist[0],'amari')
#amari=[0]*len(keisuu)
# for k in range(len(keisuu)):
# amari[len(keisuu)-k-1]=-amarilist[0]/kakuteilist[i]**(k+1)
# keisuu=keisan.tasizan(keisuu,amari)
#wari=np.linalg.norm(keisuu)
#keisuu,dummy1,dummy2=keisan.warizan(keisuu,[wari],len(keisuu)-1,0)
#print(enz.dainyu2(kakuteilist[i],xkeisuu))
break
#print(kakutei2,'kakutei')
if prelength==len(kakuteilist):
wari=np.linalg.norm(keisuu)
keisuu,dummy1,dummy2=keisan.warizan(keisuu,[wari],len(keisuu)-1,0)
print('organizing end')
prelength=len(kakuteilist)
return kakutei2,keisuu,prelength,kakutei3,xkeisuu
def kyuukai(x0keisuu):
bunkatusuu=60
maxtrytimes=100
keisan=Keisan()
keisuu=np.array(x0keisuu,dtype='complex128')
rank=len(keisuu)-1
kakuteilist=[]
enz=ButterEnzan(rank)
tanenz=ButterEnzan(rank)
kenenz=ButterEnzan(rank)
seirienz=ButterEnzan(rank)
prelength=0
kakutei2=[]
if keisuu[rank]==0:
while keisuu[rank]==0:
warilist=[1,0]#warizan
keisuu,dummy1,dummy2=keisan.warizan(keisuu,warilist,len(keisuu)-1,1)
kakuteilist.append(0)
keisuu,dummy1,dummy2=keisan.warizan(keisuu,[keisuu[len(keisuu)-1]],len(keisuu)-1,0)
xkeisuu=copy.copy(keisuu)
keisuukotaeawase=copy.copy(keisuu)
flag=0
kotaelist=[]
r=0
flag=0
tmp=1
#print(xkeisuu[0:10],keisuu[0:10])
kotaelist=patrol_tansaku(xkeisuu,0,bunkatusuu,tanenz)
#print(xkeisuu[0:10],keisuu[0:10])
kakuteilist1=[]
keisuu1=copy.copy(keisuu)
kakuteilist=[]
future_list=[]
#kakuteilist1,keisuu1=kenshou_kotae(kotaelist,keisuu,xkeisuu,kakuteilist,kenenz)
for tr in range(maxtrytimes):
kakuteilist1,keisuu1=kenshou_kotae(kotaelist,keisuu,xkeisuu,kakuteilist,kenenz)
#kakuteilist,keisuu,prelength=seiri_kakutei(kakuteilist,keisuu,xkeisuu,seirienz,prelength)
with concurrent.futures.ThreadPoolExecutor() as executor:
future1=executor.submit(patrol_tansaku,keisuu,tr/maxtrytimes,bunkatusuu,tanenz)
#future2=executor.submit(kenshou_kotae,kotaelist,keisuu,xkeisuu,kakuteilist,kenenz)
future3=executor.submit(seiri_kakutei,kakuteilist1,keisuu1,xkeisuu,seirienz,prelength)
future_list=[future1.result(),0,future3.result()]
print(xkeisuu[0:10],keisuu[0:10])
kotaelist=future_list[0]
#kakuteilist1,keisuu1=future_list[1]
kakuteilist,keisuu,prelength,tmpkakutei,xkeisuu=future_list[2]
if len(tmpkakutei)>0:
kakutei2.extend(tmpkakutei)
print(len(kakutei2),len(xkeisuu))
#print(kakuteilist,kotaelist)kakuteilist1,keisuu1=kenshou_kotae(kotaelist,keisuu,xkeisuu,kakuteilist,kenenz)
print(tr,' st try')
#print(xkeisuu[0:10],keisuu[0:10],keisuu[len(keisuu)-1])
if len(kakutei2)>=rank:
break
if len(kotaelist)>0:
#keisuu,dummy1,dummy2=warizan(keisuu,[complex(np.sin(1),np.cos(1))],len(keisuu)-1,0)
print(enz.dainyuandbibun(kotaelist[0],keisuu),'dainyu')
#print(kotaelist[0],enz.dainyu2(kotaelist[0],keisuu),enz.dainyu2(kotaelist[0],xkeisuu),enz.bibun2(kotaelist[0],keisuu),enz.bibun2(kotaelist[0],xkeisuu),len(kakuteilist),'. . ',tr)
else:
print('blanklist')
norm=0
for i in range(len(kakutei2)):
norm+=abs(enz.dainyu2(kakutei2[i],keisuukotaeawase))
print(enz.dainyu2(kakutei2[i],keisuukotaeawase))
print(norm,'norm')
return kakutei2
kyuukai([1,3,2])ちなみに、本物のプログラムで固有値や多項式の解を解く解き方は、1000項目くらいまでは行列の固有値問題を早く解く方法が確立されており、それを使うようです。それより大きい方程式は、ブロック対角化やコンパニオン行列を使った計算方法を使うようです。1000項目くらいまでは、数秒で解けるからそっちの方が早い。。