primenums_pc_bitvec.c

/* -----------------------------------------------------------------------------
 *  Copyright (C) 2019-2021 daiteq s.r.o.                http://www.daiteq.com
 *
 *  This program is distributed WITHOUT ANY WARRANTY; without even
 *  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 *  PURPOSE.
 *
 * -----------------------------------------------------------------------------
 *  Filename    : primenums_pd.c
 *  Authors     : Martin Danek
 *  Description : Sieve of Eratosthenes - reference desktop implementation
 *  Release     :
 *  Version     : 1.0
 *  Date        : 27.4.2021
 * -----------------------------------------------------------------------------
 */
#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>


void SWAR_SET_ELM(void *swar, int idx, int val, int bitsz)
{
  *((unsigned int *)swar) = ( (*((unsigned int *)swar)) & ~(((1<<bitsz)-1)<<(bitsz*idx))) | ((val & ((1<<bitsz)-1))<<(bitsz*idx));
}

//#define SWAR_GET_ELM(tp,svar,index)
int SWAR_GET_ELM(void *swar, int idx, int bitsz)
{
  return ((*((unsigned int *)swar) >> (bitsz*idx)) & ((1<<bitsz)-1));
}


//#define SWAR_SET_ARR_ELM(tp,arry,index,value)
void SWAR_SET_ARR_ELM(void *arry, int idx, int val, int bitsz)
{
  int pack = 32/bitsz;
  int ofs = idx / pack;
  int elidx = idx - ofs*pack;
  unsigned int *pw = ((unsigned int *)arry) + ofs;
  if (pack==1)
    *pw=val;
  else
    SWAR_SET_ELM(pw, elidx, val, bitsz);
}

//#define SWAR_GET_ARR_ELM(tp,arry,index)
int SWAR_GET_ARR_ELM(void *arry, int idx, int bitsz)
{
  int pack = 32/bitsz;
  int ofs = idx / pack;
  int elidx = idx - ofs*pack;
  unsigned int *pw = ((unsigned int *)arry) + ofs;
  if (pack==1)
    return *pw;
  else
    return SWAR_GET_ELM(pw, elidx, bitsz);
}

// -----------------------------------------------------------------------------

// -----------------------------------------------------------------------------
// BIT - macros


#define BIT_SET_ELM(word, pos, val) \
  ((val)?(word | (1<<pos)):(word & (~(1<<pos))))

#define BIT_GET_ELM(word, pos) \
  (word & (1<<pos))



#define VERBOSE /* */

#define PACKING 32 /* packing factor */

#define MAXNUM 1300000
#define SQRTMAXNUM 1141 // ceil(sqrt(1300000))
// #define MAXNUM 8192 //1300000
// #define SQRTMAXNUM 92 // 1141 // ceil(sqrt(1300000))
#define MAXSLOTS ((MAXNUM)/PACKING+1)
#define MAXMASKS MAXSLOTS // (((SQRTMAXNUM-1)/PACKING+1)<<1)


#define NL printf("\n")

void i_print(unsigned i) {
  unsigned char buf[22];
  unsigned idx=20;
  buf[21]='\0';
  //  buf[20]='\n';
  if (i==0) {
    buf[idx--]='0';
  }
  else {
    while ((i>0)&&(idx>0)) {
      buf[idx--]='0'+(i%10);
      i=i/10;
    }
  }
  buf[idx]=' ';
  printf("%s",&buf[idx]);
}


void i_hexprint(unsigned i) {
  unsigned char buf[22];
  unsigned idx=20;
  buf[21]='\0';
  //  buf[20]='\n';
  if (i==0) {
    buf[idx--]='0';
  }
  else {
    while ((i>0)&&(idx>0)) {
      if (i<10)
        buf[idx--]='0'+(i%16);
      else
        buf[idx--]='A'+(i%16)-10;
      i=i/16;
    }
  }
  buf[idx]=' ';
  printf("%s",&buf[idx]);
}



void main(void) {
int i,j,k,l,m;
unsigned sz;
unsigned top;
unsigned slot,slotreg;
unsigned prime, first, last;
unsigned fin, cont;

unsigned ticks_lo, ticks_hi, insns_lo, insns_hi;
unsigned sqr, val;

unsigned numbers[MAXSLOTS];
unsigned masks[MAXMASKS];
unsigned roots[MAXNUM];



top = MAXNUM-1;
sz = sizeof(unsigned);

#ifdef VERBOSE
printf("\nprimenums: MAXNUM is ");i_print(top);NL;
printf("\nprimenums: sizeof(int) is ");i_print(sz);NL;

printf("primenums: initializing array\n"); NL;
#endif

numbers[i]=0xfffffffe;
for (i=1;i<(MAXSLOTS)+1;i++) {
  numbers[i]=0xffffffff;
}


#ifdef PRECOMP_SQRT
k=0;
roots[0]=1;
for (i=1;i<top;i++) {
  roots[i]=i*i;
}

k=0;
while (roots[k]<top)  {
  k++;
}

#else

k=SQRTMAXNUM;

#endif


printf("primenums: sqrt("); i_print(top); printf(") is "); i_print(k);NL;



#ifdef VERBOSE
printf("primenums: eliminating composites\n");
#endif
j=0;
l=0;
prime=0;
cont=1;
while(cont) {
// Get the next prime
  while(numbers[j]==0) {
    printf("Skipping j=%d\n",j);
    j++;
  }
  val=numbers[j];
//   printf("numbers[%d]=%08X\n",j,numbers[j]);
  while (!((val>>l)&1))
    l++;
  prime=PACKING*j+l+1; // the actual prime
//   printf("Current prime is %d = %d * %d + %d + 1\n",prime,PACKING,j,l);

// Compute masks, set the number of active masks
  last=0;
  first=1;
  m=prime; // <<1;
  val=0;
  slot=0;
  fin=0;
  slotreg=0;
  while ((slot<MAXSLOTS)&&(fin<2)) {
    if (prime<=PACKING) {
      if (!first)
        val|=(1<<(m-1));
      else
        first=0;
//       printf("m %d val %08X\n",m,val);
      if (val&(1<<(PACKING-1))) {
        fin++;
        if (!slotreg) {
          slotreg=slot+1;
//           printf("slotreg is %d\n",slotreg);
        }
      }
      m+=prime;
      if (m>PACKING) {
//         printf("masks[%d]=%08X\n",slot,val);
        masks[slot++]=val;
        m=m%PACKING;
        val=0;
      }
    }
    else {
      if ((m-1-last)<(PACKING)) {
        // the only non-zero word
        val=(1<<(m-last-1));
        last=PACKING-m+last;
      }
      else {
        val=0;
        last+=PACKING;
      }
      if (val&(1<<(PACKING-1))) {
        fin++;
        if (!slotreg) {
          slotreg=slot+1;
//           printf("slotreg is %d\n",slotreg);
        }
      }
//       printf("last %d m %d m-last %d\n",last,m,m-last);
//       printf("masks[%d]=%08X\n",slot,val);
      if (first)
        masks[slot++]=0;
      else
        masks[slot++]=val;
      if ((first)&&(val))
        first=0;

    }
  }

  if (slot>=MAXMASKS) {
    printf("\nMAXIMUM NUMBER OF MASKS EXCEEDED - SLOT %d PRIME %d\n",slot,prime);
    exit(1);
  }

//   printf("slot %d\n",slot);
//   printf("prime %d\n",prime);
//   for (i=0;i<slot;i++) {
//     printf("masks[%d]=%08X\n",i,masks[i]);
//   }
  if (!slotreg) slotreg=MAXSLOTS;

  for (i=0;i<slotreg;i++) {
//     printf("slot %d numbers[%d]=%08X ~masks[%d]=%08X\n",slot, i, numbers[i],i%slot, ~masks[i%slot]);
    numbers[i]&=(~masks[i%slot]);
//     printf("numbers[%d]=%08X\n",i,numbers[i]);
  }
  for (i=slotreg;i<MAXSLOTS;i++) {
//     printf("slot %d numbers[%d]=%08X ~masks[%d]=%08X\n",slot, i, numbers[i],i%slotreg+slotreg, ~masks[i%slotreg+slotreg]);
    numbers[i]&=(~masks[i%slotreg+slotreg]);
//     printf("numbers[%d]=%08X\n",i,numbers[i]);
  }

  l++;
  if (l>=PACKING) {
    l=0;
    j++;
  }
  if (prime>k) {
//     printf("last prime analysed is %d\n",prime);
    cont=0;
  }
}
numbers[0]|=1;

printf("\n\n\n");

k=0;
l=0;
for (i=0;i<MAXSLOTS;i++) {
  for (j=0;j<PACKING;j++) {
    k++;
    if ((numbers[i]&(1<<j))&&(k<MAXNUM)) {
      printf("%d\n",k);
      l++;
    }
  }
}
printf("\n");

printf("primenums: primes found: %d\n\n",l);


// #ifdef VERBOSE
// printf("primenums: printing primes\n");
// #endif
// for (i=0;i<top;i++) {
// #ifdef VERBOSE
// //   val=a[i];
// //   val = SWAR_GET_ARR_ELM(a, i, BITSZ);
//   if (SWAR_GET_ARR_ELM(a, i, BITSZ)) {
// //   if ((int)(a[i/PACKING][i%PACKING])==1) {
//       i_print(i);
//       printf(" ");NL;
//   }
// #else
// //   if (a[i]==1)
//   if (SWAR_GET_ARR_ELM(a, i, BITSZ))
// //   if ((int)a[i/PACKING][i%PACKING])
// //   if (*(a+i)==1)
//     printf("1");
//   else
//     printf("0");
//
// #endif
// }
// printf("\n");

printf("primenums phase1 done\n");




printf("\n*** END OF TEST ***\n");


}