JMK no matter what

Numerical integration and formula parsing. At first, I started implementing Trapezoidal Rule; but it required an insane amount of formula evaluations; so I wrote code to generate parse tree to speed up evaluation. However, in the end, Simpson's Rule saved the day - it requires less than 50 intervals!

But, the amount of effort I threw into this code makes this worth posting anyway.

lang:cpp
#include<iostream>
#include<algorithm>
#include<string>
#include<cmath>
#include<cstdio>
using namespace std;

#define all(x) (x).begin(),(x).end()
#define CLEAR(x,with) memset(x,with,sizeof(x))

// parse tree node
struct expr
{
    /* fields */
    enum { FIXED, VARIABLE, EXPRESSION } type;
    double value;
    char op;
    expr *left, *right;

    /* constructors and destructors */
    expr(double value) : type(FIXED), value(value) {}
    expr(char op, expr* left, expr* right) : type(EXPRESSION), op(op), left(left), right(right)
    {
        // compile time optimization kkkkkkkkkkkkk
        if(left->type == FIXED && right->type == FIXED)
        {
            type = FIXED;
            value = calc(op, left->value, right->value);
            delete left;
            delete right;
        }
    }
    expr() : type(VARIABLE) {}
    ~expr() { if(type == EXPRESSION) { delete left; delete right; } }

    /* member functions */
    double calc(char op, double a, double b)
    {
        if(op == '^') return pow(a, b);
        if(op == '+') return a+b;
        if(op == '-') return a-b;
        if(op == '*') return a*b;
        return a/b;
    }

    double eval(double x)
    {
        if(type == VARIABLE) return x;
        if(type == FIXED) return value;
        return calc(op, left->eval(x), right->eval(x));
    }
};

// compile an expr* from the given string
expr* parse(const string& exp, int left, int right)
{
    int opened = 0;
    // it's an expression?
    for(int i = right; i >= left; --i)
        if(exp[i] == ')') ++opened;
        else if(exp[i] == '(') --opened;
        else if(opened == 0 && (exp[i] == '+' || exp[i] == '-'))
            return new expr(exp[i], parse(exp, left, i-1), parse(exp, i+1, right));
    // term?
    for(int i = right; i >= left; --i)
        if(exp[i] == ')') ++opened;
        else if(exp[i] == '(') --opened;
        else if(opened == 0 && (exp[i] == '*' || exp[i] == '/' || exp[i] == '^'))
            return new expr(exp[i], parse(exp, left, i-1), parse(exp, i+1, right));
    // factor!
    if(exp[left] == 'x') return new expr();
    if(exp[left] == '(') return parse(exp, left+1, right-1);
    return new expr(atof(exp.c_str() + left));
}

// helper for parse() above
expr* parse(const string& exp) { return parse(exp, 0, exp.size()-1); }

double h, b;

// Composite Simpson's Rule
double integrate(double lo, double hi, expr* f)
{
    const int SPLIT = 100;
    double h = (hi - lo) / SPLIT;
    double ret = f->eval(lo) + f->eval(hi);
    for(int i = 1; i < SPLIT; ++i)
    {
        double add = f->eval(lo + h * i);
        ret += (i % 2 ? 4 : 2) * add;
    }
    return ret * h / 3;
}

// returns a pair of <y coordinate of the center, total mass>
pair<double,double> calc(double a, double b, expr* f)
{
    double mass = integrate(a, b, f);

    // f2 = x * f(x)
    expr f2('*', new expr(), f->clone());
    double y = integrate(a, b, &f2) / mass;

    return make_pair(y, mass);
}

pair<double,double> mix(const pair<double,double>& a, const pair<double,double>& b)
{
    double mass = a.second + b.second;
    return make_pair( (a.first * a.second + b.first * b.second) / mass, mass );
}

const double PI = 2.0 * acos(0.0);

int main()
{
    while(cin >> h >> b)
    {
        if(h == 0 && b == 0) break;

        string R, T;
        cin >> R >> T;

        // water area = PI * (R-T)^2. we omit PI
        string wat = "((" + R + ")-(" + T + "))^2";
        // base area = PI*R^2
        string bas = "(" + R + ")^2";
        // wall area = baseArea - waterArea
        string wal = "(" + bas + ")-(" + wat + ")";

        // generate parse trees
        expr* waterArea = parse(wat);
        expr* baseArea = parse(bas);
        expr* wallArea = parse(wal);

        // calc the (unweighted) center of mass, and the mass of 2 parts of the cup
        pair<double,double> base = calc(0, b, baseArea);
        pair<double,double> wall = calc(b, h, wallArea);

        // the cup weighes 2.5 times more than water
        base.second *= 2.5;
        wall.second *= 2.5;

        // the cup as a whole
        pair<double,double> cup = mix(base, wall);

        // binary search on the height of the water
        double lo = b, hi = h;
        for(int iter = 0; iter < 30; ++iter)
        {
            double mid = (lo+hi) / 2;
            pair<double,double> water = calc(b, mid, waterArea);
            pair<double,double> cupAndWater = mix(cup, water);
            // if the top of the water is lower than the center of gravity, we can pour more
            if(mid < cupAndWater.first)
                lo = mid;
            else
                hi = mid;
        }

        // multiply omitted PI
        double amt = integrate(b, hi, waterArea) * PI;
        printf("Pour %.3lf litres / %.3lf cm of water.\n", amt / 1000.0, hi - b);

        delete waterArea;
        delete baseArea;
        delete wallArea;

    }
}