/* CRCW Quickhull algorithm, DC strategy                     *
 * computes convex hull of N points using p processors       *
 * in expected time O((N/p)log N).      C.W. Kessler 10/96   */

/* FOR TIME MEASUREMENT, DELETE THE CALLS TO line() ---
 * THEY ACCOUNT FOR MORE THAN 95 % OF THE RUN TIME!  */

#include <fork.h>
#include <assert.h>
#include <io.h>
#include <stdlib.h>
#include <math.h>
#include <graphic.h>

sh int N;            /* the number of points */
sh double *x, *y;    /* x and y point coordinates */
sh int *belongs_to_hull;
               /* the result vector; 
                  indicates whether a point belongs to the hull*/
sh simple_lock screen = 0;
#define PIXELS 512      /*extents of generated X-Pixmap picture*/
sh double factor = 0.00000024;  /*fuer 512*/

/* output routines: */

async void set_point( int x, int y, int width, int color )
{
 pr int i, j;
 for (i=1-width; i<width; i++)
   for (j=1-width; j<width; j++)
     set_pixel( x+i, /*PIXELS -*/ y+j, color );
}


async void print_array( pr int *pt, pr int n )
{
 pr int j;
 simple_lockup( &screen );
 printf("Array with %d points:\n", n );
 for (j=0; j<n; j++)
   printf("point %d = (%d,%d)\n", pt[j], (int)x[pt[j]], (int)y[pt[j]] );
 prS("\n");
 simple_unlock( &screen );
}


async void print_hull( void )
{
 pr int j;
 prS("Points of the convex hull:\n");
 for (j=0; j<N; j++) 
   if (belongs_to_hull[j]) {
      printf(" (%d,%d)", (int)x[j], (int)y[j] );
      set_point( x[j], y[j], 2, 2 );
   }
}


/* quickhull algorithm, applied to set S of n points (x[],y[]) */

#define cross(p,a,b) \
 /* return crossproduct of vectors p-a and b-a */ \
 ((x[a]-x[p])*(y[b]-y[p]) - (y[a]-y[p])*(x[b]-x[p]))
 /* see G. Blelloch's paper in CACM 39(3) p.95, 1996 */

#define leftturn(a,b,c) \
 /*true iff point c is lefthand of vector through points (a,b) */ \
 (cross(c,a,b)>0.0)

#define est_work(n) \
 /*float value estimating work to be done for size-n problem*/\
 ((float)((n) * (ilog2(n)+1)))

sync int *delete_right( sh int *pt, sh int *num, sh int p1, sh int p2 )
{
 pr int j;
 sh int leftcnt;
 sh int *left;
 sh int n = *num;
 sh int p = 0;
 
 $ = mpadd( &p, 1 );   /* renumber and determine number of processors */

 farm pprintf("delete_right P%d,P%d, n=%d\n", p1, p2, n );

 /* delete all points in pt located right from line p1p2: */
 farm if ($==0) left = (int *) shmalloc( n );
 left[0]=p1; left[1]=p2;
 leftcnt = 2;
 for (j=$; j<n; j+=p)
   farm
    if (!(pt[j]==p1 || pt[j]==p2))   /*p1 and p2 already in left[]*/
      if (leftturn(p1,p2,pt[j])) /* point j is lefthand to vector p1p2 */
        left[mpadd(&leftcnt,1)] = pt[j];
 *num = leftcnt;
 return left;
}


sync void inithull( sh int *pt, int n, sh int *minx, sh int *maxx )
{
 sh int p1, p2;
 sh int i;

 /* determine points p1,p2 with minimal and maximal x coordinate: */
 p1 = p2 = pt[0];   /* init. p1,p2 to first point */
 for (i=1; i<n; i++) {   /* seq. search for minimum / maximum */
   if (x[pt[i]] < x[p1])   p1 = i;
   if (x[pt[i]] > x[p2])   p2 = i;
 }
 belongs_to_hull[p1] = 1;
 belongs_to_hull[p2] = 1;
 *minx = p1; *maxx = p2;
} 


sync int seq_pivotize( sh int *pt, sh int n )    /* n>=3 is assumed */
 /* as pivot, select the point in pt[]-{p1,p2} with maximal
  * cross_prod( pivot-p1, p2-p1 )
  */
{
 sh int i, p1=pt[0], p2=pt[1];
 sh int pivotpos = 2;
 sh double maxcross = cross(pt[2], p1, p2);
 for (i=3; i<n; i++) {        /* sequential maximization */
    sh double newcross = cross(pt[i], p1, p2);
    if (newcross > maxcross) { maxcross = newcross; pivotpos = i; }
 }
 return pt[pivotpos];
}


sync int par_pivotize( sh int *pt, sh int n )   /* n>=p+3 is assumed */
{
 sh int p=0;
 sh int p1=pt[0], p2=pt[1];
 pr int mypivotpos;
 pr double mymaxcross;
 sh double *cval;                
 sh int *cidx;
 sh int s;
 pr int i;

 $ = mpadd( &p, 1 );     /* renumber processors within the group */
 mypivotpos = $+2;
 mymaxcross = cross(pt[$+2], p1, p2);
 for (i=p+2+$; i<n; i+=p) {  /* local maximizations */
    pr double mynewcross = cross(pt[i], p1, p2);
    if (mynewcross > mymaxcross)
      { mymaxcross = mynewcross; mypivotpos = i; }
 }
 cval = (double *) shalloc(p*sizeof(double)); /* temp. array for this fn.*/
 cidx = (int *) shalloc(p*sizeof(int));       /* temp. array for this fn. */
 cval[$] = mymaxcross; cidx[$] = mypivotpos;
 for (s=1; s<p; s=s*2)   /* compute max(cval[s],cval[$+s]) in step s */
   farm if (!( $ & s )) {
     if (cval[$] < cval[$+s])
       { cval[$] = cval[$+s]; cidx[$] = cidx[$+s]; }
   }
 shallfree();
 return pt[cidx[0]]; 
}


sync void seq_qh( sh int *pt, sh int n )
{
 /* DC step: select pivot point from pt */
 sh int pivotpos, pivot;
 sh int p1=pt[0], p2=pt[1];
 sh int *left1, *left2, leftcnt1, leftcnt2;

 /* DC step: select any pivot point from pt. We have p1==pt[0],p2==pt[1]*/
 if (n==2) {
    seq seq_line (x[p1],y[p1],x[p2],y[p2],0x00ff00,2);
    return;
 }
 if (n==3) {  /* one point (beyond old p1,p2) must belong to hull */
    belongs_to_hull[pt[2]] = 1;      /* saves a recursive call */
    seq seq_line (x[p1],y[p1],x[pt[2]],y[pt[2]],0xff00ff,2);
    seq seq_line (x[pt[2]],y[pt[2]],x[p2],y[p2],0xff00ff,2);
    return;
 }
 pivot = seq_pivotize( pt, n );
 belongs_to_hull[pivot] = 1;
 seq seq_line (x[p1],y[p1],x[pivot],y[pivot],0xff00ff,2);
 seq seq_line (x[pivot],y[pivot],x[p2],y[p2],0xff00ff,2);
 leftcnt1 = n;
 left1 = delete_right( pt, &leftcnt1, p1, pivot );
 seq_qh( left1, leftcnt1 );
 leftcnt2 = n;
 left2 = delete_right( pt, &leftcnt2, pivot, p2 );
 seq_qh( left2, leftcnt2 ); 
}


sync void qh( sh int *pt, sh int n )
{
 /* DC step: select any pivot point from pt. We have p1==pt[0],p2==pt[1]*/
 sh int pivot;
 sh int p = 0;                /* number of available procs */
 sh int p1=pt[0], p2=pt[1];   /* the line (p1,p2) */
 sh int *S1, *S2, n1, n2;     /* to hold subproblems */
 sh int firstrightproc;

 farm pprintf("qh(%d)\n", n);

 $ = mpadd( &p, 1 );   /*renumber $ and determine p*/
 if (n==2) {
    line (x[p1],y[p1],x[p2],y[p2],0xff,2);  /* draw hull edge */
    return;
 }
 if (n==3) {  /* one point (beyond old p1,p2) must belong to hull */
    belongs_to_hull[pt[2]] = 1;
    line (x[p1],y[p1],x[pt[2]],y[pt[2]],0xff00ff,2); /* draw hull edge */
    line (x[pt[2]],y[pt[2]],x[p2],y[p2],0xff00ff,2); /* draw hull edge */
    return;
 }
 if (p==1) { seq_qh( pt, n ); return; }

 /* as pivot, select the point in pt[]-{p1,p2} with maximal
  * cross_prod( pivot-p1, p2-p1 )
  */
 if (n>p+2)  pivot = par_pivotize( pt, n );
 else        pivot = seq_pivotize( pt, n );
 belongs_to_hull[pivot] = 1;
 line (x[p1],y[p1],x[pivot],y[pivot],0x00ff00,2); /* draw edge (p1,pivot) */
 line (x[pivot],y[pivot],x[p2],y[p2],0x00ff00,2); /* draw edge (pivot,p2) */
 /* determine subproblems: */
 n1 = n2 = n;
 S1 = delete_right( pt, &n1, p1, pivot );
 S2 = delete_right( pt, &n2, pivot, p2 );

 /* determine amount of processors dedicated to each sub-problem:*/ 
 firstrightproc = (int) ((float)p * est_work(n1) / est_work(n) ); 
 if (firstrightproc == 0 && n1 > 0) firstrightproc = 1;

 if ($ < firstrightproc)   /*split group of processors: */
    qh( S1, n1 );   else     qh( S2, n2 ); 
}


sh unsigned int starttime, stoptime;

void main( void )
{
 if (__PROC_NR__==0) {
    printf("Enter N = ");
    scanf("%d", &N );
    printf("\ngenerate %d random points.\n", N);
 }
 srand( 1717171 );    /* seed random generator */
 start {
   sh int p = groupsize();
   sh int *pt = (int *) shalloc( N * sizeof(int) );
   sh int minxpt, maxxpt;
   sh int *upper = (int *) shalloc( N * sizeof(int) );
   sh int *lower = (int *) shalloc( N * sizeof(int) );
   sh int uppercnt, lowercnt;
   pr int j;
   sh double xj;

   init_pict( PIXELS+1, PIXELS+1 );
   clear_pixels( 1 );

   belongs_to_hull = (int *) shalloc( N );
   x = (double *) shalloc( N * sizeof(double) );
   y = (double *) shalloc( N * sizeof(double) );
      
   seq {
     /* sequential to get equal problems for varying p */
     for (j=0; j<N; j++) {
#define ABS(x) (((x)>=0)?(x):(-(x)))
       xj = ((double) rand()) * factor;
       x[j] = ABS(xj);
       xj = ((double) rand()) * factor;
       y[j] = ABS(xj);
       pt[j] = j;
       set_point( (int) x[j], (int) y[j], 2, 3 );
     }
     starttime = getct();
   }
   for (j=$; j<N; j+=p) 
       belongs_to_hull[j] = 0;            /* initialization */
   inithull( pt, N, &minxpt, &maxxpt );

   /* now split original set of points:
    * into upper[] (all nodes left of (p1,p2), including p1,p2
    * and  lower[] (all points right of (p1,p2), including p1,p2
    */
   uppercnt = lowercnt = N;
   upper = delete_right( pt, &uppercnt, minxpt, maxxpt );
   lower = delete_right( pt, &lowercnt, maxxpt, minxpt );
   line (x[minxpt],y[minxpt],x[maxxpt],y[maxxpt],0xffff,2);
   qh( upper, uppercnt );
   qh( lower, lowercnt );

   seq {
      stoptime = getct();
      print_hull();
      printf("\nTime (including line()): %d PRAM CPU Cycles\nWait for picture ...\n",
             stoptime - starttime );  
   }
   write_pixmap( "HULL" );
 }
 barrier; exit(0);
}