From d5e6bacced476a9c37f9cc3d5b01eefde7694af8 Mon Sep 17 00:00:00 2001
From: ketrptr <ketrptr@cocaine.ninja>
Date: Mon, 2 Dec 2024 20:42:10 +0100
Subject: [PATCH] 	new file:   gitignore 	new file:   primes/src/sive.c

---
 gitignore         |   6 ++
 primes/Makefile   |   8 ++
 primes/src/sive.c | 189 ++++++++++++++++++++++++++++++++++++++++++++++
 ub10/a.out        | Bin 0 -> 15792 bytes
 4 files changed, 203 insertions(+)
 create mode 100644 gitignore
 create mode 100644 primes/Makefile
 create mode 100644 primes/src/sive.c
 create mode 100755 ub10/a.out

diff --git a/gitignore b/gitignore
new file mode 100644
index 0000000..86a61fb
--- /dev/null
+++ b/gitignore
@@ -0,0 +1,6 @@
+*
+!**/*.c
+!**/*.h
+# !**/**/*.h
+# !**/**/*.c
+!**/Makefile
diff --git a/primes/Makefile b/primes/Makefile
new file mode 100644
index 0000000..d8e1eb9
--- /dev/null
+++ b/primes/Makefile
@@ -0,0 +1,8 @@
+CFALGSD=-Wall -Wextra -g -ggdb 
+CFALGSF=-Wall -Wextra -Ofast -march=native -mtune=native -funroll-all-loops
+CC = gcc  
+
+fast: src/sive.c 
+	$(CC) $(CFALGSF) src/sive.c -o fastOut
+debug: src/sive.c 
+	$(CC) $(CFALGSD) src/sive.c -o debugOut
diff --git a/primes/src/sive.c b/primes/src/sive.c
new file mode 100644
index 0000000..9278af3
--- /dev/null
+++ b/primes/src/sive.c
@@ -0,0 +1,189 @@
+#include <assert.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+
+typedef uint64_t WORD;
+typedef uint8_t BOOL;
+const BOOL false = 0;
+const BOOL true = ~false;
+
+// const WORD limit = 542;
+// const WORD limit = 1e2;
+const WORD limit = 1e9;
+const WORD nr = 78498;
+
+// Sive o Atkins hit wheel
+const WORD s[] = {1, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49, 53, 59};
+const WORD wheelLen = 16;
+
+int main() {
+  BOOL *isPrime = calloc(limit, sizeof(WORD));
+  assert(isPrime != NULL);
+
+  /*
+  // Initialize the sieve with enough wheels to include limit:
+for n ← 60 × w + x where w ∈ {0,1,...,limit ÷ 60}, x ∈ s:
+    is_prime(n) ← false
+    */
+  for (WORD w = 0; w < limit / 60; w++) {
+    for (WORD i = 0; i < wheelLen; i++) {
+      WORD n = 60 * w + s[i];
+      if (n >= limit)
+        break;
+      isPrime[n] = false;
+    }
+  }
+
+  /*
+   * // Put in candidate primes:
+//   integers which have an odd number of
+//   representations by certain quadratic forms.
+// Algorithm step 3.1:
+for n ≤ limit, n ← 4x²+y² where x ∈ {1,2,...} and y ∈ {1,3,...} // all x's odd
+y's if n mod 60 ∈ {1,13,17,29,37,41,49,53}: is_prime(n) ← ¬is_prime(n)   //
+toggle state
+    */
+  WORD mask1[] = {1, 13, 17, 29, 37, 41, 49, 53};
+  WORD maskLen1 = 8;
+  for (WORD x = 1; 4 * x * x < limit; x += 1) {
+    for (WORD y = 1; 1; y += 2) {
+      WORD n = 4 * x * x + y * y;
+      if (n >= limit) {
+        break;
+      }
+      WORD mod = n % 60;
+      for (WORD i = 0; i < maskLen1; i++) {
+        if (mod == mask1[i]) {
+          isPrime[n] = ~isPrime[n];
+          // printf("3.1 %lu is %s prime\n", n, (isPrime[n]) ? "" : "not");
+          break;
+        }
+      }
+    }
+  }
+
+  /*
+  // Algorithm step 3.2:
+for n ≤ limit, n ← 3x²+y² where x ∈ {1,3,...} and y ∈ {2,4,...} // only odd x's
+    if n mod 60 ∈ {7,19,31,43}:                                 // and even y's
+        is_prime(n) ← ¬is_prime(n)   // toggle state
+    */
+  WORD mask2[] = {7, 19, 31, 43};
+  WORD maskLen2 = 4;
+  for (WORD x = 1; x * x * 3 < limit; x += 2) {
+    for (WORD y = 2; y * y < limit; y += 2) {
+
+      WORD n = 3 * x * x + y * y;
+      if (n >= limit) {
+        break;
+      }
+      WORD mod = n % 60;
+      for (WORD i = 0; i < maskLen2; i++) {
+        if (mod == mask2[i]) {
+          isPrime[n] = ~isPrime[n];
+          // printf("3.2 %lu is %s prime\n", n, (isPrime[n]) ? "" : "not");
+          break;
+        }
+      }
+    }
+  }
+
+  /*
+  // Algorithm step 3.3:
+  for n ≤ limit, n ← 3x²-y² where x ∈ {2,3,...} and y ∈ {x-1,x-3,...,1} //all
+  even/odd if n mod 60 ∈ {11,23,47,59}:                                   //
+  odd/even combos is_prime(n) ← ¬is_prime(n)   // toggle state
+  */
+  WORD mask3[] = {11, 23, 47, 59};
+  WORD maskLen3 = 4;
+  for (WORD x = 2; (3 * x * x) - ((x - 1) * (x - 1)) < limit; x += 1) {
+    WORD y = x - 1;
+    while (y >= 1) {
+      // printf("x: %lu, y:%lu\n", x, y);
+      // sleep(1);
+
+      WORD n = 3 * x * x - y * y;
+
+      if (n >= limit) {
+        // breakOuter = true;
+        if (y < 2) {
+          break;
+        }
+        y -= 2;
+        continue;
+      }
+      WORD mod = n % 60;
+      for (WORD i = 0; i < maskLen3; i++) {
+        if (mod == mask3[i]) {
+          isPrime[n] = ~isPrime[n];
+          // printf("3.3 %lu is %s prime\n", n, (isPrime[n]) ? "" : "not");
+          break;
+        }
+      }
+      y = (y > 2) ? y - 2 : 0;
+    }
+  }
+
+  /*
+   *
+// Eliminate composites by sieving, only for those occurrences on the wheel:
+for n² ≤ limit, n ← 60 × w + x where w ∈ {0,1,...}, x ∈ s, n ≥ 7:
+    if is_prime(n):
+        // n is prime, omit multiples of its square; this is sufficient
+        // because square-free composites can't get on this list
+        for c ≤ limit, c ← n² × (60 × w + x) where w ∈ {0,1,...}, x ∈ s:
+            is_prime(c) ← false
+
+    */
+
+  for (WORD w = 0; w / 60 < limit; w++) {
+    WORD n = ~0;
+    for (WORD i = 0; i < wheelLen; i++) {
+      n = 60 * w + s[i];
+      if (n * n >= limit) {
+        break;
+      }
+      // or break;?
+
+      if (n < limit && isPrime[n]) {
+        for (WORD ww = 0; n * n * (60 * ww) < limit; ww++) {
+          WORD c;
+          for (WORD ii = 0; ii < wheelLen; ii++) {
+            c = n * n * (60 * ww + s[ii]);
+            if (c < limit) {
+              isPrime[c] = false;
+              // printf("elim: %lu is %s prime\n", c, (isPrime[c]) ? "" :
+              // "not");
+            }
+          }
+        }
+      }
+    }
+  }
+
+  // isPrime[2] = true;
+  // isPrime[3] = true;
+  // isPrime[5] = true;
+  printf("2 3 5");
+  WORD sum = 3;
+  for (WORD w = 0; w / 60 < limit; w++) {
+    WORD n = ~0;
+    for (WORD i = 0; i < wheelLen; i++) {
+      n = 60 * w + s[i];
+      if (n < 7)
+        continue;
+      if (n >= limit)
+        break;
+      if (isPrime[n]) {
+        // printf(" %lu ", n);
+        sum += 1;
+      }
+    }
+  }
+  printf("\n");
+  printf("%lu\n", sum);
+
+  return 0;
+}
diff --git a/ub10/a.out b/ub10/a.out
new file mode 100755
index 0000000000000000000000000000000000000000..72e439dd9403f5c4e1ae44a35baa3c1a26250fb0
GIT binary patch
literal 15792
zcmeHOe{37o9e;M)H3ib7UCUBh=nW_<+u|l^ngvS0ZQ}F-Z7C()4+&<R*h#!Pwr4vf
zsbC8$L9?=MWfFr))dtd}F(f9m@dsmw89HG)ZBx;tshZFh7-){lSlPl9Mpy6ieee6d
zTy6<MOqw?7o^<zl-yfgvd+&Sa^SiruUykhV4h8~B$`W<GB5AzdXu@o0YU2tZtlCvA
z_Lr$k)rTN2)jDGj8-Uf6)4};z!}=va;x$lZ0o-8JN>F-8NW8^L`*LFuR6GviEh7`l
zUc5ycl@b)Y)1F0S6ie?Nyo(%A1l>H!&8bFhR;N1V;}GX1%P2e&i8sZ(DdsuI?M*RC
zJ}D=(wR1kb7>12XSnrkc@<y5G#V+<2q&BGB`g$S$LE78Lyz49-MtM=%5tR0B0gv+X
zev2ID{oLMS$DxVyFUV<CTT;osw$?4Fcw;J=9vN*MZEtIAYYpWyp%&c*>Q@ssI&F6C
zJ)nZ)szf^B2UK959Ap_Tx{&W0Y?OB5w>`Tq`{7MrA9=TD=i_Y+7kuO0H}Crz`Or8h
zCizey9c@y`CY$Kq4}Jtrs~2#*ViElz=&c3(sqKa0(6L%b-?fPTo6uX>yoG-M3&nW^
zdduq!c^kIN7tpKhNFlGH(R?A+e=yoVbTB#)OQzI7E|E~fu~aJ44<n9wKY8Q|(cxG!
z4Nd=OEIN=($5P3|36;$y(}e*=UUVX~sa?B!c6LTvLak11>s+-Zv`s~OdiO-*iCkhZ
znJ*-Ay?Z)SnRKEz)|W~k_~39R&CNxP>w=CLO%Ufv0Q;rb2x&7V0Lx9N%K5$~nXIKV
zgw7VyJ$e7Pa443jD9a<4Bv(+pFF>Xu*Cn|gWLFveX}<1=-eP@I*w{rca_I)6r{gYt
zBh8C0Jzrw9<0+S(FWFju+NDQFoigRpU*Kp_2^)c8(77n3CTwb?cmH1MU3$LcYsUtc
zp2i@>3Y!|mkui(SM?V982K)^88SpdUXTZ<E|6>OJR(IJy>?3c~*pthi+oF^`Q7Qzh
zr|cs?uQ{QIQ{DCpsH>ZPiGAJ1Flr<Zk!^aain8f>l4<HSeM-wulT1^<>5`UzL^4hN
zrk~LAV<gj*Z2Doy)PD1ydx|f8-Y&jsA9?NUzTQZ4srjURv|~GhRM)PD`^-S7ZsWb6
zpo2#{?t#WWy6iB)HM51a==%AP>H3Q5RNcmL@;ku@|68>Gwi^kqK5rM#+RwdtgZ<p>
z5<BpM{nGiu8U*+i2dJq|4b*MKO^WcTf8!ms;Hr_W2kavqkHJ&f#n%h9_GHJK5UOh*
z!oJdv!Mw2SHB<w)pP1|KiGjH~<&WY`55o5VVj%IBUA(*UGn%PaUO8W_+QrDMedO+0
zHM0CE=+v3@_9Qfu2P$^4URTRuJq@M$$rCyV0@%foDQJL^*)tazO1MS`({$7)N+W+V
zh67V}abL}`+YOt$fP}t{8wa~su3Sz1c<*47Cg30ES}Ipw{imkQ+T6oq+&ylf9GRg3
zIe|M2S51anF+2=$wOxFv_`F^GHM%ltx}sB1`3V9aiOdA%8j~D7GZ~pFz5xH2+rJ2Z
zY5blyns92y?a8H^HdE`!PjktB3iosU3CK_DE)^#f4v1ZRpoXA4Q4i3D@<hX2)takr
z=ITIsqDfcF6RkkHrv^QYYG0<0sx~{Kx^J%f*}3W+bJZcPB4}1u%M%6QWjbhNd^nan
zXqO{jH64o_V|<A5Va7)Z2g;ErlqyG_Y=Qt%3EgP{zQ|6e7*8`!F;)mcK1<}8>mb-i
zIXC6XO>}&V6Q}9yII0IpC(AaYozmLKIH@i-s)|vSNL6Q4XZ4w7LKG)vwe*0dP8zBB
zTk^h&yeA*f0p=NMje)~xzG2LqBL>GPSMJq*j<{yNH~{fXoy*<L<}>>hmwCmz)v7)d
z={$|wpID<$rNXjGJq#ytnRKq9&Q+80=odc&eg^yu_!;ms;Ag<kfS&<B1OIy&pf^eJ
zL@H58rU$J8y%O7$x30A|#jUYKESJbznN-}e8n3r3bpe#=OdLAA#2Um)n85la+rt<V
zy&HM$FV$)R*mbsAy$#s$R<(K@_!ZzYz!LCf;79*jt)2sZ2e<*ZX7mn}O5pH*6&S4#
ztiNb^%`w(*fPM;buOf%8<@E$tVr#`V1pOUAwX(i@Wy6hi7v5Pju5P$w`^UFlc{!<x
zdmFZIqdyT+>{wZUU$C?GL%}dwF!sBkeiruY-S#nU*8JnpzYKfKZ4ch&@SlPHQjF{F
zdG=Aq{xbB3Vb@PXr2QrRLDK?}{NiW8&w!r+KLdUS{0#UR@H60Nz|Vl6fe)SmS-&Ui
z^>8bt3#}2Pa*+XN!606o>C&~<;5wFNUEq3_am%F3C5*B@a0AO?fBSqjL#h%l_LDVs
z6?)@~B5UiW*e+}Nern9>LdFJ0?iNdCy(|e8UJpnwnyHK$aPBARo|Fo`lcQ4NhD9$O
zk$#u7mr}CkP};{U9z7nIwJY9Pnep+!@P<yyl4sm5YWaI?e}`orF1PaE2U$PvZMjX~
zRV&ra_-V!g#zTz4ll7jvIy<kmHe*%t)mC$Bs3p{FwKTP~HnlfxQ;V3cX>O7r3n1lh
zsd`JlclDYV($jN;fLg1(<GNInS2^p-B^I?Ez}X<@jp(~kSuKJKKeLFQo^e<%KdJs+
z>l@VOuxJb(11uE(5$NkxgK7^OKiNNu%IaX9lHa4~e@J=}TzGO3Jw1D%_}>0ZLvLa8
zmRTrP1sjwcU$HO4fmuZ>i7lhypRrm4{X+fu81(gOoeGB~g27h-tJO;9yOOKe{7D^W
zrFxdv7wcy+uziZ-c=OZA{=ehK<vR)+dg{N|zZZIo+j2vi<DBDmV{AMGy?%h?c1Hbv
z@(<$biB!1DgB)j!<LGBiaQH6kZ|8pMXGBEI<3U?$Y<Dt!FcII4RpB8lJn#D4-W_{-
zI<X`iKd5NLj&|?gu_qF>yY{Q-uH859+_5`)b9eWxk=|(Uj-9(BQN3(EEUVA;LUZrh
zb;pVo>F>1)otDKrzB*DO9xKEYt$I(QL-W_mJD#n)7~WBsP7u|sX#cn6?a_EX6CH}B
z<Fxd>=Vo}tlj-P4J`v~j_wz#~(s32bW&wNg$Pv-{@_CNTOY^aqJ{sBE6@9Oz^LVn^
zOV2dtFXW%!62@6j>BOD79lgY#4vSKu{Mc|I)(0%)3=T;Rp8+ItSrtlW3W?BQdL)$1
zWwMD}Va%oK8%g5hkz|}zJ9hRo7Gi@+yAQ?kLn;&>OQQ*cg`BY*O62m%OnROWg)Nsz
z#fZSQY^tC_y5pfjVia|Kv_iRzK6atRP;?*{8%{)r;^2y87|~cR7ehyhA$0_f4JZ52
zY6h{e3+V<-FeF}u`ZL4Bi8NjO{=-7^3H&siZv-0nk5N$O6VeR%y2}+1=s#Tk_uH8-
z7<N0OF8}YMrq9^i{88o$((`3^k?Wo4dQfrmk27D;5=TZapZ-_=jyZw!U*<`IGJcsS
zdHWv&pXNxymwA+64JVBHNLAs>d?pDUby@f_FB6n`7+iGm=KqV((ws;5GS3mbn0eBF
z@e}+K>@<fGxy17%!7%f^{b#j0iVEda_%e?ZEU}}s?;Zbx%-_j=GA|UAc_i^^?&#$|
z0v+W<;yeEb{J$XcWw_|#J^sg`b@MIe3Ep97&R+g*1i$a$kF!DWc~67FlYozV_%hEB
zl=}hk_k>eFX1<eui#-{)vnb%@)A>c?7yfmY@f7SbqO<q%X};s;k6Xr7@SMXm6*+#Q
z{tYTCvGLE=`Tj}15A^zrEcgd->s)-f4?5ilC5reNX-*dX85ro?k@jW&>U{4>gF~V7
zpNjD0IL<&vF2a}Zu1Z_TTdB<koW1Y`-+{%=m-%po`C*4?D#8;Cg6!tY{Bm=X@t1x|
z7~wfG$fT40%Y3b!`Cfl&%Ux=@eRg+D!gqb9tZl+~qnvvV8I`OA<c`sGn2wY5N3NeM
jpp(9%>XhtOQ~iN7#1@Z;x%9&Pb@bsP%EcanM^yg;#n9Uy

literal 0
HcmV?d00001