//
//////////////////////////////////////////////
// d8.c                      Dedekind numbers.
// Implementation of Doug H. Wiedemann Algorithm.
// C, OpenMP. On 8 11g cores: 15'. 800 lines.
// MIT Licensed
//
// J. M. Aranda. Madrid. 2021. 
//////////////////////////////
// References
// - A computation of the eighth Dedekind number, 
// Doug Wiedemann , Order,  volume 8, 1991.
// - Enumerating Free Distributive Lettices, 
// George Markowsky, 1989 , Report No. 89-10.
// - Cardinalities of finite distributive lattices, 1976. 
// J. Berman, P. Koehler.
/////////////////////////////////////////////////////////

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>
#include <string.H>
#include <assert.h>
#include <time.h>
#include <omp.h> 

/****** MAIN PARAMETERS ************/
#define SETCMAX  8000000
#define RSETMAX  18000
#define EXPONMAX 6
#define EFACTMAX 720
#define EXPOW2MAX 64
#define QFLMAX 65
#define CONST64 64
#define PWDIM 72
#define RULEDIM 70
#define STRINGMAX 32
//////////////////////////
typedef uint_fast16_t U16;
typedef int_fast32_t  I32;
typedef uint_fast32_t U32;
typedef uint_fast64_t U64;
typedef unsigned __int128 VAL128;
typedef uint_fast64_t VAL64;
typedef uint_fast64_t ITEM64;
//
#pragma pack(16)
static I32 CSetCard=0;
static I32 EXPONFACT,RSETCARD=0;
static I32 QORBDUA=0;
static I32 QBASE=0;
static U16 globit=0,EXPONGLOBAL=0,EXPOW;
static bool ESPOT=false;
#define DIMPAR 8
static int NUMTHR=DIMPAR;
/************************************************************************ OMP ****/
static VAL64 POWERS[PWDIM]__attribute__((aligned));
static VAL64 APORTA[SETCMAX];
static ITEM64 CSET[SETCMAX]__attribute__((aligned));
static bool ISDUPLI[SETCMAX]__attribute__((aligned));
static I32 RSET[RSETMAX]__attribute__((aligned));
static I32 RBASE[RSETMAX];
//
static I32 INDEXQ[QFLMAX]__attribute__((aligned));
static I32 INDEXF[QFLMAX];
static I32 INDEXL[QFLMAX];
static I32 CENTER[SETCMAX]__attribute__((aligned));
static I32 ORBITSZ[SETCMAX];
static I32 VALDUAL[SETCMAX]__attribute__((aligned));
static I32 ORBITDUAL[SETCMAX];
static I32 LESSCOUNT[SETCMAX]__attribute__((aligned));
static I32 GREATCOUNT[SETCMAX];

/************* OMP ************/
static ITEM64 WTOFIND[DIMPAR][SETCMAX];
static I32 ISWHERE[DIMPAR][SETCMAX];
static VAL64 JKPRODUCT[DIMPAR][SETCMAX];
/************************************************************** OMP ************/
static U64 RULENUMERIC[RULEDIM];
static char RULESTEXT[RULEDIM][RULEDIM];
static U16 PERMsmall[EFACTMAX][EXPONMAX];
static U16 PERMbig[EFACTMAX][EXPOW2MAX];
static char COMONBUFFER[STRINGMAX];
/**************************************/

static inline void clearLine(void)
__attribute__((cold));
static inline void clearLine(void)
{
printf("\r                           ");
printf("\r");
}

static inline U16 istrlen(const char *p)
__attribute__((hot));
static inline U16 istrlen(const char *p)
{
U16 r=strlen(p);
return r;
}

static  void BUFROM128(const VAL128 p)
__attribute__((cold));
static  void BUFROM128(const VAL128 p)
{
const char CHZ='0';VAL128 n=p;
for(int i=0;i<STRINGMAX;i++)COMONBUFFER[i]=' ';
int k=(STRINGMAX-1);
while(n != 0)
  {
  assert(k >= 1);
  int m=n%10;
  //if(k%5 ==0)COMONBUFFER[k--]='.';
  COMONBUFFER[k--]=(CHZ+m);
  n=(n-m);n=n/10;
  }
assert(k >=0);
for(int i=0;i<k;i++)COMONBUFFER[i]=' ';
k=0;
for(int i=0;i<STRINGMAX;i++)
  if(COMONBUFFER[i]!=' ')
    COMONBUFFER[k++]=COMONBUFFER[i];
COMONBUFFER[k]=0;
}

/************************************************************** OMP  ********/
static void configOMP(void)__attribute__((cold));
static void configOMP(void)
{
int nuc;
nuc=omp_get_num_procs();
if(nuc  >=  DIMPAR)NUMTHR=DIMPAR;
else NUMTHR=nuc;
omp_set_num_threads(NUMTHR);
printf("Detected Cores: %4d \n",nuc);
/***********************************/
omp_set_dynamic(0);omp_set_nested(0);
omp_set_max_active_levels( 1 );
omp_set_schedule(omp_sched_static,200);
}

static void genrules(const int pn)
__attribute__((cold));
static void genrules(const int pn)
{
const int dkn=pn;char l[dkn+4];int a,b,k;
assert(dkn  <=  RULEDIM);
a=3;b=3;
for(int r=0;r< dkn;r++)
  {
  if(r  ==  0){a=3;b=3;}
  for(int i=0;i<dkn;i++)l[i]='n';
  l[dkn]='f';l[dkn+1]=0;k=dkn-1-r;
  for(int i=a;i<=b;i++)
    {
    if(a  ==  (a&i))l[k]='1';else l[k]='x';
    k++;
    }
    l[dkn-1-r]='e';
    l[dkn+1]=0;
    strcpy(RULESTEXT[r],l);
    RULESTEXT[0][dkn-1]='e';
    {
    ITEM64 rvn;char ch;
    int cox,uder;

    strcpy(l,RULESTEXT[r]);uder=0;
    for(int b=dkn-1;
      (l[b]  ==  '1') && (b>=0);b--)uder++;
    rvn=0;cox=0;
    for(int b=dkn;b>=0;b--)
      {
      assert(b  <=  RULEDIM);ch=l[b];
      if(ch  ==  '1')rvn++;
      else if(ch  ==  'x')cox++;
      else if(ch  ==  'e'){rvn++;break;}
      rvn<<=1;
      }
    RULENUMERIC[r]=rvn;
    }
  a--;if(a<0){a=b;b=2*b+1;}
  }
}

static bool createset(const int pn)
__attribute__((cold));
static bool createset(const int pn)
{
ITEM64 tmp,rnm,pwm,mit;
I32 grovI,kC;
const int nk=pn;

assert(nk >= 1);
CSET[0]=0;CSET[1]=1;grovI=2;kC=grovI;
globit=1;
for(int k=globit;k<nk;k++)
  {
  globit++;pwm=POWERS[k];rnm=RULENUMERIC[k];
  for(int i=1;i<grovI;i++)
    {
    tmp=CSET[i];tmp=tmp|pwm;mit=tmp&rnm;
    if(mit != rnm)continue;
    if(kC >= SETCMAX)return false;
    CSET[kC++]=tmp;
    }
  grovI=kC;
  }
CSetCard=grovI;
for(int i=0;i <= (CSetCard-2);i++)
  assert(CSET[i] < CSET[i+1]);
/* Checked set m.ascending */
return true;
}

/******** hot functions **************/
static inline bool nbiti(
const ITEM64 n,const U16 globit)__attribute__((hot));
static inline bool nbiti(
const ITEM64 n,const U16 globit)
{
bool v;ITEM64 w;
v=false;w=n&POWERS[globit];if(w>0)v=true;
return v;
}

static int factorial(const int e)
__attribute__((cold));
static int factorial(const int e)
{
int f=1;
for(int i=2;i<=e;i++)f*=i;
return f;
}

/* modo 64b */
#define TSUXEN 4
typedef struct
  {U16 E[TSUXEN];
  } TSUXE;
typedef union
  {ITEM64 L;TSUXE P;  
  } TU64;
/*** mode ITEM64 64bit ***/
static inline bool ITEMLEQITEM(const I32 pa,const I32 pb)
__attribute__((hot)) __attribute__((pure));
static inline bool ITEMLEQITEM(const I32 pa,const I32 pb)
{
TU64 ua,ub,tmp;
/* 64b mode leq binary */
ua.L=CSET[pa];ub.L=CSET[pb];
for(int i=0;i < TSUXEN;i++)
  {
  tmp.P.E[i]=(ua.P.E[i] & ub.P.E[i]);
  if(tmp.P.E[i] != ua.P.E[i])return false;
  }
return true;
}

static void fillQless(
const I32 ini,const I32 las)
__attribute__((hot));
static void fillQless(
const I32 ini,const I32 las)
{
for(I32 ii=ini;ii<las;ii++)
  {
  const I32 igk=RSET[ii];I32 cc=0;
  cc=0;
  for(I32 k=0;k  <=  igk;k++)if(ITEMLEQITEM(k,igk))cc++;
  LESSCOUNT[igk]=cc;
  if(ii%1000  ==0)printf("\rL %09d",ii);
  }
clearLine();
}

static void fillQgreat(
const I32 ini,const I32 las)
__attribute__((hot));
static void fillQgreat(
const I32 ini,const I32 las)
{
const I32 CSetCardk=CSetCard;
for(I32 ii=ini;ii<las;ii++)
  {
  const I32 igk=RSET[ii];I32 cc=0;
  cc=0;
  for(I32 k=igk;k  <  CSetCardk;k++)if(ITEMLEQITEM(igk,k))cc++;
  GREATCOUNT[igk]=cc;
  if((ii%1000)  ==  0)printf("\rG %09d",ii);
  }
clearLine();
}

static void fillQQ(void)__attribute__((hot));
static void fillQQ(void)
{
const I32 CSetCardk=CSetCard;
const I32 RVALk=RSETCARD;
I32 q1,q2,q3;
/* Q precalculated.  ********************/
for(I32 i=0;i<SETCMAX;i++)
  {LESSCOUNT[i]=0;GREATCOUNT[i]=0;}
q1=RVALk/4;q2=2*q1;q3=3*q1;
/***************************** OMP OMP OMP *****************/
#pragma omp parallel sections
{
#pragma omp section
  fillQless(0,q1);
#pragma omp section
  fillQless(q1,q2);
#pragma omp section
  fillQless(q2,q3);
#pragma omp section
  fillQless(q3,RVALk);
#pragma omp section
  fillQgreat(0,q1);
#pragma omp section
  fillQgreat(q1,q2);
#pragma omp section
  fillQgreat(q2,q3);
#pragma omp section
  fillQgreat(q3,RVALk);
}

/* Proved checked */
for(I32 i=0;i<CSetCardk;i++)
  {
  if(ISDUPLI[i])
    {
    LESSCOUNT[i]=LESSCOUNT[CENTER[i]];
    GREATCOUNT[i]=GREATCOUNT[CENTER[i]];
    }
  else
    {
    assert(LESSCOUNT[i] >=1);
    assert(GREATCOUNT[i] >=1);
    }
  }
}

static inline I32 binsearchCset(
const int capoin,const bool mbe,const I32 ini,
const I32 fin,const U64 WTOFIND)
__attribute__((hot))__attribute__((pure));
static inline I32 binsearchCset(
const int capoin,const bool mbe,const I32 ini,
const I32 fin,const ITEM64 WTOFIND)
{
const ITEM64 wtfk=WTOFIND;const int k64=CONST64;
I32 fi,la,mi,di;int lc;
/* binary search */
if(wtfk  ==  0)lc=k64;
else lc=__builtin_clzll(wtfk);
fi=INDEXF[lc];if(ini  >  fi)fi=ini;
la=INDEXL[lc];if(fin  <  la)la=fin;
di=la-fi;
while(di  >=  0) 
  {
  ITEM64 Cmi;
  di>>=1;mi=fi+di;Cmi=CSET[mi];
  if((wtfk < Cmi))la=mi-1;
  else if((wtfk > Cmi))fi=mi+1;    
  else return mi;
  di=(la-fi);
  }
if(mbe)assert(mbe == false);
if(capoin);
return -1;
}

static inline VAL64 sumaporta(const I32 ii)
__attribute__((hot));
static inline VAL64 sumaporta(const I32 ii)
{
const I32 CSetCardk=CSetCard;const I32 CSetCardk1=CSetCard-1;
const int ntk=omp_get_thread_num();
const I32 iik=ii;
const ITEM64 Cik=CSET[ii];
const ITEM64 Cik4[4]={Cik,Cik,Cik,Cik};
VAL64 sapo;const I32 L4=(CSetCard-(CSetCard%4));

for(I32 h=0;h<CSetCardk;h++)ISWHERE[ntk][h]=-1;
for(I32 h=0;h<L4;h+=4)
  {
  const I32 hk=h;
// ******** simd omp *******
#pragma omp simd simdlen(4)
  for(I32 i=0;i< 4;i++)
    {WTOFIND[ntk][hk+i]=(Cik4[i]&CSET[hk+i]);}
  }
for(I32 h=L4;h<CSetCardk;h++)
  {WTOFIND[ntk][h]=Cik&CSET[h];
  }
for(I32 h=0;h<CSetCardk;h++)
  {
  I32 pista;
  pista=iik;if(h  <  pista)pista=h;
  ISWHERE[ntk][h]
  =binsearchCset(1,true,0,pista,WTOFIND[ntk][h]);
  }
for(I32 h=0;h<CSetCardk;h++)
  JKPRODUCT[ntk][h]= LESSCOUNT[ISWHERE[ntk][h]];
for(I32 h=0;h<CSetCardk;h++)ISWHERE[ntk][h]=-1;
for(I32 h=0;h<L4;h+=4)
  {
  const I32 hk=h;
// ******** simd omp ********
#pragma omp simd simdlen(4)
  for(I32 i=0;i< 4;i++)
    {WTOFIND[ntk][hk+i]=(Cik4[i]|CSET[hk+i]);}
  }
for(I32 h=L4;h<CSetCardk;h++)
  {WTOFIND[ntk][h]=Cik|CSET[h];
  }
for(I32 h=0;h<CSetCardk;h++)
  {
  I32 pista;
  pista=iik;if(h  >  pista)pista=h;
  ISWHERE[ntk][h]=
  binsearchCset(2,true,pista,CSetCardk1,WTOFIND[ntk][h]);
  }
for(I32 h=0;h<CSetCardk;h++)
  {JKPRODUCT[ntk][h]*=GREATCOUNT[ISWHERE[ntk][h]];}
sapo=0;
for(I32 h=0;h<CSetCardk;h++)sapo+=JKPRODUCT[ntk][h];
assert(sapo  >  0);
return sapo;
}

static VAL128 setMultiplier(void)__attribute__((hot));
static VAL128 setMultiplier(void)
{
const I32 CSetCardk=CSetCard;
const I32 RVALk=RSETCARD;
VAL128 apo8;I32 lcrum;
for(I32 i=0;i<SETCMAX;i++)APORTA[i]=0;
for(I32 i=0;i<RSETMAX;i++)RBASE[i]=0;
lcrum=0;
for(I32 j=0;j<RVALk;j++)
  {
  I32 ii=RSET[j];
  if(ORBITDUAL[ii]  <  0)RBASE[lcrum++]=ii;
  }
QBASE=lcrum;
//printf("DBG:  x4 ESPOT:%d QBASE= %d \n",ESPOT,QBASE);
if(!ESPOT)assert(QBASE == RSETCARD);
/* The HOT code */
  {
  const I32 QBASEk=QBASE;

// ******* main point of OMP ********
#pragma omp parallel for schedule(static,200) 
  for(I32 j=0;j<QBASEk;j++)
    {
    const I32 jk=j;const I32 iik=RBASE[jk];
    APORTA[iik]=sumaporta(iik);
#pragma omp critical
    if((jk%100)==0)printf("\rx4 %06d/%06d",jk,QBASEk);
    }
  }
clearLine();

if(!ESPOT)assert(QORBDUA  ==  0);
// ***** omp ******
#pragma omp parallel for
for(I32 j=0;j<RVALk;j++)
  {
  I32 ii=RSET[j];
  if(ORBITDUAL[ii]  >=  0)
    APORTA[ii]=APORTA[ORBITDUAL[ii]];
  }
if(ESPOT)
  {
  for(int i=0;i<CSetCardk;i++)
    {
    const I32 ik=i;const VAL64 APOk=APORTA[i];
    if(APOk  ==  0)continue;
/*************** omp *********/
#pragma omp parallel for 
    for(int j=ik+1;j<CSetCardk;j++)
      {
      if(APORTA[j]  !=  APOk)continue;
      if(ORBITDUAL[j]  ==  ik)continue;
      printf("INTERNAL-ERROR APOeq: %d=%d %I64d \n"
        ,ik,j,APOk);
      exit(0);
      }
    printf("\r %d ",i);
    }
  clearLine();
  }

apo8=0;
for(I32 j=0;j<RVALk;j++)
  {
  VAL128 poa;I32 ig=RSET[j];
  assert(APORTA[ig]  >  0);
  poa=(VAL128)(ORBITSZ[ig]*APORTA[ig]);
  apo8+=poa;
  }
return apo8;
}

static void fillINDEXQFL(void)__attribute__((hot));
static void fillINDEXQFL(void)
{
const I32 CSetCardk=CSetCard;
const I32 CSetCardk1=CSetCard-1;
const int k64=CONST64;
for(int i=0;i  <  QFLMAX;i++)
  {INDEXQ[i]=0;INDEXF[i]=0;INDEXL[i]=0;
  }
for(I32 i=0;i  <  CSetCardk;i++)
  {
  ITEM64 Ci=CSET[i];
  int lc;
  if(Ci  ==  0)lc=k64;
  else lc=__builtin_clzll(CSET[i]);
  assert(lc  >=  0);assert(lc  <=  k64);
  if(INDEXQ[lc]  ==  0)INDEXF[lc]=i;
  INDEXQ[lc]++;
  }
for(int i=0;i<QFLMAX;i++){INDEXL[i]=INDEXF[i]+INDEXQ[i]-1;}
  {
  I32 s=0;
  for(int i=0;i<QFLMAX;i++)s+=INDEXQ[i];
  assert(s  ==  CSetCardk);
  /* Check index integrity */
  s=globit;
  for(I32 i=0;i<CSetCardk;i+=s)
    assert(binsearchCset(0,false,0,CSetCardk1,CSET[i])  ==  i);
  }
}

static void fillVALDUAL(void)__attribute__((hot));
static void fillVALDUAL(void)
{
const int CSetCardk=CSetCard;
const int CSetCardk1=CSetCard-1;
const int biTk=globit;
const int biTk1=globit-1;
bool ba[CONST64];

for(int i=0;i<CSetCardk;i++)VALDUAL[i]=-1;
if(!ESPOT)return;
for(int i=0;i<RSETCARD;i++)
  {
  ITEM64 w;I32 j=RSET[i];if(ISDUPLI[j])continue;
  w=CSET[j];
  for(int b=0;b<CONST64;b++)ba[b]=0;
  for(int b=0;b<biTk;b++)
    {
    ITEM64 f;bool e;
    e=false;f=w&POWERS[b];if(f>0)e=true;
    ba[biTk1-b]=!e;
    }
  w=0;
  for(int b=0;b<biTk;b++)w+=ba[b]*POWERS[b];
  VALDUAL[j]=binsearchCset(0,false,0,CSetCardk1,w);
  printf("\rVALDUAL %9d",j);
  }
clearLine();
}

static void fillORBITDUAL(void)
__attribute__((hot));
static void fillORBITDUAL(void)
{
const I32 CSetCardk=CSetCard;
QORBDUA=0;
for(I32 i=0;i<CSetCardk;i++)ORBITDUAL[i]=-1;
if(!ESPOT)return;
for(I32 i=0;i<CSetCardk;i++)
  {
  I32 di;
  if(ISDUPLI[i])continue;
  if(ORBITDUAL[i]  !=  -1)continue;
  di=VALDUAL[i];if(di  ==  i)continue;
  if(CENTER[di]  ==  i)continue;
  if(ORBITSZ[di]  !=  ORBITSZ[i])continue;
  if(ORBITDUAL[CENTER[di]]  !=  -1)continue;
  ORBITDUAL[CENTER[di]]=i;QORBDUA++;
  }
}

static bool fillRSET(const int pn)
__attribute__((cold));
static bool fillRSET(const int pn)
{
I32  GLBPPX[EFACTMAX];
const I32 CSetCardk=CSetCard;
const I32 CSetCardk1=CSetCard-1;
const int pnk=pn;
const int PEfk=EXPONFACT;
const int PBXk=EXPOW;
int impa,npermut=0;
bool BITPBX[EXPOW2MAX];
bool GLBNUL[EFACTMAX];
bool PEVALE[EFACTMAX];

impa=pn;ESPOT=false;
while(impa>1 && (impa%2)  ==  0)impa/=2;
if(impa  ==  1)ESPOT=true;

for(I32 i=0;i<SETCMAX;i++)
  {ISDUPLI[i]=false;ORBITSZ[i]=-1;CENTER[i]=-1;
  }
for(U16 i=0;i<PEfk;i++)PEVALE[i]=true;
if(!ESPOT)
 for(U16 p=0;p<PEfk;p++)
  {
  bool vale=true;
  for(U16 b=0;b<pn;b++)
    {U16 d;d=PERMbig[p][b];
    if(d  >=  pnk)vale=false;
    }
  PEVALE[p]=vale;
  }
npermut=0;
for(U16 p=0;p<PEfk;p++)if(PEVALE[p])npermut++;
//printf("DBG: npermut %d \n",npermut);
assert(npermut  >  0);
/* Big for fillRSET */
for(I32 iii=0;iii<CSetCardk;iii++)
  {
  I32 sze;ITEM64 Ciii;
  if(ISDUPLI[iii])continue;
  CENTER[iii]=iii;Ciii=CSET[iii];
  for(U16 i=0;i<PEfk;i++)
    {GLBPPX[i]=-1;GLBNUL[i]=false;
    }
/* for fill vectrs */
  for(U16 p=0;p<PEfk;p++)
    {
    ITEM64 Cp=0;

    if(!ESPOT && !PEVALE[p])
      {
      GLBNUL[p]=true;GLBPPX[p]=-1;
      continue;
      }

    for(U16 b=0;b<PBXk;b++)
      {U16 d;
      d=PERMbig[p][b];
      BITPBX[d]=nbiti(Ciii,b);
      }
    for(U16 b=0;b<PBXk;b++)Cp+=BITPBX[b]*POWERS[b];
    GLBPPX[p]=binsearchCset(0,false,0,CSetCardk1,Cp);
    if(GLBPPX[p]  <  0)GLBNUL[p]=true;
    } 
  sze=0;
  for(U16 p=0;p<PEfk;p++)
    {
    I32 Cp;
    if(GLBNUL[p])continue;
    Cp=GLBPPX[p];sze++;
    for(U16 q=(p+1);q<PEfk;q++)
      if(GLBPPX[q]  ==  Cp)GLBNUL[q]=true;
    }

  assert(sze  >=  1);assert(EXPONFACT%sze  ==  0);
  ORBITSZ[iii]=sze;
  for(I32 p=0;p<PEfk;p++)
    {
    I32 j;
    if(GLBNUL[p])continue;
    j=GLBPPX[p];if(j  ==  iii)continue;
    ISDUPLI[j]=true;CENTER[j]=iii;
    ORBITSZ[j]=sze;
    }
  if((iii%1000)  ==  0)printf("\rR:%9d",iii);  
  }
  clearLine();
  {
  I32 sorb,ctd;
  for(I32 i=0;i<RSETMAX;i++)RSET[i]=0;
  ctd=0;
  for(I32 i=0;i<CSetCardk;i++)
    {if(!ISDUPLI[i])
      {
      if(ctd  >=  RSETMAX)return false;  
        // RSETMAX too small.
      RSET[ctd++]=i;
      }
    }
  RSETCARD=ctd;sorb=0;
  for(I32 i=0;i<CSetCardk;i++)
    {if(!ISDUPLI[i])sorb+=ORBITSZ[i];}
  assert(sorb  ==  CSetCardk);
  }
return true;
}

static void fillPermutations(const int e)
__attribute__((cold));
static void fillPermutations(const int e)
{
U16 PPcifra[EXPONMAX],BITPBX[EXPOW2MAX];
U16 d,dig;U16 perm;
bool ok;I32 s,t;

EXPONGLOBAL=e;EXPONFACT=factorial(e);EXPOW=POWERS[e];
assert(globit  <=  EXPOW);assert(EXPOW  <=  EXPOW2MAX);
assert(EXPONGLOBAL  <=  EXPONMAX);
assert(EXPONFACT  <=  EFACTMAX);
perm=0;
for(I32 n=0;perm  <  EXPONFACT;n++)
  {
  t=n;ok=true;
  for(U16 c=0;c<EXPONGLOBAL;c++)
    {dig=t%10;
    if(dig  >=  EXPONGLOBAL){ok=false;break;}
    PPcifra[c]=dig;t/=10;
    }
  if(!ok)continue;

  for(U16 c=0;c<EXPONGLOBAL;c++)
   for(U16 f=(c+1);f<EXPONGLOBAL;f++)
     if(PPcifra[c]  ==  PPcifra[f])
       {ok=false;break;
       }
  if(!ok)continue;

  for(U16 c=0;c<EXPONGLOBAL;c++)PERMsmall[perm][c]=PPcifra[c];
  perm++;
  }
assert(perm  ==  EXPONFACT);
for(U16 p=0;p<EXPONFACT;p++)
  {
  for(U16 n=0;n<EXPOW;n++)
    {PERMbig[p][n]=0;
    for(U16 b=0;b<EXPONGLOBAL;b++)
      {d=PERMsmall[p][b];BITPBX[d]=nbiti(n,b);
      }
    s=0;
    for(U16 b=0;b<EXPONGLOBAL;b++)if(BITPBX[b])s+=POWERS[b];
    PERMbig[p][n]=s;
    }
  }
}    

int main(void)
{
configOMP();
/* just 64bits wide. 64 x4 -> 256.*/
for(int e=0;e<64;e++)
  {
/* dont touch POWERS dim */
  POWERS[e]=1;
  for(int i=0;i<e;i++)POWERS[e]*=2;;
  }
/* just 2^6 bits wide */
int e= 6;   
/* On 11gen Intel i7, Dedekind number 8,
the last 2021 known, after 15 minutes.
*************************************/

printf("A000372(0)= 2  \n");
printf("A000372(1)= 3  \n");
int dednum=1;
for(int n=1;n<=64;n*=2)
  {
  dednum++;
  assert((unsigned)n  <=  POWERS[e]);
  // printf("==== n %d, n4 %d ====\n",n,4*n);
  // printf("PGCQFL -- n %2d -------\n",n);
  fillPermutations(e);genrules(n);
  bool st=createset(n);if(!st)continue;
  fillINDEXQFL();
  // printf("RDO -- n %2d globit %2d CSetCard %8d \n"
  // ,n,globit,CSetCard);
  st=fillRSET(n);if(!st)continue; // RSETMAX small
  fillVALDUAL();
  fillORBITDUAL();

  // printf(
  // "Q -- n %d globit %d CSetCard %d RSETCARD %d DUAL %d:%d \n"
  // ,n,globit,CSetCard,RSETCARD,ESPOT,QORBDUA);

#pragma omp flush
  fillQQ();
#pragma omp flush

//printf(
//"x4 -- n %d globit %d CSetCard %d RSETCARD %d QDUA %d n4 %d \n"
//  ,n,globit,CSetCard,RSETCARD,QORBDUA,n4);

  VAL128 FF=setMultiplier();
#pragma omp flush

//  printf(
// "ESPOT %d CSetCard %d RSETCARD %d QORBDUA %d QBASE %d \n"
// ,ESPOT,CSetCard,RSETCARD,QORBDUA,QBASE);

  BUFROM128(FF);
  printf("A000372(%d)= %s \n",dednum,COMONBUFFER);
  }
return 1;
}
//
//////////////////////
//