Development of
Algorithmic Constructions

21:14:28

19.Apr 2024

Sieve of Ulam (vertical II)

mathematical description

Instead of sieving all odd numbers by the order of rising numbers,
the numbers are sorted in columns and each colomn is sieved:

1	3	5	7	9	11	13	15	17	19	21	23	25	27	29
31	33	35	37	39	41	43	45	47	49	51	53	55	57	59
61	63	65	67	69	71	73	75	77	79	81	83	85	87	89
91	93	95	97	99	101	103	105	107	109	111	113	115	117	119
121	123	125	127	129	131	133	135	137	139	141	143	145	147	149
151	153	155	157	159	161	163	165	167	169	171	173	175	177	179
181	183	185	187	189	191	193	195	197	199	201	203	205	207	209
211	213	215	217	219	221	223	225	227	229	231	233	235	237	239
241	243	245	247	249	251	253	255	257	259	261	263	265	267	269
271	273	275	277	279	281	283	285	287	289	291	293	295	297	299
301	303	305	307	309	311	313	315	317	319	321	323	325	327	329
331	333	335	337	339	341	343	345	347	349	351	353	355	357	359
361	363	365	367	369	371	373	375	377	379	381	383	385	387	389
391	393	395	397	399	401	403	405	407	409	411	413	415	417	419
421	423	425	427	429	431	433	435	437	439	441	443	445	447	449
451	453	455	457	459	461	463	465	467	469	471	473	475	477	479
481	483	485	487	489	491	493	495	497	499	501	503	505	507	509
511	513	515	517	519	521	523	525	527	529	531	533	535	537	539
541	543	545	547	549	551	553	555	557	559	561	563	565	567	569
571	573	575	577	579	581	583	585	587	589	591	593	595	597	599
601	603	605	607	609	611	613	615	617	619	621	623	625	627	629
631	633	635	637	639	641	643	645	647	649	651	653	655	657	659
661	663	665	667	669	671	673	675	677	679	681	683	685	687	689
691	693	695	697	699	701	703	705	707	709	711	713	715	717	719
721	723	725	727	729	731	733	735	737	739	741	743	745	747	749
751	753	755	757	759	761	763	765	767	769	771	773	775	777	779
781	783	785	787	789	791	793	795	797	799	801	803	805	807	809
811	813	815	817	819	821	823	825	827	829	831	833	835	837	839
841	843	845	847	849	851	853	855	857	859	861	863	865	867	869
871	873	875	877	879	881	883	885	887	889	891	893	895	897	899

The columns started with 3, 5, 9, 15, 21, 25, 27 can be disregarded because gcd (column, 30)>1

For this example all numbers are smaller than 30*30.
Therefore after sieving the coloumns by the primes < 30 only the primes remains.

Regarding a column where the first number is r.
The first occurance of a prime pi > 5 in this column has to be of the form
pi*t = 1 mod 30
t = pi^(-1) mod 30

The inverse of pi mod 30 can be calculated one time in the beginning of the algorithm.

The first occurance of pi in the row r can be calculated by
y=r*[(pi*pi^(-1)-1) / 30] mod pi

As the value of y.th row is divisible by pi, all values y+pi | pi

All numbers which are divisible by pi are sieved out.

The first column is preset by 0, the numbers are not calculated.
The sieving vertical is made by flipping one bit from 0 to 1 from the starting point up to a.

The primes which remain on the column are representing by the zeros,
All non primes which are representant by 1 are reset to zero for the next column.

The number of zeros are counted.

You get the number of all primes concerning each colonm without the sieved primes.

Results:

Based on a Implementation in Cuda on a graphiccard Geforce 650 with 192 Cuda cores

a) a=2*3*5*7*11*13*17=510.510
numbers below a*a = 260620460100 ~ 2,6 * 10^11
columns which have to be examined 18.1 % = 92.159
Runtime : 17,51 min

b) a=2*3*5*7*11*13*17*19 = 9.699.690
numbers below a*a = 94083986096100 ~ 9,4 * 10^13
columns which have to be examined 17.1 % = 1.658.879
Runtime : ~ 5752 min

For comparison :

a) My best implementation on a single processor needs
228 min for the calculation up to 10^12
(Athlon 64 4000+ 128kbyte 1. Level Cache)
The Cuda implementation is therefore more than three time faster.

Implementation

for Cuda
in Mupad

m:=30;
p:={7,11,13,17,19,23,29};
for i from 1 to 7 do
   p_inv:=p[i]^(-1) mod  m;
   p_start[i]:=(p[i]*p_inv-1)/m;
end_for;

for j from 1 to m step 2 do
   if gcd (j, m)=1 then

      // Preset coloum with 0
      for i from 0 to m-1 do
         liste[i]:=0;
      end_for;

      // Sieving
      for i from 1 to nops (p) do
        start:=p_start[i]*j mod p[i];
        print (j, p[i], start);
        while start<m do
            liste[start]:=1;
            start:=start+p[i]
        end_while;
      end_for;

      // Result
      for i from 0 to m-1 do
         if liste [i]=0 then
            print (i*m+j, isprime (i*m+j));
         end_if;
      end_for;

   end_if;
end_for;

1	3	5	7	9	11	13	15	17	19	21	23	25	27	29
31	33	35	37	39	41	43	45	47	49	51	53	55	57	59
61	63	65	67	69	71	73	75	77	79	81	83	85	87	89
91	93	95	97	99	101	103	105	107	109	111	113	115	117	119
121	123	125	127	129	131	133	135	137	139	141	143	145	147	149
151	153	155	157	159	161	163	165	167	169	171	173	175	177	179
181	183	185	187	189	191	193	195	197	199	201	203	205	207	209
211	213	215	217	219	221	223	225	227	229	231	233	235	237	239
241	243	245	247	249	251	253	255	257	259	261	263	265	267	269
271	273	275	277	279	281	283	285	287	289	291	293	295	297	299
301	303	305	307	309	311	313	315	317	319	321	323	325	327	329
331	333	335	337	339	341	343	345	347	349	351	353	355	357	359
361	363	365	367	369	371	373	375	377	379	381	383	385	387	389
391	393	395	397	399	401	403	405	407	409	411	413	415	417	419
421	423	425	427	429	431	433	435	437	439	441	443	445	447	449
451	453	455	457	459	461	463	465	467	469	471	473	475	477	479
481	483	485	487	489	491	493	495	497	499	501	503	505	507	509
511	513	515	517	519	521	523	525	527	529	531	533	535	537	539
541	543	545	547	549	551	553	555	557	559	561	563	565	567	569
571	573	575	577	579	581	583	585	587	589	591	593	595	597	599
601	603	605	607	609	611	613	615	617	619	621	623	625	627	629
631	633	635	637	639	641	643	645	647	649	651	653	655	657	659
661	663	665	667	669	671	673	675	677	679	681	683	685	687	689
691	693	695	697	699	701	703	705	707	709	711	713	715	717	719
721	723	725	727	729	731	733	735	737	739	741	743	745	747	749
751	753	755	757	759	761	763	765	767	769	771	773	775	777	779
781	783	785	787	789	791	793	795	797	799	801	803	805	807	809
811	813	815	817	819	821	823	825	827	829	831	833	835	837	839
841	843	845	847	849	851	853	855	857	859	861	863	865	867	869
871	873	875	877	879	881	883	885	887	889	891	893	895	897	899

1	3	5	7	9	11	13	15	17	19	21	23	25	27	29
31	33	35	37	39	41	43	45	47	49	51	53	55	57	59
61	63	65	67	69	71	73	75	77	79	81	83	85	87	89
91	93	95	97	99	101	103	105	107	109	111	113	115	117	119
121	123	125	127	129	131	133	135	137	139	141	143	145	147	149
151	153	155	157	159	161	163	165	167	169	171	173	175	177	179
181	183	185	187	189	191	193	195	197	199	201	203	205	207	209
211	213	215	217	219	221	223	225	227	229	231	233	235	237	239
241	243	245	247	249	251	253	255	257	259	261	263	265	267	269
271	273	275	277	279	281	283	285	287	289	291	293	295	297	299
301	303	305	307	309	311	313	315	317	319	321	323	325	327	329
331	333	335	337	339	341	343	345	347	349	351	353	355	357	359
361	363	365	367	369	371	373	375	377	379	381	383	385	387	389
391	393	395	397	399	401	403	405	407	409	411	413	415	417	419
421	423	425	427	429	431	433	435	437	439	441	443	445	447	449
451	453	455	457	459	461	463	465	467	469	471	473	475	477	479
481	483	485	487	489	491	493	495	497	499	501	503	505	507	509
511	513	515	517	519	521	523	525	527	529	531	533	535	537	539
541	543	545	547	549	551	553	555	557	559	561	563	565	567	569
571	573	575	577	579	581	583	585	587	589	591	593	595	597	599
601	603	605	607	609	611	613	615	617	619	621	623	625	627	629
631	633	635	637	639	641	643	645	647	649	651	653	655	657	659
661	663	665	667	669	671	673	675	677	679	681	683	685	687	689
691	693	695	697	699	701	703	705	707	709	711	713	715	717	719
721	723	725	727	729	731	733	735	737	739	741	743	745	747	749
751	753	755	757	759	761	763	765	767	769	771	773	775	777	779
781	783	785	787	789	791	793	795	797	799	801	803	805	807	809
811	813	815	817	819	821	823	825	827	829	831	833	835	837	839
841	843	845	847	849	851	853	855	857	859	861	863	865	867	869
871	873	875	877	879	881	883	885	887	889	891	893	895	897	899

Development ofAlgorithmic Constructions

Sieve of Ulam (vertical II)

mathematical description

Results:

Implementation

Development of
Algorithmic Constructions

1	3	5	7	9	11	13	15	17	19	21	23	25	27	29
31	33	35	37	39	41	43	45	47	49	51	53	55	57	59
61	63	65	67	69	71	73	75	77	79	81	83	85	87	89
91	93	95	97	99	101	103	105	107	109	111	113	115	117	119
121	123	125	127	129	131	133	135	137	139	141	143	145	147	149
151	153	155	157	159	161	163	165	167	169	171	173	175	177	179
181	183	185	187	189	191	193	195	197	199	201	203	205	207	209
211	213	215	217	219	221	223	225	227	229	231	233	235	237	239
241	243	245	247	249	251	253	255	257	259	261	263	265	267	269
271	273	275	277	279	281	283	285	287	289	291	293	295	297	299
301	303	305	307	309	311	313	315	317	319	321	323	325	327	329
331	333	335	337	339	341	343	345	347	349	351	353	355	357	359
361	363	365	367	369	371	373	375	377	379	381	383	385	387	389
391	393	395	397	399	401	403	405	407	409	411	413	415	417	419
421	423	425	427	429	431	433	435	437	439	441	443	445	447	449
451	453	455	457	459	461	463	465	467	469	471	473	475	477	479
481	483	485	487	489	491	493	495	497	499	501	503	505	507	509
511	513	515	517	519	521	523	525	527	529	531	533	535	537	539
541	543	545	547	549	551	553	555	557	559	561	563	565	567	569
571	573	575	577	579	581	583	585	587	589	591	593	595	597	599
601	603	605	607	609	611	613	615	617	619	621	623	625	627	629
631	633	635	637	639	641	643	645	647	649	651	653	655	657	659
661	663	665	667	669	671	673	675	677	679	681	683	685	687	689
691	693	695	697	699	701	703	705	707	709	711	713	715	717	719
721	723	725	727	729	731	733	735	737	739	741	743	745	747	749
751	753	755	757	759	761	763	765	767	769	771	773	775	777	779
781	783	785	787	789	791	793	795	797	799	801	803	805	807	809
811	813	815	817	819	821	823	825	827	829	831	833	835	837	839
841	843	845	847	849	851	853	855	857	859	861	863	865	867	869
871	873	875	877	879	881	883	885	887	889	891	893	895	897	899