;-------------------------------------------------------------
;+
; NAME:
;     PRIME
; PURPOSE:
;     Return an array with the specified number of prime numbers.
; EXPLANATATION:
;     This procedure is similar to PRIMES in the standard IDL distribution,
;     but stores results in a common block, and so is much faster 
;
; CALLING SEQUENCE:
;       p = prime(n)
; INPUTS:
;       n = desired number of primes, scalar positive integer
; OUTPUTS:
;       p = resulting array of primes, vector of positive integers
; COMMON BLOCKS:
;       prime_com
; NOTES:
;       Note: Primes that have been found in previous calls are
;         remembered and are not regenerated.
; MODIFICATION HISTORY:
;       R. Sterner  17 Oct, 1985.
;       R. Sterner,  5 Feb, 1993 --- fixed a bug that missed a few primes.
;       Converted to IDL V5          March 1999
;
; Copyright (C) 1985, Johns Hopkins University/Applied Physics Laboratory
; This software may be used, copied, or redistributed as long as it is not
; sold and this copyright notice is reproduced on each copy made.  This
; routine is provided as is without any express or implied warranties
; whatsoever.  Other limitations apply as described in the file disclaimer.txt.
;-
;-------------------------------------------------------------
 
	function prime,n, help=hlp
 
	common prime_com, max, pmax
 
	if (n_params(0) lt 1) or keyword_set(hlp) then begin
	  print,' Return an array with the specified number of prime numbers.'
	  print,' p = prime(n)'
	  print,'   n = desired number of primes.    in'
	  print,'   p = resulting array of primes.   out'
	  print,' Note: Primes that have been found in previous calls are'
	  print,'   remembered and are not regenerated.'
	  return, -1
	endif
 
	if n_elements(max) eq 0 then max = 0	; Make MAX defined.
	if n le max then return, pmax[0:n-1]	; Enough primes in memory.
	p = lonarr(n)				; Need to find primes.
	if max eq 0 then begin			; Have none now. Start with 8.
	  p[0] = [2,3,5,7,11,13,17,19]
	  if n le 8 then return, p[0:n-1]	; Need 8 or less.
	  i = 8					; Need more than 8.
	  t = 19L				; Search start value.
	endif else begin			; Start with old primes.
	  p[0] = pmax				; Move old primes into big arr.
	  i = max				; Current prime count.
	  t = p[max-1]				; Biggest prime so far.
	endelse
 
loop:	if i eq n then begin			; Have enough primes.
	  max = n				; Count.
	  pmax = p				; Array of primes.
	  return, p				; Return primes.
	endif
loop2:	t = t + 2				; Next test value, t.
	it = 1					; Start testing with 1st prime.
loop3:	pr = p[it]				; Pick next test prime.
	pr2 = pr*pr				; Square it.
	if pr2 gt t then begin			; Selected prime > sqrt(t)?
	  i = i + 1				; Yes, count
	  p[i-1] = t				; and store new prime.
	  goto, loop				; Go check if done.
	endif
	if pr2 eq t then goto, loop2		; Test number, t, was a square.
	if (t mod pr) eq 0 then goto, loop2	; Curr prime divides t.
	it = it + 1				; Check next prime.
	goto, loop3
	end