!*****************************************************
!*  Program to demonstrate Chebyshev economization   *
!* ------------------------------------------------- *
!* Reference: BASIC Scientific Subroutines, Vol. II  *
!* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. *
!*                                                   *
!*              F90 Version by J.-P. Moreau, Paris.  *
!*                      (www.jpmoreau.fr)            *
!* ------------------------------------------------- *
!* SAMPLE RUN:                                       *
!*                                                   *
!* What is the degree of the input polynomial: ? 15  *
!*                                                   *
!* What is the degree of the desired economized      *
!* polynomial: ? 9                                   *
!*                                                   *
!* What is the range of the input polynomial: ? 1.57 *
!*                                                   *
!* Input the coefficients:                           *
!* C( 0) = ? 0                                       *
!* C( 1) = ? 1                                       *
!* C( 2) = ? 0                                       *
!* C( 3) = ? -.166666666                             *
!* C( 4) = ? 0                                       *
!* C( 5) = ? .00833333333                            *
!* C( 6) = ? 0                                       *
!* C( 7) = ? -.0001984127                            *
!* C( 8) = ? 0                                       *
!* C( 9) = ? .000002755732                           *
!* C( 10) = ? 0                                      *
!* C( 11) = ? -.000000025052109                      *
!* C( 12) = ? 0                                      *
!* C( 13) = ? .00000000016059045                     *
!* C( 14) = ? 0                                      *
!* C( 15) = ? -.00000000000076471635                 *
!*                                                   *
!* The Chebyshev series coefficients are:            *
!*                                                   *
!* A( 0) =  0.0000000000                             *
!* A( 1) =  1.1334708982                             *
!* A( 2) =  0.0000000000                             *
!* A( 3) = -0.1378841454                             *
!* A( 4) =  0.0000000000                             *
!* A( 5) =  0.0044798168                             *
!* A( 6) =  0.0000000000                             *
!* A( 7) = -0.0000674667                             *
!* A( 8) =  0.0000000000                             *
!* A( 9) =  0.0000005865                             *
!* A(10) =  0.0000000000                             *
!* A(11) = -0.0000000033                             *
!* A(12) =  0.0000000000                             *
!* A(13) =  0.0000000000                             *
!* A(14) =  0.0000000000                             *
!* A(15) =  0.0000000000                             *
!*                                                   *
!* The economized polynomial coefficients are:       *
!*                                                   *
!* C( 0) =  0.0000000000                             *
!* C( 1) =  0.9999999767                             *
!* C( 2) =  0.0000000000                             *
!* C( 3) = -1.6666647620                             *
!* C( 4) =  0.0000000000                             *
!* C( 5) =  0.0083329009                             *
!* C( 6) =  0.0000000000                             *
!* C( 7) = -0.0001980098                             *
!* C( 8) =  0.0000000000                             *
!* C( 9) =  0.0000025906                             *
!*                                                   *
!*****************************************************
PROGRAM CHEBECON

integer i,m,m1

real*8, allocatable :: B(:,:)
real*8, allocatable :: A(:), C(:)

real*8  x0

  print *,' '
  write(*,'(" What is the degree of the input polynomial: ")',advance='no')
  read *, m
  write(*,'(" What is the degree of the desired economized polynomial: ")',advance='no')
  read *, m1

  allocate (B(m+1,m+1))
  allocate (A(m+1), C(m+1))

  write(*,'(" What is the range of input polynomial: ")',advance='no')
  read *, x0
  print *,' '
  print *, 'Input the coefficients:'
  print *,' '
  do i = 0, m
    write(*,25,advance='no') i
    read *, C(i)
  end do
  print *,' '

  call Cheby_Econ(m,m1,x0,A,B,C)

  print *,'The Chebyshev series coefficients are:'
  print *,' '

  do i = 0, m
    if (dabs(A(i)) < 1.d-50) A(i) = 0.d0
	write(*,50) i, A(i)
  end do
  print *,' '
  
  pause '<Enter> to continue...'

  print *,'The economized polynomial coefficients are:'
  print *,' '

  do i = 0, m1
    if (dabs(C(i)) < 1.d-50) C(i) = 0.d0
    write(*,51) i, C(i)   
  end do
  print *,' '
  stop

25 format(' C(',i2,') = ')
50 format(' A(',i2,') = ',f13.10)
51 format(' C(',i2,') = ',f13.10)
end


!********************************************************
!* Chebyshev series coefficients evaluation subroutine  *
!* The order of the polynomial is n. The coefficients   *
!* are returned in the array B[i,j), i is the degree of *
!* the polynomial, j is the coefficient order.          *
!********************************************************
Subroutine Cheby_Ser(m,n,B)
  parameter(SIZE=10)  
  integer m,n,i,j
  real*8 B(m+1,m+1)
  ! Establish t0 and t1 coefficients
  B(0,0) = 1; B(1,0) = 0; B(1,1) = 1
  ! Return if order is less than two
  if (n < 2) return
  do i = 2, n
    do j = 1, i
      ! Basic recursion relation
      B(i,j) = 2.d0 * B(i-1,j-1) - B(i-2,j)
    end do
    B(i,0) = -B(i-2,0)
  end do
  return
end


!************************************************************
!* Chebyshev economization subroutine. The program takes    *
!* the input polynomial coefficients, C(i), and returns the *
!* Chebyshev series coefficients, A(i). The degree of the   *
!* series passed to the routine is m. The degree of the     *
!* series returned is m1. The maximum range of x is x0 used *
!* for scaling. Note that the input series coefficients are *
!* nulled during the process, and then set equal to the     *
!* economized series coefficients.                          *
!************************************************************
Subroutine Cheby_Econ(m,m1,x0,A,B,C)
  integer i,n,m,m1
  real*8  bb,x0;  
  real*8  B(m+1,m+1)
  real*8  A(m+1),C(m+1)
  ! Start by scaling the input coefficients according to C(i)
  bb = x0
  do i = 1, m
    C(i) = C(i) * bb; bb = bb * x0
  end do
  ! Call Chebyshev series coefficients subroutine
  do n = m, 0, -1
    call Cheby_Ser(m,n,B)
    A(n) = C(n) / B(n,n)
    do l = 0, n
      !Chebyshev series of order l is substracted out of the polynomial
      C(l) = C(l) - A(n) * B(n,l)
    end do
  end do
  ! Perform truncation
  do i = 0, m1
    do j = 0, i
      C(j) = C(j) + A(i) * B(i,j)
    end do
  end do 
  ! Convert back to the integererval X0
  bb = 1.d0 / x0
  do i = 1, m1
    C(i) = C(i) * bb
    bb = bb / x0
  end do
  return
end

! End of file chebecon.f90