/*
 * C version of prime.f, a program to calculate 2^859433 - 1.
 * The algorithm is credited to Slowinski of CRI.
 *
 * compile with:
 *
 * cc prime.c -o prime -lm
 *
 * if FTIME is defined, ftime() is used to obtaining timing data.
 * otherwise, gettimeofday() is used.
 */

#include <stdio.h>
#include <math.h>
#include <malloc.h>
#include <sys/types.h>

#ifdef FTIME
#include <sys/timeb.h>
#else
#include <sys/time.h>
#endif

/* #define EXP 859433*/ /* will calculate 2^859433 - 1, the 33rd Mersenne prime */
/* #define EXP 216091 */ /* will calculate 2^216091 - 1, the 31st Mersenne prime */
#define EXP 132049 /* will calculate 2^132049 - 1, the 30th Mersenne prime */

main(argc, argv)
int argc;
char *argv[];
{
#ifdef FTIME
   struct timeb tp0, tp1;
#else
   struct timeval tp0, tp1;
#endif
   double t, ftp0, ftp1;
   double c0 = 0, c1 = 0, c2 = 0;
   int n;
   int *a, *b, *c;
   int ddigit = 4;     /* number decimal digits per part */
   int vdig = (int)pow((double)10, (double)ddigit); /* decimal value of part */
   int i, j, k, maxpart, mx;
   unsigned int size;
   char format[64];

   n = EXP;
   if (argc == 2)  /* command-line param for new n? */
      n = (atoi(argv[1]) > 0) ? atoi(argv[1]) : EXP;
   printf("calculating 2^%d - 1\n", n);

   size = (n+1)*sizeof(int);
   printf ("allocating %d bytes for working arrays\n", 3*size);

   if ((a = (int *)malloc(size)) == NULL) {
      fprintf(stderr, "cannot allocate %d bytes for a\n", size);
      exit(1);
   }
   if ((b = (int *)malloc(size)) == NULL) {
      fprintf(stderr, "cannot allocate %d bytes for b\n", size);
      exit(1);
   }
   if ((c = (int *)malloc(size)) == NULL) {
      fprintf(stderr, "cannot allocate %d bytes for c\n", size);
      exit(1);
   }

   for (i = 1; i <= n; i++) {
      a[i] = 0; /* set results array to null */
      b[i] = 2; /* load up multiplicand array with n 2s */
      c[i] = 0; /* set carry array to null */
   }

   a[1] = 1; /* set initial intermediate result to 1 */
   mx = 1;

#ifdef FTIME
   ftime(&tp0);
#else
   gettimeofday (&tp0, (struct timezone *)NULL);
#endif

#pragma _CRI parallel shared(n, a, b, c, mx, vdig) private(i, j)
   for (j = 1; j <= n; j++) {

#pragma _CRI taskloop vector
      for (i = 1; i <= mx; i++) {
         a[i] *= b[j];
      }

#pragma _CRI taskloop vector
      for (i = 2; i <= mx+1; i++)
         c[i] = a[i-1]/vdig;

#pragma _CRI guard
      if (c[mx+1] > 0)
         mx++;
#pragma _CRI endguard

#pragma _CRI taskloop vector
      for (i = 1; i <= mx; i++) {
         a[i] = (a[i] % vdig) + c[i];
      }
   }
#pragma _CRI endparallel

   /* propogate carry through one last time */
   j = 0;
   for (i = 1; i <= n; i++) {
      k = a[i] + j;
      a[i] = k % vdig;
      j = k/vdig;
   }

#ifdef FTIME
   ftime(&tp1);
   ftp0 = (double)tp0.time+(double)tp0.millitm/(double)1000000;
   ftp1 = (double)tp1.time+(double)tp1.millitm/(double)1000000;
#else
   gettimeofday (&tp1, (struct timezone *)NULL);
   ftp0 = (double)tp0.tv_sec+(double)tp0.tv_usec/(double)1000000;
   ftp1 = (double)tp1.tv_sec+(double)tp1.tv_usec/(double)1000000;
#endif

   t = (ftp1-ftp0 == 0) ? 0.000001 : ftp1-ftp0;
   printf("wall time for main loop:         %.5lf seconds\n", t);

   maxpart = 0;

   for (i = 0; i <= n; i++)
      if (a[i] != 0)
         maxpart = i;

   printf ("max parts used: %d (each part is %d decimal digits)\n",
      maxpart, ddigit);

   /* subtract one from answer and see if we match Slow's number */

   a[1]--;

   if (a[1] < 0) {

      a[1] = vdig - 1;

      for (i = 2; i <= maxpart; i++) {
         a[i]--;
         if (a[i] < 0)
            a[i] = vdig - 1;
         else
            break;
      }
   }

   printf("result:\n");

   sprintf(format, "%s%d%s", "%0", ddigit, "d");
   j = 0;
   for (i = maxpart; i >= 1; i--) {
      printf (format, a[i]);
      j++;
      if (j % 15 == 0)
         printf ("\n");
   }
   printf ("\n");

   return (0);
}