#include <stdio.h> /* (C) Achim Flammenkamp, University Bielefeld, 2026-01-18 */
#include <math.h>  /* extern long double expl(long double),logl(long double); */
#include <stdlib.h>/* the  0th SHCN = 1st HCN = 1  ,  1st SHCN = 2nd HCN = 2  */
#define  version "2026-01-25 19:41"/* Intel Xeon E3-1275v5@3.7GHz & gcc 15.2.1*/
/* time hcn_bene 4198890 4540250489 4198891  -->  ERROR:  MPE/N too small     */
/* time hcn_bene       0          1   10000  --> 41 sec   (2255135 many HCNs) */
/* this code can also read an upperbound file containing 1 byte for each SHCN */
/* time hcn_bene       0          1   10000  --> 17.4 sec (2255135 many HCNs) */
/* time hcn_bene 4198890 4540250489 4198891  -->  9.5 sec (1292+2 many HCNs)  */

#ifndef  uns64
#define  uns64  unsigned long
#endif
#ifndef  uns32
#define  uns32  unsigned
#endif
#ifndef  int32
#define  int32  int
#endif
#ifndef  int16
#define  int16  short
#endif
#ifndef  double80
#define  double80  long double
#endif
#ifndef  uns8
#define  uns8  unsigned char
#endif
#ifndef  bool 
#define  bool   uns8
#endif

#ifndef  N
#define  N  (1<<22)  // maximal number of used prime numbers
#endif 
#ifndef  ANZ
#define  ANZ  16  // maximal number of different prime exponents >= 0
#endif
#ifndef  MPE
#define  MPE  40  // maximal (two) power exponent
#endif
#ifndef  EPS 
#define  EPS  8e-14   // numerical absolute precision to distinguish double80s
#endif
#define  MHSS (1<<11)   // maximal number of hcns between successive shcns
#define  UPPER_BOUND  "max_x2b_upper_bound"  // default upper bound x2bene file
#ifndef  MAX_TABLE
#define  MAX_TABLE 4200000 // maximal index of SHCN which we read an upper bound
#endif
#define  MUB   3.5  // highest upper bound of max_n benefit_x(n)*x^2 ever needed

static uns32     anz;    // total number of computed prime numbers
static double80  eps, max_double_value;
static uns32     p[N];  // list of prime numbers
static double80  lnp[N];  // list of logarithmns of prime numbers
static double80  slnp[N+1]; // sum of logarithmns of prime numbers
static double80  d_incr[MPE];  // list of logarithmns of  1 + 1/(arg+1)
static double80  ln1n[MPE];  // list of logarithmns of  1+arg
static int16     sgn[N];  // list of up and downs of each prime number 
static uns32     max_two_exp, c;  // counter for used entries in hcn[]
typedef  struct { uns32 a, ev[ANZ], ea[ANZ]; double80 lv, ld, x2b; } HCN;
static HCN       hcn[MHSS];    // list of sorted by magnitude possible hcns
static FILE      *fp;   // for file to read upper bounds of maxub of shcn


uns32 next_prime()
{  // trial division by already calculated primes is fast enough for us
   static uns32  h=1;
   uns32  j, pri;
   if (anz <= 1)
      return  2 + anz;
   for (pri=p[anz-1]+2;;pri+=2)
   {  h += (p[h]*p[h] <= pri);
      for (j=1;j<h;j++)
         if (!(pri%p[j]))
            goto over;
      return  pri;
      over: ;
   }
}


void update_next_prime()
{  static double80  integer= 0.0;
   static double80  fraction= 0.0;
   double80  ln, floorln;
   if (anz == N)
   {  fprintf(stderr,"N too small\n");
      exit(6);
   }
   p[anz]= next_prime();
   ln= logl((double80)p[anz]);
   floorln= (uns32)ln;
   lnp[anz]= ln;
   sgn[anz]= 0;
   integer += floorln;
   fraction += ln - floorln;
   if (fraction >= 1.0)
   {  fraction -= 1.0;
      integer += 1.0;
   }
   slnp[++anz]= integer + fraction;
}


void list_expo_hcn(HCN next)
{  uns32  i;
   for (i=0;i<next.a;i++)
      if (next.ea[i])
      {  printf(" %u",next.ev[i]);
         if (next.ea[i] > 1)
            printf("^%u",next.ea[i]);
      }
}


void check_and_store(HCN next, double80 val)
{  uns32  h, i, j;
   for (h=0,i=c;h<i;)
   {  j= h + (i-h)/2;
      if (next.lv <= hcn[j].lv+eps && next.lv >= hcn[j].lv-eps)
      {  fprintf(stderr,"eps %Le too large for sorting\n",eps);
         fprintf(stderr,"delta=%Le\n",next.lv-hcn[j].lv);
         exit(3);
      }
      if (next.lv < hcn[j].lv)
         i= j;
      else
         h= j+1;
   };
   j= h;
   if (j && next.ld <= hcn[j-1].ld+eps)
      return;
   if (j==c || next.lv < hcn[j].lv)
   {  if (j==c && next.ld > hcn[j-1].ld+eps || j<c && next.ld+eps < hcn[j].ld)
      {  if (c == MHSS)
         {  fprintf(stderr,"MHSS too small\n");
            exit(7);  
         }
         for (h=j,j=c;j>h;j--)
            hcn[j]= hcn[j-1];
         hcn[h]= next;
         c++;
      }
      else if (j < c)
      {  hcn[j++]= next;
         for (h=j;h<c;h++)
            if (hcn[h].ld > next.ld+eps)
               break;
         if (h != j)
         {  for (;h<c;h++,j++)
               hcn[j]= hcn[h];
            c= j;
         }
      }
   }
   return;
}


uns32 exp_down(HCN *next, uns32 i)
{  uns32  j, k;
   next->ea[i]--;
   if (next->ev[i+1] == next->ev[i]-1)
   {  next->ea[i+1]++;
      k=1;
      if (!next->ea[i])
      {  for (j=i;j<next->a;j++)
         {  next->ev[j]= next->ev[j+1];
            next->ea[j]= next->ea[j+1];
         }
         next->a--; 
         k=3;
      }
   }
   else
   {  k=0;
      if (next->ea[i])
      {  next->a++;
         if (next->a == ANZ)
         {  fprintf(stderr,"ANZ too small\n");
            exit(4);
         }
         for (j= next->a;j>i;j--)
         {  next->ev[j]= next->ev[j-1];
            next->ea[j]= next->ea[j-1];
         }
         k=2;
         i++;  // tricky split
      }
      next->ev[i]--;
      next->ea[i]= 1;
   }
   return k;
}


uns32 exp_up(HCN *next, uns32 i)
{  uns32  j, k;
   next->ea[i]--;
   if (i && next->ev[i-1] == next->ev[i]+1)
   {  next->ea[i-1]++;
      k=1;
      if (!next->ea[i])
      {  for (j=i;j<next->a;j++)
         {  next->ev[j]= next->ev[j+1];
            next->ea[j]= next->ea[j+1];
         }
         next->a--; 
         k=3;
      }
   }
   else
   {  k=0;
      if (next->ea[i])
      {  next->a++;
         if (next->a == ANZ)
         {  fprintf(stderr,"ANZ too small\n");
            exit(4);
         }
         for (j=next->a;j>i;j--)
         {  next->ev[j]= next->ev[j-1];
            next->ea[j]= next->ea[j-1];
         }
         k=2;
      }
      next->ev[i]++;
      next->ea[i]= 1;
   }
   if (!i && next->ev[0] > max_two_exp)
   {  if (max_two_exp == MPE-1)
      {  fprintf(stderr,"MPE too small\n");
         exit(5);
      }
      d_incr[next->ev[0]]= logl(1.0+1.0/(double80)(next->ev[0]+1));
      ln1n[next->ev[0]]= logl((double80)(next->ev[0]+1));
      max_two_exp++;
   }
   return k;
}


void hcn_updown(HCN, uns32, uns32, uns32, uns32, bool,
                double80, double80, double80);

void hcn_down(HCN next, uns32 n0, uns32 h0, uns32 n1, uns32 h1, bool flag,
              double80 val, double80 x, double80 maxub)
{  uns32  i, j, k;

   for (i=h1;(int)i>=0;i--)
   { 
      if (sgn[n1] <= 0)
      {  double80  nval= next.lv, nd= next.ld, nb= next.x2b;
         sgn[n1]--;
         next.lv -= lnp[n1];
         next.ld -= d_incr[next.ev[i]-1];
         next.x2b += x*(x*d_incr[next.ev[i]-1]-lnp[n1]);
         if (next.x2b < maxub)
         {  k= exp_down(&next,i);
            if (next.lv >= hcn[0].lv-eps && next.lv < val-eps)
               check_and_store(next,val);
            if (i+(k==2) < next.a && n1 >= (k&1) && i >= (k==3))
            {  uns32  i1=i-(k==3)+(k==2), n1k= n1-(k&1);
               if (k==2 && h0 > i)
                  hcn_updown(next,n0,h0+1,n1k,i1,flag,val,x,maxub);
               else if (k==3 && h0 > i+1)
                  hcn_updown(next,n0,h0-1,n1k,i1,flag,val,x,maxub);
               else if ((k&1) && h0 == i+1)
                  hcn_updown(next,n0-1+next.ea[h0-(k>>1)],h0-(k>>1)+1,n1k,i1,
                             flag,val,x,maxub);
               else
                  hcn_updown(next,n0,h0,n1k,i1,flag,val,x,maxub);
            }
            exp_up(&next,i+((k^k>>1)&1));
         }
         next.x2b= nb;
         next.ld= nd;
         next.lv= nval;
         sgn[n1]++;
      }
      n1 -= next.ea[i];
   }
}


void hcn_up(HCN next, uns32 n0, uns32 h0, uns32 n1, uns32 h1, bool flag,
            double80 val, double80 x, double80 maxub)
{  uns32  i, j, k;
   for (i=h0;i<=next.a;i++)
   { 
      if (sgn[n0] >= 0)
      {  double80  nval= next.lv, nd= next.ld, nb= next.x2b;
         sgn[n0]++;
         next.lv += lnp[n0];
         if (next.lv > max_double_value)
         {  fprintf(stderr,"double80 precision insufficient for ln(val)\n");
            exit(3);
         }
         next.ld += d_incr[next.ev[i]];
         if (next.ld > max_double_value)
         {  fprintf(stderr,"double80 precision insufficient for ln(d)\n");
            exit(3);
         }
         next.x2b -= x*(x*d_incr[next.ev[i]]-lnp[n0]);
         if (next.x2b < maxub)
         {  k= exp_up(&next,i);
            if (i == next.a-(k==2))
               if (n0+1 == anz)
                  update_next_prime();
            if (next.lv >= hcn[0].lv-eps && next.lv < val-eps)
               check_and_store(next,val);
            if (k==2 && (int)h1 >= (int)i)
               hcn_updown(next,n0+(k&1),i,n1,h1+1,flag,val,x,maxub);
            else if (k==3 && (int)h1 >= (int)i)
               hcn_updown(next,n0+(k&1),i,n1,h1-1,flag,val,x,maxub);
            else if ((k&1) && h1==i-1)
               hcn_updown(next,n0+(k&1),i,n1+1-next.ea[h1],h1-1,flag,val,x,maxub);
            else
               hcn_updown(next,n0+(k&1),i,n1,h1,flag,val,x,maxub);
            exp_down(&next,i-(k&1));
         }
         next.x2b= nb;
         next.ld= nd;
         next.lv= nval;
         sgn[n0]--;
      }
      n0 += next.ea[i];
   }
   n0 -= next.ea[next.a];
}


void hcn_updown(HCN next, uns32 n0, uns32 h0, uns32 n1, uns32 h1, bool flag,
                double80 val, double80 x, double80 maxub)
{  
   if (next.lv >= val-eps)
      hcn_down(next,n0,h0,n1,h1,flag,val,x,maxub);
   else if (next.lv > hcn[0].lv+eps || !flag)
   {  if (!flag)
      {  hcn_up(next,n0,h0,n1,h1,flag,val,x,maxub);
         flag= 1;
      }
      hcn_down(next,n0,h0,n1,h1,flag,val,x,maxub);
   }
}


void hcn_preample(uns64 k, HCN next, uns32 n, uns8 z)
{
   printf("%10luth hcn%c lnv=%21.12Lf  lnd=%20.12Lf  ",k,z,next.lv,next.ld);
   printf("ndpe=%2u  ndp=%7u",next.a,n);
}


void print_shcn(uns32 i, uns64 k, HCN next, uns32 nn)
{  // i is the index of the SHCN
   hcn_preample(k,next,nn,'S');
   list_expo_hcn(next);
   printf("\n");
}  


void print_hcn(uns64 k, HCN next)
{  uns32  i, n;
   for (n=i=0;i<next.a;i++)
      n += next.ea[i];
   hcn_preample(k,next,n,' ');
   list_expo_hcn(next);
   printf("\n");
}


void recompute_lv_ld(HCN *next)  // to keep rounding errors of lv and ld low
{  uns32  i, m;
   next->lv= 0.0;
   m= 0;
   next->ld= 0.0;
   for (i=0;i<next->a;i++)
   {  m += next->ea[i];
      next->lv += (next->ev[i]-next->ev[i+1])*slnp[m];
      next->ld += next->ea[i]*ln1n[next->ev[i]];
   }
}


uns32 shcn(uns32 first, uns64 offset, uns32 last)
{
   uns32  h, i, j, k, n, m, kanz;
   uns64  total;
   double80  val, d, maxub, x, dx, xprev, xnew, maxx2b;
   uns8  *maxub_table;

   setlinebuf(stdout);
   fprintf(stderr,"size of (double80) = %lu bytes and ",sizeof(double80));
   for (x=2.0,dx=0.5,j=0;x>1.0;x-=dx,j++,dx*=0.5);
   fprintf(stderr,"maximal mantissa precision %u bits\n",j-1);
   max_double_value= eps*0.25/dx;
   anz= 0;
   slnp[0]= 0.0;
   update_next_prime();
   max_two_exp= 0;
   ln1n[0]= 0.0;
   d_incr[0]= logl(1.0+1.0/1.0);
   x= 0.0;
   val= 0.0;
   d= 0.0;
   hcn[0].a= 0;
   hcn[0].ev[0]= 0;
   hcn[0].ea[0]= ~0U - 1;
   hcn[0].lv= val;
   hcn[0].ld= d;
   hcn[0].x2b= 0.0;
   n= m= 0;
   total= offset;  // we count normally the shcns form 0, but the hcns from 1
   kanz= 0; // number of SHCNs (starting from 0) we know an upper bound of maxx2b
   maxub_table= (uns8*) 0;
   if (fp)
   {  if (!(maxub_table= malloc(MAX_TABLE)))
         return  fprintf(stderr,"insufficient memory due to MAX_TABLE\n"), 8;
      kanz= fread(maxub_table,1,MAX_TABLE,fp);
      fprintf(stderr,"%u max_x2bene upper bound values read\n",kanz);
      fclose(fp);
   }
   else
      fprintf(stderr,"no upper bounds of max_x2b values accessable\n");
   for (k=0;k<last /* && k<kanz */ ;k++)    // loop for shcns
   {
      if (first <= k)
      {  recompute_lv_ld(&hcn[0]);
         print_shcn(k, total, hcn[0], n);
      }
      xprev=x;
      i= hcn[0].a;
      m= n;
      x= lnp[n]/d_incr[hcn[0].ev[i]];
      for (n=j=0;j<hcn[0].a;j++)
      {  xnew= lnp[n]/d_incr[hcn[0].ev[j]];
         if (xnew < x)
         {  x= xnew;
            i= j;
            m= n;
         }
         n += hcn[0].ea[j];
      }
      val += lnp[m];     // ln( next shcn )
      d += d_incr[hcn[0].ev[i]];  // ln( d(next shcn) )
      // here we decode the upper bound of max_x2b for the next SHCN
      // it must always be > max_{n\in{HCNs}} benefit_eps(n)/eps^2
      maxub= (k < kanz ? (double)maxub_table[k]/256.0 : 1.0)*MUB;
      if (first <= k)
      {  double80  oldlv, oldld;
         c= 1;
         oldlv= hcn[0].lv;
         hcn[0].lv= 0.0;
         oldld= hcn[0].ld;
         hcn[0].ld= 0.0;
// compute and collect all hcns from current up to next SHCN in hcn[c++]
         hcn_updown(hcn[0],0,0,n-1,hcn[0].a-1,0,lnp[m],xprev,maxub);
         hcn[0].lv= oldlv;
         hcn[0].ld= oldld;
         maxx2b= 0.0;
         for (j=1;j<c;j++)
         {  if (hcn[j].x2b > maxx2b)
               maxx2b= hcn[j].x2b;
            hcn[j].lv += oldlv;
            hcn[j].ld += oldld;
            recompute_lv_ld(&hcn[j]);
            print_hcn(total+j,hcn[j]);
         }
         total += c;
      }
      hcn[0].lv= val;
      hcn[0].ld= d;
      j= exp_up(&hcn[0], i);
      if (i == hcn[0].a-(j==2))
         if (++n == anz)
            update_next_prime();
   }
   recompute_lv_ld(&hcn[0]);
   print_shcn(k, total, hcn[0], n);
   if (maxub_table)
      free(maxub_table);
return 0;
}


int32 main(uns32 argc, char *argv[])
{
   int32  first, last, i;
   uns64  offset;
   char  *upper_bound_file;
   fprintf(stderr,"%s  N=%u  ANZ=%u  MPE=%u  MHSS=%u  EPS=%.1e\n",version,
           N,ANZ,MPE,MHSS,EPS);
   if (argc <= 1)
   {  fprintf(stderr,"usage: %s  index_of_shcn [hcn_offset [last_shcn "
                     "[upper_bound_file [eps]]]]\n",argv[0]);
      return 1;
   }
   first= 0;
   if (argc > 1 && (1 != sscanf(argv[1],"%d",&first) || first < 0))
   {  fprintf(stderr,"index_of_shcn must be >= 0\n");
      return 2;
   }
   offset= 1;
   if (argc > 2 && (1 != sscanf(argv[2],"%ld",&offset) || offset <= 0))
   {  fprintf(stderr,"hcn_offset must be > 0\n");
      return 2;
   }
   last= ~0;
   if (argc > 3 && (1 != sscanf(argv[3],"%d",&last) || last <= first))
   {  fprintf(stderr,"last_shcn must be > index_of_shcn\n");
      return 2;
   }
   upper_bound_file=  argc > 4 ? argv[4] : UPPER_BOUND ;
   if (!(fp=fopen(upper_bound_file,"r")) )
   {  if (argc > 4)
      {  fprintf(stderr,"Error: can not access %s to read upper bounds\n",
                 upper_bound_file);
         return 2;
      }
      else
         fprintf(stderr,"Warning: can not access %s to read upper bounds\n",
                 upper_bound_file);
   }
   eps= EPS;
   if (argc > 5 && (1 != sscanf(argv[5],"%Lf",&eps) || eps <= 0.0 ||
                    eps > 0.69314718))
   {  fprintf(stderr,"eps must be > 0 and < ln(2)\n");
      return 2;
   }
   fprintf(stderr,"first_shcn=%u  hcn_offset=%lu  last_shcn=%u  "
           "eps=%Le\n",first,offset,last,eps);
   return  shcn(first,offset,last);
}