I'm implementing Newton's method for finding roots in a precompiled library. The common case will be to work with functions from Float64
to Float64
, and I want an optimized compiled version of it to exist in the library. I will, of course, implement a type generic version, too, but Julia will need some way of differentiating the methods by signature so it knows which one to call at runtime. Presently, the implementation is:
#Float64 optimized version
function findrootNewton( func, funcder, guess::Float64,
rtol::Float64=1e-12, abstol::Float64=1e-12, maxiter::Int=100 )
#make sure tolerances are positive
if rtol <= 0.0
error("findrootNewton: rtol must be a positive number")
end
if abstol <= 0.0
error("findrootNewton: abstol must be a positive number")
end
if maxiter <= 0.0
error("findrootNewton: maxiter must be a positive integer")
end
converged::Bool = false
oldx::Float64 = guess
newx::Float64 = oldx - ((func(oldx) / funcder(oldx))::Float64)
absdiff = abs(oldx - newx)
iter = 2
while (absdiff < abstol || absdiff < rtol * abs(newx)) && iter <= maxiter
oldx = newx
newx = oldx - func(oldx) / funcder(oldx)
absdiff = abs(oldx - newx)
iter += 1
end #while (absdiff < abstol || absdiff < rtol * abs(newx)) && newxiter <= maxiter
if iter <= maxiter
converged = true
end
return (newx, converged)
end #findzeroNewton
#Generic version
function findrootNewton( func, funcder, guess::Number,
rtol::Real=1e-12, abstol::Real=1e-12, maxiter::Int=100 )
#make sure tolerances are positive
if rtol <= 0
error("findrootNewton: rtol must be a positive number")
end
if abstol <= 0
error("findrootNewton: abstol must be a positive number")
end
if maxiter <= 0
error("findrootNewton: maxiter must be a positive integer")
end
converged::Bool = false
newx = oldx - func(oldx) / funcder(oldx)
oldx = convert(typeof(newx), guess)
absdiff = abs(oldx - newx)
iter = 2
while (absdiff < abstol || absdiff < rtol * abs(newx)) && iter <= maxiter
oldx = newx
newx = oldx - func(oldx) / funcder(oldx)
absdiff = abs(oldx - newx)
iter += 1
end #while (absdiff < abstol || absdiff < rtol * abs(newx)) && newxiter <= maxiter
if iter <= maxiter
converged = true
end
return (newx, converged)
end #findzeroNewton
This has not been used/debugged, yet. For example, I'm not checking for the derivative hitting zero because the particular case I'm coding this for does not need it.
Notice that if guess
, rtol
, and abstol
arguments are given as Float64
, but the functions don't return Float64
but BigFloat
, then the code will fail at the point of the type assertion in the definition of newx
, even though there is a method appropriate to generic functions available. Is there a way to avoid this problem?
Edit: I can specify the type of data a variable stores in Julia, for example:
x::Float64 = 2.5
Is it possible to similarly specify the signature of a variable that can store (pointers to) functions?
Float64(... expression returning BigFloat ...)
. Also, in Julia it is often recommended to start off with less typing, as the typing often flows from compile time type inference and dynamically at run-time. – Boundsconvert(...)
or constructors attempt to convert one type into another. I think you may want the latter here, but perhaps I am misunderstanding the question. – Boundsfunction f(); 2; end::Int
– Moorfowl