素数的sas实现

1.DATA STEP

data tmp(drop=i);
 n=2;
 output;
 do n=3 to 10000 by 2;
  do i=2 to n-1;
   if mod(n,i)=0 and i^=n-1 then leave;
   if i=n-1 then output;
  end;
 end;
run;
proc print;run;

2.macro from SAS_L

option mprint;
%macro prime(n);
data prime;
do i=1 to &n;
    prime=1;
    do j=2 to ceil(i/2);
        if mod(i,j)=0 then prime=0;
    end;
    if prime=1 then output;
end;
run;
%mend;
%prime(10000);

3. Rick Aster

%LET TOTAL = 2000; * Limits the number of prime numbers generated ;
%LET DIM = 1000; * The size of the sieve arrays used ;
%LET TIME = 30; * 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 10000; * Is X prime? ;
      IF DATETIME() >= TIMEOUT THEN STOP; * Time limit reached ;
      IF X = SQUARE THEN DO; * Extend sieve;
         IMAX + 1;
         IF IMAX >= &DIM THEN STOP; * 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 STOP; * Output maximum reached. ;
      ELSE IF N <= &DIM THEN DO; * Add prime number to sieve. ;
         P{N} = X;
         M{N} = X*X;
         END;
      END;
   STOP;
RUN;

4. Rick Wicklin 's algorithm

/** The Sieve of Eratosthenes **/
proc iml;
max = 10000; 
list = 2:max;                 /** find prime numbers in [2, max] **/
primes = j(1, ncol(list), .); /** allocate space to store primes **/

numPrimes = 0;
p = list[1];                   /** 2 is the first prime **/
do while (p##2 <= max);        /** search through sqrt(max) **/
   idx = loc( mod(list, p)=0 );/** find all numbers divisible by p **/
   list = remove(list, idx);   /** remove them. They are not prime. **/
   numPrimes = numPrimes + 1;
   primes[numPrimes] = p;      /** include p in the list of primes **/
   p = list[1];                /** next prime number to sift **/
end;

/** include remaining numbers greater than sqrt(max) **/
k = numPrimes;
n = k + ncol(list);
primes[k+1:n] = list;          /** put them at the end of the list **/
primes = primes[ ,1:n];        /** get rid of excess storage space **/

print primes;
quit;

5.discussions

http://communities.sas.com/message/107293