Do loop and mod fuction

[ChuksManuel]

Hello,

I've been trying to write a code to answer this practice question in SAS but i haven't had success.

Use a do loop, the mod function and an if-statement to put out all prime numbers between 4 and 300.  Consider a prime number to be any number greater than 3 that is not divisible by 2, 3, 5 or 7.    

Here are my codes

data one;
Do K=4 to 300;
if mod(k,2)=0 then output;
If mod(k,3)=0 then output;
if mod(k,5)=0 then output;
if mod(k,7)=0 then output;
End;
Proc print; run;

I had a problem in the results. Please can anyone review these to let me know how best i can modify or get the prime numbers between 4 and 300.

Thanks

[heffo]

What problem are you having? It works for me, except that you get multiple hits on the same number, but you just have to add else to the if cases. 

data want;
	Do K=4 to 300;
		if mod(k,2)=0 then output;
		else If mod(k,3)=0 then output;
		else if mod(k,5)=0 then output;
		else if mod(k,7)=0 then output;
	End;
run;

Proc print data=want;
run;

If it is actually prime number (and not just as an exercise) that you are after, I found this code: http://www.globalstatements.com/shortcuts/88a.html 

[ChuksManuel]

Thanks.

I was trying to get the prime numbers between 4 and 300. The code wasn't giving me that.

[Astounding]

You overlooked the word NOT:

"not divisible by 2, 3, 5, or 7."

[ChuksManuel]

That's the part i've been struggling with.

Any help is appreciated.

[Patrick · Accepted solution]

@ChuksManuel 

Here one way to get the logic right. You could now try to also write a second loop for your list of prime numbers instead of spelling out everything using AND conditions.

data want;
  do K=4 to 300;
    if mod(k,2) ne 0 
      and mod(k,3) ne 0 
      and mod(k,5) ne 0 
      and mod(k,7) ne 0 then
      output;
  end;
run;

proc print data=want;
run;

And here a more general approach which returns actual prime numbers up to the &upper_bound you define.

%let upper_bound=300;
data _null_;
  length prim_num 8;
  call missing(prim_num);
  dcl hash h1(ordered:'y');
  dcl hiter iter1('h1');
  h1.defineKey('prim_num');
  h1.defineData('prim_num');
  h1.defineDone();

  h1.add(key:2,data:2);
  h1.add(key:3,data:3);

  /* loop only over odd numbers */
  do k=5 to &upper_bound by 2;
    prim_flg='1';
    rc=iter1.first();
    do while(rc eq 0);
      if mod(k,prim_num)=0 then 
        do;
          /* K not a prime number: Try next value of K */
          prim_flg='0';
          leave;
        end;
      rc=iter1.next();
    end;
    /* add new prime number to hash */
    if prim_flg='1' then h1.add(key:k,data:k);
  end;

  /* write list of found prime numbers to data set WANT */
  h1.output(dataset:'want');
  stop;
run;

proc print data=want;
run;

[ChuksManuel]

Thank you that worked.

[Tom]

Hash is both more complex and takes more storage than as simple array would.

%let max_primes=1000;
%let max_number=300;
data want;
  array primes (&max_primes) _temporary_;
  do prime=2,3 by 2 to &max_number until(last>=dim(primes));
    notprime=0;
    do i=1 to last while(not notprime);
      notprime= not mod(prime,primes[i]) ;
    end;
    if not notprime then do;
      last+1;
      primes[last]=prime;
      output;
      put prime @;
    end;
  end;
  keep prime;
run;

2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101
103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193
197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293

NOTE: The data set WORK.WANT has 62 observations and 1 variables.

[Patrick]

@Tom 

"Hash is both more complex and takes more storage than as simple array would"

More complex: In my opinion not really once you get the hang of it.

Unlike an array is flexible and doesn't need the elements predefined so also not sure that the memory argument is really true once the numbers get bigger. So yes, the hash requires some initial memory just to define the structure but once done it then only requires more memory for actually identified primary numbers. 

[FreelanceReinh]

@Patrick: Nice application of the hash object.

A minor change improves the performance substantially (for larger values of upper_bound😞

do while(rc eq 0 & prim_num**2<=k);

Run times on my computer: 203 s --> 1 s for upper_bound=1000000.

Similarly for @Tom's solution, which is then even faster (run time 0.5 s in the above example):

do i=1 to last while(not notprime & primes[i]**2<=prime);

[Patrick]

@FreelanceReinh 

Nice!

I get it that we don't have to test with all already found primary numbers. I don't get it (mathematically) why it's prim_num**2<=k . Can you please enlighten me? 

[Tom]

---

@Patrick wrote:

@FreelanceReinh 

Nice!

I get it that we don't have to test with all already found primary numbers. I don't get it (mathematically) why it's prim_num**2<=k . Can you please enlighten me? 

---

You do not need to test whether the current number can be divided by a prime that is larger than its square root.  If said prime could divide the number being tested the result would have been a prime number less than the square root.  But you already tested all of those primes.

[Patrick]

@Tom 

I do understand the consequences of the formula BUT I don't understand its reason. 

"If said prime could divide the number being tested the result would have been a prime number less than the square root."

Oh! Yes, of course. I had to read this like 10 times but I got it now. Thanks!

[FreelanceReinh]

Thanks @Tom. More precisely: "... the result would have been divisible by a prime number less than the square root. (...)".

Example: sqrt(6853)=82.7..., hence it's sufficient to divide by the primes up to 79 (next prime 83 is already greater than that square root). We don't need to find out that mod(6853, 89)=0 in order to prove that 6853 is not prime, because 6853/89=77 must necessarily have prime factors <=77<=sqrt(6853), indeed: 7 and 11. So, we got this information already when we saw mod(6853, 7)=0.

@Patrick: My first idea was to write prim_num<=sqrt(k), but I figured that the equivalent criterion prim_num**2<=k (or even prim_num*prim_num<=k) can be evaluated slightly faster as it involves only integer multiplication rather than the computationally more expensive SQRT function.

[Ksharp]

The fastest algorithm is Prime Number Sieve :

data want;
array x{300} _temporary_;
do i=1 to 300;
  x{i}=i;
end;

do i=2 to int(sqrt(300));
  if not missing(x{i}) then do;
     do j=i+1 to dim(x);
       if not missing(x{j}) then do; if mod(x{j},x{i})=0 then call missing(x{j});end; 
     end;
  end;
end;

do k=4 to 300;
  if not missing(x{k}) then do;prime_number=x{k};output;end;
end;
keep prime_number;
run;