# Maclaurin series of Black Box Function

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)

• I'm not clear what the function is supposed to do in languages without a REPL. Does it print to stdout, return a string, or whichever is more convenient for me? Also, is there an omitted special case whereby if the parameters are all 0 the output should not exclude the constant term? Oct 26, 2011 at 20:44
• Output method is not specified, as long as it's formatted according to the spec given. I'm a fan of avoiding stdin/stdout disadvantages. As far as when all output is 0, it would be fine to generate blank output or 0 either one. I'll make those both clear, thanks! Oct 26, 2011 at 20:53
• One of your examples is sqrt(1+x), but you allowed us to assume that ff is continuous and everywhere differentiable. Oct 26, 2011 at 22:38
• Updated: everywhere/infinitely differentiable for x > -1. My goal was that caring about discontinuities would not be part of the problem, but in the case of sqrt(1+x), I had just lifted off of the samples on the Wiki page without thinking. I believe you can find the answer as long as ff is infinitely differentiable in at least small range centered around 0. Oct 27, 2011 at 2:00

## 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>


### 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>

• I think it's less that the approach is numerically unstable and more that the problem is ill-conditioned. As a rough estimate, based on extrapolating the condition numbers of the inverse Vandermonde matrices for the simple polynomial fit approach, 8 significant decimal digits is as much as it's reasonable to hope for when n=8 even assuming infinite precision for the intermediate results. Oct 28, 2011 at 22:01

## 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

• While correct, this answer involves knowledge of the function itself to generate the series. The idea behind the question is assuming the only thing you can know are input/output pairs of the function, and not what the original function actually is. Apr 24, 2012 at 16:33
• @mellamokb No, your interpretation is not correct. #1 is the function (black box, as requested), and is a parameter. See usage example above and compare with the first example in the question. Apr 24, 2012 at 17:15
• I understand that. I mean the Mathematica engine could make use of symbolic math to generate the series. I can't tell from your example whether it is making limit approximations, or just calculating the derivatives directly from the input. In any case, Mathematica is probably too high-level for this problem, since it's involving direct calculation of limits and calculating the series completely manually. Apr 24, 2012 at 18:02
• @mellamokb calculating the series completely manually is not stated in the Q. Anyway, the best suited the language for the problem, the least fun (but the better golf). Apr 24, 2012 at 18:09
• That's the whole point / challenge of the question. It's implicit simply because in most programming languages, Sin[x] is not understood in any semantic sense by the language other than knowing how to perform the calculation on a given input. Since Mathematica can interpret Sin[x] at a higher-level than just numeric calculation, that implicit expectation breaks down. I would feel the same way, for example, to a solution given in TI-89 Basic as Taylor(ff, x, 1, 8). So ff is passed in instead of directly put in as an argument; it's still no challenge. Anyway +1 for a correct answer :) Apr 24, 2012 at 18:18

## 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.