polyroots-fortran で多項式を解く

多項式のお手軽なソルバー。
こういった複素数関連のソルバを自力で実装すると複素数の丸め誤差が出てしまうのでかなりありがたいところです

setup

前準備として fpm をインストールしておく (ココを参照).

githubからソースコードを持ってきてビルド

git clone https://github.com/jacobwilliams/polyroots-fortran.git
cd polyroots-fortran/
fpm build --profile release
fpm install
sudo cp ~/.local/lib/libpolyroots-fortran.a /usr/local/lib/
sudo cp ~/.local/include/polyroots_module.mod /usr/local/include/
ls /usr/local/lib
ls /usr/local/include

計算

yahoo知恵袋に 超難問 自作問題 があるのでこれを解いてみる。

32x^5 + 16x^4 -32x^3 -12x^2 + 6x + 1 = 0

program test
use polyroots_module
implicit none

integer, parameter :: deg = 5
real(8) :: p(deg+1) = [32,16,-32,-12,6,1], zr(deg), zi(deg)
integer :: i, istat

call polyroots( deg, p, zr, zi, istat )

do i = 1, deg
    if(istat==0) write(*,'(2ES16.5)') zr(i), zi(i)
end do

end program test

• deg で次数を指定
• zrはリアルパート、ziはイマジナリーパート
• p に係数を入力
• call polyroots で計算
• istat = 0 で成功

status output:

• 0 : success
• -1 : leading coefficient is zero
• -2 : no roots found
• >0 : the number of zeros found

gfortran test.f90 -L/usr/local/lib -lpolyroots-fortran -llapack -lblas -I/usr/local/include

lapackとblasもリンクしておく (toml参照)

$ ./a.out
     8.41254E-01     0.00000E+00
    -9.59493E-01     0.00000E+00
    -6.54861E-01     0.00000E+00
     4.15415E-01     0.00000E+00
    -1.42315E-01     0.00000E+00

どうやら全て実数解になるらしい。 解は知恵袋の回答によると

x 〜 -0.959493
x 〜 -0.654861
x 〜 -0.142315
x 〜 0.415415
x 〜 0.841254

おそらくいい感じに計算できた。