5
\$\begingroup\$

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

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>
Improve this answer
\$\endgroup\$
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ – Peter Taylor Oct 28 '11 at 22:01
5
\$\begingroup\$

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>
Improve this answer
\$\endgroup\$
2
\$\begingroup\$

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
Improve this answer
\$\endgroup\$
  • \$\begingroup\$ 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. \$\endgroup\$ – mellamokb Apr 24 '12 at 16:33
  • \$\begingroup\$ @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. \$\endgroup\$ – Dr. belisarius Apr 24 '12 at 17:15
  • \$\begingroup\$ 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. \$\endgroup\$ – mellamokb Apr 24 '12 at 18:02
  • \$\begingroup\$ @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). \$\endgroup\$ – Dr. belisarius Apr 24 '12 at 18:09
  • \$\begingroup\$ 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 :) \$\endgroup\$ – mellamokb Apr 24 '12 at 18:18
2
\$\begingroup\$

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.

Improve this answer
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.