c   8.Newton$BK!$GO"N)J}Dx<0$r2r$/(B
c     $B?9@n(B $BLwBg(B , 22980138
c
      program ex8
c
      implicit none
      integer :: n,i,j
      real :: a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,
     -        b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,
     -        x0,y0,xk1,yk1,epsilon
      epsilon = 1.0e-4    ! $B<}B+@:EY(B

      write(6,*)'Newton$BK!$GO"N)J}Dx<0$r2r$/(B'
      write(6,300) epsilon
      do i=1,2
       if (i==1) then
        write(6,*)''
        write(6,*)'(1)'
        write(6,*)'$BO"N)J}Dx<0(B x**2+y**2-1=0'
        write(6,*)'           x**3-y=0     $B$N2r(B'

        a0= 0.0e0           !$BBh(B1$B<0$N78?t(B
        a1= 0.0e0
        a2= 0.0e0
        a3= 0.0e0
        a4= 1.0e0
        a5= 0.0e0
        a6= 1.0e0
        a7= 0.0e0
        a8= 0.0e0
        a9=-1.0e0
       
        b0= 1.0e0           !$BBh(B2$B<0$N78?t(B
        b1= 0.0e0
        b2= 0.0e0
        b3= 0.0e0
        b4= 0.0e0
        b5= 0.0e0
        b6= 0.0e0
        b7= 0.0e0
        b8=-1.0e0
        b9= 0.0e0
       
       else if (i==2) then
        write(6,*)''
        write(6,*)'(2)'
        write(6,*)'$BO"N)J}Dx<0(B x**2+y**2-1=0'
        write(6,*)'           x+y**2=0     $B$N2r(B'

        a0= 0.0e0           !$BBh(B1$B<0$N78?t(B
        a1= 0.0e0
        a2= 0.0e0
        a3= 0.0e0
        a4= 1.0e0
        a5= 0.0e0
        a6= 1.0e0
        a7= 0.0e0
        a8= 0.0e0
        a9=-1.0e0
       
        b0= 0.0e0           !$BBh(B2$B<0$N78?t(B
        b1= 0.0e0
        b2= 0.0e0
        b3= 0.0e0
        b4= 0.0e0
        b5= 0.0e0
        b6= 1.0e0
        b7= 1.0e0
        b8= 0.0e0
        b9= 0.0e0
       end if     
       do j=1,2
        if  (j==1)  then
         x0=-0.1e0           ! $B=i4|CM(B
         y0=-1.0e0
        else if (j==2) then
         x0= 0.1e0
         y0= 1.0e0
        end if
        write(6,400) x0,y0

        call newton(a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,
     -              b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,
     -              x0,y0,xk1,yk1,n,epsilon)
 
        write(6,200) n      
        write(6,100)'x=',xk1
        write(6,100)'y=',yk1
        write(6,*)''
       enddo
      enddo
      stop
  100 format(a3,e18.9)
  200 format(1x,'$BH?I|2s?t(B',i5) 
  300 format(1x,'$B<}B+@:EY(B',1x,e18.9)
  400 format(1x,'$B=i4|CM(B(x,y)=(',f4.1,',',f4.1,')') 
      end

c-----Newton$BK!$G(B2$B$D$N(Bx,y$B$N(B3$B<!J}Dx<0$N(B
c-----$B2r$HH?I|2s?t$r5a$a$k%5%V%k!<%A%s(B
      subroutine newton(a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,
     -                  b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,
     -                  x0,y0,xk1,yk1,n,epsilon)
       implicit none
       integer :: n
       real :: f,g,dfx,dgx,dfy,dgy,fs,gs,dfxs,dgxs,dfys,dgys,
     -         a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,
     -         b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,
     -         x0,y0,xk,yk,xk1,yk1,epsilon
      
       n=1
       xk=x0
       yk=y0
       fs  =  f(a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,xk,yk)
       gs  =  g(b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,xk,yk) 
       dfxs=dfx(a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,xk,yk)
       dgxs=dgx(b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,xk,yk)
       dfys=dfy(a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,xk,yk)
       dgys=dgy(b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,xk,yk)

       xk1=xk-( fs*dgys-gs*dfys)/(dfxs*dgys-dfys*dgxs)
       yk1=yk-(-fs*dgxs+gs*dfxs)/(dfxs*dgys-dfys*dgxs)
  
        do while (sqrt((xk1-xk)**2+(yk1-yk)**2) >= epsilon)
          n=n+1
          xk=xk1
          yk=yk1

          fs  =  f(a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,xk,yk)
          gs  =  g(b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,xk,yk) 
          dfxs=dfx(a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,xk,yk)
          dgxs=dgx(b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,xk,yk)
          dfys=dfy(a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,xk,yk)
          dgys=dgy(b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,xk,yk)

          xk1=xk-( fs*dgys-gs*dfys)/(dfxs*dgys-dfys*dgxs)
          yk1=yk-(-fs*dgxs+gs*dfxs)/(dfxs*dgys-dfys*dgxs)
        enddo

      return
      end

c-----$B4X?t(B f(x,y)
      function f(a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,xk,yk)
      implicit none
      real :: f,xk,yk,a0,a1,a2,a3,a4,a5,a6,a7,a8,a9
      
      f= a0*xk**3 +a1*xk**2*yk +a2*xk*yk**2 +a3*yk**3 +a4*xk**2
     -  +a5*xk*yk +a6*yk**2    +a7*xk       +a8*yk    +a9

      return
      end

c-----$B4X?t(B g(x,y)
      function g(b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,xk,yk)
      implicit none
      real :: g,xk,yk,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9
      
      g= b0*xk**3 +b1*xk**2*yk +b2*xk*yk**2 +b3*yk**3 +b4*xk**2
     -  +b5*xk*yk +b6*yk**2    +b7*xk       +b8*yk    +b9

      return
      end

c-----$B4X?t(B df(x,y)/dx
      function dfx(a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,xk,yk)
      implicit none
      real :: dfx,xk,yk,a0,a1,a2,a3,a4,a5,a6,a7,a8,a9

      dfx= 3.0e0*a0*xk**2 +2.0e0*a1*xk*yk +a2*yk**2 +a3*0.0e0
     -    +2.0e0*a4*xk +a5*yk +a6*0.0e0 +a7 +a8*0.0e0 +a9*0.0e0

      return
      end

c-----$B4X?t(B dg(x,y)/dx
      function dgx(b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,xk,yk)
      implicit none
      real :: dgx,xk,yk,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9

      dgx= 3.0e0*b0*xk**2 +2.0e0*b1*xk*yk +b2*yk**2 +b3*0.0e0
     -    +2.0e0*b4*xk +b5*yk +b6*0.0e0 +b7 +b8*0.0e0 +b9*0.0e0

      return
      end

c-----$B4X?t(B df(x,y)/dy
      function dfy(a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,xk,yk)
      implicit none
      real :: dfy,xk,yk,a0,a1,a2,a3,a4,a5,a6,a7,a8,a9
      
      dfy= a0*0.0e0 +a1*xk**2 +2.0e0*a2*xk*yk +3.0e0*a3*yk**2 
     -    +a4*0.0e0 +a5*xk +2.0e0*a6*yk +a7*0.0e0 +a8 +a9*0.0e0

      return
      end

c-----$B4X?t(B dg(x,y)/dy
      function dgy(b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,xk,yk)
      implicit none
      real :: dgy,xk,yk,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9
      
      dgy= b0*0.0e0 +b1*xk**2 +2.0e0*b2*xk*yk +3.0e0*b3*yk**2 
     -    +b4*0.0e0 +b5*xk +2.0e0*b6*yk +b7*0.0e0 +b8 +b9*0.0e0

      return
      end
