#include "maplec.h"

typedef struct {

   /* The MKernelVector is needed in the tasks */

   MKernelVector kv;

   /* The continuation function sums the two values computed by its

       two children */

   ALGEB left, right;

} AddContArg;

/* Mark the values from the left and right children, if they have

   been computed */

void MarkAddContFunction( void *a )

{

   AddContArg *args;

   args = (AddContArg*)a;

   if ( args->left != NULL )

       MapleGcMark( args->kv, args->left );

   if ( args->right != NULL )

       MapleGcMark( args->kv, args->right );

}

/* The continuation function, sum the values from the two children */

int AddTaskContFunction( void *p, int arg_number, void *a )

{

   AddContArg *args;

   args = (AddContArg*)a;

   switch( arg_number )

   {

       case MAPLE_TASK_ROOT:

           /* If arg_number is 0, then parent is the value passed

               into MapleStartRootTask */

           *(ALGEB*)p = MapleNumericAdd( args->kv, args->left, args->right );

           break;

       case 1:

           /* If arg_number is 1, then we update the left field of

               our parent struct. */

           ((AddContArg*)p)->left = MapleNumericAdd( args->kv, args->left, args->right );

           break;

       case 2:

           /* If arg_number is 2, then we update the right field of

               our parent struct. */

           ((AddContArg*)p)->right = MapleNumericAdd( args->kv, args->left, args->right );

           break;

   }

   free( args );

   return 1;

}

/* The worker task take a range as the argument */

typedef struct {

   MKernelVector kv;

   M_INT start, end;

} AddTaskArg;

/* The main task function */

int AddTaskFunction( void *p, int arg_number, void *a )

{

   AddTaskArg *args;

   args = (AddTaskArg*)a;

   /* check the size of the range, if it is big, divide range in half

       and create tasks for the two new ranges. */

   if ( args->end - args->start > 1000 )

   {

       M_INT mid;

       AddTaskArg *newargs;

       AddContArg *cont;

       mid = (args->end - args->start)/2 + args->start;

       newargs = (AddTaskArg*)malloc( sizeof( AddTaskArg ) );

       cont = (AddContArg*)malloc( sizeof( AddContArg ) );

       newargs->start = args->start;

       newargs->end = mid;

       args->start = mid + 1;

       cont->kv = newargs->kv = args->kv;

       cont->left = cont->right = NULL;

       /* create the continuation task, to sum the results of the

           two child tasks. */

       MapleCreateContinuationTask( args->kv, AddTaskContFunction, cont, MarkAddContFunction );

       /* create the child tasks, there are not garbage collectable

           bits in the task struct, so no Mark function is required */

       MapleStartChildTask( args->kv, 1, AddTaskFunction, newargs, NULL );

       MapleStartChildTask( args->kv, 2, AddTaskFunction, args, NULL );

   }

   else

   {

       M_INT i, t;

       /* compute the sum over this range */

       t = args->start;

       for ( i = args->start+1; i <= args->end; i++ )

       {

           t += i;

       }

       switch ( arg_number )

       {

           case MAPLE_TASK_ROOT:

               /* if this is the root task, then p is the value

                   passed into MapleStartRootTask  */

               *(ALGEB*)p = ToMapleInteger( args->kv, t );

               break;

           case 1:

               /* if arg_number is 1, them we assign the computed

                   value to the left field of the parent struct */

               ((AddContArg*)p)->left = ToMapleInteger( args->kv, t );

               break;

           case 2:

               /* if arg_number is 2, them we assign the computed

                   value to the right field of the parent struct */

               ((AddContArg*)p)->right = ToMapleInteger( args->kv, t );

               break;

       }

       free( args );

   }

   return 1;

}

/* The main entry point for this algorithm */

EXP_DECL ALGEB LANG_STDCALL StartAddTask( MKernelVector kv, ALGEB args )

{

   AddTaskArg *arg;

   ALGEB ret;

   /* create the root task */

   arg = (AddTaskArg*)malloc( sizeof(AddTaskArg) );

   arg->kv = kv;

   arg->start = MapleToM_INT( kv, MapleExpseqSelect( kv, args, 1 ) );

   arg->end = MapleToM_INT( kv, MapleExpseqSelect( kv, args, 2 ) );

   MapleStartRootTask( kv, &ret, AddTaskFunction, arg, NULL, NULL, 0 );

   return ret;

}

struct MandelStruct {

   /* The MKernelVector is needed in the tasks */

   MKernelVector kv;

   /* X and Y are shared among tasks, but used as read only */

   double *X, *Y;

   /* output is shared among tasks, but areas of the table are

     only written to by one task */

   ALGEB output;

   /* the 2D boundaries of the area this task is

     responsible for */

   M_INT xLow, xHigh;

   M_INT yLow, yHigh;

   /* some parameters used in the final computation */

   M_INT iter;

   double bailout;

};

/* need to mark the output array */

void MarkMandelStruct( void *a )

{

   struct MandelStruct *arg = (struct MandelStruct*)a;

   MapleGcMark( arg->kv, arg->output );

}

/* a utility function used for creating the task structures */

struct MandelStruct* NewMandelStruct( struct MandelStruct *old,

       M_INT xL, M_INT xH, M_INT yL, M_INT yH )

{

   struct MandelStruct *newArgs;

   newArgs = (struct MandelStruct*)malloc( sizeof( struct MandelStruct ) );

   newArgs->X = old->X;

   newArgs->Y = old->Y;

   newArgs->output = old->output;

   newArgs->kv = old->kv;

   newArgs->bailout = old->bailout;

   newArgs->iter = old->iter;

   newArgs->xLow = xL;

   newArgs->xHigh = xH;

   newArgs->yLow = yL;

   newArgs->yHigh = yH;

   return newArgs;

}

/* the main task function */

int MandelTaskFunction( void *parent, int arg_number, void *self )

{

   M_INT w,h;

   struct MandelStruct *args;

   /* unused */

   parent = NULL;

   arg_number = 0;

   args = (struct MandelStruct*)self;

   w = args->xHigh - args->xLow;

   h = args->yHigh - args->yLow;

   /* check the area of the region, if too large divide it up */

   if ( w * h > 100 )

   {

       struct MandelStruct *newArgs1,*newArgs2,*newArgs3;

       /* create the task structures for the child tasks */

       newArgs1 = NewMandelStruct( args, args->xLow, w/2+args->xLow, args->yLow, h/2+args->yLow );

       /* update the current task structure so it can be reused*/

       args->xLow = newArgs1->xHigh+1;

       args->yLow = newArgs1->yHigh+1;

       newArgs2 = NewMandelStruct( args, newArgs1->xLow, newArgs1->xHigh, args->yLow, args->yHigh );

       newArgs3 = NewMandelStruct( args, args->xLow, args->xHigh, newArgs1->yLow, newArgs1->yHigh );

       /* create a continuation task that does not do anything */

       MapleCreateContinuationTask( args->kv, NULL, NULL, NULL );

       /* start the child tasks  */

       MapleStartChildTask( args->kv, 1, MandelTaskFunction, newArgs1, MarkMandelStruct );

       MapleStartChildTask( args->kv, 2, MandelTaskFunction, newArgs2, MarkMandelStruct );

       MapleStartChildTask( args->kv, 3, MandelTaskFunction, newArgs3, MarkMandelStruct );

       MapleStartChildTask( args->kv, 4, MandelTaskFunction, args, MarkMandelStruct );

   }

   else

   {

       /* compute the Mandelbrot set for the given region */

       double Xtemp, Ytemp, Xc, Yc, Xold, Yold, tmp;

       M_INT k, index[3];

       RTableData val;

       for ( w = args->xLow; w <= args->xHigh; w++ )

       {

           for ( h = args->yLow; h <= args->yHigh; h++ )

           {

               Xtemp = args->X[w-1];

               Ytemp = args->Y[h-1];

               Xc = Xtemp;

               Yc = Ytemp;

               k = 0;

               while ( k < args->iter )

               {

                   Xold = Xtemp;

                   Yold = Ytemp;

                   Xtemp = Xold*Xold-Yold*Yold+Xc;

                   Ytemp = 2*Xold*Yold+Yc;

                   tmp = Xtemp*Xtemp+Ytemp*Ytemp;

                   if ( tmp >= args->bailout )

                   {

                       /* update the output rtable structure */

                       index[0] = h;

                       index[1] = w;

                       index[2] = 1;

                       val.float64 = k;

                       RTableAssign( args->kv, args->output, index, val );

                       index[2] = 2;

                       val.float64 = sqrt(tmp);

                       RTableAssign( args->kv, args->output, index, val );

                       break;

                   }

                   ++k;

               }

           }

       }

       free( args );

   }

   return 1;

}

/* the main Mandelbrot entry point */

EXP_DECL ALGEB LANG_STDCALL Mandelbrot( MKernelVector kv, ALGEB args )

{

   M_INT i;

   M_INT w, h, iter;

   M_INT bounds[6];

   double x1, x2, y1, y2, bailout, *X, *Y;

   struct MandelStruct *arg;

   ALGEB output;

   RTableSettings settings;

   /* decode the arguments */

   w = MapleToM_INT( kv, MapleExpseqSelect( kv, args, 1 ) );

   h = MapleToM_INT( kv, MapleExpseqSelect( kv, args, 2 ) );

   iter = MapleToM_INT( kv, MapleExpseqSelect( kv, args, 3 ) );

   x1 = MapleToFloat64( kv, MapleExpseqSelect( kv, args, 4 ) );

   x2 = MapleToFloat64( kv, MapleExpseqSelect( kv, args, 5 ) );

   y1 = MapleToFloat64( kv, MapleExpseqSelect( kv, args, 6 ) );

   y2 = MapleToFloat64( kv, MapleExpseqSelect( kv, args, 7 ) );

   bailout = MapleToFloat64( kv, MapleExpseqSelect( kv, args, 8 ) );

   /* pre-compute the points */

   X = (double *)malloc( sizeof( double )*w );

   for ( i = 0; i < w; i++ )

   {

       X[i] = x1 + (x2-x1)*(i-1)/(w-1);

   }

   Y = (double *)malloc( sizeof( double )*h );

   for ( i = 0; i < h; i++ )

   {

       Y[i] = y1 + (y2-y1)*(i-1)/(h-1);

   }

   /* create the root task */

   arg = (struct MandelStruct*)malloc( sizeof( struct MandelStruct ) );

   arg->X = X;

   arg->Y = Y;

   /* create the output rtable */

   RTableGetDefaults( kv, &settings );

   settings.data_type = RTABLE_FLOAT64;

   settings.num_dimensions = 3;

   bounds[0] = 1;

   bounds[1] = h;

   bounds[2] = 1;

   bounds[3] = w;

   bounds[4] = 1;

   bounds[5] = 2;

   arg->output = output = RTableCreate( kv, &settings, NULL, bounds );

   arg->xLow = 1;

   arg->xHigh = w;

   arg->yLow = 1;

   arg->yHigh = h;

   arg->iter = iter;

   arg->bailout = bailout;

   arg->kv = kv;

   /* start the root task */

   MapleStartRootTask( kv, NULL, MandelTaskFunction, arg, MarkMandelStruct, NULL, 0 );

   free( X );

   free( Y );

   return output;

}