'*********************************************************
'*  Calculate a limited development of a real function   *
'*  f(x) at point x0 with step h (to order 5).           *
'* ----------------------------------------------------- *
'* SAMPLE RUN:                                           *
'*                                                       *
'* Limited development of f(x) at x0:                    *
'*                                                       *
'* Function to develop: 1/(1-x)                          *
'*                                                       *
'* At point x0= 0                                        *
'* Value of step h= 0.1                                  *
'*                                                       *
'* a1 =  1.0000416250                                    *
'* a2 =  1.0000416250                                    *
'* a3 =  1.0011207928                                    *
'* a4 =  1.0011207928                                    *
'* a5 =  1.0136468470                                    *
'*                                                       *
'* (The correct answer is a1=a2=a3=a4=a5=1)              *     
'*                                                       *
'* Function to develop: Ln(x)                            *
'*                                                       *
'* At point x0= 1                                        *
'* Value of step h= 0.1                                  *
'*                                                       *
'* a1 =  1.0000057552                                    *
'* a2 = -0.5000050520                                    *
'* a3 =  0.3334508142                                    *
'* a4 = -0.2501062334                                    *
'* a5 =  0.2011301593                                    *
'*                                                       *
'* (the correct answer is a1=1, a2=-0.5, a3=1/3,         *
'*  a4=-1/4 and a5=1/5)                                  *
'* ----------------------------------------------------- *
'* Ref.: "Mathematiques en Turbo Pascal By Alain         *
'*        Reverchon and Marc Ducamp, Armand Colin        *
'*        Editeur, Paris, 1991" [BIBLI 03].              *
'* ----------------------------------------------------- *
'* Subroutines used:                                     *
'*                                                       *
'*   500: Sample function f(x)                           *
'*  1000: Nth derivative of a discrete function          *
'*        (N=1 to 5).                                    *
'*                                                       *  
'*                          Basic Version By J-P Moreau. *
'*                               (www.jpmoreau.fr)       *
'*********************************************************
'Note: as you can see this direct method is not suitable to
'      calculate a limited development in an accurate way.
'      The next programs (ltddev1 to ltddev3) will allow us
'      to calculate limited developments at x=0 up to a high
'      order with a good accuracy for functions f*g, f/g or
'      gof, knowing the limited developments of f(x) and g(x)
'      at x=0.
'-----------------------------------------------------------}
'Program Ltddev
DEFDBL A-H, O-Z
DEFINT I-N

DIM itable(5, 7)

DATA 3,45,-9,1,0,0,60,3,270,-27,2,0,0,180,4,-488,338,-72,7,0,240
DATA 4,-1952,676,-96,7,0,240,5,1938,-1872,783,-152,13,288
 
  'initialize table of coefficients
  FOR i = 1 TO 5
    FOR j = 1 TO 7
      READ itable(i, j)
    NEXT j
  NEXT i

  CLS
  PRINT
  PRINT " Limited development of a real function f(x):"
  PRINT
  PRINT " Function to develop: 1/(1-x)"
  PRINT
  INPUT " At point x0= ", x0
  INPUT " Value of step h= ", h
  PRINT
 
  s = 1#

  FOR n = 1 TO 5
    s = s * n
    GOSUB 1000  'r = arndederiv(x0,n,h)
    PRINT "  a"; n; " = ";
    PRINT USING "##.##########"; R / s
  NEXT n
  PRINT

END 'of main program

'sample function 1 f(xx)=1/(1-xx)
500 IF ABS(xx - 1) > 1E-12 THEN
    f = 1# / (1# - xx)
  ELSE
    f = 1E+12
  END IF
RETURN

'sample function 2 f(xx)=Log(xx)
'500 if xx >1e-16 then
'    f = LOG(x)
'  else
'    f = -1e16
'  END IF
'return                      }

'--------------------------------------------------
'returns the value of nth derivative of function
'f(x) at point x0 with increment step h (n=1 to 5).
'Note that for polynomial functions h must not be
'too small. Here h=1 gives best results than h=0.01
'--------------------------------------------------
1000 'real*8 Function arnderiv(x0,n,h)

  IF (n > 5) THEN
    R = 0#
    RETURN
  END IF
 
  xx = x0: GOSUB 500: fa = f
  t = 0#
  FOR i = 1 TO itable(n, 1)
    xx = x0 - i * h: GOSUB 500: fb = f
    xx = x0 + i * h: GOSUB 500: fc = f
    IF (n MOD 2) <> 0 THEN
      t = t + itable(n, i + 1) * (fc - fb)
    ELSE
      t = t + itable(n, i + 1) * ((fc - fa) + (fb - fa))
    END IF
  NEXT i
  t = t / itable(n, 7)
  FOR i = 1 TO n
    t = t / h
  NEXT i
  R = t
RETURN

'end of file ltddev.bas