/* pprime.c :   gcc -O3 -o pprime pprime.c
   recherche d'un polynôme f(x) = a*x^2 + b*x + c
   tel que f(0), ... , f(n) soient premiers pour n aussi grand que
   possible.
   prendre trois nombres premiers, en déduire a, b et c
   f(0) = c              = i  est premier
   f(1) = a + b + c      = j  est premier
   f(2) = 4*a + 2*b + c  = k  est premier
   a = (i+k)/2  -  j;
   b = 2*j - (3*i+k)/2;
   c = i;
   d1 = f(1) - f(0), d2 = f(2) -2*f(1) + f(0)

   Jean-Paul Davalan <jpdvl@wanadoo.fr>

   programme
   http://perso.wanadoo.fr/jean-paul.davalan/arit/pprime/pprime.c
   
   page web de présentation
   http://perso.wanadoo.fr/jean-paul.davalan/arit/pprime/index.html
 */
 
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
 
#define MINLEN         30
#define isprime(z) ((cr[(z)>>5] & (1<<((z) & 31)))==0)
typedef struct {
  int     a, b, c, l;
} poly;
poly    poltab[100000];
int     curtab = 0, tab[100], mtab[100];
void    polyprint(FILE * f, int a, int b, int c);              
 
int     main(int argc, char *argv[])
{
  int    *cr, ii, rg, r, i1, j1, k1, a0, b0, c0, a1, b1, c1, i, j,
          k, jabs, kabs, a, b, c, n, u, l, Max, max, max1, max2, record = 0,
          d1, d2, t1, t2, N, len, licite, new, tbegin, tlast;
  float   pcent = 10.0, pcentcur;
  char    st[200], *p, filename[] = "results.out";
  FILE   *fileout = fopen(filename, "w");
  tbegin = time(NULL);
  if (fileout == NULL)
    fileout = stderr;
/* Lecture de la ligne de commande
   -------------------------------
 */
  if (argc < 3) {
    fprintf(stderr, "usage   : N n1 [n2 n3]\n"
       "  N premiers (crible d'Erathosthène)\n"
       "  f(0) <= n1, -n2 <= f(1) <= n2, -n3 <= f(2) <= n3\n"
       "exemple : %s 5000000 3000\n", argv[0]);
    exit(0);
  }
  Max = atoi(argv[1]);
  max = max1 = max2 = atoi(argv[2]);
  if (argc > 3)
    max1 = max2 = atoi(argv[3]);
  if (argc > 4)
    max2 = atoi(argv[4]);
  max = (max % 2 == 0) ? max - 1 : max;
  max1 = (max1 % 2 == 0) ? max1 - 1 : max1;
  max2 = (max2 % 2 == 0) ? max2 - 1 : max2;
/* écriture de la ligne de commande dans le fichier
   ------------------------------------------------                         
 */
  fprintf(fileout, "# ");
  for (ii = 0; ii < argc; ii++)
    fprintf(fileout, " %s", argv[ii]);
  fprintf(fileout, "\n"
     "# ligne longueur \"a*x^2 + b*x + c\"\n");
/* construction du crible d'Érathostène
   ------------------------------------
   élimine 0, 1 et 2 pour n'avoir que des premiers impairs
 */
  Max = ((Max + 31) >> 5) << 5;
  cr = (int *)malloc((Max >> 5) * sizeof(int));
  memset(cr, 0, Max >> 5);      /* vrai = 0 = premier */
  j = 2;
  do {
    for (i = j * j; i < Max; i += j) {
      cr[i >> 5] |= (1 << (i & 31));
    }
    while ((cr[++j >> 5] & (1 << (j & 31))) != 0);
  } while (j * j < Max);
  cr[0] |= 7;                   // 0, 1, 2
/* parcours i = f(0) < max, j = f(1) < max1, k = f(2) < max2
   ---------------------------------------------------------
   f(0) ... sont sauvegardés dans tab[]
 */
  for (i = 3; i <= max; i += 2) {
    if (!isprime(i))
      continue;
    pcentcur = (i - 3.0) * 100.0 / (max - 3.0);
    if (pcentcur >= pcent) {
      printf(" %.1f %%", pcentcur);
      pcent += 10.0;
    } else
      printf(";");
    tab[0] = i;
    for (j = -max1; j <= max1; j += 2) {                               
      jabs = (j < 0) ? -j : j;
      if (j == i || !isprime(jabs))
        continue;
      tab[1] = j;
      for (k=2*j-i+2; k <= max2; k += 2) {
        kabs = (k < 0) ? -k : k;
        if (k == i || k == j || !isprime(kabs))
          continue;
        tab[2] = k;
        d2 = k - 2 * j + i;
        d1 = k - j;
        n = k;
        l = 3;
        licite = 1;
        a1 = a = (i + k) / 2 - j;
        b1 = b = 2 * j - (3 * i + k) / 2;
        c1 = c = i;
    /* la variable croît 3, 4, ... */
        do {
          d1 += d2;
          n += d1;
          N = (n < 0) ? -n : n;
          tab[l] = n;
      /* vérification de la valeur obtenue dans la 1ère table */
          licite = (N > 2 && N < Max && isprime(N));
          for (t1 = 0; licite && t1 < l; t1++)  /* vérif. dans la table */
            licite = (tab[t1] != n && tab[t1] != -n);
          if (licite)
            l++;
        } while (licite);
    /* la variable décroît -1, -2 ... les valeurs sont mises dans mtab[] */
        for (rg = 0, ii = 1, licite = 1, d1 = j - i, u = c; licite; ii++) {
          d1 -= d2;             /* variante (différences finies) */
          u -= d1;
          N = (u < 0) ? -u : u;                          
          mtab[ii] = u;
      /* vérification de la valeur obtenue, dans les deux tables */
          licite = (N < Max && isprime(N));
          for (t1 = 0; licite && t1 < l; t1++)
            licite = (tab[t1] != u && tab[t1] != -u);
          for (t1 = 1; licite && t1 < ii; t1++)
            licite = (mtab[t1] != u && mtab[t1] != -u);
          if (licite)
            rg++;
        }
        len = l + rg;
        if (len > record || len >= MINLEN) {
          p = &st[0];
          st[0] = '\x0';
          new = 1;
          if (rg > 0) {
            r = -rg;
            i1 = (a * r + b) * r + c;
            j1 = (a * ++r + b) * r + c;
            k1 = (a * ++r + b) * r + c;
            a1 = (i1 + k1) / 2 - j1;
            b1 = 2 * j1 - (3 * i1 + k1) / 2;
            c1 = i1;
          }
      /* mise dans la table des polynômes */
          for (u = 0, licite = 1; licite && u < curtab; u++) {
            if (c1 < 0) {
              a0 = -a1, b0 = -b1, c0 = -c1;
            } else {
              a0 = a1, b0 = b1, c0 = c1;
            }
            licite = (a0 != poltab[u].a || b0 != poltab[u].b || c0 != poltab[u].c);
          }
          if (licite && u == curtab) {
            poltab[curtab].a = a0;
            poltab[curtab].b = b0;                                        
            poltab[curtab].c = c0;
            poltab[curtab].l = len;
            curtab++;
          }
          printf(".");
          if (len > record) {
            record = len;
          }
          fflush(stdout);
        }
      }
    }
  }
  printf(" 100.0 %%");
  fprintf(fileout,
     "# ============================================\n");
  for (t1 = 1, t2 = 0; t1 <= record; t1++) {
    for (u = 1; u < curtab; u++) {
      if (poltab[u].l == t1) {
        t2++;
        fprintf(fileout, "%3d %3d\t \"", t2, poltab[u].l);
        polyprint(fileout, poltab[u].a, poltab[u].b, poltab[u].c);
        fprintf(fileout, "\"\n");
      }
    }
  }
  fprintf(fileout,
     "# ============================================\n");
  tlast = time(NULL) - tbegin;
  fprintf(fileout,
     "# durée totale (en temps réel) %5d min %2d s\n", tlast / 60, tlast % 60);
  printf(
     "\n# durée totale (en temps réel) %5d min %2d s\n", tlast / 60, tlast % 60);
  printf(
     "# %d polynômes ont été retenus\n"
     "# les résultats sont écrits dans le fichier %s\n", t2, filename);             
  if (fileout != stderr)
    fclose(fileout);
  return 0;
}
 
/* polyprint
   mise en forme et affichage sur f du polynôme a*x^2 + b*x + c
 */
void    polyprint(FILE * f, int a, int b, int c)
{
  char    str[100], *p;
  int     v;
  p = str;
  *p = '\x0';
  v = (a < 0) ? -a : a;
  if (a < 0)
    strcat(p, "- ");
  p = p + strlen(p);
  if (v > 1)
    sprintf(p, "%d*x^2 ", v);
  else
    sprintf(p, "x^2 ");
  if (b != 0) {
    v = (b < 0) ? -b : b;
    if (b < 0)
      strcat(p, "- ");
    else
      strcat(p, "+ ");
    p = p + strlen(p);
    if (v > 1)
      sprintf(p, "%d*x ", v);
    else
      sprintf(p, "x ");
  } 
  v = (c < 0) ? -c : c;
  if (c < 0)
    strcat(p, "- ");
  else
    strcat(p, "+ ");
  p = p + strlen(p);
  sprintf(p, "%d", v);
  fprintf(f, "%s", str);
}