topic Another math puzzle in SAS Programming

Another math puzzle

FriedEgg — Thu, 06 Oct 2011 19:55:51 GMT

Hi,

Here is another good problem I found to help test you solver skills. This question is curtosy of Google. Specifically the Google Treasure Hunt 2008. Which is available for anyone to use here: http://treasurehunt.appspot.com

Since it varys the question I will paste the version I receieved here so everyone can work on the same set problem:

Question:

Find the smallest number that can be expressed as
the sum of 19 consecutive prime numbers,
the sum of 21 consecutive prime numbers,
the sum of 405 consecutive prime numbers,
the sum of 781 consecutive prime numbers,
and is itself a prime number.

For example, 41 is the smallest prime number that can be expressed as
the sum of 3 consecutive primes (11 + 13 + 17 = 41) and
the sum of 6 consecutive primes (2 + 3 + 5 + 7 + 11 + 13 = 41).

Another math puzzle

Ksharp — Fri, 07 Oct 2011 09:23:50 GMT

It looks like SAS can not offer me a prime numbers generator.

Does anyone tell me how to call prime numbers generator in SAS.

The generator I make spends me lots of time to get prime numbers, So I only give the answer for 19 and 21.

prim_sum=857
p[5]=11
p[6]=13
p[7]=17
p[8]=19
p[9]=23
p[10]=29
p[11]=31
p[12]=37
p[13]=41
p[14]=43
p[15]=47
p[16]=53
p[17]=59
p[18]=61
p[19]=67
p[20]=71
p[21]=73
p[22]=79
p[23]=83




prim_sum=953
p[4]=7
p[5]=11
p[6]=13
p[7]=17
p[8]=19
p[9]=23
p[10]=29
p[11]=31
p[12]=37
p[13]=41
p[14]=43
p[15]=47
p[16]=53
p[17]=59
p[18]=61
p[19]=67
p[20]=71
p[21]=73
p[22]=79
p[23]=83
p[24]=89









data _null_;
array p{4000} _temporary_;
 do i=2 to 20000;
  count=0;
  do j=1 to i;
   if mod(i,j)=0 then count+1;
  end;
  if count=2 then do;x+1;p{x}=i;end;
 end;
*get 19 consecutive prim number;
 do m=1 to x-18;
  prim_sum=0;
  do n=m to m+18;
  prim_sum+p{n};
  end;
  if prim_sum in p then do;
   put prim_sum= ;
   do _n=m to m+18;
    put p{_n}=;
   end;
   stop;
  end;
 end;
run;









data _null_;
array p{4000} _temporary_;
 do i=2 to 20000;
  count=0;
  do j=1 to i;
   if mod(i,j)=0 then count+1;
  end;
  if count=2 then do;x+1;p{x}=i;end;
 end;
*get 21 consecutive prim number;
 do m=1 to x-20;
  prim_sum=0;
  do n=m to m+20;
  prim_sum+p{n};
  end;
  if prim_sum in p then do;
   put prim_sum= ;
   do _n=m to m+20;
    put p{_n}=;
   end;
   stop;
  end;
 end;
run;

Ksharp

Another math puzzle

DF — Fri, 07 Oct 2011 10:53:00 GMT

I haven't done much on this one, but the below is the method I was planning to use to get a list of primes. It's written in 9.1.3, so you might have to alter the hash functions for the newer versions. Also I think it can be made to run a little better - I'm sure it's been proven that you can stop testing for factors when you pass sqrt(integer), although I'd need to check.

For reference, the algorithm used in the sieve of erastothenes. Rather than dividing by all integers, you only need to divide by known primes (all non-primes have prime factors).

On my system running for 20000 integers takes a few seconds.

options compress=no;

data z;

format integer 8.;

format prime 8.;

if _n_ = 1 then do;

declare hash h(ordered:'a');

h.defineKey('prime');

h.defineData('prime');

h.defineDone();

declare hiter hi('h');

end;

do i = 2 to 20000;

integer = i;

divisor_found = 0;

rc = hi.first();

do while (rc = 0);

if mod(integer,prime) = 0 then do;

divisor_found = 1;

leave;

end;

rc = hi.next();

end;

if divisor_found = 0 then do;

putlog _all_;

prime = integer;

rc = h.add();

end;

h.output(dataset:'work.primes');

run;

Another math puzzle

Ksharp — Fri, 07 Oct 2011 14:01:35 GMT

Thanks. DF

I also notice this. So recode to optimize it.

But it still looks like that I need super computer to get the 405 781 .

data _null_;
array p{100000} _temporary_;
retain x 1;
p{1}=2;
 do i=2 to 200000;
  prime=1;_k=1;
  do while(not missing(p{_k}));
  if mod(i,p{_k})=0 then prime=0;
  _k+1;
  end;
  if prime then do; x+1;p{x}=i;end;
 end;


*get 405 consecutive prim number;
 do m=1 to x-404;
  prim_sum=0;
  do n=m to m+404;
  prim_sum+p{n};
  end;
  if prim_sum in p then do;
   put prim_sum= ;
   do _n=m to m+404;
    put p{_n}=;
   end;
   stop;
  end;
 end;
 run;

Ksharp

Another math puzzle

Ksharp — Fri, 07 Oct 2011 14:11:42 GMT

But I have gotten 100 consecutive prime number.

prim_sum=2413

p[1]=2

p[2]=3

p[3]=5

p[4]=7

p[5]=11

p[6]=13

p[7]=17

p[8]=19

p[9]=23

p[10]=29

p[11]=31

p[12]=37

p[13]=41

p[14]=43

p[15]=47

p[16]=53

p[17]=59

p[18]=61

p[19]=67

p[20]=71

p[21]=73

p[22]=79

p[23]=83

p[24]=89

p[25]=97

p[26]=101

p[27]=103

p[28]=107

p[29]=109

p[30]=113

p[31]=127

p[32]=131

p[33]=137

p[34]=139

p[35]=149

p[36]=151

p[37]=157

p[38]=163

p[39]=167

p[40]=173

p[41]=179

p[42]=181

p[43]=191

p[44]=193

p[45]=197

p[46]=199

p[47]=211

p[48]=223

p[49]=227

p[50]=229

p[51]=233

p[52]=239

p[53]=241

p[54]=251

p[55]=257

p[56]=263

p[57]=269

p[58]=271

p[59]=277

p[60]=281

p[61]=283

p[62]=293

p[63]=307

p[64]=311

p[65]=313

p[66]=317

p[67]=331

p[68]=337

p[69]=347

p[70]=349

p[71]=353

p[72]=359

p[73]=367

p[74]=373

p[75]=379

p[76]=383

p[77]=389

p[78]=397

p[79]=401

p[80]=409

p[81]=419

p[82]=421

p[83]=431

p[84]=433

p[85]=439

p[86]=443

p[87]=449

p[88]=457

p[89]=461

p[90]=463

p[91]=467

p[92]=479

p[93]=487

p[94]=491

p[95]=499

p[96]=503

p[97]=509

p[98]=521

p[99]=523

p[100]=541

NOTE: DATA st

real ti

cpu tim

Ksharp

Another math puzzle

FriedEgg — Fri, 07 Oct 2011 16:19:48 GMT

So, the first trick to this puzzle is figuring out the best and most efficient way to calculate prime numbers. I will give you simple way to generate a list of primes.

data prime;

do n=1 to 10**8;

m=n+int((n-2)/8)+int((n-4)/8);

p=m+int((m-3)/2);

a=p+abs(p-3/2)+1/2;

output;

end;

keep a;

run;

With this method I can calculate 100,000,000 prime numbers in a matter of secods...

Another math puzzle

DF — Fri, 07 Oct 2011 19:14:00 GMT

I don't think primes can be calculating using such a method, without checking for factors. For example, for n=23, you get:

n=23 m=27 p=39 a=77 _ERROR_=0 _N_=1

77 is divisible by 7 and 11, and therefore is not prime.

Another math puzzle

FriedEgg — Fri, 07 Oct 2011 19:21:18 GMT

I got this method from the Encyclopedia of Integer Sequences, just assumed it worked... I will check to make sure I implemented it correctly, as 77 should not be a number the formula calculates.

Another math puzzle

FriedEgg — Fri, 07 Oct 2011 19:33:12 GMT

After further reading this method is only valid for the first 15 prime numbers The formula breaks down on the 16th iteration (a(16)=49).

Another math puzzle

Ksharp — Sat, 08 Oct 2011 03:11:26 GMT

Haha.

Finally I got 405 .

prim_sum=541283
p[8]=19
p[9]=23
p[10]=29
p[11]=31
p[12]=37
p[13]=41
p[14]=43
p[15]=47
p[16]=53
p[17]=59
p[18]=61
p[19]=67
p[20]=71
p[21]=73
p[22]=79
p[23]=83
p[24]=89
p[25]=97
p[26]=101
p[27]=103
p[28]=107
p[29]=109
p[30]=113
p[31]=127
p[32]=131
p[33]=137
p[34]=139
p[35]=149
p[36]=151
p[37]=157
p[38]=163
p[39]=167
p[40]=173
p[41]=179
p[42]=181
p[43]=191
p[44]=193
p[45]=197
p[46]=199
p[47]=211
p[48]=223
p[49]=227
p[50]=229
p[51]=233
p[52]=239
p[53]=241
p[54]=251
p[55]=257
p[56]=263
p[57]=269
p[58]=271
p[59]=277
p[60]=281
p[61]=283
p[62]=293
p[63]=307
p[64]=311
p[65]=313
p[66]=317
p[67]=331
p[68]=337
p[69]=347
p[70]=349
p[71]=353
p[72]=359
p[73]=367
p[74]=373
p[75]=379
p[76]=383
p[77]=389
p[78]=397
p[79]=401
p[80]=409
p[81]=419
p[82]=421
p[83]=431
p[84]=433
p[85]=439
p[86]=443
p[87]=449
p[88]=457
p[89]=461
p[90]=463
p[91]=467
p[92]=479
p[93]=487
p[94]=491
p[95]=499
p[96]=503
p[97]=509
p[98]=521
p[99]=523
p[100]=541
p[101]=547
p[102]=557
p[103]=563
p[104]=569
p[105]=571
p[106]=577
p[107]=587
p[108]=593
p[109]=599
p[110]=601
p[111]=607
p[112]=613
p[113]=617
p[114]=619
p[115]=631
p[116]=641
p[117]=643
p[118]=647
p[119]=653
p[120]=659
p[121]=661
p[122]=673
p[123]=677
p[124]=683
p[125]=691
p[126]=701
p[127]=709
p[128]=719
p[129]=727
p[130]=733
p[131]=739
p[132]=743
p[133]=751
p[134]=757
p[135]=761
p[136]=769
p[137]=773
p[138]=787
p[139]=797
p[140]=809
p[141]=811
p[142]=821
p[143]=823
p[144]=827
p[145]=829
p[146]=839
p[147]=853
p[148]=857
p[149]=859
p[150]=863
p[151]=877
p[152]=881
p[153]=883
p[154]=887
p[155]=907
p[156]=911
p[157]=919
p[158]=929
p[159]=937
p[160]=941
p[161]=947
p[162]=953
p[163]=967
p[164]=971
p[165]=977
p[166]=983
p[167]=991
p[168]=997
p[169]=1009
p[170]=1013
p[171]=1019
p[172]=1021
p[173]=1031
p[174]=1033
p[175]=1039
p[176]=1049
p[177]=1051
p[178]=1061
p[179]=1063
p[180]=1069
p[181]=1087
p[182]=1091
p[183]=1093
p[184]=1097
p[185]=1103
p[186]=1109
p[187]=1117
p[188]=1123
p[189]=1129
p[190]=1151
p[191]=1153
p[192]=1163
p[193]=1171
p[194]=1181
p[195]=1187
p[196]=1193
p[197]=1201
p[198]=1213
p[199]=1217
p[200]=1223
p[201]=1229
p[202]=1231
p[203]=1237
p[204]=1249
p[205]=1259
p[206]=1277
p[207]=1279
p[208]=1283
p[209]=1289
p[210]=1291
p[211]=1297
p[212]=1301
p[213]=1303
p[214]=1307
p[215]=1319
p[216]=1321
p[217]=1327
p[218]=1361
p[219]=1367
p[220]=1373
p[221]=1381
p[222]=1399
p[223]=1409
p[224]=1423
p[225]=1427
p[226]=1429
p[227]=1433
p[228]=1439
p[229]=1447
p[230]=1451
p[231]=1453
p[232]=1459
p[233]=1471
p[234]=1481
p[235]=1483
p[236]=1487
p[237]=1489
p[238]=1493
p[239]=1499
p[240]=1511
p[241]=1523
p[242]=1531
p[243]=1543
p[244]=1549
p[245]=1553
p[246]=1559
p[247]=1567
p[248]=1571
p[249]=1579
p[250]=1583
p[251]=1597
p[252]=1601
p[253]=1607
p[254]=1609
p[255]=1613
p[256]=1619
p[257]=1621
p[258]=1627
p[259]=1637
p[260]=1657
p[261]=1663
p[262]=1667
p[263]=1669
p[264]=1693
p[265]=1697
p[266]=1699
p[267]=1709
p[268]=1721
p[269]=1723
p[270]=1733
p[271]=1741
p[272]=1747
p[273]=1753
p[274]=1759
p[275]=1777
p[276]=1783
p[277]=1787
p[278]=1789
p[279]=1801
p[280]=1811
p[281]=1823
p[282]=1831
p[283]=1847
p[284]=1861
p[285]=1867
p[286]=1871
p[287]=1873
p[288]=1877
p[289]=1879
p[290]=1889
p[291]=1901
p[292]=1907
p[293]=1913
p[294]=1931
p[295]=1933
p[296]=1949
p[297]=1951
p[298]=1973
p[299]=1979
p[300]=1987
p[301]=1993
p[302]=1997
p[303]=1999
p[304]=2003
p[305]=2011
p[306]=2017
p[307]=2027
p[308]=2029
p[309]=2039
p[310]=2053
p[311]=2063
p[312]=2069
p[313]=2081
p[314]=2083
p[315]=2087
p[316]=2089
p[317]=2099
p[318]=2111
p[319]=2113
p[320]=2129
p[321]=2131
p[322]=2137
p[323]=2141
p[324]=2143
p[325]=2153
p[326]=2161
p[327]=2179
p[328]=2203
p[329]=2207
p[330]=2213
p[331]=2221
p[332]=2237
p[333]=2239
p[334]=2243
p[335]=2251
p[336]=2267
p[337]=2269
p[338]=2273
p[339]=2281
p[340]=2287
p[341]=2293
p[342]=2297
p[343]=2309
p[344]=2311
p[345]=2333
p[346]=2339
p[347]=2341
p[348]=2347
p[349]=2351
p[350]=2357
p[351]=2371
p[352]=2377
p[353]=2381
p[354]=2383
p[355]=2389
p[356]=2393
p[357]=2399
p[358]=2411
p[359]=2417
p[360]=2423
p[361]=2437
p[362]=2441
p[363]=2447
p[364]=2459
p[365]=2467
p[366]=2473
p[367]=2477
p[368]=2503
p[369]=2521
p[370]=2531
p[371]=2539
p[372]=2543
p[373]=2549
p[374]=2551
p[375]=2557
p[376]=2579
p[377]=2591
p[378]=2593
p[379]=2609
p[380]=2617
p[381]=2621
p[382]=2633
p[383]=2647
p[384]=2657
p[385]=2659
p[386]=2663
p[387]=2671
p[388]=2677
p[389]=2683
p[390]=2687
p[391]=2689
p[392]=2693
p[393]=2699
p[394]=2707
p[395]=2711
p[396]=2713
p[397]=2719
p[398]=2729
p[399]=2731
p[400]=2741
p[401]=2749
p[402]=2753
p[403]=2767
p[404]=2777
p[405]=2789
p[406]=2791
p[407]=2797
p[408]=2801
p[409]=2803
p[410]=2819
p[411]=2833
p[412]=2837
NOTE: DATA statement used (Total process time):
      real time           13:35.23
      cpu time            11:59.45











data _null_;
array p{500000} _temporary_;
retain x 1;
p{1}=2;
 do i=2 to 1000000;
  prime=1;_k=1;
  do while(not missing(p{_k}));
  if mod(i,p{_k})=0 then do;prime=0;leave;end;
  _k+1;
  end;
  if prime then do; x+1;p{x}=i;end;
 end;


*get 405 consecutive prim number;
 do m=1 to x-404;
  prim_sum=0;
  do n=m to m+404;
  prim_sum+p{n};
  end;
  if prim_sum in p then do;
   put prim_sum= ;
   do _n=m to m+404;
    put p{_n}=;
   end;
   stop;
  end;
 end;
 run;

Ksharp

Another math puzzle

Ksharp — Sat, 08 Oct 2011 07:53:04 GMT

I got 781. With the combiantion of Rick Aster's prime number generator. and also thanks Art297.

prim_sum=2174729
P[3]=5
P[4]=7
P[5]=11
P[6]=13
P[7]=17
P[8]=19
P[9]=23
P[10]=29
P[11]=31
P[12]=37
P[13]=41
P[14]=43
P[15]=47
P[16]=53
P[17]=59
P[18]=61
P[19]=67
P[20]=71
P[21]=73
P[22]=79
P[23]=83
P[24]=89
P[25]=97
P[26]=101
P[27]=103
P[28]=107
P[29]=109
P[30]=113
P[31]=127
P[32]=131
P[33]=137
P[34]=139
P[35]=149
P[36]=151
P[37]=157
P[38]=163
P[39]=167
P[40]=173
P[41]=179
P[42]=181
P[43]=191
P[44]=193
P[45]=197
P[46]=199
P[47]=211
P[48]=223
P[49]=227
P[50]=229
P[51]=233
P[52]=239
P[53]=241
P[54]=251
P[55]=257
P[56]=263
P[57]=269
P[58]=271
P[59]=277
P[60]=281
P[61]=283
P[62]=293
P[63]=307
P[64]=311
P[65]=313
P[66]=317
P[67]=331
P[68]=337
P[69]=347
P[70]=349
P[71]=353
P[72]=359
P[73]=367
P[74]=373
P[75]=379
P[76]=383
P[77]=389
P[78]=397
P[79]=401
P[80]=409
P[81]=419
P[82]=421
P[83]=431
P[84]=433
P[85]=439
P[86]=443
P[87]=449
P[88]=457
P[89]=461
P[90]=463
P[91]=467
P[92]=479
P[93]=487
P[94]=491
P[95]=499
P[96]=503
P[97]=509
P[98]=521
P[99]=523
P[100]=541
P[101]=547
P[102]=557
P[103]=563
P[104]=569
P[105]=571
P[106]=577
P[107]=587
P[108]=593
P[109]=599
P[110]=601
P[111]=607
P[112]=613
P[113]=617
P[114]=619
P[115]=631
P[116]=641
P[117]=643
P[118]=647
P[119]=653
P[120]=659
P[121]=661
P[122]=673
P[123]=677
P[124]=683
P[125]=691
P[126]=701
P[127]=709
P[128]=719
P[129]=727
P[130]=733
P[131]=739
P[132]=743
P[133]=751
P[134]=757
P[135]=761
P[136]=769
P[137]=773
P[138]=787
P[139]=797
P[140]=809
P[141]=811
P[142]=821
P[143]=823
P[144]=827
P[145]=829
P[146]=839
P[147]=853
P[148]=857
P[149]=859
P[150]=863
P[151]=877
P[152]=881
P[153]=883
P[154]=887
P[155]=907
P[156]=911
P[157]=919
P[158]=929
P[159]=937
P[160]=941
P[161]=947
P[162]=953
P[163]=967
P[164]=971
P[165]=977
P[166]=983
P[167]=991
P[168]=997
P[169]=1009
P[170]=1013
P[171]=1019
P[172]=1021
P[173]=1031
P[174]=1033
P[175]=1039
P[176]=1049
P[177]=1051
P[178]=1061
P[179]=1063
P[180]=1069
P[181]=1087
P[182]=1091
P[183]=1093
P[184]=1097
P[185]=1103
P[186]=1109
P[187]=1117
P[188]=1123
P[189]=1129
P[190]=1151
P[191]=1153
P[192]=1163
P[193]=1171
P[194]=1181
P[195]=1187
P[196]=1193
P[197]=1201
P[198]=1213
P[199]=1217
P[200]=1223
P[201]=1229
P[202]=1231
P[203]=1237
P[204]=1249
P[205]=1259
P[206]=1277
P[207]=1279
P[208]=1283
P[209]=1289
P[210]=1291
P[211]=1297
P[212]=1301
P[213]=1303
P[214]=1307
P[215]=1319
P[216]=1321
P[217]=1327
P[218]=1361
P[219]=1367
P[220]=1373
P[221]=1381
P[222]=1399
P[223]=1409
P[224]=1423
P[225]=1427
P[226]=1429
P[227]=1433
P[228]=1439
P[229]=1447
P[230]=1451
P[231]=1453
P[232]=1459
P[233]=1471
P[234]=1481
P[235]=1483
P[236]=1487
P[237]=1489
P[238]=1493
P[239]=1499
P[240]=1511
P[241]=1523
P[242]=1531
P[243]=1543
P[244]=1549
P[245]=1553
P[246]=1559
P[247]=1567
P[248]=1571
P[249]=1579
P[250]=1583
P[251]=1597
P[252]=1601
P[253]=1607
P[254]=1609
P[255]=1613
P[256]=1619
P[257]=1621
P[258]=1627
P[259]=1637
P[260]=1657
P[261]=1663
P[262]=1667
P[263]=1669
P[264]=1693
P[265]=1697
P[266]=1699
P[267]=1709
P[268]=1721
P[269]=1723
P[270]=1733
P[271]=1741
P[272]=1747
P[273]=1753
P[274]=1759
P[275]=1777
P[276]=1783
P[277]=1787
P[278]=1789
P[279]=1801
P[280]=1811
P[281]=1823
P[282]=1831
P[283]=1847
P[284]=1861
P[285]=1867
P[286]=1871
P[287]=1873
P[288]=1877
P[289]=1879
P[290]=1889
P[291]=1901
P[292]=1907
P[293]=1913
P[294]=1931
P[295]=1933
P[296]=1949
P[297]=1951
P[298]=1973
P[299]=1979
P[300]=1987
P[301]=1993
P[302]=1997
P[303]=1999
P[304]=2003
P[305]=2011
P[306]=2017
P[307]=2027
P[308]=2029
P[309]=2039
P[310]=2053
P[311]=2063
P[312]=2069
P[313]=2081
P[314]=2083
P[315]=2087
P[316]=2089
P[317]=2099
P[318]=2111
P[319]=2113
P[320]=2129
P[321]=2131
P[322]=2137
P[323]=2141
P[324]=2143
P[325]=2153
P[326]=2161
P[327]=2179
P[328]=2203
P[329]=2207
P[330]=2213
P[331]=2221
P[332]=2237
P[333]=2239
P[334]=2243
P[335]=2251
P[336]=2267
P[337]=2269
P[338]=2273
P[339]=2281
P[340]=2287
P[341]=2293
P[342]=2297
P[343]=2309
P[344]=2311
P[345]=2333
P[346]=2339
P[347]=2341
P[348]=2347
P[349]=2351
P[350]=2357
P[351]=2371
P[352]=2377
P[353]=2381
P[354]=2383
P[355]=2389
P[356]=2393
P[357]=2399
P[358]=2411
P[359]=2417
P[360]=2423
P[361]=2437
P[362]=2441
P[363]=2447
P[364]=2459
P[365]=2467
P[366]=2473
P[367]=2477
P[368]=2503
P[369]=2521
P[370]=2531
P[371]=2539
P[372]=2543
P[373]=2549
P[374]=2551
P[375]=2557
P[376]=2579
P[377]=2591
P[378]=2593
P[379]=2609
P[380]=2617
P[381]=2621
P[382]=2633
P[383]=2647
P[384]=2657
P[385]=2659
P[386]=2663
P[387]=2671
P[388]=2677
P[389]=2683
P[390]=2687
P[391]=2689
P[392]=2693
P[393]=2699
P[394]=2707
P[395]=2711
P[396]=2713
P[397]=2719
P[398]=2729
P[399]=2731
P[400]=2741
P[401]=2749
P[402]=2753
P[403]=2767
P[404]=2777
P[405]=2789
P[406]=2791
P[407]=2797
P[408]=2801
P[409]=2803
P[410]=2819
P[411]=2833
P[412]=2837
P[413]=2843
P[414]=2851
P[415]=2857
P[416]=2861
P[417]=2879
P[418]=2887
P[419]=2897
P[420]=2903
P[421]=2909
P[422]=2917
P[423]=2927
P[424]=2939
P[425]=2953
P[426]=2957
P[427]=2963
P[428]=2969
P[429]=2971
P[430]=2999
P[431]=3001
P[432]=3011
P[433]=3019
P[434]=3023
P[435]=3037
P[436]=3041
P[437]=3049
P[438]=3061
P[439]=3067
P[440]=3079
P[441]=3083
P[442]=3089
P[443]=3109
P[444]=3119
P[445]=3121
P[446]=3137
P[447]=3163
P[448]=3167
P[449]=3169
P[450]=3181
P[451]=3187
P[452]=3191
P[453]=3203
P[454]=3209
P[455]=3217
P[456]=3221
P[457]=3229
P[458]=3251
P[459]=3253
P[460]=3257
P[461]=3259
P[462]=3271
P[463]=3299
P[464]=3301
P[465]=3307
P[466]=3313
P[467]=3319
P[468]=3323
P[469]=3329
P[470]=3331
P[471]=3343
P[472]=3347
P[473]=3359
P[474]=3361
P[475]=3371
P[476]=3373
P[477]=3389
P[478]=3391
P[479]=3407
P[480]=3413
P[481]=3433
P[482]=3449
P[483]=3457
P[484]=3461
P[485]=3463
P[486]=3467
P[487]=3469
P[488]=3491
P[489]=3499
P[490]=3511
P[491]=3517
P[492]=3527
P[493]=3529
P[494]=3533
P[495]=3539
P[496]=3541
P[497]=3547
P[498]=3557
P[499]=3559
P[500]=3571
P[501]=3581
P[502]=3583
P[503]=3593
P[504]=3607
P[505]=3613
P[506]=3617
P[507]=3623
P[508]=3631
P[509]=3637
P[510]=3643
P[511]=3659
P[512]=3671
P[513]=3673
P[514]=3677
P[515]=3691
P[516]=3697
P[517]=3701
P[518]=3709
P[519]=3719
P[520]=3727
P[521]=3733
P[522]=3739
P[523]=3761
P[524]=3767
P[525]=3769
P[526]=3779
P[527]=3793
P[528]=3797
P[529]=3803
P[530]=3821
P[531]=3823
P[532]=3833
P[533]=3847
P[534]=3851
P[535]=3853
P[536]=3863
P[537]=3877
P[538]=3881
P[539]=3889
P[540]=3907
P[541]=3911
P[542]=3917
P[543]=3919
P[544]=3923
P[545]=3929
P[546]=3931
P[547]=3943
P[548]=3947
P[549]=3967
P[550]=3989
P[551]=4001
P[552]=4003
P[553]=4007
P[554]=4013
P[555]=4019
P[556]=4021
P[557]=4027
P[558]=4049
P[559]=4051
P[560]=4057
P[561]=4073
P[562]=4079
P[563]=4091
P[564]=4093
P[565]=4099
P[566]=4111
P[567]=4127
P[568]=4129
P[569]=4133
P[570]=4139
P[571]=4153
P[572]=4157
P[573]=4159
P[574]=4177
P[575]=4201
P[576]=4211
P[577]=4217
P[578]=4219
P[579]=4229
P[580]=4231
P[581]=4241
P[582]=4243
P[583]=4253
P[584]=4259
P[585]=4261
P[586]=4271
P[587]=4273
P[588]=4283
P[589]=4289
P[590]=4297
P[591]=4327
P[592]=4337
P[593]=4339
P[594]=4349
P[595]=4357
P[596]=4363
P[597]=4373
P[598]=4391
P[599]=4397
P[600]=4409
P[601]=4421
P[602]=4423
P[603]=4441
P[604]=4447
P[605]=4451
P[606]=4457
P[607]=4463
P[608]=4481
P[609]=4483
P[610]=4493
P[611]=4507
P[612]=4513
P[613]=4517
P[614]=4519
P[615]=4523
P[616]=4547
P[617]=4549
P[618]=4561
P[619]=4567
P[620]=4583
P[621]=4591
P[622]=4597
P[623]=4603
P[624]=4621
P[625]=4637
P[626]=4639
P[627]=4643
P[628]=4649
P[629]=4651
P[630]=4657
P[631]=4663
P[632]=4673
P[633]=4679
P[634]=4691
P[635]=4703
P[636]=4721
P[637]=4723
P[638]=4729
P[639]=4733
P[640]=4751
P[641]=4759
P[642]=4783
P[643]=4787
P[644]=4789
P[645]=4793
P[646]=4799
P[647]=4801
P[648]=4813
P[649]=4817
P[650]=4831
P[651]=4861
P[652]=4871
P[653]=4877
P[654]=4889
P[655]=4903
P[656]=4909
P[657]=4919
P[658]=4931
P[659]=4933
P[660]=4937
P[661]=4943
P[662]=4951
P[663]=4957
P[664]=4967
P[665]=4969
P[666]=4973
P[667]=4987
P[668]=4993
P[669]=4999
P[670]=5003
P[671]=5009
P[672]=5011
P[673]=5021
P[674]=5023
P[675]=5039
P[676]=5051
P[677]=5059
P[678]=5077
P[679]=5081
P[680]=5087
P[681]=5099
P[682]=5101
P[683]=5107
P[684]=5113
P[685]=5119
P[686]=5147
P[687]=5153
P[688]=5167
P[689]=5171
P[690]=5179
P[691]=5189
P[692]=5197
P[693]=5209
P[694]=5227
P[695]=5231
P[696]=5233
P[697]=5237
P[698]=5261
P[699]=5273
P[700]=5279
P[701]=5281
P[702]=5297
P[703]=5303
P[704]=5309
P[705]=5323
P[706]=5333
P[707]=5347
P[708]=5351
P[709]=5381
P[710]=5387
P[711]=5393
P[712]=5399
P[713]=5407
P[714]=5413
P[715]=5417
P[716]=5419
P[717]=5431
P[718]=5437
P[719]=5441
P[720]=5443
P[721]=5449
P[722]=5471
P[723]=5477
P[724]=5479
P[725]=5483
P[726]=5501
P[727]=5503
P[728]=5507
P[729]=5519
P[730]=5521
P[731]=5527
P[732]=5531
P[733]=5557
P[734]=5563
P[735]=5569
P[736]=5573
P[737]=5581
P[738]=5591
P[739]=5623
P[740]=5639
P[741]=5641
P[742]=5647
P[743]=5651
P[744]=5653
P[745]=5657
P[746]=5659
P[747]=5669
P[748]=5683
P[749]=5689
P[750]=5693
P[751]=5701
P[752]=5711
P[753]=5717
P[754]=5737
P[755]=5741
P[756]=5743
P[757]=5749
P[758]=5779
P[759]=5783
P[760]=5791
P[761]=5801
P[762]=5807
P[763]=5813
P[764]=5821
P[765]=5827
P[766]=5839
P[767]=5843
P[768]=5849
P[769]=5851
P[770]=5857
P[771]=5861
P[772]=5867
P[773]=5869
P[774]=5879
P[775]=5881
P[776]=5897
P[777]=5903
P[778]=5923
P[779]=5927
P[780]=5939
P[781]=5953
P[782]=5981
P[783]=5987
NOTE: DATA statement used (Total process time):
      real time           1:09.87
      cpu time            1:06.65




%LET TOTAL = 1000000; * Limits the number of prime numbers generated ;
%LET DIM = 1000000; * The size of the sieve arrays used ;
%LET TIME = 2400; * Time limit in seconds ;

DATA _NULL_;
   ARRAY P{&DIM} _TEMPORARY_; * Prime numbers;
   ARRAY M{&DIM} _TEMPORARY_; * Multiples of prime numbers;
   TIMEOUT = DATETIME() + &TIME; * Time limit;
   *FILE PRINT NOTITLES;
   SQUARE = 4;

   DO X = 2 TO CONSTANT('EXACTINT'); * Is X prime? ;
      IF DATETIME() >= TIMEOUT THEN STOP; * Time limit reached ;
      IF X = SQUARE THEN DO; * Extend sieve;
         IMAX + 1;
         IF IMAX >= &DIM THEN leave; * Sieve size limit reached. ;
         SQUARE = M{IMAX + 1}; 
         CONTINUE;
         END;
      * Find least prime factor (LPF). ;
      LPF = 0;
      DO I = 1 TO IMAX UNTIL (LPF);
         DO WHILE (M{I} < X); * Update sieve with new multiple. ;
            M{I} + P{I};
            END;
         IF M{I} = X THEN LPF = P{I};
         END;
      IF LPF THEN CONTINUE; * Composite number found. ;
      *PUT X @; * Write prime number in output. ;
      N + 1;
      IF N >= &TOTAL THEN leave; * Output maximum reached. ;
      ELSE IF N <= &DIM THEN DO; * Add prime number to sieve. ;
         P{N} = X;
         M{N} = X*X;
         END;
      END;






*get 781 consecutive prim number;
 do _m=1 to &dim -780;
  prim_sum=0;
  do _n=_m to _m+780;
  prim_sum+p{_n};
  end;
 
  if prim_sum in p then do;
   put prim_sum= ;
   do _n=_m to _m+780;
    put p{_n}=;
   end;          stop;          
                 end; 

end;
 run;

Ksharp

Another math puzzle

Peter_C — Sun, 09 Oct 2011 18:16:00 GMT

I assume a list of the first 100000+ primes can be found on internet

Another math puzzle

Ksharp — Mon, 10 Oct 2011 02:33:25 GMT

Peter.

But if want 781, then at least need to generate one million prime numbers.

Another math puzzle

DF — Mon, 10 Oct 2011 13:01:58 GMT

You can download the first 50 million here: http://primes.utm.edu/lists/small/millions/ .

However, as the author points out it should be possible to generate these yourself faster than downloading them. Certainly if you were writing in C++ etc. it would be.

Another math puzzle

Peter_C — Mon, 10 Oct 2011 15:42:02 GMT

if primes are important to you then downloading the collection once might be worth the time

Another math puzzle

Peter_C — Mon, 10 Oct 2011 15:45:09 GMT

downloading the first two million primes took less than 10 secs. I might take longer to find the c-compiler :smileyconfused:

Another math puzzle

FriedEgg — Mon, 10 Oct 2011 21:42:09 GMT

Hi,

I have spent a minute on another method to calculate prime numbers, it could gain effeciency by have a better method to reduce loops for higher numbers. You could write a method in C using proc proto, I have begun on something.

data test;

i=2; output; i=3; output; i=5;

do while(i<10**7);

n=0;

j=3;

do while(j<=int(sqrt(i)) and n=0);

if mod(i,j)=0 then n+1;

j+2;

end;

if n=0 then output;

i+2;

end;

run;

I avoid even numbers because the only even prime number is 2 and the products of even numbers are also even so I avoid testing even numbers as factors. With the code above I find 664,579 prime numbers in a little under a minute on a server with load. I confirm output set against the dataset from website discussed here and have 100% match. As the number of tests increases the time per test increases logaritmically. Better logic for reducing the number of factors to test would dramatically improve performance for testing to more prime numbers.

Another math puzzle

FriedEgg — Mon, 10 Oct 2011 22:32:31 GMT

Here is a macro to generate prime numbers also.

%macro get_primes(from, to);

data pnbrs;

%if &from<=2 %then %do; i=2; output; %end;

do i=&from to &to;

if mod(i,2) ne 0 then output;

end;

run;

sasfile work.pnbrs open;

%let nbr=&from;

data pnbrs;

modify pnbrs;

if i>&nbr and mod(i,&nbr)=0 then remove;

run;

%do %while(%eval(&nbr**2)<=%eval(&to));

proc sql noprint;

select min(i) into :nbr from pnbrs where i>&nbr;

quit;

data pnbrs;

modify pnbrs;

if i>&nbr and mod(i,&nbr)=0 then remove;

run;

%end;

sasfile work.pnbrs close;

%mend;

%get_primes(2,10**7);

Another math puzzle

Ksharp — Tue, 11 Oct 2011 05:03:46 GMT

FriedEgg

However. I do not think your method is as far efficient as Rick Wicklin 's algorithm.

Another math puzzle

FriedEgg — Tue, 11 Oct 2011 05:19:56 GMT

No it is not very efficient, in my opinion there are a number of clear and obvious flaws. That being said, his method uses SAS/IML, which I do not have...

Here is a pretty fast version, written in C that I will try later to implement similar functionality via proc proto:

http://www.fpx.de/fp/Software/sieve.c

#include <stdio.h>

#include <malloc.h>

#include <time.h>

#define TEST(f,x) (*(f+(x)/16)&(1<<(((x)%16L)/2)))

#define SET(f,x) *(f+(x)/16)|=1<<(((x)%16L)/2)

int

main(int argc, char *argv[])

{

unsigned char *feld=NULL, *zzz;

unsigned long teste=1, max, mom, hits=1, count, alloc, s=0, e=1;

time_t begin;

if (argc > 1)

max = atol (argv[1]) + 10000;

else

max = 14010000L;

while (feld==NULL)

zzz = feld = malloc (alloc=(((max-=10000L)>>4)+1L));

for (count=0; count<alloc; count++) *zzz++ = 0x00;

printf ("Searching prime numbers to : %ld\n", max);

begin = time (NULL);

while ((teste+=2) < max)

if (!TEST(feld, teste)) {

if (++hits%2000L==0) {printf (" %ld. prime number\x0d", hits); fflush(stdout);}

for (mom=3L*teste; mom<max; mom+=teste<<1) SET (feld, mom);

}

printf (" %ld prime numbers foundn %ld secs.\n\nShow prime numbers",

hits, time(NULL)-begin);

while (s<e) {

printf ("\n\nStart of Area : "); fflush (stdout); scanf ("%ld", &s);

printf ("End of Area : "); fflush (stdout); scanf ("%ld", &e);

count=s-2; if (s%2==0) count++;

while ((count+=2)<e) if (!TEST(feld,count)) printf ("%ld\t", count);

}

free (feld);

return 0;

}