#include <exec/types.h>
#include <math.h>
#include "ppc.h"

#ifdef _M68000
#include <proto/exec.h>
#else
#include <PowerUP/GCCLib/powerup_protos.h>
#endif

ULONG   errno = 0;

//#define NO_SIGNALS

/* Signal stuff */
#define SIGBREAKB_CTRL_C  12                            /* from <dos/dos.h> */
#define SIGBREAKF_CTRL_C  (1 << SIGBREAKB_CTRL_C)       /* from <dos/dos.h> */
ULONG   SignalBreak;

#ifdef _M68000                                          /* M680x0 */
#define AllocBreak      SignalBreak = AllocSignal(SIGBREAKB_CTRL_C);
#define FreeBreak       FreeSignal(SignalBreak);
#define CheckBreak      if (SetSignal(0, SIGBREAKF_CTRL_C) & SIGBREAKF_CTRL_C) break;

#else                                                   /* PowerPC */
#define AllocBreak      SignalBreak = PPCAllocSignal(SIGBREAKB_CTRL_C);
#define FreeBreak       PPCFreeSignal(SignalBreak);
#define CheckBreak      if (PPCSetSignal(0, SIGBREAKF_CTRL_C) & SIGBREAKF_CTRL_C) break;
#endif

/* Color defines */
#define MIN_COL 0
#define MAX_COL 255
#define STEP    3

/* function prototypes */
VOID mandelbrot         (struct ppcIO *ppcIO);
VOID mandcubed          (struct ppcIO *ppcIO);
VOID mandfourth         (struct ppcIO *ppcIO);
VOID mandexpo           (struct ppcIO *ppcIO);
VOID julia              (struct ppcIO *ppcIO);
VOID juliacubed         (struct ppcIO *ppcIO);
VOID juliafourth        (struct ppcIO *ppcIO);
VOID juliaexpo          (struct ppcIO *ppcIO);
VOID dragon             (struct ppcIO *ppcIO);
VOID salamander         (struct ppcIO *ppcIO);
VOID meteors            (struct ppcIO *ppcIO);
VOID pearls             (struct ppcIO *ppcIO);
VOID formula            (struct ppcIO *ppcIO);

/* main() routine for PPC */
ULONG calcmain(struct ppcIO *ppcIO)
{
        ULONG   i;

#ifndef NO_SIGNALS
        /* Allocate Break Signal */
        AllocBreak;
#endif

        /* Clear image buffer */
        for (i = 0; i < (ppcIO->width * ppcIO->height); i++) {
                ppcIO->buffer[i] = 0;
        }

        switch (ppcIO->type) {
                case TYPE_MANDELBROT:
                        switch (ppcIO->exponent) {
                                case 0:
                                case 1:
                                        break;
                                case 2:
                                        mandelbrot(ppcIO);
                                        break;
                                case 3:
                                        mandcubed(ppcIO);
                                        break;
                                case 4:
                                        mandfourth(ppcIO);
                                        break;
                                default:
                                        mandexpo(ppcIO);
                                        break;
                        }
                        break;
                case TYPE_JULIA:
                        switch (ppcIO->exponent) {
                                case 0:
                                case 1:
                                        break;
                                case 2:
                                        julia(ppcIO);
                                        break;
                                case 3:
                                        juliacubed(ppcIO);
                                        break;
                                case 4:
                                        juliafourth(ppcIO);
                                        break;
                                default:
                                        juliaexpo(ppcIO);
                                        break;
                        }
                        break;
                case TYPE_DRAGON:
                        dragon(ppcIO);
                        break;
                case TYPE_SALAMANDER:
                        salamander(ppcIO);
                        break;
                case TYPE_METEORS:
                        meteors(ppcIO);
                        break;
                case TYPE_PEARLS:
                        pearls(ppcIO);
                        break;
                case TYPE_FORMULA:
                        formula(ppcIO);
                        break;
        }

#ifndef NO_SIGNALS
        /* Free Break Signal */
        FreeBreak;
#endif

        return 0;
}

#define plot_pixel(x, y, iter, ppcIO)   \
{                                       \
        ULONG   value, phase;           \
                                        \
        if (iter >= ppcIO->maxiter)     \
                ppcIO->buffer[y * ppcIO->width + x] = 0;        \
        else {                                                  \
                phase = iter * 6 / ppcIO->repeat;               \
                value = (ULONG) (((double) iter * 6.0 / (double) ppcIO->repeat - (double) phase) * 255.0);      \
                                                                        \
                if ((ppcIO->evencol > 0) && ((iter & 1) == 0)) {        \
                        phase += ppcIO->evencol;                        \
                        if (phase >= 6)                                 \
                                phase -= 6;                             \
                }                                                       \
                if ((ppcIO->oddcol > 0) && ((iter & 1) == 1)) {         \
                        phase += ppcIO->oddcol;                         \
                        if (phase >= 6)                                 \
                                phase -= 6;                             \
                }                                                       \
                                                                        \
                switch (phase) {                                        \
                        case 0:                                         \
                                ppcIO->buffer[y * ppcIO->width + x] = (MAX_COL << 16) | (value << 8);                   \
                                break;                                  \
                        case 1:                                         \
                                ppcIO->buffer[y * ppcIO->width + x] = ((MAX_COL - value) << 16) | (MAX_COL << 8);       \
                                break;                                  \
                        case 2:                                         \
                                ppcIO->buffer[y * ppcIO->width + x] = (MAX_COL << 8) | value;                           \
                                break;                                  \
                        case 3:                                         \
                                ppcIO->buffer[y * ppcIO->width + x] = ((MAX_COL - value) << 8) | MAX_COL;               \
                                break;                                  \
                        case 4:                                         \
                                ppcIO->buffer[y * ppcIO->width + x] = (value << 16) | MAX_COL;                          \
                                break;                                  \
                        case 5:                                         \
                                ppcIO->buffer[y * ppcIO->width + x] = (MAX_COL << 16) | (MAX_COL - value);              \
                                break;                                  \
                }                                                       \
        }                                                               \
}

double  xa, ya, da, v;

#define cadd(x1, y1, x, y)                      \
{                                               \
        x += x1;                                \
        y += y1;                                \
}

#define csub(x1, y1, x, y)                      \
{                                               \
        x -= x1;                                \
        y -= y1;                                \
}

#define cmul(x1, y1, x, y)                      \
{                                               \
        xa = x1 * x - y1 * y;                   \
        ya = x1 * y + x * y1;                   \
        x = xa;                                 \
        y = ya;                                 \
}

#define cdiv(x1, y1, x, y)                      \
{                                               \
        da = x1 * x1 + y1 * y1;                 \
        xa = (x * x1 + y * y1) / da;            \
        ya = (y * x1 - x * y1) / da;            \
        x = xa;                                 \
        y = ya;                                 \
}

#define cimul(x1, y1, x, y)                     \
{                                               \
        x = -y1;                                \
        y = x1;                                 \
}

#define csqr(x1, y1, x, y)                      \
{                                               \
        x = x1 * x1 - y1 * y1;                  \
        y = 2.0 * x1 * y1;                      \
}

#define ccube(x1, y1, x, y)                     \
{                                               \
        x = x1 * x1 * x1 - 3.0 * x1 * y1 * y1;  \
        y = 3.0 * x1 * x1 * y1 - y1 * y1 * y1;  \
}

#define cpow(p, x, y)                           \
{                                               \
        da = pow(x * x + y * y, p / 2.0);       \
        if (x == 0.0)                           \
                if (y >= 0.0)                   \
                        v = PID2;               \
                else                            \
                        v = -PID2;              \
        else                                    \
                v = atan2(y, x);                \
        x = da * cos(v * p);                    \
        y = da * sin(v * p);                    \
}



VOID mandelbrot(struct ppcIO *ppcIO)
{
        ULONG  i, j, iter;
        double x1, y1, xs, ys, x, y, xstep, ystep, count;

        x     = ppcIO->xmin;
        y     = ppcIO->ymin;
        xstep = (ppcIO->xmax - x) / (double) (ppcIO->width);
        ystep = (ppcIO->ymax - y) / (double) (ppcIO->height);
        count = 0.0;
        for (j = 0; j < ppcIO->height; j++) {
                for (i = 0; i < ppcIO->width; i++) {
                        x1 = x;
                        y1 = y;
                        iter = 0;
                        xs = x1 * x1;
                        ys = y1 * y1;
                        while ((iter++ < ppcIO->maxiter) && (fabs(xs + ys) < 4.0)) {
                                y1 = 2.0 * x1 * y1 + y;
                                x1 = xs - ys + x;
                                xs = x1 * x1;
                                ys = y1 * y1;
                        }
                        plot_pixel(i, j, iter, ppcIO);
                        count += (double) iter;
                        x += xstep;
                }
                x = ppcIO->xmin;
                y += ystep;

#ifndef NO_SIGNALS
                /* check if someone shot us down */
                CheckBreak;
#endif
        }
        ppcIO->average = count / (double) (ppcIO->width * ppcIO->height);
}

VOID mandcubed(struct ppcIO *ppcIO)
{
        ULONG   i, j, iter;
        double  x1, y1, x2, xs, ys, x, y, xstep, ystep, count;

        x     = ppcIO->xmin;
        y     = ppcIO->ymin;
        xstep = (ppcIO->xmax - ppcIO->xmin) / (double) (ppcIO->width);
        ystep = (ppcIO->ymax - ppcIO->ymin) / (double) (ppcIO->height);
        count = 0.0;
        for (j = 0; j < ppcIO->height; j++) {
                for (i = 0; i < ppcIO->width; i++) {
                        x1 = x;
                        y1 = y;
                        iter = 0;
                        xs = x1 * x1;
                        ys = y1 * y1;
                        while ((iter++ < ppcIO->maxiter) && (fabs(xs + ys) < 4.0)) {
                                x2 = xs * x1 - 3.0 * x1 * ys + x;
                                y1 = 3.0 * xs * y1 - ys * y1 + y;
                                x1 = x2;
                                xs = x1 * x1;
                                ys = y1 * y1;
                        }
                        plot_pixel(i, j, iter, ppcIO);
                        count += (double) iter;
                        x += xstep;
                }
                x = ppcIO->xmin;
                y += ystep;

#ifndef NO_SIGNALS
                /* check if someone shot us down */
                CheckBreak;
#endif
        }
        ppcIO->average = count / (double) (ppcIO->width * ppcIO->height);
}

VOID mandfourth(struct ppcIO *ppcIO)
{
        ULONG   i, j, iter;
        double  x1, y1, x2, xs, ys, x, y, xstep, ystep, count;

        x     = ppcIO->xmin;
        y     = ppcIO->ymin;
        xstep = (ppcIO->xmax - ppcIO->xmin) / (double) (ppcIO->width);
        ystep = (ppcIO->ymax - ppcIO->ymin) / (double) (ppcIO->height);
        count = 0.0;
        for (j = 0; j < ppcIO->height; j++) {
                for (i = 0; i < ppcIO->width; i++) {
                        x1 = x;
                        y1 = y;
                        iter = 0;
                        xs = x1 * x1;
                        ys = y1 * y1;
                        while ((iter++ < ppcIO->maxiter) && (fabs(xs + ys) < 4.0)) {
                                x2 = xs * xs + ys * ys - 6.0 * xs * ys + x;
                                y1 = 4.0 * x1 * y1 * (xs - ys) + y;
                                x1 = x2;
                                xs = x1 * x1;
                                ys = y1 * y1;
                        }
                        plot_pixel(i, j, iter, ppcIO);
                        count += (double) iter;
                        x += xstep;
                }
                x = ppcIO->xmin;
                y += ystep;

#ifndef NO_SIGNALS
                /* check if someone shot us down */
                CheckBreak;
#endif
        }
        ppcIO->average = count / (double) (ppcIO->width * ppcIO->height);
}

VOID mandexpo(struct ppcIO *ppcIO)
{
        ULONG   i, j, k, iter;
        double  x1, y1, x2, y2, xs, ys, x, y, xstep, ystep, count;

        x     = ppcIO->xmin;
        y     = ppcIO->ymin;
        xstep = (ppcIO->xmax - ppcIO->xmin) / (double) (ppcIO->width);
        ystep = (ppcIO->ymax - ppcIO->ymin) / (double) (ppcIO->height);
        count = 0.0;
        for (j = 0; j < ppcIO->height; j++) {
                for (i = 0; i < ppcIO->width; i++) {
                        x1 = x;
                        y1 = y;
                        iter = 0;
                        xs = x1 * x1;
                        ys = y1 * y1;
                        while ((iter++ < ppcIO->maxiter) && (fabs(xs + ys) < 4.0)) {
                                x2 = x1 * (xs * xs - 10.0 * xs * ys + 5.0 * ys * ys);
                                y2 = y1 * (ys * ys - 10.0 * xs * ys + 5.0 * xs * xs);
                                for (k = 5; k < ppcIO->exponent; k++)
                                        cmul(x1, y1, x2, y2);

                                x1 = x2 + x;
                                y1 = y2 + y;

                                xs = x1 * x1;
                                ys = y1 * y1;
                        }
                        plot_pixel(i, j, iter, ppcIO);
                        count += (double) iter;
                        x += xstep;
                }
                x = ppcIO->xmin;
                y += ystep;

#ifndef NO_SIGNALS
                /* check if someone shot us down */
                CheckBreak;
#endif
        }
        ppcIO->average = count / (double) (ppcIO->width * ppcIO->height);
}

VOID julia(struct ppcIO *ppcIO)
{
        ULONG   i, j, iter;
        double  x1, y1, xs, ys, x, y, jx, jy, xstep, ystep, count;

        jx    = ppcIO->xseed;
        jy    = ppcIO->yseed;
        x     = ppcIO->xmin;
        y     = ppcIO->ymin;
        xstep = (ppcIO->xmax - ppcIO->xmin) / (double) (ppcIO->width);
        ystep = (ppcIO->ymax - ppcIO->ymin) / (double) (ppcIO->height);
        count = 0.0;
        for (j = 0; j < ppcIO->height; j++) {
                for (i = 0; i < ppcIO->width; i++) {
                        x1 = x;
                        y1 = y;
                        iter = 0;
                        xs = x1 * x1;
                        ys = y1 * y1;
                        while ((iter++ < ppcIO->maxiter) && (fabs(xs + ys) < 4.0)) {
                                y1 = 2.0 * x1 * y1 + jy;
                                x1 = xs - ys + jx;
                                xs = x1 * x1;
                                ys = y1 * y1;
                        }
                        plot_pixel(i, j, iter, ppcIO);
                        count += (double) iter;
                        x += xstep;
                }
                x = ppcIO->xmin;
                y += ystep;

#ifndef NO_SIGNALS
                /* check if someone shot us down */
                CheckBreak;
#endif
        }
        ppcIO->average = count / (double) (ppcIO->width * ppcIO->height);
}

VOID juliacubed(struct ppcIO *ppcIO)
{
        ULONG   i, j, iter;
        double  x1, y1, x2, xs, ys, x, y, jx, jy, xstep, ystep, count;

        jx    = ppcIO->xseed;
        jy    = ppcIO->yseed;
        x     = ppcIO->xmin;
        y     = ppcIO->ymin;
        xstep = (ppcIO->xmax - ppcIO->xmin) / (double) (ppcIO->width);
        ystep = (ppcIO->ymax - ppcIO->ymin) / (double) (ppcIO->height);
        count = 0.0;
        for (j = 0; j < ppcIO->height; j++) {
                for (i = 0; i < ppcIO->width; i++) {
                        x1 = x;
                        y1 = y;
                        iter = 0;
                        xs = x1 * x1;
                        ys = y1 * y1;
                        while ((iter++ < ppcIO->maxiter) && (fabs(xs + ys) < 4.0)) {
                                x2 = xs * x1 - 3.0 * x1 * ys + jx;
                                y1 = 3.0 * xs * y1 - ys * y1 + jy;
                                x1 = x2;
                                xs = x1 * x1;
                                ys = y1 * y1;
                        }
                        plot_pixel(i, j, iter, ppcIO);
                        count += (double) iter;
                        x += xstep;
                }
                x = ppcIO->xmin;
                y += ystep;

#ifndef NO_SIGNALS
                /* check if someone shot us down */
                CheckBreak;
#endif
        }
        ppcIO->average = count / (double) (ppcIO->width * ppcIO->height);
}

VOID juliafourth(struct ppcIO *ppcIO)
{
        ULONG   i, j, iter;
        double  x1, y1, x2, xs, ys, x, y, jx, jy, xstep, ystep, count;

        jx    = ppcIO->xseed;
        jy    = ppcIO->yseed;
        x     = ppcIO->xmin;
        y     = ppcIO->ymin;
        xstep = (ppcIO->xmax - ppcIO->xmin) / (double) (ppcIO->width);
        ystep = (ppcIO->ymax - ppcIO->ymin) / (double) (ppcIO->height);
        count = 0.0;
        for (j = 0; j < ppcIO->height; j++) {
                for (i = 0; i < ppcIO->width; i++) {
                        x1 = x;
                        y1 = y;
                        iter = 0;
                        xs = x1 * x1;
                        ys = y1 * y1;
                        while ((iter++ < ppcIO->maxiter) && (fabs(xs + ys) < 4.0)) {
                                x2 = xs * xs + ys * ys - 6.0 * xs * ys + jx;
                                y1 = 4.0 * x1 * y1 * (xs - ys) + jy;
                                x1 = x2;
                                xs = x1 * x1;
                                ys = y1 * y1;
                        }
                        plot_pixel(i, j, iter, ppcIO);
                        count += (double) iter;
                        x += xstep;
                }
                x = ppcIO->xmin;
                y += ystep;

#ifndef NO_SIGNALS
                /* check if someone shot us down */
                CheckBreak;
#endif
        }
        ppcIO->average = count / (double) (ppcIO->width * ppcIO->height);
}

VOID juliaexpo(struct ppcIO *ppcIO)
{
        ULONG   i, j, k, iter;
        double  x1, y1, x2, y2, xs, ys, x, y, jx, jy, xstep, ystep, count;

        jx    = ppcIO->xseed;
        jy    = ppcIO->yseed;
        x     = ppcIO->xmin;
        y     = ppcIO->ymin;
        xstep = (ppcIO->xmax - ppcIO->xmin) / (double) (ppcIO->width);
        ystep = (ppcIO->ymax - ppcIO->ymin) / (double) (ppcIO->height);
        count = 0.0;
        for (j = 0; j < ppcIO->height; j++) {
                for (i = 0; i < ppcIO->width; i++) {
                        x1 = x;
                        y1 = y;
                        iter = 0;
                        xs = x1 * x1;
                        ys = y1 * y1;
                        while ((iter++ < ppcIO->maxiter) && (fabs(xs + ys) < 4.0)) {
                                x2 = x1 * (xs * xs - 10.0 * xs * ys + 5.0 * ys * ys);
                                y2 = y1 * (ys * ys - 10.0 * xs * ys + 5.0 * xs * xs);
                                for (k = 5; k < ppcIO->exponent; k++)
                                        cmul(x1, y1, x2, y2);

                                x1 = x2 + jx;
                                y1 = y2 + jy;

                                xs = x1 * x1;
                                ys = y1 * y1;
                        }
                        plot_pixel(i, j, iter, ppcIO);
                        count += (double) iter;
                        x += xstep;
                }
                x = ppcIO->xmin;
                y += ystep;

#ifndef NO_SIGNALS
                /* check if someone shot us down */
                CheckBreak;
#endif
        }
        ppcIO->average = count / (double) (ppcIO->width * ppcIO->height);
}

VOID dragon(struct ppcIO *ppcIO)
{
        ULONG   i, j, iter;
        double  x1, y1, x2, xs, ys, x, y, jx, jy, xstep, ystep, count;

        jx    = ppcIO->xseed;
        jy    = ppcIO->yseed;
        x     = ppcIO->xmin;
        y     = ppcIO->ymin;
        xstep = (ppcIO->xmax - ppcIO->xmin) / (double) (ppcIO->width);
        ystep = (ppcIO->ymax - ppcIO->ymin) / (double) (ppcIO->height);
        count = 0.0;
        for (j = 0; j < ppcIO->height; j++) {
                for (i = 0; i < ppcIO->width; i++) {
                        x1 = x;
                        y1 = y;
                        iter = 0;
                        xs = x1 * x1;
                        ys = y1 * y1;
                        while ((iter++ < ppcIO->maxiter) && (fabs(xs + ys) < 4.0)) {
                                x2 = (x1 - xs + ys) * jx - (y1 - 2.0 * x1 * y1) * jy;
                                y1 = (x1 - xs + ys) * jy + (y1 - 2.0 * x1 * y1) * jx;

                                x1 = x2;
                                xs = x1 * x1;
                                ys = y1 * y1;
                        }
                        plot_pixel(i, j, iter, ppcIO);
                        count += (double) iter;
                        x += xstep;
                }
                x = ppcIO->xmin;
                y += ystep;

#ifndef NO_SIGNALS
                /* check if someone shot us down */
                CheckBreak;
#endif
        }
        ppcIO->average = count / (double) (ppcIO->width * ppcIO->height);
}

VOID salamander(struct ppcIO *ppcIO)
{
        ULONG   i, j, iter;
        double  x1, y1, x2, xs, ys, x, y, jx, jy, xstep, ystep, count;

        jx    = ppcIO->xseed;
        jy    = ppcIO->yseed;
        x     = ppcIO->xmin;
        y     = ppcIO->ymin;
        xstep = (ppcIO->xmax - ppcIO->xmin) / (double) (ppcIO->width);
        ystep = (ppcIO->ymax - ppcIO->ymin) / (double) (ppcIO->height);
        count = 0.0;
        for (j = 0; j < ppcIO->height; j++) {
                for (i = 0; i < ppcIO->width; i++) {
                        x1 = x;
                        y1 = y;
                        iter = 0;
                        xs = x1 * x1;
                        ys = y1 * y1;
                        while ((iter++ < ppcIO->maxiter) && (fabs(xs + ys) < 4.0)) {
                                x2 = (xs - ys) * jx - (2.0 * x1 * y1) * jy - 1.0;
                                y1 = (xs - ys) * jy + (2.0 * x1 * y1) * jx;

                                x1 = x2;
                                xs = x1 * x1;
                                ys = y1 * y1;
                        }
                        plot_pixel(i, j, iter, ppcIO);
                        count += (double) iter;
                        x += xstep;
                }
                x = ppcIO->xmin;
                y += ystep;

#ifndef NO_SIGNALS
                /* check if someone shot us down */
                CheckBreak;
#endif
        }
        ppcIO->average = count / (double) (ppcIO->width * ppcIO->height);
}

VOID meteors(struct ppcIO *ppcIO)
{
        ULONG   i, j, iter;
        double  x1, y1, x2, xs, ys, x, y, xstep, ystep, count;

        x     = ppcIO->xmin;
        y     = ppcIO->ymin;
        xstep = (ppcIO->xmax - ppcIO->xmin) / (double) (ppcIO->width);
        ystep = (ppcIO->ymax - ppcIO->ymin) / (double) (ppcIO->height);
        count = 0.0;
        for (j = 0; j < ppcIO->height; j++) {
                for (i = 0; i < ppcIO->width; i++) {
                        x1 = x;
                        y1 = y;
                        iter = 0;
                        xs = x1 * x1;
                        ys = y1 * y1;
                        while ((iter++ < ppcIO->maxiter) && (fabs(xs + ys) < 4.0)) {
                                x2 = x * x1 - y * y1 - (x * x1 + y * y1) / (xs + ys);
                                y1 = x * y1 + y * x1 + (x * y1 - y * x1) / (xs + ys);

                                x1 = x2;
                                xs = x1 * x1;
                                ys = y1 * y1;
                        }
                        plot_pixel(i, j, iter, ppcIO);
                        count += (double) iter;
                        x += xstep;
                }
                x = ppcIO->xmin;
                y += ystep;

#ifndef NO_SIGNALS
                /* check if someone shot us down */
                CheckBreak;
#endif
        }
        ppcIO->average = count / (double) (ppcIO->width * ppcIO->height);
}

VOID pearls(struct ppcIO *ppcIO)
{
        ULONG   i, j, iter;
        double  x1, y1, x2, xs, ys, x, y, jx, jy, xstep, ystep, count;

        jx    = ppcIO->xseed;
        jy    = ppcIO->yseed;
        x     = ppcIO->xmin;
        y     = ppcIO->ymin;
        xstep = (ppcIO->xmax - ppcIO->xmin) / (double) (ppcIO->width);
        ystep = (ppcIO->ymax - ppcIO->ymin) / (double) (ppcIO->height);
        count = 0.0;
        for (j = 0; j < ppcIO->height; j++) {
                for (i = 0; i < ppcIO->width; i++) {
                        x1 = x;
                        y1 = y;
                        iter = 0;
                        xs = x1 * x1;
                        ys = y1 * y1;
                        while ((iter++ < ppcIO->maxiter) && (fabs(xs + ys) < 4.0)) {
                                x2 = jx * x1 - jy * y1 - ((jx * x1 + jy * y1) / (xs + ys));
                                y1 = jx * y1 + jy * x1 + ((jx * y1 - jy * x1) / (xs + ys));

                                x1 = x2;
                                xs = x1 * x1;
                                ys = y1 * y1;
                        }
                        plot_pixel(i, j, iter, ppcIO);
                        count += (double) iter;
                        x += xstep;
                }
                x = ppcIO->xmin;
                y += ystep;

#ifndef NO_SIGNALS
                /* check if someone shot us down */
                CheckBreak;
#endif
        }
        ppcIO->average = count / (double) (ppcIO->width * ppcIO->height);
}

UBYTE   form[80];
double  stack[2*20];
ULONG   stackptr;

#define initstack()     \
{                       \
        stackptr = 0;   \
}

#define push(x, y)              \
{                               \
        stack[stackptr++] = x;  \
        stack[stackptr++] = y;  \
}

#define pop(x, y)               \
{                               \
        y = stack[--stackptr];  \
        x = stack[--stackptr];  \
}

typedef struct complex {
        double  re, im;
} complex;

complex perform_formula(struct ppcIO *ppcIO, double zx, double zy, double cx, double cy)
{
        ULONG   i, j, number;
        complex z;
        double  sx, sy, hx, hy, k, pnt;

        initstack();
        i = 0;
        sx = sy = 0.0;
        while (form[i] != 0) {
                k = pnt = 0.0;
                number = FALSE;
                while ((form[i] == 0x2e) || ((form[i] >= 0x30) && (form[i] <= 0x39))) {
                        number = TRUE;
                        if (form[i] == 0x2e) {
                                pnt = 0.1;
                                i++;
                        } else
                                if (pnt == 0.0)
                                        k = k * 10.0 + (form[i++] - 0x30);
                                else {
                                        k = k + pnt * (form[i++] - 0x30);
                                        pnt *= 0.1;
                                }
                }
                if (number)
                        push(k, 0.0);

                switch (form[i++]) {
                        case 0x7a:                                                      /* "z" */
                                push(zx, zy);
                                break;
                        case 0x63:                                                      /* "c" */
                                push(cx, cy);
                                break;
                        case 0x69:                                                      /* "i" */
                                push(0.0, 1.0);
                                break;
                        case 0x6e:                                                      /* "n" */
                                push((double) ppcIO->exponent, 0.0);
                                break;

                        case 0x2b:                                                      /* "+" */
                                pop(hx, hy);
                                pop(sx, sy);
                                cadd(hx, hy, sx, sy);
                                push(sx, sy);
                                break;
                        case 0x2d:                                                      /* "-" */
                                pop(hx, hy);
                                pop(sx, sy);
                                csub(hx, hy, sx, sy);
                                push(sx, sy);
                                break;
                        case 0x2a:                                                      /* "*" */
                                pop(hx, hy);
                                pop(sx, sy);
                                cmul(hx, hy, sx, sy);
                                push(sx, sy);
                                break;
                        case 0x2f:                                                      /* "/" */
                                pop(hx, hy);
                                pop(sx, sy);
                                cdiv(hx, hy, sx, sy);
                                push(sx, sy);
                                break;
                        case 0x5e:                                                      /* "^" */
                                pop(k, hy);
                                pop(sx, sy);
                                if (k == (double) ((ULONG) k)) {
                                        if (k == 2.0)
                                                csqr(sx, sy, hx, hy);
                                        if (k == 3.0)
                                                ccube(sx, sy, hx, hy);
                                        if (k >= 4.0) {
                                                hx = sx * sx * sx * sx + sy * sy * sy * sy - 6.0 * sx * sx * sy * sy;
                                                hy = 4.0 * sx * sy * (sx * sx - sy * sy);
                                                for (j = 4; j < k; j++)
                                                        cmul(sx, sy, hx, hy);
                                        }
                                        push(hx, hy);
                                } else {
                                        cpow(k, sx, sy);
                                        push(sx, sy);
                                }
                                break;

                        case 0x78:                                                      /* "x" */
                                push(zx, 0.0);
                                break;
                        case 0x79:                                                      /* "y" */
                                push(0.0, zy);
                                break;
                        case 0x61:                                                      /* "a" */
                                push(cx, 0.0);
                                break;
                        case 0x62:                                                      /* "b" */
                                push(0.0, cy);
                                break;

                        case 0xb2:                                                      /* "²" */
                                pop(hx, hy);
                                csqr(hx, hy, sx, sy);
                                push(sx, sy);
                                break;
                        case 0xb3:                                                      /* "³" */
                                pop(hx, hy);
                                ccube(hx, hy, sx, sy);
                                push(sx, sy);
                                break;

                }
        }
        pop(sx, sy);
        z.re = sx;
        z.im = sy;

        return z;
}

VOID formula_mand(struct ppcIO *ppcIO)
{
        ULONG   i, j, iter;
        complex z;
        double  x1, y1, xs, ys, x, y, xstep, ystep, count;

        x     = ppcIO->xmin;
        y     = ppcIO->ymin;
        xstep = (ppcIO->xmax - ppcIO->xmin) / (double) (ppcIO->width);
        ystep = (ppcIO->ymax - ppcIO->ymin) / (double) (ppcIO->height);
        count = 0.0;
        for (j = 0; j < ppcIO->height; j++) {
                for (i = 0; i < ppcIO->width; i++) {
                        x1 = x;
                        y1 = y;
                        iter = 0;
                        xs = x1 * x1;
                        ys = y1 * y1;
                        while ((iter++ < ppcIO->maxiter) && (fabs(xs + ys) < 4.0)) {
                                z = perform_formula(ppcIO, x1, y1, x, y);
                                x1 = z.re;
                                y1 = z.im;
                                xs = x1 * x1;
                                ys = y1 * y1;
                        }
                        plot_pixel(i, j, iter, ppcIO);
                        count += (double) iter;
                        x += xstep;
                }
                x = ppcIO->xmin;
                y += ystep;

#ifndef NO_SIGNALS
                /* check if someone shot us down */
                CheckBreak;
#endif
        }
        ppcIO->average = count / (double) (ppcIO->width * ppcIO->height);
}

VOID formula_julia(struct ppcIO *ppcIO)
{
        ULONG   i, j, iter;
        complex z;
        double  x1, y1, xs, ys, x, y, jx, jy, xstep, ystep, count;

        jx    = ppcIO->xseed;
        jy    = ppcIO->yseed;
        x     = ppcIO->xmin;
        y     = ppcIO->ymin;
        xstep = (ppcIO->xmax - ppcIO->xmin) / (double) (ppcIO->width);
        ystep = (ppcIO->ymax - ppcIO->ymin) / (double) (ppcIO->height);
        count = 0.0;
        for (j = 0; j < ppcIO->height; j++) {
                for (i = 0; i < ppcIO->width; i++) {
                        x1 = x;
                        y1 = y;
                        iter = 0;
                        xs = x1 * x1;
                        ys = y1 * y1;
                        while ((iter++ < ppcIO->maxiter) && (fabs(xs + ys) < 4.0)) {
                                z = perform_formula(ppcIO, x1, y1, jx, jy);
                                x1 = z.re;
                                y1 = z.im;
                                xs = x1 * x1;
                                ys = y1 * y1;
                        }
                        plot_pixel(i, j, iter, ppcIO);
                        count += (double) iter;
                        x += xstep;
                }
                x = ppcIO->xmin;
                y += ystep;

#ifndef NO_SIGNALS
                /* check if someone shot us down */
                CheckBreak;
#endif
        }
        ppcIO->average = count / (double) (ppcIO->width * ppcIO->height);
}

VOID formula(struct ppcIO *ppcIO)
{
        LONG    i, j, k, len, ka, kz;
        ULONG   offset[80];
        UBYTE   operator;

        i = 0;
        while ((form[i] = ppcIO->formel[i]) != 0)
                i++;
        form[i + 1] = 0;

        for (i = 0; form[i] != 0; i++) {
                j = 0;
                while ((form[i + j] == 0x2e) || ((form[i + j] >= 0x30) && (form[i + j] <= 0x39)))
                        j++;
                if (j > 0)
                        offset[i] = j;
                else
                        offset[i] = 1;
        }
        len = i;

        while (offset[0] < len) {
                i = 0;
                ka = kz = -1;
                while (i < len) {
                        if (form[i] == 0x28)                                            /* "(" */
                                ka = i;
                        if ((form[i] == 0x29) && (ka > -1)){                            /* ")" */
                                kz = i;
                                form[ka] = 0x20;
                                form[kz] = 0x20;
                                offset[ka] = kz - ka + 1;
                                j = ka + 1;
                                while ((form[j] != 0x2b) && (form[j] != 0x2d) && (form[j] != 0x2a) && (form[j] != 0x2f) && (form[j] != 0x5e))
                                        j += offset[j];
                                operator = form[j];
                                for (k = j; k < kz - 1; k++)
                                        form[k] = form[k + 1];
                                form[kz - 1] = operator;
                                ka = -1;
                        }
                        i += offset[i];
                }
        }

        for (i = 0; i < len; i++)
                while (form[i] == 0x20) {
                        for (j = i + 1; j < len; j++)
                                form[j - 1] = form[j];
                        form[len - 1] = 0;
                        len--;
                }

        switch (ppcIO->style) {
                case 0:
                        formula_mand(ppcIO);
                        break;
                case 1:
                        formula_julia(ppcIO);
                        break;
        }
}
