# 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? – Peter Taylor Oct 26 '11 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! – mellamokb Oct 26 '11 at 20:53
• One of your examples is sqrt(1+x), but you allowed us to assume that ff is continuous and everywhere differentiable. – Peter Taylor Oct 26 '11 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. – mellamokb Oct 27 '11 at 2:00

### 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  # 1.00000x-0.16666x^3+0.00833x^5-0.00020x^7

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

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

def ff(x) Math.sin(x) end
puts M  # <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. – Peter Taylor Oct 28 '11 at 22:01

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


## 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. – mellamokb Apr 24 '12 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. – Dr. belisarius Apr 24 '12 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. – mellamokb Apr 24 '12 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). – Dr. belisarius Apr 24 '12 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 :) – mellamokb Apr 24 '12 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.