/*I am interested in two functions: getPrimes(n,t) – where n is a positive integer and t can be one of two values. If it is 1 getPrimes returns the first n primes. If it is 2 getPrimes will return all primes less than or equal to n. You can assume the value of n will not exceed 1e8. isPrime(vec) – where vec is a nRow row, 1 column matrix of integer values no greater than 4e14. It is to return a matrix of 1’s and 0’s where vec[i] = 1 indicates the ith element of vec is prime and 0 it is not prime. The input test vector will contain between 100 and 1000 elements (but the function should be able to accepts larger or smaller vectors). Examples: getPrimes(10,1) would return [2,3,5,7,11,13,17,19,23,29] getPrimes(10,2) would return [2,3,5,7] isPrime([-1,0,2,4,5]) would return [0,0,1,0,1]*/ //Sieve of Atkin function MySieve = Function( {limit}, // limit = 100; sqrt_limit = SQRT(limit) + 1; sieve = 1 :: limit; sieve[1 :: limit] = 0; For(j = 1, j < sqrt_limit, j++, For( k = 1, k < sqrt_limit, k++, l = (4 * j * j) + (k * k); if(l <= limit & ( Modulo(l,12) == 1 | Modulo(l,12) == 5 ), if(sieve[l] == 0, sieve[l] = l, sieve[l] = 0)); l = (3 * j * j) + (k * k); if(l < limit & Modulo(l,12) == 7, if(sieve[l] == 0, sieve[l] = l, sieve[l] = 0)); if( j > k, l = (3 * j * j) - (k * k); if( l < limit & Modulo(l,12) == 11, if(sieve[l] == 0, sieve[l] = l, sieve[l] = 0)); ); ); ); if( limit >= 2, sieve[2] = 2); if( limit >= 3, sieve[3] = 3); For(r = 5, r <= sqrt_limit, r++, if(sieve[r] == r, q = r * r; For(i = q, i < limit, i += q, sieve[i] = 0; ); ); ); v = Loc( sieve > 0 ); Primes = sieve[v]; ); ////getPrimes(n,t) getPrimes = Function( {n,t}, Match( t, 1, limit = Ceiling( 3.067 * (n+5) ^ 1.1106 ), 2, limit = n); Ans = MySieve(limit); if(t == 1, Ans = Ans[1::n], Ans = Ans); ); ////isPrime(vec) isPrime = Function( {vec}, rows = nRows(Vec); max = Max(Vec); top = Sqrt(max) + 1; Pvec = vec; neg = Loc( Pvec < 2 ); Pvec[neg] = 0; low = Loc( Pvec == 2 | Pvec == 3 | Pvec == 5 ); Pvec[low] = 1; etc = Loc( Modulo(Pvec,2) == 0 | Modulo(Pvec,3) == 0 | Modulo(Pvec,5) == 0 | Modulo(Pvec,7) == 0 ); Pvec[etc] = 0; Key = MySieve(top); Krows = nRows(Key); For(i = 1, i <= rows, i++, m = 1; if(Pvec[i] > 1, For(m = 1, m <= Krows & Pvec[i] > Key[m], m++, if(Modulo(Pvec[i],Key[m]) == 0, Pvec[i] = 0, Pvec[i] = Pvec[i]); if(Pvec[i] == 0, Break()); ); ); if(m >= Krows & Pvec[i] > 1, Pvec[i] = 1); ); rem = Loc( Pvec > 1 ); Pvec[rem] = 1; PrimeVec = Pvec; );