## 素数的sas实现

2018年10月22日 ⁄ 综合 ⁄ 共 1859字 ⁄ 字号 评论关闭

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);
```
```%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