! Copyright 2012 Jan Skowron & Andrew Gould ! ! Licensed under the Apache License, Version 2.0 (the "License"); ! you may not use this file except in compliance with the License. ! You may obtain a copy of the License at ! ! http://www.apache.org/licenses/LICENSE-2.0 ! ! Unless required by applicable law or agreed to in writing, software ! distributed under the License is distributed on an "AS IS" BASIS, ! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. ! See the License for the specific language governing permissions and ! limitations under the License. ! !-------------------------------------------------------------------! ! ! The authors also make this file available under the terms of ! GNU Lesser General Public License version 2 or any later version. ! (text of the LGPL licence 2 in NOTICE file) ! !-------------------------------------------------------------------! ! ! A custom in the scientific comunity is (regardless of the licence ! you chose to use or distribute this software under) ! that if this code was important in the scientific process or ! for the results of your scientific work, we kindly ask you for the ! appropriate citation of the Paper (Skowron & Gould 2012), and ! we would be greatful if you pass the information about ! the proper citation to anyone whom you redistribute this software to. ! !-------------------------------------------------------------------! ! ! No Subroutine ! ! 1 cmplx_roots_gen - general polynomial solver, works for random degree, not as fast or robust as cmplx_roots_5 ! 2 cmplx_roots_5 - complex roots finding algorithm taylored for 5th order polynomial (with failsafes for polishing) ! 3 sort_5_points_by_separation - sorting of an array of 5 points, 1st most isolated, 4th and 5th - closest ! 4 sort_5_points_by_separation_i - sorting same as above, returns array of indicies rather than sorted array ! 5 find_2_closest_from_5 - finds closest pair of 5 points ! 6 cmplx_laguerre - Laguerre's method with simplified Adams' stopping criterion ! 7 cmplx_newton_spec - Newton's method with stopping criterion calculated every 10 steps ! 8 cmplx_laguerre2newton - three regime method: Laguerre's, Second-order General method and Newton's ! 9 solve_quadratic_eq - quadratic equation solver ! 10 solve_cubic_eq - cubic equation solver based on Lagrange's method ! 11 divide_poly_1 - division of the polynomial by (x-p) ! ! fortran 77 code ! ! Paper: Skowron & Gould 2012 ! "General Complex Polynomial Root Solver and Its Further Optimization for Binary Microlenses" ! ! for a full text see: ! http://www.astrouw.edu.pl/~jskowron/cmplx_roots_sg/ ! or http://arxiv.org/find/astro-ph ! or http://www.adsabs.harvard.edu/abstract_service.html ! see also file NOTICE and LICENSE ! ! ver. 2012.03.03 initial ! ver. 2014.03.12 bug fix ! ver. 2016.01.21 bug fix ! ver. 2016.04.28 bug fix ! !-------------------------------------------------------------------! ! _1_ CMPLX_ROOTS_GEN ! !-------------------------------------------------------------------! subroutine cmplx_roots_gen(roots, poly, degree, * polish_roots_after, use_roots_as_starting_points) ! This subroutine finds roots of a complex polynomial. ! It is general, however less fast or robust than cmplx_roots_5 ! which contains failsafe checks in the polishing stage, but is ! designed only for 5th order polynomials. ! It uses a new dynamic root finding algorithm (see the Paper). ! ! It can use Laguerre's method (subroutine cmplx_laguerre) ! or Laguerre->SG->Newton method (subroutine ! cmplx_laguerre2newton - this is default choice) to find ! roots. It divides polynomial one by one by found roots. At the ! end it finds last root from Viete's formula for quadratic ! equation. Finally, it polishes all found roots using a full ! polynomial and Newton's or Laguerre's method (default is ! Laguerre's - subroutine cmplx_laguerre). ! You can change default choices by commenting out and uncommenting ! certain lines in the code below. ! ! Note: ! - we solve for the last root with Viete's formula rather ! than doing full Laguerre step (which is time consuming ! and unnecessary) ! - we do not introduce any preference to real roots ! - in Laguerre implementation we omit unneccesarry calculation of ! absolute values of denominator ! - we do not sort roots. If you need to sort ! roots - we have provided sorting subroutine called: ! sort_5_points_by_separation, which sorts points from most ! isolated to most close. Algorithm in this routine can be ! easily used for number of points different than 5. ! implicit none ! roots - array which will hold all roots that had been found. ! If the flag 'use_roots_as_starting_points' is set to ! .true., then instead of point (0,0) we use value from ! this array as starting point for cmplx_laguerre ! poly - is an array of polynomial cooefs, length = degree+1, ! poly(1) is a constant term: ! 1 2 3 ! poly(1) x^0 + poly(2) x^1 + poly(3) x^2 + ... ! degree - degree of the polynomial and size of 'roots' array ! polish_roots_after - after all roots have been found by dividing ! original polynomial by each root found, ! you can opt in to polish all roots using full ! polynomial ! use_roots_as_starting_points - usually we start Laguerre's ! method from point (0,0), but you can decide to use the ! values of 'roots' array as starting point for each new ! root that is searched for. This is useful if you have ! very rough idea where some of the roots can be. ! integer degree ! intent(in) double complex poly(degree+1) ! intent(in), coeffs of the polynomial double complex roots(degree) ! intent(inout) logical polish_roots_after, use_roots_as_starting_points ! intent(in) double complex poly2(degree+1) double complex zero parameter(zero=(0d0,0d0)) integer i, n, iter logical success double complex coef, prev do i=1,degree+1 poly2(i)=poly(i) enddo ! initialize starting points if(.not.use_roots_as_starting_points) then do i=1,degree roots(i)=zero enddo endif ! skip small degree polynomials from doing Laguerre's method if(degree.le.1)then if(degree.eq.1) roots(1)=-poly(1)/poly(2) return endif do n=degree, 3, -1 ! find root with Laguerre's method !call cmplx_laguerre(poly2, n, roots(n), iter, success) ! or ! find root with (Laguerre's method -> SG method -> Newton's method) call cmplx_laguerre2newton(poly2, n, roots(n), * iter, success, 2) if(.not.success) then roots(n)=zero call cmplx_laguerre(poly2, n, roots(n), iter, success) endif ! divide the polynomial by this root coef=poly2(n+1) do i=n,1,-1 prev=poly2(i) poly2(i)=coef coef=prev+roots(n)*coef enddo ! coef now holds a remainder - should be close to 0 enddo ! find all but last root with Laguerre's method !call cmplx_laguerre(poly2, 2, roots(2), iter, success) ! or call cmplx_laguerre2newton(poly2, 2, roots(2), * iter, success, 2) if(.not.success) then call solve_quadratic_eq(roots(2),roots(1),poly2) else ! calculate last root from Viete's formula roots(1)=-(roots(2)+poly2(2)/poly2(3)) endif if(polish_roots_after)then do n=1, degree ! ! polish roots one-by-one with a full polynomial call cmplx_laguerre(poly, degree, roots(n), iter, success) ! or !call cmplx_newton_spec(poly, degree, roots(n), iter, success) enddo endif return end !-------------------------------------------------------------------! ! _2_ CMPLX_ROOTS_5 ! !-------------------------------------------------------------------! subroutine cmplx_roots_5 * (roots, first_3_roots_order_changed, poly, polish_only) implicit none ! Subroutine finds or polishes roots of a complex polynomial ! (degree=5) ! This routine is especially tailored for solving binary lens ! equation in form of 5th order polynomial. ! ! Use of this routine, in comparission to 'cmplx_roots_gen' can yield ! consideribly faster code, because it makes polishing of the roots ! (that come in as a guess from previous solutions) secure by ! implementing additional checks on the result of polishing. ! If those checks are not satisfied then routine reverts to the ! robust algorithm. These checks are designed to work for 5th order ! polynomial originated from binary lens equation. ! ! Usage: ! ! polish_only == false - I do not know the roots, routine should ! find them from scratch. At the end it ! sorts roots from the most distant to closest. ! Two last roots are the closest (in no particular ! order). ! polish_only = true - I do know the roots pretty well, for example ! I have changed the coefficiens of the polynomial ! only a bit, so the two closest roots are ! most likely still the closest ones. ! If the output flag 'first_3_roots_order_changed' ! is returned as 'false', then first 3 returned roots ! are in the same order as initialy given to the ! routine. The last two roots are the closest ones, ! but in no specific order (!). ! If 'first_3_roots_order_changed' is 'true' then ! it means that all roots had been resorted. ! Two last roots are the closest ones. First is most ! isolated one. ! ! ! If you do not know the position of the roots just use flag ! polish_only=.false. In this case routine will find the roots by ! itself. ! Returns all five roots in the 'roots' array. ! ! poly - is an array of polynomial cooefs, length = degree+1 ! poly(1) x^0 + poly(2) x^1 + poly(3) x^2 + poly(4) x^3 + ... ! roots - roots of the polynomial ('out' and optionally 'in') ! ! ! Jan Skowron 2011 ! integer degree parameter(degree=5) ! specifying this speeds up optimization a little double complex roots(degree) ! intent(inout) logical first_3_roots_order_changed ! intent(out) double complex poly(degree+1) ! intent(in) logical polish_only ! intent(in) !------------------------------------------------------------------ ! double complex remainder, roots_robust(degree) double precision d2min integer iter, m, root4, root5, kk, go_to_robust, i, i2, loops double complex poly2(degree+1) !integer sorted_ind(degree) double complex zero parameter(zero=(0d0,0d0)) logical succ do kk=1,degree roots_robust(kk) = roots(kk) enddo !--------------------------------------- go_to_robust=0 if(.not.polish_only) then do kk=1,degree roots(kk)=zero enddo go_to_robust=1 endif first_3_roots_order_changed=.false. do loops=1,3 111 continue ! ROBUST ! (we do not know the roots) if(go_to_robust>0) then if(go_to_robust>2)then ! something is wrong do kk=1,degree roots(kk)=roots_robust(kk) ! return not-polished roots, because polishing creates errors enddo return endif do kk=1,degree+1 poly2(kk)=poly(kk) ! copy coeffs enddo do m=degree,4,-1 ! find the roots one-by-one (until 3 are left to be found) call cmplx_laguerre2newton(poly2, m, roots(m), * iter, succ, 2) if(.not.succ)then roots(m)=zero call cmplx_laguerre(poly2, m, roots(m), iter, succ) endif ! divide polynomial by this root call divide_poly_1(poly2, remainder, roots(m), poly2, m) enddo ! find last 3 roots with cubic euqation solver (Lagrange's method) call solve_cubic_eq(roots(1),roots(2),roots(3),poly2) ! all roots found ! sort roots - first will be most isolated, last two will be the closest call sort_5_points_by_separation(roots) ! copy roots in case something will go wrong during polishing do kk=1,degree roots_robust(kk)=roots(kk) enddo ! set flag, that roots have been resorted first_3_roots_order_changed=.true. endif ! go_to_robust>0 ! POLISH ! (we know the roots approximately, and we guess that last two are closest) !--------------------- do kk=1,degree+1 poly2(kk)=poly(kk) ! copy coeffs enddo do m=1,degree-2 !do m=1,degree ! POWN - polish only with Newton (option) ! polish roots with full polynomial call cmplx_newton_spec * (poly2, degree, roots(m), iter, succ) if(.not.succ)then ! go back to robust go_to_robust=go_to_robust+1 do kk=1,degree ! restart roots(kk)=zero enddo goto 111 ! go back endif enddo ! m=1,degree-2 ! comment out division and quadratic if you (POWN) polish with Newton only do m=1,degree-2 call divide_poly_1(poly2, remainder, roots(m), * poly2, degree-m+1) enddo ! last two roots are found with quadratic equation solver ! (this is faster and more robust, although little less accurate) call solve_quadratic_eq(roots(degree-1), * roots(degree ),poly2) ! all roots found and polished ! TEST ORDER ! test closest roots if they are the same pair as given to polish call find_2_closest_from_5(root4,root5,d2min, roots) ! check if the closest roots are not too close, this could happen ! when using polishing with Newton only, when two roots erroneously ! colapsed to the same root. This check is not needed for polishing ! 3 roots by Newton and using quadratic for the remaining two. ! If the real roots are so close indeed (very low probability), this will just ! take more time and the unpolished result be returned at the end !if(d2min<1d-18) then ! POWN - polish only with Newton ! go_to_robust=go_to_robust+1 ! POWN - polish only with Newton !else ! POWN - polish only with Newton if((root4.lt.degree-1).or.(root5.lt.degree-1)) then ! after polishing some of the 3 far roots become one of the 2 closest ones ! go back to robust if(go_to_robust.gt.0)then ! if came from robust ! copy two most isolated roots as starting points for new robust do i=1,degree-3 roots(degree-i+1)=roots_robust(i) enddo else ! came from users initial guess ! copy first 2 roots (except the closest ones) i2=degree do i=1,degree if((i.ne.root4).and.(i.ne.root5))then roots(i2)=roots(i) i2=i2-1 endif if(i2.le.3) exit ! do not copy those that will be done by cubic in robust enddo endif go_to_robust=go_to_robust+1 else ! root4 and root5 comes from the initial last pair ! most common case return endif !endif ! POWN - polish only with Newton !--------------------- enddo ! loops return end !-------------------------------------------------------------------! ! _3_ SORT_5_POINTS_BY_SEPARATION ! !-------------------------------------------------------------------! subroutine sort_5_points_by_separation(points) ! Sort array of five points ! Most isolated point will become the first point in the array ! The closest points will be the last two points in the array ! ! Algorithm works well for all dimensions. We put n=5 as ! a hardcoded value just for optimization purposes. implicit none integer n parameter(n=5) ! works for different n as well, but is faster for n as constant (optimization) double complex points(n) ! intent(inout) integer sorted_points(n) double complex savepoints(n) integer i call sort_5_points_by_separation_i(sorted_points, points) do i=1,n savepoints(i)=points(i) enddo do i=1,n points(i)=savepoints(sorted_points(i)) enddo return end !-------------------------------------------------------------------! ! _4_ SORT_5_POINTS_BY_SEPARATION_I ! !-------------------------------------------------------------------! subroutine sort_5_points_by_separation_i(sorted_points, points) ! Return index array that sorts array of five points ! Index of the most isolated point will appear on the first place ! of the output array. ! The indices of the closest 2 points will be at the last two ! places in the 'sorted_points' array ! ! Algorithm works well for all dimensions. We put n=5 as ! a hardcoded value just for optimization purposes. implicit none integer n parameter(n=5) ! works for different n as well, but is faster for n as constant (optimization) integer sorted_points(n) ! intent(out) double complex points(n) ! intent(in) double precision dmin, d1, d2, d double precision distances2(n,n) integer ki, kj, ind2, put double precision neigh1st(n), neigh2nd(n) double complex p do kj=1, n distances2(kj,kj)=1d100 enddo dmin=1d100 do kj=1, n do ki=1, kj-1 p=points(ki)-points(kj) d=dble(conjg(p)*p) distances2(ki,kj)=d distances2(kj,ki)=d enddo enddo ! find neighbours do kj=1, n neigh1st(kj)=1d100 neigh2nd(kj)=1d100 enddo do kj=1, n do ki=1, n d=distances2(kj,ki) if(d.lt.neigh2nd(kj))then if(d.lt.neigh1st(kj))then neigh2nd(kj)=neigh1st(kj) neigh1st(kj)=d else neigh2nd(kj)=d endif endif enddo enddo ! initialize sorted_points do ki=1,n sorted_points(ki)=ki enddo ! sort the rest 1..n-2 do kj=2,n d1=neigh1st(kj) d2=neigh2nd(kj) put=1 do ki=kj-1,1,-1 ind2=sorted_points(ki) d=neigh1st(ind2) if(d.ge.d1) then if(d.eq.d1)then if(neigh2nd(ind2)>d2)then put=ki+1 exit endif else put=ki+1 exit endif endif sorted_points(ki+1)=sorted_points(ki) enddo sorted_points(put)=kj enddo return end !-------------------------------------------------------------------! ! _5_ FIND_2_CLOSEST_FROM_5 ! !-------------------------------------------------------------------! subroutine find_2_closest_from_5(i1,i2,d2min, points) ! Returns indices of the two closest points out of array of 5 implicit none integer n parameter(n=5) ! will work for other n too, but it is faster with n as constant integer i1,i2 ! intent(out) double complex points(n) ! intent(in) double precision d2min ! intent(out) double precision d2min1, d2 integer i,j double complex p d2min1=1d100 do j=1,n do i=1,j-1 p=points(i)-points(j) d2=dble(conjg(p)*p) if(d2.le.d2min1)then i1=i i2=j d2min1=d2 endif enddo enddo d2min=d2min1 end !-------------------------------------------------------------------! ! _6_ CMPLX_LAGUERRE ! !-------------------------------------------------------------------! subroutine cmplx_laguerre * (poly, degree, root, iter, success) implicit none ! Subroutine finds one root of a complex polynomial using ! Laguerre's method. In every loop it calculates simplified ! Adams' stopping criterion for the value of the polynomial. ! ! Uses 'root' value as a starting point (!!!!!) ! Remember to initialize 'root' to some initial guess or to ! point (0,0) if you have no prior knowledge. ! ! poly - is an array of polynomial cooefs ! length = degree+1, poly(1) is constant ! 1 2 3 ! poly(1) x^0 + poly(2) x^1 + poly(3) x^2 + ... ! degree - a degree of the polynomial ! root - input: guess for the value of a root ! output: a root of the polynomial ! iter - number of iterations performed (the number of polynomial ! evaluations and stopping criterion evaluation) ! success - is false if routine reaches maximum number of iterations ! ! For a summary of the method go to: ! http://en.wikipedia.org/wiki/Laguerre's_method ! integer MAX_ITERS, FRAC_JUMP_EVERY, FRAC_JUMP_LEN parameter(MAX_ITERS=200) ! Laguerre is used as a failsafe ! constants needed to break cycles in the scheme parameter(FRAC_JUMP_EVERY=10,FRAC_JUMP_LEN=10) double precision FRAC_JUMPS(FRAC_JUMP_LEN) double precision pi parameter(pi = 3.141592653589793d0) double precision faq ! jump length double precision FRAC_ERR parameter(FRAC_ERR = 2.0d-15) ! fractional error for kind=8 (see. Adams 1967 Eqs 9 and 10) ! subroutine parameters integer degree double complex poly(degree+1) ! intent(in) integer iter ! intent(out) double complex root ! intent(inout) logical success ! intent(out) double complex p ! value of polynomial double complex dp ! value of 1st derivative double complex d2p_half ! value of 2nd derivative integer i, k logical good_to_go double complex denom, denom_sqrt, dx, newroot double precision ek, absroot, abs2p double complex fac_netwon, fac_extra, F_half, c_one_nth double precision one_nth, n_1_nth, two_n_div_n_1 double complex zero, c_one parameter(zero=(0d0,0d0),c_one=(1d0,0d0)) double precision stopping_crit2 !--------------------------------------- data FRAC_JUMPS /0.64109297d0, 0.91577881d0, 0.25921289d0, * 0.50487203d0, 0.08177045d0, 0.13653241d0, * 0.306162d0 , 0.37794326d0, 0.04618805d0, * 0.75132137d0/ ! some random numbers iter=0 success=.true. ! next if-endif block is an EXTREME failsafe, not usually needed, and thus turned off in this version. if(.false.)then ! change false-->true if you would like to use caution about having first coefficient == 0 if(degree.lt.0) then write(*,*) 'Error: cmplx_laguerre: degree<0' return endif if(poly(degree+1).eq.zero) then if(degree.eq.0) return call cmplx_laguerre_wrapper * (poly, degree-1, root, iter, success) return endif if(degree.le.1)then if(degree.eq.0) then ! we know from previous check than poly(1) not equal zero success=.false. write(*,*) 'Warning: cmplx_laguerre: '// * 'degree=0 and poly(1)/=0, no roots' return else root=-poly(1)/poly(2) return endif endif endif ! end EXTREME failsafe good_to_go=.false. one_nth=1d0/degree n_1_nth=(degree-1d0)*one_nth two_n_div_n_1=2d0/n_1_nth c_one_nth=dcmplx(one_nth,0d0) do i=1,MAX_ITERS ! prepare stoping criterion ek=abs(poly(degree+1)) absroot=abs(root) ! calculate value of polynomial and its first two derivatives p =poly(degree+1) dp =zero d2p_half=zero do k=degree,1,-1 ! Horner Scheme, see for eg. Numerical Recipes Sec. 5.3 how to evaluate polynomials and derivatives d2p_half=dp + d2p_half*root dp =p + dp*root p =poly(k)+p*root ! b_k ! Adams, Duane A., 1967, "A stopping criterion for polynomial root finding", ! Communications of the ACM, Volume 10 Issue 10, Oct. 1967, p. 655 ! ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/67/55/CS-TR-67-55.pdf ! Eq 8. ek=absroot*ek+abs(p) enddo iter=iter+1 abs2p=dble(conjg(p)*p) if(abs2p.eq.0d0) return stopping_crit2=(FRAC_ERR*ek)**2 if(abs2p.lt.stopping_crit2) then ! (simplified a little Eq. 10 of Adams 1967) ! do additional iteration if we are less than 10x from stopping criterion if(abs2p.lt.0.01d0*stopping_crit2) then return ! return immediately, because we are at very good place else good_to_go=.true. ! do one iteration more endif else good_to_go=.false. ! reset if we are outside the zone of the root endif faq=1d0 denom=zero if(dp.ne.zero)then fac_netwon=p/dp fac_extra=d2p_half/dp F_half=fac_netwon*fac_extra denom_sqrt=sqrt(c_one-two_n_div_n_1*F_half) !G=dp/p ! gradient of ln(p) !G2=G*G !H=G2-2d0*d2p_half/p ! second derivative of ln(p) !denom_sqrt=sqrt( (degree-1)*(degree*H-G2) ) ! NEXT LINE PROBABLY CAN BE COMMENTED OUT if(dble(denom_sqrt).ge.0d0)then ! use realpart(), dble() or dreal() in g77 or ifort77 ! real part of a square root is positive for probably all compilers. You can ! test this on your compiler and if so, you can omit this check denom=c_one_nth+n_1_nth*denom_sqrt else denom=c_one_nth-n_1_nth*denom_sqrt endif endif if(denom.eq.zero)then !test if demoninators are > 0.0 not to divide by zero dx=(absroot+1d0)*exp(dcmplx(0d0, * FRAC_JUMPS(mod(i,FRAC_JUMP_LEN)+1)*2*pi)) ! make some random jump else dx=fac_netwon/denom !dx=degree/denom endif newroot=root-dx if(newroot.eq.root) return ! nothing changes -> return if(good_to_go)then ! this was jump already after stopping criterion was met root=newroot return endif if(mod(i,FRAC_JUMP_EVERY).eq.0) then ! decide whether to do a jump of modified length (to break cycles) faq=FRAC_JUMPS(mod(i/FRAC_JUMP_EVERY-1,FRAC_JUMP_LEN)+1) !write(*,'(3i5,f11.6)') i, i/FRAC_JUMP_EVERY, mod(i/FRAC_JUMP_EVERY-1,FRAC_JUMP_LEN)+1, faq newroot=root-faq*dx ! do jump of some semi-random length (0true if you would like to use caution about having first coefficient == 0 if(degree.lt.0) then write(*,*) 'Error: cmplx_newton_spec: degree<0' return endif if(poly(degree+1).eq.zero) then if(degree.eq.0) return call cmplx_newton_spec_wrapper * (poly, degree-1, root, iter, success) return endif if(degree.le.1)then if(degree.eq.0) then ! we know from previous check than poly(1) not equal zero success=.false. write(*,*) 'Warning: cmplx_newton_spec: degree=0 and'// * ' poly(1)/=0, no roots' return else root=-poly(1)/poly(2) return endif endif endif ! end EXTREME failsafe good_to_go=.false. stopping_crit2 = 0d0 ! value not important, will be initialized anyway on the first loop (because mod(1,10)==1) do i=1,MAX_ITERS faq=1d0 ! prepare stoping criterion ! calculate value of polynomial and its first two derivatives p =poly(degree+1) dp =zero if(mod(i,10).eq.1) then ! calculate stopping criterion every tenth iteration ek=abs(poly(degree+1)) absroot=abs(root) do k=degree,1,-1 ! Horner Scheme, see for eg. Numerical Recipes Sec. 5.3 how to evaluate polynomials and derivatives dp =p + dp*root p =poly(k)+p*root ! b_k ! Adams, Duane A., 1967, "A stopping criterion for polynomial root finding", ! Communications of the ACM, Volume 10 Issue 10, Oct. 1967, p. 655 ! ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/67/55/CS-TR-67-55.pdf ! Eq 8. ek=absroot*ek+abs(p) enddo stopping_crit2=(FRAC_ERR*ek)**2 else ! calculate just the value and derivative do k=degree,1,-1 ! Horner Scheme, see for eg. Numerical Recipes Sec. 5.3 how to evaluate polynomials and derivatives dp =p + dp*root p =poly(k)+p*root ! b_k enddo endif iter=iter+1 abs2p=dble(conjg(p)*p) if(abs2p.eq.0d0) return if(abs2p.lt.stopping_crit2) then ! (simplified a little Eq. 10 of Adams 1967) if(dp.eq.zero) return ! if we have problem with zero, but we are close to the root, just accept ! do additional iteration if we are less than 10x from stopping criterion if(abs2p.lt.0.01d0*stopping_crit2) then return ! return immediately, because we are at very good place else good_to_go=.true. ! do one iteration more endif else good_to_go=.false. ! reset if we are outside the zone of the root endif if(dp.eq.zero)then ! problem with zero dx=(abs(root)+1d0)* * exp(dcmplx(0d0,FRAC_JUMPS(mod(i,FRAC_JUMP_LEN)+1)*2*pi)) ! make some random jump else dx=p/dp ! Newton method, see http://en.wikipedia.org/wiki/Newton's_method endif newroot=root-dx if(newroot.eq.root) return ! nothing changes -> return if(good_to_go)then ! this was jump already after stopping criterion was met root=newroot return endif if(mod(i,FRAC_JUMP_EVERY).eq.0) then ! decide whether to do a jump of modified length (to break cycles) faq=FRAC_JUMPS(mod(i/FRAC_JUMP_EVERY-1,FRAC_JUMP_LEN)+1) newroot=root-faq*dx ! do jump of some semi-random length (0true if you would like to use caution about having first coefficient == 0 if(degree.lt.0) then write(*,*) 'Error: cmplx_laguerre2newton: degree<0' return endif if(poly(degree+1).eq.zero) then if(degree.eq.0) return call cmplx_laguerre2newton_wrapper * (poly, degree-1, root, iter, success, starting_mode) return endif if(degree.le.1)then if(degree.eq.0) then ! we know from previous check than poly(1) not equal zero success=.false. write(*,*) 'Warning: cmplx_laguerre2newton: '// * 'degree=0 and poly(1)/=0, no roots' return else root=-poly(1)/poly(2) return endif endif endif ! end EXTREME failsafe j=1 good_to_go=.false. mode=starting_mode ! mode=2 full laguerre, mode=1 SG, mode=0 newton do ! infinite loop, just to be able to come back from newton, if more than 10 iteration there !------------------------------------------------------------- mode 2 if(mode.ge.2) then ! LAGUERRE'S METHOD one_nth=1d0/degree n_1_nth=(degree-1d0)*one_nth two_n_div_n_1=2d0/n_1_nth c_one_nth=dcmplx(one_nth,0d0) do i=1,MAX_ITERS ! faq=1d0 ! prepare stoping criterion ek=abs(poly(degree+1)) absroot=abs(root) ! calculate value of polynomial and its first two derivatives p =poly(degree+1) dp =zero d2p_half=zero do k=degree,1,-1 ! Horner Scheme, see for eg. Numerical Recipes Sec. 5.3 how to evaluate polynomials and derivatives d2p_half=dp + d2p_half*root dp =p + dp*root p =poly(k)+p*root ! b_k ! Adams, Duane A., 1967, "A stopping criterion for polynomial root finding", ! Communications of the ACM, Volume 10 Issue 10, Oct. 1967, p. 655 ! ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/67/55/CS-TR-67-55.pdf ! Eq 8. ek=absroot*ek+abs(p) enddo abs2p=dble(conjg(p)*p) !abs(p) iter=iter+1 if(abs2p.eq.0d0) return stopping_crit2=(FRAC_ERR*ek)**2 if(abs2p.lt.stopping_crit2) then ! (simplified a little Eq. 10 of Adams 1967) ! do additional iteration if we are less than 10x from stopping criterion if(abs2p.lt.0.01d0*stopping_crit2) then ! ten times better than stopping criterion return ! return immediately, because we are at very good place else good_to_go=.true. ! do one iteration more endif else good_to_go=.false. ! reset if we are outside the zone of the root endif denom=zero if(dp.ne.zero)then fac_netwon=p/dp fac_extra=d2p_half/dp F_half=fac_netwon*fac_extra abs2_F_half=dble(conjg(F_half)*F_half) if(abs2_F_half.le.0.0625d0)then ! F<0.50, F/2<0.25 ! go to SG method if(abs2_F_half.le.0.000625d0)then ! F<0.05, F/2<0.025 mode=0 ! go to Newton's else mode=1 ! go to SG endif endif denom_sqrt=sqrt(c_one-two_n_div_n_1*F_half) ! NEXT LINE PROBABLY CAN BE COMMENTED OUT if(dble(denom_sqrt).ge.0d0)then ! use realpart(), dble() or dreal() for g77 or ifort77 ! real part of a square root is positive for probably all compilers. You can ! test this on your compiler and if so, you can omit this check denom=c_one_nth+n_1_nth*denom_sqrt else denom=c_one_nth-n_1_nth*denom_sqrt endif endif if(denom.eq.zero)then !test if demoninators are > 0.0 not to divide by zero dx=(abs(root)+1d0)*exp( * dcmplx(0d0,FRAC_JUMPS(mod(i,FRAC_JUMP_LEN)+1)*2*pi)) ! make some random jump else dx=fac_netwon/denom endif newroot=root-dx if(newroot.eq.root) return ! nothing changes -> return if(good_to_go)then ! this was jump already after stopping criterion was met root=newroot return endif if(mode.ne.2) then root=newroot j=i+1 ! remember iteration number exit ! go to Newton's or SG endif if(mod(i,FRAC_JUMP_EVERY).eq.0) then ! decide whether to do a jump of modified length (to break cycles) faq=FRAC_JUMPS(mod(i/FRAC_JUMP_EVERY-1,FRAC_JUMP_LEN)+1) newroot=root-faq*dx ! do jump of some semi-random length (0 0.0 not to divide by zero dx=(abs(root)+1d0)*exp( * dcmplx(0d0,FRAC_JUMPS(mod(i,FRAC_JUMP_LEN)+1)*2*pi)) ! make some random jump else fac_netwon=p/dp fac_extra=d2p_half/dp F_half=fac_netwon*fac_extra abs2_F_half=dble(conjg(F_half)*F_half) if(abs2_F_half.le.0.000625d0)then ! F<0.05, F/2<0.025 mode=0 ! set Newton's, go there after jump endif dx=fac_netwon*(c_one+F_half) ! SG endif newroot=root-dx if(newroot.eq.root) return ! nothing changes -> return if(good_to_go)then ! this was jump already after stopping criterion was met root=newroot return endif if(mode.ne.1) then root=newroot j=i+1 ! remember iteration number exit ! go to Newton's endif if(mod(i,FRAC_JUMP_EVERY).eq.0) then ! decide whether to do a jump of modified length (to break cycles) faq=FRAC_JUMPS(mod(i/FRAC_JUMP_EVERY-1,FRAC_JUMP_LEN)+1) newroot=root-faq*dx ! do jump of some semi-random length (0 0.0 not to divide by zero dx=(abs(root)+1d0)*exp( * dcmplx(0d0,FRAC_JUMPS(mod(i,FRAC_JUMP_LEN)+1)*2*pi)) ! make some random jump else dx=p/dp endif newroot=root-dx if(newroot.eq.root) return ! nothing changes -> return if(good_to_go)then root=newroot return endif ! this loop is done only 10 times. So skip this check !if(mod(i,FRAC_JUMP_EVERY).eq.0) then ! decide whether to do a jump of modified length (to break cycles) ! faq=FRAC_JUMPS(mod(i/FRAC_JUMP_EVERY-1,FRAC_JUMP_LEN)+1) ! newroot=root-faq*dx ! do jump of some semi-random length (0',val x=x1 val=poly(3) do i=2,1,-1 val=val*x+poly(i) enddo write(*,'(2f19.15,a3,2f19.15)') x,' ->',val endif end !-------------------------------------------------------------------! ! _10_ SOLVE_CUBIC_EQ ! !-------------------------------------------------------------------! subroutine solve_cubic_eq(x0,x1,x2,poly) ! Cubic equation solver for complex polynomial (degree=3) ! http://en.wikipedia.org/wiki/Cubic_function Lagrange's method implicit none double complex x0, x1, x2 ! intent(out) double complex poly(*) ! intent(in) ! coeffs of the polynomial ! poly - is an array of polynomial cooefs, length = degree+1, poly(1) is constant ! 1 2 3 4 ! poly(1) x^0 + poly(2) x^1 + poly(3) x^2 + poly(4) x^3 double complex zeta, zeta2, zero double precision third parameter(zero=(0d0,0d0)) parameter(zeta =(-0.5d0, 0.8660254037844386d0)) ! sqrt3(1) parameter(zeta2=(-0.5d0,-0.8660254037844386d0)) ! sqrt3(1)**2 parameter(third=0.3333333333333333d0) ! 1/3 double complex s0, s1, s2 double complex E1 ! x0+x1+x2 double complex E2 ! x0x1+x1x2+x2x0 double complex E3 ! x0x1x2 double complex A, B, a_1, E12 double complex delta, A2 double complex val, x integer i a_1=poly(4)**(-1) E1=-poly(3)*a_1 E2=poly(2)*a_1 E3=-poly(1)*a_1 s0=E1 E12=E1*E1 A=2d0*E1*E12-9d0*E1*E2+27d0*E3 ! = s1^3 + s2^3 B=E12-3d0*E2 ! = s1 s2 ! quadratic equation: z^2-Az+B^3=0 where roots are equal to s1^3 and s2^3 A2=A*A delta=sqrt(A2-4d0*(B*B*B)) if( dble(conjg(A)*delta).ge.0d0 )then ! scallar product to decide the sign yielding bigger magnitude ! use dble(), dreal() or realpart() instead of real() not to loose precision in ifort77 or g77 s1=(0.5d0*(A+delta))**third else s1=(0.5d0*(A-delta))**third endif if(s1.eq.zero)then s2=zero else s2=B/s1 endif x0=third*(s0+s1+s2) x1=third*(s0+s1*zeta2+s2*zeta ) x2=third*(s0+s1*zeta +s2*zeta2) if(.false.)then ! print the results x=x0 val=poly(4) do i=3,1,-1 val=val*x+poly(i) enddo write(*,'(2f19.15,a3,2f19.15)') x,' ->',val x=x1 val=poly(4) do i=3,1,-1 val=val*x+poly(i) enddo write(*,'(2f19.15,a3,2f19.15)') x,' ->',val x=x2 val=poly(4) do i=3,1,-1 val=val*x+poly(i) enddo write(*,'(2f19.15,a3,2f19.15)') x,' ->',val endif end !-------------------------------------------------------------------! ! _11_ DIVIDE_POLY_1 ! !-------------------------------------------------------------------! subroutine divide_poly_1(polyout, remainder, p, polyin, degree) ! Subroutine will divide complex polynomial 'polyin' by (x-p) ! results will be returned in polynomial 'polyout' of degree-1 ! The remainder of the division will be returned in 'remainder' ! ! You can provide same array as 'polyin' and 'polyout' - this ! routine will work fine, though it will not set to zero the ! unused, highest coefficient in the output array. You just have ! remember the proper degree of a polynomial. implicit none integer degree ! intent(in) double complex polyout(degree) ! intent(out) double complex remainder ! intent(out) double complex p ! intent(in) double complex polyin(degree+1) ! intent(in) ! coeffs of the polynomial ! poly - is an array of polynomial cooefs, length = degree+1, poly(1) is constant ! 1 2 3 ! poly(1) x^0 + poly(2) x^1 + poly(3) x^2 + ... integer i double complex coef, prev coef=polyin(degree+1) do i=1,degree polyout(i)=polyin(i) enddo do i=degree,1,-1 prev=polyout(i) polyout(i)=coef coef=prev+p*coef enddo remainder=coef return end !-------------------------------------------------------------------! ! _12_ EVAL_POLY ! !-------------------------------------------------------------------! double complex function eval_poly(x, poly, degree, errk) ! Evaluation of the complex polynomial 'poly' of a given degree ! at the point 'x'. This routine calculates also the simplified ! Adams' (1967) stopping criterion. ('errk' should be multiplied ! by 2d-15 for double precision, real*8, arithmetic) implicit none double complex x ! intent(in) integer degree ! intent(in) double precision errk ! intent(out) double complex poly(degree+1) ! intent(in) ! coeffs of the polynomial ! poly - is an array of polynomial cooefs, length = degree+1, poly(1) is constant ! 1 2 3 ! poly(1) x^0 + poly(2) x^1 + poly(3) x^2 + ... integer i double complex val double precision absx ! prepare stoping criterion errk=abs(poly(degree+1)) val=poly(degree+1) absx=abs(x) do i=degree,1,-1 ! Horner Scheme, see for eg. Numerical Recipes Sec. 5.3 how to evaluate polynomials and derivatives val=val*x+poly(i) ! Adams, Duane A., 1967, "A stopping criterion for polynomial root finding", ! Communications of the ACM, Volume 10 Issue 10, Oct. 1967, p. 655 ! ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/67/55/CS-TR-67-55.pdf ! Eq 8. errk=errk*absx+abs(val) enddo eval_poly=val ! if(abs(val)<2d-15*errk) return ! (simplified a little Eq. 10 of Adams 1967) return end !-------------------------------------------------------------------! ! _13_ MULTIPLY_POLY_1 ! !-------------------------------------------------------------------! subroutine multiply_poly_1(polyout, p, polyin, degree) ! Subroutine will multiply polynomial 'polyin' by (x-p) ! results will be returned in polynomial 'polyout' of degree+1 ! ! You can provide same array as 'polyin' and 'polyout' - this ! routine will work fine. implicit none integer degree ! intent(in) ! OLD degree, new will be +1 double complex polyout(degree+2) ! intent(out) double complex p ! intent(in) double complex polyin(degree+1) ! intent(in) ! coeffs of the polynomial ! poly - is an array of polynomial cooefs, length = degree+1, poly(1) is constant ! 1 2 3 ! poly(1) x^0 + poly(2) x^1 + poly(3) x^2 + ... integer i do i=1,degree+1 ! copy polyout(i)=polyin(i) enddo polyout(degree+2)=polyout(degree+1) do i=degree+1,2,-1 polyout(i)=polyout(i-1)-polyout(i)*p enddo polyout(1)=-polyout(1)*p return end !-------------------------------------------------------------------! ! _14_ CREATE_POLY_FROM_ROOTS ! !-------------------------------------------------------------------! subroutine create_poly_from_roots(poly, degree, a, roots) ! Routine will build polynomial from a set of points given in ! the array 'roots'. These points will be zeros of the resulting ! polynomial. ! ! poly - is an array of polynomial cooefs, length = degree+1, poly(1) is constant ! 1 2 3 ! poly(1) x^0 + poly(2) x^1 + poly(3) x^2 + ... ! degree - is and integer denoting size of the 'roots' array ! a - gives the leading coefficient of the resutling polynomial ! roots - input array of points, size=degree ! ! This subroutine works, but it is not optimal - it will work ! up to a polynomial of degree~50, if you would like to have ! more robust routine, up to a degree ~ 2000, you should ! split your factors into a binary tree, and then multiply ! leafs level by level from down up with subroutine like: ! multiply_poly for arbitraty polynomial multiplications ! not multiply_poly_1. ! implicit none integer degree ! intent(in) double complex poly(degree+1) ! intent(out) double complex a ! intent(in) double complex roots(degree) ! intent(in) ! integer i double complex zero parameter(zero=(0d0,0d0)) ! do i=1,degree+1 poly(i)=zero enddo poly(1)=a ! leading coeff of the polynomial do i=1,degree call multiply_poly_1(poly, roots(i), poly, i-1) enddo end !-------------------------------------------------------------------! ! _15_ CMPLX_LAGUERRE_WRAPPER ! !-------------------------------------------------------------------! subroutine cmplx_laguerre_wrapper * (poly, degree, root, iter, success) implicit none ! subroutine parameters integer degree double complex poly(degree+1) ! intent(in) integer iter ! intent(out) double complex root ! intent(inout) logical success ! intent(out) call cmplx_laguerre(poly, degree, root, iter, success) end !-------------------------------------------------------------------! ! _16_ CMPLX_NEWTON_SPEC_WRAPPER ! !-------------------------------------------------------------------! subroutine cmplx_newton_spec_wrapper * (poly, degree, root, iter, success) implicit none integer degree double complex poly(degree+1) ! intent(in) integer iter ! intent(out) double complex root ! intent(inout) logical success ! intent(out) call cmplx_newton_spec(poly, degree, root, iter, success) end !-------------------------------------------------------------------! ! _17_ CMPLX_LAGUERRE2NEWTON_WRAPPER ! !-------------------------------------------------------------------! subroutine cmplx_laguerre2newton_wrapper * (poly, degree, root, iter, success, starting_mode) implicit none ! subroutine parameters integer degree double complex poly(degree+1) ! intent(in) integer iter ! intent(out) double complex root ! intent(inout) logical success ! intent(out) integer starting_mode ! intent(in) call cmplx_laguerre2newton * (poly, degree, root, iter, success, starting_mode) end !-------------------------------------------------------------------!