WIP: Support 0.6

.travis.yml

@@ -3,7 +3,6 @@ os:
   - linux
   - osx
 julia:
-  - 0.4
   - 0.5
   - nightly
 notifications:

REQUIRE

@@ -1,2 +1,3 @@
-julia 0.4
+julia 0.5
 Compat 0.8.8
+SpecialFunctions 0.1

src/FastGaussQuadrature.jl

@@ -1,6 +1,6 @@
 module FastGaussQuadrature
 
-using Compat
+using Compat, SpecialFunctions
 
 export gausslegendre
 export gausschebyshev
@@ -11,6 +11,8 @@ export gausslobatto
 export gaussradau
 export besselroots
 
+import SpecialFunctions: besselj, airyai, airyaiprime
+
 include("gausslegendre.jl")
 include("gausschebyshev.jl")
 include("gausslaguerre.jl")

src/besselroots.jl

@@ -1,4 +1,4 @@
-function besselroots(nu::Float64, n::Int64) 
+function besselroots(nu::Float64, n::Int64)
 #BESSELROOTS    The first N roots of the function J_v(x)
 
 # DEVELOPERS NOTES:
@@ -10,17 +10,17 @@ function besselroots(nu::Float64, n::Int64)
 #   V > 5 --> moderately accurate for the 6 first zeros and good
 #     approximations for the others (McMahon's expansion)
 
-# This code was originally written by L. L. Peixoto in MATLAB. 
+# This code was originally written by L. L. Peixoto in MATLAB.
 # Later modified by A. Townsend to work in Julia
-    
+
 if ( n == 0 )
     x = []
 elseif( n > 0 && nu == 0 )
     x = McMahon( nu, n )
     correctFirstFew = Tabulate( )
     x[1:min(n,20)] = correctFirstFew[1:min(n,20)]
     x
-elseif ( n>0 && nu >= -1 && nu <= 5 ) 
+elseif ( n>0 && nu >= -1 && nu <= 5 )
     x = McMahon( nu, n )
     correctFirstFew = Piessens( nu )
     x[1:min(n,6)] = correctFirstFew[1:min(n,6)]
@@ -33,55 +33,55 @@ end
 
 
 function McMahon( nu::Float64, n::Int64 )
-# McMahon's expansion. This expansion gives very accurate approximation 
-# for the sth zero (s >= 7) in the whole region V >=- 1, and moderate
-# approximation in other cases.
-b = Array(Float64,n)
-mu = 4*nu^2
-a1 = 1. / 8.
-a3 = (7.*mu-31.) / 384.
-a5 = 4.*(3779.+mu*(-982.+83.*mu)) / 61440. # Evaluate via Horner's method.
-a7 = 6.*(-6277237.+mu*(1585743.+mu*(-153855.+6949.*mu))) / 20643840.;
-a9 = 144.*(2092163573.+mu*(-512062548.+mu*(48010494.+mu*(-2479316.+70197.*mu)))) / 11890851840.
-a11 = 720.*(-8249725736393.+mu*(1982611456181.+mu*(-179289628602.+mu*(8903961290. + 
-      mu*(-287149133.+5592657.*mu))))) / 10463949619200.
-a13 = 576.*(423748443625564327. + mu*(-100847472093088506.+mu*(8929489333108377. + 
-    mu*(-426353946885548.+mu*(13172003634537.+mu*(-291245357370. + mu*4148944183.)))))) / 13059009124761600.
-for k=1:n b[k] = .25*(2*nu+4*k-1)*pi end # beta
-# Evaluate using Horner's scheme: 
-x = b - (mu-1)*( ((((((a13./b.^2 + a11)./b.^2 + a9)./b.^2 + a7)./b.^2 + a5)./b.^2 + a3)./b.^2 + a1)./b)    
+    # McMahon's expansion. This expansion gives very accurate approximation
+    # for the sth zero (s >= 7) in the whole region V >=- 1, and moderate
+    # approximation in other cases.
+    b = Array{Float64}(n)
+    mu = 4*nu^2
+    a1 = 1. / 8.
+    a3 = (7.*mu-31.) / 384.
+    a5 = 4.*(3779.+mu*(-982.+83.*mu)) / 61440. # Evaluate via Horner's method.
+    a7 = 6.*(-6277237.+mu*(1585743.+mu*(-153855.+6949.*mu))) / 20643840.;
+    a9 = 144.*(2092163573.+mu*(-512062548.+mu*(48010494.+mu*(-2479316.+70197.*mu)))) / 11890851840.
+    a11 = 720.*(-8249725736393.+mu*(1982611456181.+mu*(-179289628602.+mu*(8903961290. +
+          mu*(-287149133.+5592657.*mu))))) / 10463949619200.
+    a13 = 576.*(423748443625564327. + mu*(-100847472093088506.+mu*(8929489333108377. +
+        mu*(-426353946885548.+mu*(13172003634537.+mu*(-291245357370. + mu*4148944183.)))))) / 13059009124761600.
+    for k=1:n b[k] = .25*(2*nu+4*k-1)*pi end # beta
+    # Evaluate using Horner's scheme:
+    x = b - (mu-1)*( ((((((a13./b.^2 + a11)./b.^2 + a9)./b.^2 + a7)./b.^2 + a5)./b.^2 + a3)./b.^2 + a1)./b)
 end
 
 
 function Tabulate( )
-# First 20 roots of J0(x) are precomputed (using Wolfram Alpha):
-y = Array(Float64,20)
-y = [   2.4048255576957728
-        5.5200781102863106
-        8.6537279129110122
-        11.791534439014281
-        14.930917708487785
-        18.071063967910922
-        21.211636629879258
-        24.352471530749302
-        27.493479132040254
-        30.634606468431975
-        33.775820213573568
-        36.917098353664044
-        40.058425764628239
-        43.199791713176730
-        46.341188371661814
-        49.482609897397817
-        52.624051841114996
-        55.765510755019979
-        58.906983926080942
-        62.048469190227170  ]
+    # First 20 roots of J0(x) are precomputed (using Wolfram Alpha):
+    y = Array{Float64}(20)
+    y = [   2.4048255576957728
+            5.5200781102863106
+            8.6537279129110122
+            11.791534439014281
+            14.930917708487785
+            18.071063967910922
+            21.211636629879258
+            24.352471530749302
+            27.493479132040254
+            30.634606468431975
+            33.775820213573568
+            36.917098353664044
+            40.058425764628239
+            43.199791713176730
+            46.341188371661814
+            49.482609897397817
+            52.624051841114996
+            55.765510755019979
+            58.906983926080942
+            62.048469190227170  ]
 end
 
-function Piessens( nu::Float64 ) 
-# Piessens's Chebyshev series approximations (1984). Calculates the 6 first
-# zeros to at least 12 decimal figures in region -1 <= V <= 5:
-    y = Array(Float64,6);
+function Piessens( nu::Float64 )
+    # Piessens's Chebyshev series approximations (1984). Calculates the 6 first
+    # zeros to at least 12 decimal figures in region -1 <= V <= 5:
+    y = Array{Float64}(6);
 
     C = [
        2.883975316228  8.263194332307 11.493871452173 14.689036505931 17.866882871378 21.034784308088
@@ -114,8 +114,8 @@ function Piessens( nu::Float64 )
        0.000000000008  0.000000000011               0.               0.               0.               0.
       -0.000000000003 -0.000000000005               0.               0.               0.               0.
        0.000000000001  0.000000000002               0.               0.               0.               0.]
-    
-    T = Array(Float64,size(C,1))
+
+    T = Array{Float64}(size(C,1))
     pt = (nu-2)/3
     T[1], T[2] = 1., pt
     for k = 2:size(C,1)-1
@@ -124,4 +124,4 @@ function Piessens( nu::Float64 )
     y[1:6] = T'*C;
     y[1] *= sqrt(nu+1);                  # Scale the first root.
     y
-end
+end

src/gausshermite.jl

@@ -39,20 +39,20 @@ function hermpts_asy( n::Int )
 
 x0 = HermiteInitialGuesses( n ) # get initial guesses
 t0 = x0./sqrt(2n+1)
-theta0 = acos(t0)               # convert to theta-variable
+theta0 = acos.(t0)               # convert to theta-variable
 val = x0;
 for k = 1:20
     val = hermpoly_asy_airy(n, theta0);
-    dt = -val[1]./(sqrt(2)*sqrt(2n+1)*val[2].*sin(theta0))
-    theta0 = theta0 - dt;                        # Newton update
+    dt = -val[1]./(sqrt(2).*sqrt(2n+1).*val[2].*sin.(theta0))
+    theta0 .-= dt;                        # Newton update
     if norm(dt,Inf) < sqrt(eps(Float64))/10
        break
     end
 end
-t0 = cos(theta0)
+t0 = cos.(theta0)
 x = sqrt(2n+1)*t0                          #back to x-variable
-ders = x.*val[1] + sqrt(2)*val[2]
-w = (exp(-x.^2)./ders.^2)';            # quadrature weights
+ders = x.*val[1] .+ sqrt(2).*val[2]
+w = exp.(-x.^2)./ders.^2;            # quadrature weights
 
 x = (x, w)
 end
@@ -66,14 +66,14 @@ val = x0
 for kk = 1:10
     val = hermpoly_rec(n, x0)
     dx = val[1]./val[2]
-    dx[ isnan( dx ) ] = 0
+    dx[ isnan.( dx ) ] = 0
     x0 = x0 - dx
     if norm(dx, Inf)<sqrt(eps(Float64))
         break
     end
 end
 x = x0/sqrt(2)
-w = exp(-x.^2)./val[2].^2           # quadrature weights
+w = exp.((-).(x.^2))./val[2].^2           # quadrature weights
 
 x = (x, w)
 end
@@ -82,8 +82,8 @@ function hermpoly_rec( n::Int, x0)
 # HERMPOLY_rec evaluation of scaled Hermite poly using recurrence
 
 # evaluate:
-Hold = exp(-x0.^2/4)
-H = x0.*exp(-x0.^2/4)
+Hold = exp.(x0.^2./(-4))
+H = x0.*exp.(x0.^2./(-4))
 for k = 1:n-1
     Hold, H = H, (x0.*H./sqrt(k+1) - Hold./sqrt(1+1/k))
 end
@@ -97,15 +97,15 @@ function hermpoly_asy_airy(n::Int, theta)
 # theta-space.
 
 musq = 2n+1;
-cosT = cos(theta)
-sinT = sin(theta)
-sin2T = 2*cosT.*sinT
-eta = .5*theta - .25*sin2T
+cosT = cos.(theta)
+sinT = sin.(theta)
+sin2T = 2.*cosT.*sinT
+eta = 0.5.*theta .- 0.25.*sin2T
 chi = -(3*eta/2).^(2/3)
 phi = (-chi./sinT.^2).^(1/4)
 C = 2*sqrt(pi)*musq^(1/6)*phi
-Airy0 = real(airy(musq.^(2/3)*chi))
-Airy1 = real(airy(1,musq.^(2/3)*chi))
+Airy0 = real.(airyai.(musq.^(2/3).*chi))
+Airy1 = real.(airyaiprime.(musq.^(2/3).*chi))
 
 # Terms in (12.10.43):
 a0 = 1; b0 = 1
@@ -179,9 +179,9 @@ function HermiteInitialGuesses( n::Int )
 # [2] F. G. Tricomi, Sugli zeri delle funzioni di cui si conosce una
 # rappresentazione asintotica, Ann. Mat. Pura Appl. 26 (1947), pp. 283-300.
 
-# Error if n < 20 because initial guesses are based on asymptotic expansions: 
-@assert n>=20    
-    
+# Error if n < 20 because initial guesses are based on asymptotic expansions:
+@assert n>=20
+
 # Gatteschi formula involving airy roots [1].
 # These initial guess are good near x = sqrt(n+1/2);
 if mod(n,2) == 1
@@ -209,9 +209,9 @@ airyrts_exact = [-2.338107410459762           # Exact Airy roots.
     -12.828776752865757]
 airyrts[1:10] = airyrts_exact  # correct first 10.
 
-x_init = sqrt(abs(nu + 2^(2/3)*airyrts*nu^(1/3) + 1/5*2^(4/3)*airyrts.^2*nu^(-1/3) +
-    (11/35-a^2-12/175*airyrts.^3)/nu + (16/1575*airyrts+92/7875*airyrts.^4)*2^(2/3)*nu^(-5/3) -
-    (15152/3031875*airyrts.^5+1088/121275*airyrts.^2)*2^(1/3)*nu^(-7/3)))
+x_init = sqrt.(abs.(nu .+ (2^(2/3)).*airyrts.*nu^(1/3) .+ (1/5*2^(4/3)).*airyrts.^2.*nu^(-1/3) .+
+    (11/35-a^2-12/175).*airyrts.^3./nu .+ ((16/1575).*airyrts.+(92/7875).*airyrts.^4).*2^(2/3).*nu^(-5/3) .-
+    ((15152/3031875).*airyrts.^5.+(1088/121275).*airyrts.^2).*2^(1/3).*nu^(-7/3)))
 x_init_airy = real( flipdim(x_init,1) )
 
 # Tricomi initial guesses. Equation (2.1) in [1]. Originally in [2].
@@ -223,14 +223,14 @@ nu = (4*m+2*a+2)
 rhs = (4*m-4*collect(1:m)+3)./nu*pi
 
 for k = 1:7
-    val = Tnk0 - sin(Tnk0) - rhs
-    dval = 1 - cos(Tnk0)
+    val = Tnk0 .- sin.(Tnk0) .- rhs
+    dval = 1 .- cos.(Tnk0)
     dTnk0 = val./dval
-    Tnk0 = Tnk0 - dTnk0
+    Tnk0 = Tnk0 .- dTnk0
 end
 
-tnk = cos(Tnk0/2).^2
-x_init_sin = sqrt(nu*tnk - (5./(4*(1-tnk).^2) - 1./(1-tnk)-1+3*a^2)/3/nu)
+tnk = cos.(Tnk0./2).^2
+x_init_sin = sqrt.(nu*tnk .- (5./(4.*(1-tnk).^2) .- 1./(1.-tnk).-1.+3*a^2)./3./nu)
 
 # Patch together
 p = 0.4985+eps(Float64)
@@ -251,7 +251,7 @@ end
 function hermpts_gw( n::Int )
 # Golub--Welsch algorithm. Used here for n<=20.
 
-    beta = sqrt(.5*(1:n-1))              # 3-term recurrence coeffs
+    beta = sqrt.(0.5.*(1:n-1))              # 3-term recurrence coeffs
     T = diagm(beta, 1) + diagm(beta, -1)   # Jacobi matrix
     (D, V) = eig(T)                      # Eigenvalue decomposition
     indx = sortperm(D)                  # Hermite points
@@ -263,4 +263,4 @@ function hermpts_gw( n::Int )
     x = x[ii]
     w = w[ii]
     return (x,w)
-end
+end