As part of a numerical algorithm, I want to compute x/\sqrt{y} for real numbers x and y. In general is it preferable to compute x/sqrt(y) or x*sqrt(1/y) or x*1/sqrt(y)?
As a slightly less trivial example, is there an optimal way to compute x/\sqrt{y^3/3}?
You can test it against the results of BigFloats:
d(f, x, y) = f(x, y) - f(big(x), big(y))
I got deviations in the order of 10^-17 for both methods and Float64.
x*1/sqrt(y) is by the way identical to x/sqrt(y).
x/sqrt(y) requires the fewest operations and characters, and it should be at least as good as the other versions in terms of accuracy.
I get absolutely identical results in all cases:
julia> d(f,x,y) = f(x,y) - f(big(x),big(y))
julia> f1(x,y) = x/sqrt(y);
julia> f2(x,y) = x*sqrt(1/y);
julia> f3(x,y) = x*(1/sqrt(y));
julia> f4(x,y) = sqrt(x^2/y);
julia> x = rand(); y = rand();
julia> d(f1,x,y) == d(f2,x,y) == d(f3,x,y) == d(f4,x,y)
true
julia> d(f1,x,y)
-4.407656372008221707500556004683899793700299580520962546832391917766617849819689e-18
Does that makes sense or are we tricked by something? The lowered codes are not identical.
If we define x and y as BigFloats and convert them back instead in d, then we get a difference:
julia> x = rand(BigFloat)
0.9793180540915932364829244784180117806999674427403793663059170823324311264381625
julia> y = rand(BigFloat)
0.3233069578761485009887950171010027218596435775913905575889005425608602252122216
julia> d(f,x,y) = f(Float64(x),Float64(y)) - f(big(x),big(y))
d (generic function with 1 method)
julia> d(f1,x,y)
-9.759393809419778848327739693818875724638531567376830791086116873272363299254071e-17
julia> d(f2,x,y)
1.2445066683083351960144893667997530525361468432623169208913881399493925681857e-16
julia> d(f3,x,y)
-9.759393809419778848327739693818875724638531567376830791086116873272363299254071e-17
julia> d(f4,x,y)
-9.759393809419778848327739693818875724638531567376830791086118600506074318142996e-17
If you wanted to do better, you could work out the Taylor series for this specific function and implement it that way, but that seems unlikely to be worth it. This also typically requires splitting the domain into separate regions and doing different series centered around different points and piecing them together in order to get optimal accuracy and speed.
Thanks all for your replies. With @lungben’s BigFloat method I can now check the accuracy of different functions myself.
Edit: my second question didn’t make sense…
I’m not sure what you mean here — this function does not have a Taylor series in y. It has a multivariate Puiseux series with exactly one term: x/\sqrt{y}. So, that road leads in a circle.
(Unless you mean computing the function near some particular nonzero value of y, and doing a series expansion there? It seems like it would be hard to get a performance boost that way unless you are looking at a very small domain.)
Sure, that 
I’d test it around the domain you’re interested in. For example:
xs = randn(1000)
ys = abs.(randn(1000))
d(f, x, y) = (f(x, y) - f(big(x), big(y)))/eps(f(x,y)))
histogram(vec(d.((x,y)->x/sqrt(y), xs, ys')), label="x/sqrt(y)", xlabel="ulps")
histogram!(vec(d.((x,y)->x*sqrt(1/y), xs, ys')), label="x*sqrt(1/y)", alpha=0.4)
Or:
julia> extrema(((x,y),)->d((x,y)->x/sqrt(y),x,y), Iterators.product(xs, ys))
(-1.445641472157348115162593165051000134893206566567640842822595139202883960327676, 1.431037408382898990935666476069582924726289593427886810488154045658482591603165)
julia> extrema(((x,y),)->d((x,y)->x*sqrt(1/y),x,y), Iterators.product(xs, ys))
(-1.786730856066382917222061843607036797336143841622304373374413114890563230145751, 1.747706262134240363279541031882310620807941690399233701771142762467278005276572)
There should be slight differences for certain inputs, e.g.:
julia> x, y = 0.5164620277491543, 0.8164777692882677;
julia> x/sqrt(y)
0.5715657726427575
julia> x*sqrt(1/y)
0.5715657726427577
julia> exact = big(x)/sqrt(big(y))
0.5715657726427576028583756984664763036711584135326458415886212315578020257028185
julia> Float64(x/sqrt(y) - exact)
-6.150312393821535e-17
julia> Float64(x*sqrt(1/y) - exact)
4.95191785243003e-17
Hehe, indeed. But it does happen that for some random numbers there is no difference, and I was not lucky in the first random choice: 
julia> x = rand(); y = rand();
julia> d(f1,x,y) == d(f2,x,y) == d(f3,x,y) == d(f4,x,y)
false
julia> x = rand(); y = rand();
julia> d(f1,x,y) == d(f2,x,y) == d(f3,x,y) == d(f4,x,y)
false
julia> x = rand(); y = rand();
julia> d(f1,x,y) == d(f2,x,y) == d(f3,x,y) == d(f4,x,y)
true
julia> x
0.7188671601251424
julia> y
0.6747963454138877
Herbie seems to recommend the x * y^(-0.5) variant, although the (very slight) increase in accuracy detected on the 256-point training set does not seem to be reliably found in the 8000-point test set. I would say that this tends to corroborate the general trend of answers here, that there is no easy way to rewrite x/\sqrt y:
Herbie 1.4 with seed 172253712
Find help on https://herbie.uwplse.org/, exit with Ctrl-D
herbie> (FPCore (x y) :name "x / sqrt(y)" (/ x (sqrt y)))
(FPCore
(x y)
:herbie-status
ex-start
:herbie-time
6874.320068359375
:herbie-error-input
((256 0.27058823529411763) (8000 0.23190398799849982))
:herbie-error-output
((256 0.23921568627450981) (8000 0.24678084760595073))
:name
"x / sqrt(y)"
:precision
binary64
(* x (pow y -0.5)))
However, it looks like your second formula (x/\sqrt{y^3/3}) is already quite good on average, but could still be rewritten significantly more accurately in the (arguably less elegant and undoubtedly more expensive) way described below:
herbie> (FPCore (x y) :name "x/sqrt(y^3/3)" (/ x (/ (pow y 3) 3)))
(FPCore
(x y)
:herbie-status
imp-start
:herbie-time
12470.798095703125
:herbie-error-input
((256 7.9874078337477235) (8000 7.963538107495137))
:herbie-error-output
((256 0.9712115675397124) (8000 1.0766656794517426))
:name
"x/sqrt(y^3/3)"
:precision
binary64
(*
(/
(/
(* (* (cbrt x) (cbrt x)) (cbrt 3.0))
(pow (pow (cbrt y) 2.0) (/ 3.0 2.0)))
(/ (pow (* (cbrt y) (cbrt y)) (/ 3.0 2.0)) (cbrt 3.0)))
(/ (cbrt x) (/ (pow y 1.0) (cbrt 3.0)))))