C ニュートン法による sin x / x の逆関数
      PROGRAM REP5_3
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION INVFNC
      WRITE(*,*)'X='
      READ(*,*)X
      S=INVFNC(X)
      WRITE(*,'(A23,F22.18)')'Inverse Cardinal =',S
      IF (S.EQ.0.0D0) THEN
         Y=1.0D0  
      ELSE
         Y=DSIN(S)/S
      END IF
      WRITE(*,'(A23,F22.18)')'Final Error      =',X-Y
      END
C
      DOUBLE PRECISION FUNCTION INVFNC(X)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IF (DABS(X).EQ.1.0D0) THEN
          INVFNC=0.0D0
          RETURN
      ELSE IF (DABS(X).GT.0.999999) THEN
          IP=DLOG(1-X)/DLOG(2.0D0)/2
C         WRITE(*,*)IP,DSQRT(6.0D0),2.0D0**IP,(1-X)*2.0D0**(-IP*2),
C     +              (1+3*(1-X)/20)
          INVFNC=DSQRT(6.0D0)*2.0D0**IP*DSQRT((1-X)*2.0D0**(-IP*2))
     +           *(1+3*(1-X)/20)
          RETURN
      END IF
      INVFNC=1.0D0             ! 近似解の初期値
      EPSLON=1.0D-15
      DO 100 I=1,80
          Z=F(X,INVFNC)/G(X,INVFNC)
          INVFNC=INVFNC-Z
          WRITE(*,300)I,'-th Iteration: ',INVFNC,'; Error: ',Z
          IF (DABS(Z).LT.EPSLON) RETURN
  100 CONTINUE
  300 FORMAT(1H ,I2,A15,F18.15,A9,F18.15)
      END
C
      DOUBLE PRECISION FUNCTION F(X,Y)
      DOUBLE PRECISION X,Y
C     F=X*Y-DSIN(Y)
      IF (DABS(Y).GT.0.000001) THEN
          F=X-DSIN(Y)/Y
      ELSE
          F=X-1+Y*Y/6
      END IF
      RETURN
      END
C
      DOUBLE PRECISION FUNCTION G(X,Y)
      DOUBLE PRECISION X,Y
C     G=X-DCOS(Y)
      IF (DABS(Y).GT.0.000001) THEN
          G=(-DCOS(Y)+DSIN(Y)/Y)/Y
      ELSE
C         G=-Y/3+Y*Y*Y/30
          G=0.001
      END IF
      RETURN
      END