Mandelbrat Generator in Python

Here's some code I wrote a while back. As the title suggests, it generates the mandelbrat set using some methods from the excellent NumPy library. I've included a few photos to give you an idea of its "capabilities".

# This code generates the Mandelbrot series fractal.
# Author: Cooper Bell


import math
import numpy
import matplotlib.pyplot as plt

xmin, xmax = -2, 1		# extent of x range
ymin, ymax = -1, 1		# extent of y range
max_iteration = 200		# max number of iteration for unbound series
delta = .001			# increment for points
inner = 512				# color inside bulb
mand = []				# points in the mandelbrot set

x0 = [i for i in numpy.arange(xmin, xmax, delta)]
y0 = [i for i in numpy.arange(ymin, ymax, delta)]

for y in y0:
	row = []
	for x in x0:
	
		z, c = 0, x + 1j*y		# iteration for determining set
		iteration = 0			# initialize iteration number
		
		p = math.sqrt((x - 0.25)**2 + y**2)
		
		if (x < p - 2*p**2 + 0.25):			  # checks if point is in the cardioid
			color = inner
			
		elif (x + 1)**2 + y**2 < float(1/16): # checks if point is in 2-period bulb
			color = inner
			
		else:
			while (numpy.abs(z) < 2) and (iteration < max_iteration):
		
				z = z*z + c
				iteration += 1
			
				if numpy.abs(z) >= 2:	# check if z is not set member
					color = (float(iteration)/max_iteration)*512
				else:
					color = inner
				
		row.append(color)
	mand.append(row)

plt.figure()
plt.imshow(numpy.array(mand), cmap='flag', extent=(xmin, xmax, ymin, ymax))
plt.axis('off')
plt.show()

The Sieve of Atkin in Common Lisp [UPDATE]

So I was able to retool the functions used in the previous code using iteration instead of recursion in order to evaluate them. Now I'm able to generate as many primes as I need, however it did take it a few seconds to calculate all of the primes under 1,000,000. I'm not sure if it will ever be able to hang with C++ when it comes to speed, but it seems to perform better than python anyway.

;;;; Generates a list of primes using the "Sieve of Atkin."
;;;; More on this algorithm can be found at:
;;;; http://en.wikipedia.org/wiki/Sieve_of_Atkin
;;;; Author: Cooper Bell

(setf *limit* 50000)    ; limit for maximum prime in list


;;; generate list of "nil" of size lim
(defun nullness (lim)
  (let ((lst '()))
    (do ((x 1 (+ x 1)))
        ((> x lim) lst)
      (setf lst (cons nil lst)))))


;;; sieve the primes
(defun sieve (lim)
  (setf is_prime (nullness lim))
  (do ((x 1 (+ x 1)))
      ((> x (+ (round (sqrt lim)) 1)) is_prime)
    (do ((y 1 (+ y 1)))
        ((> y (+ (round (sqrt lim)) 1)) is_prime)

      (setf n (+ (* 4 x x) (* y y)))
      (if (and (<= n lim) (or (eql (mod n 12) 1) (eql (mod n 12) 5)))
        (setf (nth n is_prime) (not (nth n is_prime))))

      (setf n (+ (* 3 x x) (* y y)))
      (if (and (<= n lim) (eql (mod n 12) 7))
        (setf (nth n is_prime) (not (nth n is_prime))))

      (setf n (- (* 3 x x) (* y y)))
      (if (and (> x y) (<= n lim) (eql (mod n 12) 11))
        (setf (nth n is_prime) (not (nth n is_prime)))))))


;;; set the square multiples of primes to false
(defun square-sieve (lst lim)
  (do ((n 0 (+ n 1)))
      ((> n (+ (round (sqrt lim))) 1) lst)
    (if (eql (nth n lst) 't)
      (do ((k (* n n) (+ k (* n n))))
          ((> k (- lim 1)) lst)
        (setf (nth k lst) nil)))))

;;; print the primes
(defun prime-print (primes)
  (let ((lst '())
        (max-length (list-length primes)))
    (do ((n 1 (+ n 1)))
        ((> n max-length) (cons 2 (cons 3 (reverse lst))))
      (if (eql (nth n primes) 't)
        (setf lst (cons n lst))))))


;;; produce the list of primes (better formatting)
(defun prime-list (lim)
  (prime-print (square-sieve (sieve lim) lim)))


(print (prime-list *limit*))

The Sieve of Atkin in Common Lisp

So since I've been learning a bit of common lisp lately, I've decided to implement one of my favorite go-to algorithms for generating a list or primes: the sieve of Atkin. I won't go to in depth on the underlying mathematics of the algorithm (Wikipedia does a pretty good job as it is), but it basically sieves primes based on the solutions to a few specific quadratic forms. As of right now, the code doesn't work as well as I'd like. The speed seems fine, but generating a list of primes over 5000 gives me a stack overflow (whereas the C++ version I implemented earlier can generate primes under 1,000,000+). I'm thinking the recursive functions are the main culprits, even though I tried to make them tail-recursive in order to avoid this very issue . . . Anyway, if I find a way to optimize the code any better I'll post it here later.

;;;; Generates a list of primes using the "Sieve of Atkin."
;;;; More on this algorithm can be found at:
;;;; http://en.wikipedia.org/wiki/Sieve_of_Atkin

(setf *limit* 5000)    ; limit for maximum prime in list


;;; generate list of "nil" of size lim
(defun null-list (lim) 
  (nullness lim '()))

(defun nullness (lim lst)
  (if (eql lim 0)
    lst
    (nullness (- lim 1) (cons nil lst))))


;;; sieve the primes
(defun sieve (lim)
  (setf is_prime (null-list lim))
  (do ((x 1 (+ x 1)))
      ((> x (+ (round (sqrt lim)) 1)) is_prime)
    (do ((y 1 (+ y 1)))
        ((> y (+ (round (sqrt lim)) 1)) is_prime)

      (setf n (+ (* 4 x x) (* y y)))
      (if (and (<= n lim) (or (eql (mod n 12) 1) (eql (mod n 12) 5)))
        (setf (nth n is_prime) (not (nth n is_prime))))

      (setf n (+ (* 3 x x) (* y y)))
      (if (and (<= n lim) (eql (mod n 12) 7))
        (setf (nth n is_prime) (not (nth n is_prime))))

      (setf n (- (* 3 x x) (* y y)))
      (if (and (> x y) (<= n lim) (eql (mod n 12) 11))
        (setf (nth n is_prime) (not (nth n is_prime)))))))


;;; set the square multiples of primes to false
(defun square-sieve (lst lim)
  (do ((n 0 (+ n 1)))
      ((> n (+ (round (sqrt lim))) 1) lst)
    (if (eql (nth n lst) 't)
      (do ((k (* n n) (+ k (* n n))))
          ((> k (- lim 1)) lst)
        (setf (nth k lst) nil)))))


;;; print the primes
(defun print-primes (lst lim)
  (prime-print-print lst 0 '()))

(defun prime-print-print (lst n primes)
  (if (eql n (list-length lst))
    (cons 2 (cons 3 (reverse primes)))
    (if (eql (nth n lst) 't)
      (prime-print-print lst (+ n 1) (setf primes (cons n primes)))
      (prime-print-print lst (+ n 1) primes))))


;;; print the primes for the above limit
(print (print-primes (square-sieve (sieve *limit*) *limit*) *limit*))

Instantiation

Hola, I'm starting this blog in order to document my meager coding attempts and perhaps post any cool *nix stuff I find as well. Right now I'm a CS student, fighting my way through coursework and the procrastination that comes with it, but I hope to at least post here weekly. Anyway, enough blabbing, let's get to the code!