Implement the GTH algorithm

This is to import [gth_solve.py](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/gth_solve.py) to the Julia part (I also imported the tests as well).
(I have had a brief discussion with @ZacCranko on this matter.)

1.
This is why I think `gth_solve` is necessary also for the Julia side:

I guess that the source of the problem is the same as in NumPy and SciPy, but it seems the routine (probably `lufact`) from the Julia standard library loses precision for some inputs with very small float numbers.

Here's an example:

```julia
for i in 12:18
    a = 1e-1 ^ i
    P = [1-a a; a/2 1-a/2]
    mc = MarkovChain(P)
    print("a = 1e-", i, "\n")
    print(mc_compute_stationary(mc), "\n")
end
```

The output is:

```julia
a = 1e-12
[0.3333530884328846,0.6666469115671155]
a = 1e-13
[0.33326424872384786,0.6667357512761523]
a = 1e-14
[0.33315566870167473,0.6668443312983253]
a = 1e-15
[0.3569592970777264,0.6430407029222737]
a = 1e-16
[0.3105159921038858,0.6894840078961141]
a = 1e-17
[0.0,1.0]
a = 1e-18
[0.0,1.0]
```

The correct answer is [1/3, 2/3] for any value of `a`.

A problem occurs also for "KMR Markov matrices" with small epsilon.

The "GTH" algorithm is robust to this kind of inputs, as it does not involve subtraction. This is why I included "gth_solve.py" into the python side, and I believe this is also necessary for the julia side.

2.
I de-vectorized the python version following the general advice for Julia. Surprisingly, this de-vectorized version is 10x faster than the direct translation of the python version using array indexing for some random 1000x1000 matrices.

3.
Docstring is yet to be imported.

4.
I tried to follow a suggestion by @ZacCranko and manipulated the code, and then it appears the routine now returns an output with rational number elements for a matrix input with rational number elements:

```julia
P = [2//3 1//3; 1//4 3//4]
gth_solve(P)
```

Output:

```julia
2-element Array{Rational{Int64},1}:
 3//7
 4//7
```

But I don't think I really understand what is going on..


Anyway, this is effectively my first code in Julia, and definitely needs reviews by experienced Julia programmers... Thanks!

@oyamad I think you are right about needing to have floats. 

Here is a "julian" way of doing this:

```
# T<: Real catches any kind of real matrix, converts it to float64, passes it to main method
gth_solve{T <:  Real}(A::Matrix{T}) = gth_solve(float64(A))

function gth_solve{T <: FloatingPoint}(A::AbstractMatrix{T})
    if size(A, 1) != size(A, 2)
        throw(ArgumentError("matrix must be square"))
    end

    A1 = copy(A)
    T_x = Float64

    n = size(A1, 1)

    x = zeros(T, n)
    # .... continue with your original code
```