多項式の解の推定のためのプログラムをpythonで作る。
引き続き多項式の解を求めるためのプログラムを書き続けていました。ある程度の数の解を求めるためのプログラムができたので、書き出しておきます。
import numpy as np
import scipy.linalg as linalg
import random
import math
def warizan(keisuu,bunbo,krank,brank):
p=keisuu
b=bunbo
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(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 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
class ButterEnzan():
def __init__(self,length):
self.value=length
self.keta=length
self.log2bord=math.ceil(np.log2(length))
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)]
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))])
def dainyu2(self,xval,keisuu0):#maxlength<=self.keta
x=xval
pkei=keisuu0
maxketa=len(keisuu0)-1
maxlog2=math.ceil(np.log2(maxketa))
if maxketa>self.keta:
print('error')
ruilist=[0]*self.yousosu
kekka=pkei[maxketa]
for i in range(self.log2bord):
k=self.nii[i]
ruilist[k]=x**(2**(self.log2bord-i-1))
kekka+=pkei[maxketa-self.bangou[k]]*ruilist[k]
for j in range(k):
if i<self.log2bord-1 and self.bangou[k+1+j]<=maxketa:
ruilist[k+1+j]=ruilist[k]*ruilist[j]
kekka+=pkei[maxketa-self.bangou[k+1+j]]*ruilist[k+1+j]
elif self.bangou[k+1+j]<=maxketa and i==self.log2bord-1:
kekka+=pkei[maxketa-self.bangou[k+1+j]]*ruilist[k]*ruilist[j]
return kekka
def bibun2(self,xval,keisuu0):
x=xval
pkei=keisuu0
maxketa=len(keisuu0)-1
maxlog2=math.ceil(np.log2(maxketa))
if maxketa>self.keta:
print('error')
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:
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):
if i<self.log2bord-1 and self.bangou[k+1+j]<=maxketa-1:
ruilist[k+1+j]=ruilist[k]*ruilist[j]
kekka+=pkei[maxketa-self.bangou[k+1+j]-1]*ruilist[k+1+j]*(self.bangou[k+1+j]+1)
elif self.bangou[k+1+j]<=maxketa-1 and i==self.log2bord-1:
kekka+=pkei[maxketa-self.bangou[k+1+j]-1]*ruilist[k]*ruilist[j] *(self.bangou[k+1+j]+1)
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 tasizan(keisuu1,keisuu2):
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
def tasizan(keisuu1,keisuu2):
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
def kyuukai(x0keisuu):
bunkatusuu=30
maxtrytimes=100
xkeisuu=np.array(x0keisuu,dtype='complex128')
keisuu=xkeisuu
rank=len(xkeisuu)-1
xkai=[random.random() for i in range(rank)]
kyori=np.zeros((rank,rank),dtype=complex)
kakuteilist=[]
enz=ButterEnzan(rank)
if xkeisuu[rank]==0:
while xkeisuu[rank]==0:
warilist=[1,0]#warizan
keisuu,dummy1,dummy2=warizan(keisuu,warilist,len(keisuu)-1,1)
kakuteilist.append(0)
keisuu,dummy1,dummy2=warizan(keisuu,[keisuu[len(keisuu)-1]],len(keisuu)-1,0)
flag=0
bkeisuu=bibunkeisuu(keisuu)
kotaelist=[]
r=0
flag=0
for tr in range(maxtrytimes):
kotaelist=[]
if len(keisuu)<=1:
print('keisuu ended')
break
for j in range(bunkatusuu+1):
flag=0
rad=j/bunkatusuu*2*np.pi+np.random.uniform(-np.pi/bunkatusuu,np.pi/bunkatusuu)
kakudop=complex(np.cos(rad),np.sin(rad))
kakudon=complex(np.cos(rad),-np.sin(rad))
x=0
for i in range(100):
# fx=dainyu(x,keisuu)
fx=enz.dainyu2(x,keisuu)
# bx=bibun(x,keisuu)
bx=enz.bibun2(x,keisuu)
sabun=fx/bx*kakudon
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-5:
kotaelist.append(x)
#print('appended')
flag=1
break
else:
for k in range(100):
sabun=fx/bx
x-=sabun
# fx=dainyu(x,keisuu)
# bx=bibun(x,keisuu)
fx=enz.dainyu2(x,keisuu)
bx=enz.bibun2(x,keisuu)
#sabun=fx/bx
#print(x,fx,bx,sabun,'sabun')
if abs(fx)<1.0e-5:
kotaelist.append(x)
#print('appended')
flag=1
break
if flag==1:
flag=0
break
i=0
j=0
k=0
if len(kotaelist)>=0:
#keisuu,dummy1,dummy2=warizan(keisuu,[complex(np.sin(1),np.cos(1))],len(keisuu)-1,0)
print(dainyu(kotaelist[0],keisuu),dainyu(kotaelist[0],xkeisuu),'dainyu')
for i in range(len(kotaelist)):
for j in range(100):
# fx0=dainyu(kotaelist[i],keisuu)
# fx1=dainyu(kotaelist[i],xkeisuu)
fx0=enz.dainyu2(kotaelist[i],keisuu)
fx1=enz.dainyu2(kotaelist[i],xkeisuu)
fx0abs=abs(fx0)
fx1abs=abs(fx1)
# fx0b=bibun(kotaelist[i],keisuu)
# fx1b=bibun(kotaelist[i],xkeisuu)
fx0b=enz.bibun2(kotaelist[i],keisuu)
fx1b=enz.bibun2(kotaelist[i],xkeisuu)
#if tr>0:
#print(kotaelist[i],fx0,fx1,'hantei',fx0abs<1.0e-7,fx1abs<1.0e-13,len(kakuteilist),tr,fx0/fx0b,fx1/fx1b)
if fx0abs<1.0e-11 and fx1abs<1.0e-13:
kakuteilist.append(kotaelist[i])
keisuu,dummy1,dummy2=warizan(keisuu,[1,-kotaelist[i]],len(keisuu)-1,1)
keisuu,dummy1,dummy2=warizan(keisuu,[keisuu[len(keisuu)-1]],len(keisuu)-1,0)
print(kotaelist[i],fx0,fx1,fx0abs<1.0e-11,fx1abs<1.0e-13,len(kakuteilist),'. succeed. ',tr)
break
if np.isnan(kotaelist[i]) or np.isinf(kotaelist[i]):
break
if fx0abs<1 and fx0abs>1.0e-11:
#if fx0abs/abs(fx0b)<abs(kotaelist[i])*9.9e-12:
#kotaelist[i]-=fx0/fx0b/fx0abs*abs(fx0b)*abs(kotaelist[i])*9.9e-12
kotaelist[i]-=fx0/fx0b
continue
if fx1abs<1 and fx1abs>1.0e-30:
#if fx1abs/abs(fx1b)<9.9e-13*abs(kotaelist[i]):
#kotaelist[i]-=fx1/fx1b/ fx1abs*abs(fx1b)*9.9e-13*abs(kotaelist[i])
kotaelist[i]-=fx1/fx1b
continue
if fx0abs>100 and fx1abs>100 and j>20:
break
if fx0abs>1 and fx1abs>1 and j>50:
break
#if fx0abs>1.0e-10 and j>50:
#break
return kakuteilist
x=[np.random.uniform(-1,1)*3]*500
x[0]=1.0
import time
start=time.time()
lists=kyuukai(x)
end=time.time()
print('process time is ',end-start)
kotae=[dainyu(lists[i],x) for i in range(len(lists))]
print(kotae[20:30],kotae[100:110],lists[0:10])
print('value length is ',len(lists))
想定した通りに解を高い精度で求めることができているようです。
今回は複素数で表される解xの複素数的なところでいう角度(偏角というらしいです。)を固定して、x＝0のところからニュートン法で近似を続けることによって、一度に5〜6個の解を探索することができるようになりました。
500係数がある式で、1600秒かけて150くらいの解しか見つけられなかったのは残念。
この記事が気に入ったらサポートをしてみませんか？