Chebyshev coefficients

Chebyshev coefficients is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Chebyshev coefficients are the basis of polynomial approximations of functions. Write a program to generate Chebyshev coefficients.

Calculate coefficients: cosine function, 10 coefficients, interval 0 1

C

<lang C>// Program to calculate Chebyshev coefficients // Code taken from Numerical Recipes in C 1/e

1. include <math.h>
2. define PI 3.141592653589793

void chebft(float a, float b, float c[], int n, float (*func)(float)) {

     int k,j;
     float fac,bpa,bma,f[300];

     bma = 0.5 * (b-a);
     bpa = 0.5 * (b+a);
     for(k = 0;k<n;k++) {
           float y = cos(PI*(k+0.5)/n);
           f[k] = (*func)(y*bma+bpa);
     }
     fac = 2.0/n;
     for (j = 0;j<n;j++) {
           double sum = 0.0;
           for(k = 0;k<n;k++)
                 sum += f[k] * cos(PI*j*(k+0.5)/n);
           c[j] = fac*sum;
     }

}</lang>

J

From 'J for C Programmers: Calculating Chebyshev Coefficients [[1]] <lang J> chebft =: adverb define

f =. u 0.5 * (+/y) - (-/y) * 2 o. o. (0.5 + i. x) % x

  (2 % x) * +/ f * 2 o. o. (0.5 + i. x) *"0 1 (i. x) % x

) </lang> Calculate coefficients: <lang J>

     10 (2&o.) chebft 0 1

1.64717 _0.232299 _0.0537151 0.00245824 0.000282119 _7.72223e_6 _5.89856e_7 1.15214e_8 6.59629e_10 _1.00227e_11 </lang>

Perl 6

Translation of: C

<lang perl6>sub chebft ( Code $func, Real $a, Real $b, Int $n ) {

   my $bma = 0.5 * ( $b - $a );
   my $bpa = 0.5 * ( $b + $a );

   my @pi_n = ( (^$n).list »+» 0.5 ) »*» ( pi / $n );
   my @f    = ( @pi_n».cos »*» $bma »+» $bpa )».$func;
   my @sums = map { [+] @f »*« ( @pi_n »*» $_ )».cos }, ^$n;

   return @sums »*» ( 2 / $n );

}

say .fmt('%+13.7e') for chebft &cos, 0, 1, 10;</lang>

Output:

+1.6471695e+00
-2.3229937e-01
-5.3715115e-02
+2.4582353e-03
+2.8211906e-04
-7.7222292e-06
-5.8985565e-07
+1.1521427e-08
+6.5962992e-10
-1.0021994e-11

Racket

Translation of: C

<lang racket>#lang typed/racket (: chebft (Real Real Nonnegative-Integer (Real -> Real) -> (Vectorof Real))) (define (chebft a b n func)

 (define b-a/2 (/ (- b a) 2))
 (define b+a/2 (/ (+ b a) 2))
 (define pi/n (/ pi n))
 (define fac (/ 2 n))

 (define f (for/vector : (Vectorof Real)
             ((k : Nonnegative-Integer (in-range n)))
             (define y (cos (* pi/n (+ k 1/2))))
             (func (+ (* y b-a/2) b+a/2))))

 (for/vector : (Vectorof Real)
   ((j : Nonnegative-Integer (in-range n)))
   (define s (for/sum : Real
               ((k : Nonnegative-Integer (in-range n)))
               (* (vector-ref f k)
                  (cos (* pi/n j (+ k 1/2))))))
   (* fac s)))

(module+ test

 (chebft 0 1 10 cos))

Tim Brown 2015</lang>

Output:

'#(1.6471694753903137
   -0.2322993716151719
   -0.05371511462204768
   0.0024582352669816343
   0.0002821190574339161
   -7.722229155637806e-006
   -5.898556451056081e-007
   1.1521427500937876e-008
   6.596299173544651e-010
   -1.0022016549982027e-011)

REXX

Translation of: C

<lang rexx>/*REXX program calculates ten Chebyshev coefficients for the range 0 ──► 1.*/ /*Pafnuty Lvovich Chebysheff: Chebysheff [English transliteration] */ /*─────────────────────────── Chebyshov [ " " */ /*─────────────────────────── Tchebychev [French " */ /*─────────────────────────── Tchebysheff [ " " */ /*─────────────────────────── Tschebyschow [German " */ /*─────────────────────────── Tschebyschev [ " " */ /*─────────────────────────── Tschebyschef [ " " */ /*─────────────────────────── Tschebyscheff [ " " */ numeric digits length(pi()) /*DIGITS default is 9, but use more. */ parse arg a b n . /*obtain optional arguments from the CL.*/ if a== | a==',' then a=0 /*A not specified? Then use the default*/ if b== | b==',' then b=1 /*B " " " " " " */ if n== | n==',' then n=10 /*N " " " " " " */ bma=(b-a)/2 /*calculate one-half of the difference. */ bpa=(b+a)/2 /* " " " " sum. */

           do k=0  for n
           y=cos(pi*(k+.5)/n)
           f.k=y*bma+bpa
           end   /*k*/

fac=2/n

           do j=0  for n;  $=0
           piJ@n=pi*j/n
                          do m=0  for n
                          $=$+f.m*cos(piJ@n*(m+.5))
                          end   /*m*/
           cheby.j=fac*$
           end  /*j*/

exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────subroutines───────────────────────────────────────────────*/ cos: procedure; parse arg x; x=r2r(x); a=abs(x); numeric fuzz 5

               if a=pi then return -1;    if a=pi*.5 | a=pi*2 then return 0
               pi3=pi/3; if a=pi3  then return .5; if a=2*pi3 then return -.5
               z=1;  _=1;    q=x*x
                 do k=2 by 2 until p=z; p=z; _=-_*q/(k*(k-1)); z=z+_; end /*k*/
               return z

pi: pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862; return pi r2r: return arg(1) // (pi()*2) /*normalize radians ──► a unit circle*/</lang>