Evil Science A whole load of stuff

31Mar/110

ASieve of Atkin prime number generator using TSQL

I worked out an implementation of the Sieve of Atkin for use in solving Project Euler problems. It's probably not as efficient as it could be, but it does the trick for me.

set nocount on

drop table #atkin
create table #atkin
(
	num float
	, isPrime bit
)
create clustered index primenums
	on #atkin(num)

declare @limit as bigint
declare @sqrt as float
declare @x as int
declare @y as int
declare @n as int

set @limit = 1000000
set @sqrt = sqrt(@limit)
set @x = 1

print 'begin outer loop'

/*
	put in candidate primes:
	integers which have an odd number of representations by certain quadratic forms
*/

while (@x <= @sqrt)
	begin

		set @y = 1

		while (@y <=@sqrt)
			begin

				set @n = (4 * (@x * @x)) + (@y * @y)
				if ((select num from #atkin where num =@n) is null)
					insert into #atkin (num, isPrime) values (@n,0)

				if (@n <= @limit AND (@n % 12 = 1 or @n % 12 =5))
					update #atkin set isPrime = isPrime ^ 1 where num = @n

				set @n = (3 * (@x * @x)) + (@y * @y)
				if ((select num from #atkin where num =@n) is null)
					insert into #atkin (num, isPrime) values (@n,0)

				if (@n <= @limit AND @n % 12 = 7)
					update #atkin set isPrime = isPrime ^ 1 where num = @n

				set @n = (3 * (@x * @x)) - (@y * @y)
				if ((select num from #atkin where num =@n) is null)
					insert into #atkin (num, isPrime) values (@n,0)

				if (@x >= @y AND @n <= @limit  AND @n % 12 = 11)
					update #atkin set isPrime = isPrime ^ 1 where num = @n

				set @y = @y + 1

			end

		set @x = @x + 1

	end 

declare @nsqr as bigint
declare @k as bigint

/*
	eliminate composites by sieving
*/
set @n = 5
while ( @n <= @sqrt)

	begin

		print @n

		if ((select isprime from #atkin where num = @n)=1)
			begin

				set @nsqr = @n * @n

				set @k = @nsqr

				while (@k <= @limit)

					begin

						update #atkin set isPrime = 0 where num = @k

						set @k = @k + @nsqr
						print '@k increment ' + str(@k)

					end

			end

		set @n = @n + 1

		print '@n increment ' + str(@n)

	end

	insert into #atkin (num, isPrime) values (2,1)
	insert into #atkin (num, isPrime) values (3,1)

delete from	#atkin where isprime <> 1