Split-op in Haskell

CONTRIBUTORS.md

@@ -65,5 +65,4 @@ Trashtalk
 Cyrus Burt
 <br>
 Patrik Tesarik
-<br>
-
+<br>

book.json

@@ -138,7 +138,7 @@
  },
  {
  "lang": "f90",
- "name": "Fortran"
+ "name": "Fortran90"
  }
  ]
  }

contents/cooley_tukey/code/julia/fft.jl

@@ -1,3 +1,5 @@
+using FFTW
+
 #simple DFT function
 function DFT(x)
  N = length(x)
@@ -24,8 +26,8 @@ function cooley_tukey(x)
  n = 0:N-1
  half = div(N,2)
  factor = exp.(-2im*pi*n/N)
- return vcat(x_odd + x_even .* factor[1:half],
- x_odd - x_even .* factor[1:half])
+ return vcat(x_odd .+ x_even .* factor[1:half],
+ x_odd .- x_even .* factor[1:half])
 
 end
 
@@ -37,7 +39,7 @@ function bitreverse(a::Array)
 
  bit_indices = []
  for i = 1:length(indices)
- push!(bit_indices, bits(indices[i]))
+ push!(bit_indices, bitstring(indices[i]))
  end
 
  # Now stripping the unnecessary numbers
@@ -50,11 +52,11 @@ function bitreverse(a::Array)
  bit_indices[i] = reverse(bit_indices[i])
  end
 
- #replacing indices
+ # Replacing indices
  for i = 1:length(indices)
  indices[i] = 0
  for j = 1:digits
- indices[i] += 2^(j-1) * parse(string(bit_indices[i][end-j]))
+ indices[i] += 2^(j-1) * parse(Int, string(bit_indices[i][end-j]))
  end
  indices[i] += 1
  end

contents/cooley_tukey/cooley_tukey.md

@@ -70,17 +70,17 @@ For some reason, though, putting code to this transformation really helped me fi
 
 {% method %}
 {% sample lang="jl" %}
-[import:2-11, lang:"julia"](code/julia/fft.jl)
+[import:4-13, lang:"julia"](code/julia/fft.jl)
 {% sample lang="c" %}
 [import:8-19, lang:"c_cpp"](code/c/fft.c)
 {% sample lang="cpp" %}
 [import:23-33, lang:"c_cpp"](code/c++/fft.cpp)
 {% sample lang="hs" %}
-[import:2-11, lang:"julia"](code/julia/fft.jl)
+[import:4-13, lang:"julia"](code/julia/fft.jl)
 {% sample lang="py" %}
 [import:5-11, lang:"python"](code/python/fft.py)
 {% sample lang="scratch" %}
-[import:2-11, lang:"julia"](code/julia/fft.jl)
+[import:4-13, lang:"julia"](code/julia/fft.jl)
 {% endmethod %}
 
 In this function, we define `n` to be a set of integers from $$0 \rightarrow N-1$$ and arrange them to be a column.
@@ -115,7 +115,7 @@ With recursion, we can reduce the complexity to $$\sim \mathcal{O}(n \log n)$$,
 In the end, the code looks like:
 {% method %}
 {% sample lang="jl" %}
-[import:14-31, lang:"julia"](code/julia/fft.jl)
+[import:16-32, lang:"julia"](code/julia/fft.jl)
 {% sample lang="c" %}
 [import:20-39, lang:"c_cpp"](code/c/fft.c)
 {% sample lang="cpp" %}
@@ -125,7 +125,7 @@ In the end, the code looks like:
 {% sample lang="py" %}
 [import:13-24, lang:"python"](code/python/fft.py)
 {% sample lang="scratch" %}
-[import:14-31, lang:"julia"](code/julia/fft.jl)
+[import:16-32, lang:"julia"](code/julia/fft.jl)
 {% endmethod %}
 
 As a side note, we are enforcing that the array must be a power of 2 for the operation to work.

contents/euclidean_algorithm/code/fortran/euclidean.f90

@@ -0,0 +1,49 @@
+INTEGER FUNCTION euclid_sub(a, b)
+ IMPLICIT NONE
+ INTEGER, INTENT(INOUT) :: a, b
+
+ a = ABS(a)
+ b = ABS(b)
+
+ DO WHILE (a /= b)
+
+ IF (a > b) THEN
+ a = a - b
+ ELSE
+ b = b - a
+ END IF
+ END DO
+
+ euclid_sub = a
+
+END FUNCTION euclid_sub 
+
+INTEGER FUNCTION euclid_mod(a, b)
+ IMPLICIT NONE
+ INTEGER, INTENT(INOUT) :: a, b
+ INTEGER :: temp
+
+ DO WHILE (b > 0)
+ temp = b
+ b = MODULO(a,b)
+ a = temp
+ END DO
+
+ euclid_mod = a
+
+END FUNCTION euclid_mod
+
+PROGRAM euclidean
+
+ IMPLICIT NONE
+ INTEGER :: a, b, euclid_sub, euclid_mod
+
+ a = 24
+ b = 27
+ WRITE(*,*) 'Subtraction method: GCD is: ', euclid_sub(a, b)
+
+ a = 24
+ b = 27
+ WRITE(*,*) 'Modulus method: GCD is: ', euclid_mod(a, b)
+
+END PROGRAM euclidean 

contents/euclidean_algorithm/code/fortran/interactive_euclidean.f90

@@ -0,0 +1,60 @@
+INTEGER FUNCTION euclid_sub(a, b)
+ IMPLICIT NONE
+ INTEGER, INTENT(INOUT) :: a, b
+
+
+ a = ABS(a)
+ b = ABS(b)
+
+ DO WHILE (a /= b)
+
+ IF (a > b) THEN
+ a = a - b
+ ELSE
+ b = b - a
+ END IF
+ END DO
+
+ euclid_sub = a
+END FUNCTION euclid_sub 
+
+INTEGER FUNCTION euclid_mod(a, b)
+ IMPLICIT NONE
+ INTEGER, INTENT(INOUT) :: a, b
+ INTEGER :: temp
+
+
+ DO WHILE (b > 0)
+ temp = b
+ b = MODULO(a,b)
+ a = temp
+ END DO
+
+ euclid_mod = a
+
+END FUNCTION euclid_mod
+
+PROGRAM euclidean
+
+ IMPLICIT NONE
+ INTEGER :: a, b, euclid_sub, euclid_mod, temp_a, temp_b, ioerror
+
+ DO
+ WRITE(*,*) 'Calculate greatest common divisor. Give two integers:'
+ READ(*, '(i10)', iostat=ioerror) temp_a, temp_b
+
+ IF (ioerror == 0) THEN
+ EXIT
+ END IF
+
+ WRITE(*,*) 'Entered numbers are not integers. Try again.'
+ END DO
+
+ a = temp_a
+ b = temp_b
+ WRITE(*,*) 'Subtraction method: GCD is: ', euclid_sub(a, b)
+
+ a = temp_a 
+ b = temp_b 
+ WRITE(*,*) 'Modulus method: GCD is: ', euclid_mod(a, b)
+END PROGRAM euclidean 

contents/euclidean_algorithm/euclidean_algorithm.md

@@ -39,6 +39,8 @@ The algorithm is a simple way to find the *greatest common divisor* (GCD) of two
 [import:12-25, lang="julia"](code/julia/euclidean.jl)
 {% sample lang="nim" %}
 [import:13-24, lang="nim"](code/nim/euclid_algorithm.nim)
+{% sample lang="f90" %}
+[import:1-19, lang="Fortran"](code/fortran/euclidean.f90)
 {% endmethod %}
 
 Here, we simply line the two numbers up every step and subtract the lower value from the higher one every timestep. Once the two values are equal, we call that value the greatest common divisor. A graph of `a` and `b` as they change every step would look something like this:
@@ -84,6 +86,8 @@ Modern implementations, though, often use the modulus operator (%) like so
 [import:1-10, lang="julia"](code/julia/euclidean.jl)
 {% sample lang="nim" %}
 [import:1-11, lang="nim"](code/nim/euclid_algorithm.nim)
+{% sample lang="f90" %}
+[import:21-34, lang="Fortran"](code/fortran/euclidean.f90)
 {% endmethod %}
 
 Here, we set `b` to be the remainder of `a%b` and `a` to be whatever `b` was last timestep. Because of how the modulus operator works, this will provide the same information as the subtraction-based implementation, but when we show `a` and `b` as they change with time, we can see that it might take many fewer steps:
@@ -134,6 +138,8 @@ The Euclidean Algorithm is truly fundamental to many other algorithms throughout
 [import, lang="julia"](code/julia/euclidean.jl)
 {% sample lang="nim" %}
 [import, lang="nim" %](code/nim/euclid_algorithm.nim)
+{% sample lang="f90" %}
+[import, lang="Fortran"](code/fortran/euclidean.f90)
 {% endmethod %}
 
 

contents/gaussian_elimination/code/julia/gaussian_elimination.jl

@@ -10,7 +10,7 @@ function gaussian_elimination(A::Array{Float64,2})
  for col = 1:(cols-1)
 
  # Step 1: finding the maximum element for each column
- max_index = indmax(abs.(A[row:end,col])) + row-1
+ max_index = argmax(abs.(A[row:end,col])) + row-1
 
  # Check to make sure matrix is good!
  if (A[max_index, col] == 0)
@@ -50,7 +50,7 @@ function back_substitution(A::Array{Float64,2})
  cols = size(A,2)
 
  # Creating the solution Vector
- soln = Vector{Float64}(rows)
+ soln = zeros(rows)
 
  for i = rows:-1:1
  sum = 0.0

contents/gaussian_elimination/gaussian_elimination.md

@@ -420,7 +420,7 @@ In code, this involves keeping a rolling sum of all the values we substitute in
 
 {% method %}
 {% sample lang="jl" %}
-[import:47-67, lang:"julia"](code/julia/gaussian_elimination.jl)
+[import:47-64, lang:"julia"](code/julia/gaussian_elimination.jl)
 {% sample lang="c" %}
 [import:50-62, lang:"c_cpp"](code/c/gaussian_elimination.c)
 {% sample lang="rs" %}

contents/graham_scan/code/julia/graham.jl

@@ -1,4 +1,4 @@
-type Point
+struct Point
  x::Float64
  y::Float64
 end
@@ -14,8 +14,8 @@ function graham_scan(points::Vector{Point})
  sort!(points, by = item -> item.y)
 
  # Sort all other points according to angle with that point
- other_points = sort(points[2:end], by = item -> atan2(item.y - points[1].y,
-  item.x - points[1].x))
+ other_points = sort(points[2:end], by = item -> atan(item.y - points[1].y,
+ item.x - points[1].x))
 
  # Place points sorted by angle back into points vector
  for i in 1:length(other_points)
@@ -46,7 +46,11 @@ end
 
 function main()
  # This hull is just a simple test so we know what the output should be
- points = [Point(2,1.9), Point(1, 1), Point(2, 4), Point(3, 1), Point(2, 0)]
+ points = [
+ Point(-5.,2), Point(5,7), Point(-6,-12), Point(-14,-14), Point(9,9),
+ Point(-1,-1), Point(-10,11), Point(-6,15), Point(-6,-8), Point(15,-9),
+ Point(7,-7), Point(-2,-9), Point(6,-5), Point(0,14), Point(2,8)
+ ]
  hull = graham_scan(points)
  println(hull)
 end

contents/quantum_systems/code/haskell/Energy.hs

@@ -0,0 +1,14 @@
+import Data.Array.CArray
+import Data.Complex
+import Math.FFT (dft, idft) -- Binding to fftw
+
+type Vector = CArray Int (Complex Double)
+
+calculateEnergy :: Double -> Vector -> Vector -> Vector -> Double
+calculateEnergy dx kin pot wfc = (* dx) . sum . map realPart $ elems total
+ where
+ total = liftArray2 (+) kineticE potentialE
+ potentialE = wfcConj .* pot .* wfc
+ kineticE = wfcConj .* idft (kin .* dft wfc)
+ wfcConj = liftArray conjugate wfc
+ a .* b = liftArray2 (*) a b

contents/quantum_systems/quantum_systems.md

@@ -227,6 +227,8 @@ This ultimately looks like this:
 {% method %}
 {% sample lang="jl" %}
 [import, lang:"julia"](code/julia/energy.jl)
+{% sample lang="hs" %}
+[import, lang:"haskell"](code/haskell/Energy.hs)
 {% sample lang="c" %}
 [import:28-, lang:"c_cpp"](code/c/energy.c)
 {% sample lang="py" %}

contents/split-operator_method/code/gnuplot/plot_output.plt

@@ -0,0 +1,24 @@
+# use like: gnuplot -e "folder='/path/to/data/directory'" plot_output.plt
+
+set terminal gif animate delay 10
+set output folder.'/wavefunction.gif'
+
+set key off
+
+set xrange [-10:10] 
+set yrange [0:1]
+set y2range [0:50]
+
+set ytics nomirror autofreq tc lt 1
+set ylabel 'Psi(x)' tc lt 1
+
+set y2tics nomirror autofreq tc lt 2
+set y2label 'V(x)' tc lt 2
+
+list = system('ls '.folder.'/output*.dat')
+
+do for [file in list] {
+ plot file u 1:2 smooth csplines, \
+ file u 1:3 smooth csplines axes x1y2
+}
+