Saturday, December 7, 2013

Nand2Tetris

So I've found this course online that intends to teach you how to build a complete computer from first principles (i.e. from the digital logic level and up). I think I'm going to try and put a dent in it over the christmas break, so I'll update as I go through it. Since I'm about to complete a course on computer architecture (using this book), I'm excited to "fill in the blanks" with what I've undoubtedly missed in my studies.

Friday, November 22, 2013

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*))

Wednesday, November 20, 2013

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*))

Tuesday, November 19, 2013

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!