Maclaurin series of Black Box Function

Question

Given a Black Box function ff that takes a single floating-point parameter, write a function M to generate the Maclaurin series for that function, according to the following specifications:

• take a parameter n, and generate all terms from x^0 to x^n where 0 <= n <= 10
• coefficients for each term should be given to at least 5 decimal places for non-integral values; underflows can simply be ignored or replaced with 0 and left out if you wish, or you can use scientific notation and more decimal places
• assume ff is a continuous function and everywhere/infinitely differentiable for x > -1, and returns values for input x in the range: -10^9 <= ff(x) <= 10^9
• programming languages with functions as first-class data types can take ff as a parameter; all other languages assume a function called ff has already been defined.
• output should be formatted like # + #x + #x^2 + #x^5 - #x^8 - #x^9 ..., i.e., terms having 0 coefficient should be left out, and negative coefficients should be written - #x^i instead of + -#x^i; if the output terms are all zero, you can output either 0 or null output
• output can be generated either to stdout or as a return value of type string from the function

Examples (passing ff as a first parameter):

> M(x => sin(x), 8)
x - 0.16667x^3 + 0.00833x^5 - 0.00020x^7
> M(x => e^x, 4)
1 + x + 0.5x^2 + 0.16667x^3 + 0.04167x^4
> M(x => sqrt(1+x), 4)
1 + 0.5x - 0.125x^2 + 0.0625x^3 - 0.03906x^4
> M(x => sin(x), 0)
0 (or blank is also accepted)

Answer 1

Ruby, 173 chars

M=->n{r=(-n..n).map{|x|ff(x/2e2)};((0..n).map{|t|u=r[n-t];t<n&&r=r[2..-1].zip(r).map{|x,y|(x-y)*1e2/(t+1)};u.abs<1e-5?'':'%+.5f'%u+(t<2??x[0,t]:'x^%d'%t)}*'').sub(/^\+/,'')}

Unfortunately this approach is similar to Keith's solution and thus also numerically unstable for large n.

def ff(x) Math.sin(x) end
puts M[8]  # 1.00000x-0.16666x^3+0.00833x^5-0.00020x^7

def ff(x) Math.exp(x) end
puts M[4]  # 1.00000+1.00000x+0.50000x^2+0.16667x^3+0.04167x^4

def ff(x) Math.sqrt(1+x) end
puts M[4]  # 1.00000+0.50000x-0.12500x^2+0.06251x^3-0.03907x^4

def ff(x) Math.sin(x) end
puts M[0]  # <empty line>

Answer 2

Python, 188 chars

D=lambda f:lambda x:100*(f(x+.005)-f(x-.005))
def M(f,n):
 s="";d=1
 for i in xrange(n+1):
  c=f(0)/d;f=D(f);d*=i+1
  if c:s+='%+.5f'%c+['','x','x^%d'%i][(i>0)+(i>1)]
 print s[s[:1]=='+':]

This code takes lots of liberties with accuracy and numerical stability, but it generates something close to the right answer for the given examples.

Running it:

import math
M(lambda x:math.sin(x),8)
M(lambda x:math.exp(x),4)
M(lambda x:math.sqrt(x+1),4)
M(lambda x:math.sin(x),0)

gives:

1.00000x-0.16666x^3+0.00833x^5-0.00020x^7
1.00000+1.00000x+0.50000x^2+0.16667x^3+0.04167x^4
1.00000+0.50000x-0.12500x^2+0.06251x^3-0.03907x^4
<a blank line>

Answer 3

Mathematica, 39 chars

InputForm@N@Normal@Series[#1,{x,0,#2}]&

Usage

%[Sin[x], 8]

Output

x - 0.16666666666666666*x^3 + 0.008333333333333333*x^5 - 0.0001984126984126984*x^7

Answer 4

GCC, 549 chars

#define f __float128
#define D long double
#define L(v)for(v=0;v<n;v++)
void G(f a[],int n,double d){f m[n][n],s;int
i,j,k;L(i)L(j)m[i][j]=j?m[i][j-1]*d*(i-n/2):1;L(i){s=m[i][i];L(j)m[i][j]/=s;a[i]/=s;L(j)if(j^i){s=m[j][i];L(k)m[j][k]-=s*m[i][k];a[j]-=s*a[i];}}}
#define F(a)L(i)a[i]=(f)ff(d*(i-n/2));G(a,n,d);
void M(int n){f p[++n],q[n];double d=1;int
i,j=1;F(p)while(j){d/=2;F(q)j=0;L(i){j|=(p[i]-q[i])>1E-7;p[i]=q[i];}}L(i)if(p[i]<-5E-6){printf("- %.5Lfx^%d ",(D)-p[i],i);j=1;}else if(p[i]>5E-6){printf("%s%.5Lfx^%d ",j?"+ ":"",(D)p[i],i);j=1;}}

Example usage (ungolfed):

#include <math.h>
#include <stdio.h>

double ff(double x)
{
    return sin(x);
}

#define f __float128
#define D long double
#define L(v) for(v=0;v<n;v++)
void G(f a[],int n,double d) {
    f m[n][n],s;
    int i,j,k;
    L(i) L(j) m[i][j] = j ? m[i][j-1]*d*(i-n/2) : 1;
    L(i) {
        s=m[i][i]; L(j) m[i][j]/=s; a[i]/=s;
        L(j) if(j^i) {
            s=m[j][i]; L(k) m[j][k]-=s*m[i][k]; a[j]-=s*a[i];
        }
    }
}
#define F(a) L(i) a[i]=(f)ff(d*(i-n/2)); G(a,n,d);
void M(int n) {
    f p[++n],q[n];
    double d=1;
    int i,j=1;

    F(p)
    while(j) {
        d/=2;
        F(q)
        j=0;
        L(i) {
            j|=(p[i]-q[i])>1E-7;
            p[i]=q[i];
        }
    }

    L(i)
        if (p[i]<-5E-6) { printf("- %.5Lfx^%d ",(D)-p[i],i); j=1; }
        else if (p[i]>5E-6) { printf("%s%.5Lfx^%d ",j?"+ ":"",(D)p[i],i); j=1; }
}

void main()
{
    M(8);
}

Compiles with just gcc maclaurin.c -lm

I'm sure that someone who actually knows C can shave at least 10% off, but if nothing else it gets the ball rolling.

With a bit of playing with matrix inverses I've discovered a shorter, but less numerically stable, method which exploits the particular structure of the matrices to special-case the Gaussian elimination:

Alternate, 536 chars

#define f __float128
#define D long double
#define L(v)for(v=0;v<n;v++)
void G(f a[],int n,double d){int i,t;f F[n];L(i)F[i]=i?i*F[i-1]:1;while(--n>0){f c=0;t=2*(n&1)-1;for(i=0;i<=n;i++){t=-t;c+=t*a[i]/F[i]/F[n-i];}L(i)a[i]-=c*pow(i,n);L(i)c/=d;a[n]=c;}}
#define F(a)L(i)a[i]=(f)ff(d*i);G(a,n,d);
void M(int n){f p[++n],q[n];double d=.5;int i,j=1;F(p)while(j){d/=2;F(q)j=0;L(i){j|=(p[i]-q[i])>1E-7;p[i]=q[i];}}L(i)if(p[i]<-5E-6){printf("- %.5Lfx^%d ",(D)-p[i],i);j=1;}else if(p[i]>5E-6){printf("%s%.5Lfx^%d ",j?"+ ":"",(D)p[i],i);j=1;}}

However, this really needs a better error check, so that when it hits the point at which the largest term in the polynomial starts blowing out of control it will quit the loop.