topic Re: Do loop and mod fuction in SAS Programming

Do loop and mod fuction

ChuksManuel — Fri, 24 May 2019 00:32:07 GMT

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

Re: Do loop and mod fuction

heffo — Fri, 24 May 2019 00:48:07 GMT

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

Re: Do loop and mod fuction

Astounding — Fri, 24 May 2019 01:06:44 GMT

You overlooked the word NOT:

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

Re: Do loop and mod fuction

ChuksManuel — Fri, 24 May 2019 01:09:03 GMT

Thanks.

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

Re: Do loop and mod fuction

ChuksManuel — Fri, 24 May 2019 01:09:50 GMT

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

Any help is appreciated.

Re: Do loop and mod fuction

Patrick — Fri, 24 May 2019 01:51:59 GMT

@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;

Re: Do loop and mod fuction

ChuksManuel — Fri, 24 May 2019 01:46:44 GMT

Thank you that worked.

Re: Do loop and mod fuction

Tom — Fri, 24 May 2019 02:08:08 GMT

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.

Re: Do loop and mod fuction

Patrick — Fri, 24 May 2019 03:54:25 GMT

@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.

Re: Do loop and mod fuction

FreelanceReinh — Fri, 24 May 2019 08:31:02 GMT

@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);

Re: Do loop and mod fuction

Patrick — Fri, 24 May 2019 12:27:06 GMT

@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?

Re: Do loop and mod fuction

Tom — Fri, 24 May 2019 12:52:56 GMT

@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.

Re: Do loop and mod fuction

Patrick — Fri, 24 May 2019 13:07:48 GMT

@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!

Re: Do loop and mod fuction

Ksharp — Fri, 24 May 2019 13:12:13 GMT

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;

Re: Do loop and mod fuction

FreelanceReinh — Fri, 24 May 2019 14:12:26 GMT

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.

Re: Do loop and mod fuction

FreelanceReinh — Fri, 24 May 2019 18:42:55 GMT

Thanks @Ksharp for chiming in. Yes, given that it's one of the oldest algorithms (Eratosthenes, 3rd century B.C.), it's stunningly fast.

I think your current implementation is somewhat below its best, though: The idea of setting every x{i}-th array element to missing could be implemented using a DO ... BY statement rather than calling the MOD function to identify the elements to be canceled. This change reduces the run time for the primes up to 1 million from about 2 to 0.1 s (on my computer).

The version below takes only 6 seconds for the 5,761,455 primes up to 100 million.

%let max=100000000;

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

do i=2 to int(sqrt(&max));
  if not missing(x{i}) then do;
     do j=i*i to dim(x) by x{i};
       x{j}=.;
     end;
  end;
end;

do k=2 to &max;
  if not missing(x{k}) then do; prime_number=x{k}; output; end;
end;
keep prime_number;
run;

I was tempted to replace the first DO loop by specifying initial values for the array, (1:&max), but that turned out to be significantly slower.

Edit: By skipping the even numbers >2 the run time can be further reduced to <5.0 seconds.

data want;
array x{&max} _temporary_;
do i=2, 3 to &max by 2;
  x{i}=i;
end;

do i=3 to int(sqrt(&max)) by 2;
  if not missing(x{i}) then do;
     do j=i*i to dim(x) by x{i};
       x{j}=.;
     end;
  end;
end;

do k=2, 3 to &max by 2;
  if not missing(x{k}) then do; prime_number=x{k}; output; end;
end;
keep prime_number;
run;