'*********************************************************
'*  Calculate a limited development of a real function   *
'*  f(x)og(x) at x=0 up to order 25, knowing the limited *
'*  developments of f(x) and g(x).                       *
'* ----------------------------------------------------- *
'* SAMPLE RUN:                                           *
'*                                                       *
'* Limited development of a real function f(x)og(x):     *
'*                                                       *
'* Function to develop: exp(sin(x))                      *
'*                                                       *
'* Limited development of f(x)=exp(x) is:                *
'* 1 + x + x^2/2! + x^3/3! + ... + x^n/n! + ...          *
'*                                                       *
'* Limited development of g(x)=sin(x) is:                *
'* x - x^3/3! +x^5/5! - ... + (-1)^n x^2n+1/(2n+1)! + ...*
'*                                                       *
'* Order of development (max=25): 10                     *
'*                                                       *
'* Input coefficients of limited dev. of f(x):           *
'*  a0 = 1                                               *
'*  a1 = 1                                               *
'*  a2 = 0.5                                             *
'*  a3 = 0.16666667                                      *
'*  a4 = 0.04166667                                      *
'*  a5 = 0.00833333                                      *
'*  a6 = 0.00138889                                      *
'*  a7 = 0.00019841                                      *
'*  a8 = 0.00002480                                      *
'*  a9 = 0.00000276                                      *
'*  a10 = 0.000000276                                    *
'*                                                       *
'* Input coefficients of limited dev. of g(x):           *
'*  b0 = 0                                               *
'*  b1 = 1                                               *
'*  b2 = 0                                               *
'*  b3 = -0.16666667                                     *
'*  b4 = 0                                               *
'*  b5 = 0.00833333                                      *
'*  b6 = 0                                               *
'*  b7 = -0.00019841                                     *
'*  b8 = 0                                               *
'*  b9 = 0.00000276                                      *
'*  b10 = 0                                              *                        
'*                                                       *
'* The coefficients of limited dev. of f(x)og(x) are:    *
'*  c0 =  1.0000000                                      *
'*  c1 =  1.0000000                                      *
'*  c2 =  0.5000000                                      *
'*  c3 =  0.0000000                                      *
'*  c4 = -0.1250000                                      *
'*  c5 = -0.0666667                                      *
'*  c6 = -0.0041667                                      *
'*  c7 =  0.0111111                                      *
'*  c8 =  0.0053819                                      *
'*  c9 =  0.0001764                                      *
'*  c10 = -0.0008132                                     *
'*                                                       *
'* ----------------------------------------------------- *
'* Ref.: "Mathematiques en Turbo Pascal By Alain         *
'*        Reverchon and Marc Ducamp, Armand Colin        *
'*        Editeur, Paris, 1991" [BIBLI 03].              *
'*                                                       *
'*                          Basic Version By J-P Moreau. *
'*                               (www.jpmoreau.fr)       *
'*********************************************************
'PROGRAM LTDDEV3
DEFDBL A-H, O-Z
DEFINT I-N

NSIZE = 25

DIM t1(NSIZE + 1), t2(NSIZE + 1), R(NSIZE + 1)
DIM binome(NSIZE + 1, NSIZE + 1)
DIM id(NSIZE + 1)

  CLS
  PRINT
  PRINT " Limited development of a real function f(x)og(x):"
  PRINT
  PRINT " Function to develop: exp(sin(x))"
  PRINT
  PRINT " Limited development of f(x)=exp(x) is:"
  PRINT "  1 + x + x^2/2! + x^3/3! + ... + x^n/n! + ..."
  PRINT
  PRINT " Limited development of g(x)=sin(x) is:"
  PRINT "  x - x^3/3! +x^5/5! - ... + (-1)^n x^2n+1/(2n+1)! + ..."
  PRINT
  INPUT " Order of development (max=25): ", n
  PRINT
  PRINT " Input coefficients of limited dev. of f(x):"
  FOR i = 0 TO n
    PRINT "  a"; i; " = "; : INPUT "", t1(i + 1)
  NEXT i
  PRINT
  PRINT " Input coefficients of limited dev. of g(x):"
  FOR i = 0 TO n
    PRINT "  b"; i; " = "; : INPUT "", t2(i + 1)
  NEXT i
  PRINT

  GOSUB 1000  'call ar_dlcomp(n,T1,T2,R)

  PRINT " The coefficients of limited dev. of f(x)og(x) are:"
  FOR i = 0 TO n
    PRINT "  c"; i; " = "; : PRINT USING "##.#######"; R(i + 1)
  NEXT i
  PRINT
  PRINT

END

'Function ProductBin
500 u = t2(id(2) + 1)
  ip = i: iq = 1
  FOR k = 2 TO i
    u = u * t2(id(k + 1) + 1)
    IF id(k + 1) = id(k) THEN
      iq = iq + 1
    ELSE
      u = u * binome(ip + 1, iq + 1)
      ip = ip - iq
      iq = 1
    END IF
  NEXT k
  ProductBin = u
RETURN


'Subroutine ar_dlcomp(n,T1,T2,R)
1000 IF n > NSIZE THEN RETURN
  IF t2(1) <> 0# THEN RETURN
  R(1) = t1(1)
  'calculate coefficients of binome
  binome(2, 1) = 1#
  binome(2, 2) = 1#
  FOR i = 2 TO n
    binome(i + 1, 1) = 1#
    binome(i + 1, i + 1) = 1#
    FOR j = 1 TO i - 1
      binome(i + 1, j + 1) = binome(i, j) + binome(i, j + 1)
    NEXT j
  NEXT i
  'calculate coefficients of gof
  FOR m = 1 TO n
    x = t1(2) * t2(m + 1)
    FOR i = 2 TO m
      'calculate s=gamma(i,m)
      s = 0#
      FOR j = 1 TO i - 1
	id(j + 1) = 1
      NEXT j
      ipos = 1
10    FOR j = i - ipos + 1 TO i - 1
	id(j + 1) = id(i - ipos + 1)
      NEXT j
      id(i + 1) = m
      FOR j = 1 TO i - 1
	id(i + 1) = id(i + 1) - id(j + 1)
      NEXT j
20    IF id(i + 1) < id(i) THEN GOTO 30
	ipos = 1
	GOSUB 500 'calculate ProductBin
	s = s + ProductBin
	id(i + 1) = id(i + 1) - 1
	id(i) = id(i) + 1
      GOTO 20
30    ipos = ipos + 1
      IF ipos < i THEN
	id(i - ipos + 1) = id(i - ipos + 1) + 1
	GOTO 10
      END IF
      x = x + s * t1(i + 1)
    NEXT i
    R(m + 1) = x
  NEXT m
RETURN

' end of file ltddev3.bas