Mandel + Mandel custom FPU code

// test.cpp
#include <stdio.h>
#include <time.h>

int mandel(double x, double y, int max_iters)
{
    double  a, b, a2;
    int     iter;

    a = 0.0;
    b = 0.0;
    for (iter = 0; iter < max_iters; iter++)
    {
        a2 = (a*a - b*b) + x;
        b = 2.0*a*b + y;
        a = a2;
        if (a*a + b*b >= 4.0)
        {
            max_iters = iter;
            break;
        }
    }
    return max_iters;
}

int mandel_custom(double x, double y, int max_iters)
{
    double  result;
    int     iter, two;


    two = 2;
    _asm {
        // Set computation to 64-bit
        fldz                    // st0=a                                        // push a
        fldz                    // st0=b, st1=a                                 // push b
        fild    two             // st0=2, st1=b, st2=a                          // push 2
        fld     x               // st0=x, st1=2, st2=b, st3=a                   // push x
        fld     y               // st0=y, st1=x, st2=2, st3=b, st4=a            // push y

        // Compute initial a2 value (we know a and b are 0)
        // a2 = (a*a - b*b) + x;
        fld     st(1)           // st0=a2, st1=y, st2=x, st3=2, st4=b, st5=a
    }

    for (iter = 0; iter < max_iters; iter++)
    {
        // Note:  Top-level entry of this loop:  st0=y, st1=x, st2=2, st3=b, st4=a
        _asm {
            // st0=a2

            // b = 2.0*a*b + y;
            fld     st(5)           // st0=a, st1=a2,               st2=y, st3=x, st4=2, st5=b, st6=a
            fmul    st(0),st(4)     // st0=2*a, st1=a2,             st2=y, st3=x, st4=2, st5=b, st6=a
            fmul    st(0),st(5)     // st0=2*a*b, st1=a2,           st2=y, st3=x, st4=2, st5=b, st6=a
            fadd    st(0),st(2)     // st0=(2*a*b)+y, st1=a2,       st2=y, st3=x, st4=2, st5=b, st6=a
            fst     st(5)           // st0=b, st1=a2,               st2=y, st3=x, st4=2, st5=b, st6=a
            // st0=b

            // b*b
            fmul    st(0),st(0)     // st0=b*b, st1=a2,             st2=y, st3=x, st4=2, st5=b, st6=a
            fld     st(0)           // st0=b*b, st1=b*b, st2=a2,    st3=y, st4=x, st5=2, st6=b, st7=a

            // a = a2
            fxch    st(2)           // st0=a2, st1=b*b, st2=b*b,    st3=y, st4=x, st5=2, st6=b, st7=a
            fst     st(7)
            // st0=a

            // if (a*a + b*b >= 4.0)
            // a*a
            fmul    st(0),st(0)     // st0=a*a, st1=b*b, st2=b*b,           st3=y, st4=x, st5=2, st6=b, st7=a
            fsubr   st(2),st(0)     // st0=a*a, st1=b*b, st2=(a*a - b*b),   st3=y, st4=x, st5=2, st6=b, st7=a
            faddp   st(1),st(0)     // st0=(a*a + b*b), st1=(a*a - b*b),    st2=y, st3=x, st4=2, st5=b, st6=a
            fstp    result          // st0=(a*a - b*b),                     st1=y, st2=x, st3=2, st4=b, st5=a

            // Pre-compute a2 for inner loop
            // a2 = (a*a - b*b) + x;
            fadd    st(0),st(2)     // st0=(a*a - b*b) + x,     st1=y, st2=x, st3=2, st4=b, st5=a
            // st0=a2
        }
        if (result >= 4.0)
        {
            max_iters = iter;
            break;
        }
    }

    // Clean the fpu stack
    _asm {
        // st0=a2 st1=y, st2=x, st3=2, st4=b, st5=a
        // Reset FPU
        fstp    st(0)       // stx=a2, st0=y, st1=x, st2=2, st3=b, st4=a        // pop a2
        fstp    st(0)       // stx=a2, stx=y, st0=x, st1=2, st2=b, st3=a        // pop y
        fstp    st(0)       // stx=a2, stx=y, stx=x, st0=2, st1=b, st2=a        // pop x
        fstp    st(0)       // stx=a2, stx=y, stx=x, stx=2, st0=b, st1=a        // pop 2
        fstp    st(0)       // stx=a2, stx=y, stx=x, stx=2, stx=b, st0=a        // pop b
        fstp    st(0)       // stx=a2, stx=y, stx=x, stx=2, stx=b, stx=a        // pop a
    }

    return max_iters;
}

int main(int argc, char* argv[])
{
    int         lnI, max, max_custom;
    double      x, y, step_x, step_y, x_min, x_max, y_min, y_max;
    clock_t     start, end, start_custom, end_custom;


    // Initialize
    x_min   = -3.0;         x_max   = 2.0;
    y_min   = -2.0;         y_max   = 2.0;
    step_x  = +0.0005;      step_y  = +0.0005;

    // Iterate through the rectangle
    start = clock();
    for (lnI = 0; lnI < 5; lnI++)
    {
        max = 0;
        for (x = x_min; x < x_max; x += step_x)
            for (y = y_min; y < y_max; y += step_y)
                max += mandel(x, y, 40);

    }
    end = clock();

    // Iterate through the rectangle
    start_custom = clock();
    for (lnI = 0; lnI < 5; lnI++)
    {
        max_custom = 0;
        for (x = x_min; x < x_max; x += step_x)
            for (y = y_min; y < y_max; y += step_y)
                max_custom += mandel_custom(x, y, 40);

    }
    end_custom = clock();

    printf("Max       : %d, time: %lf\n", max, (double)(end - start) / (lnI * CLOCKS_PER_SEC));
    printf("Max custom: %d, time: %lf\n", max_custom, (double)(end_custom - start_custom) / (lnI * CLOCKS_PER_SEC));

    _asm int 3;     // Injected to cause it to stop even in release mode
    return 0;
}