Strange compile behavior for @turbo

I have a very long sequence of scalar mathematical operations [example of code shown at bottom of thread]. The below can be compiled but if I add one additional line:

```julia
b = sqrt(taupa^2+1)
```

I get a: LoadError: "Reduction not found." error

Now if I change this line to the equivalent:
```julia
b = (taupa^2+1)^0.5
```

then I do not get an error?

Any thoughts on what might be going on?

P.S. Thanks for the absolutely amazing package

```julia
@turbo for i = eachindex(y)
y[i] = y[i] + y0
xi = y[i] / a1
eta = x[i] / a1
xisign = sign(xi)
etasign = sign(eta)
xi = xi * xisign
eta = eta * etasign
backside = xi > p2

xi = (p - xi) * backside + (xi * !backside)

c0 = cos(twoT * xi)
ch0 = cosh(twoT * eta)
s0 = sin(twoT * xi)
sh0 = sinh(twoT * eta)

ar = twoT * c0 * ch0
ai = twoT * -s0 * sh0

# --- j = 6 ------
y1r = - bet[6]
y1i = zeroT

z1r = -twoT * 6 * bet[6]
z1i = zeroT

y0r = ar * y1r - ai * y1i - bet[6-1]
y0i = ar * y1i + ai * y1r

z0r = ar * z1r - ai * z1i - twoT * (6 - 1) * bet[6-1]
z0i = ar * z1i + ai * z1r

# --- j = 4 ------
y1r = ar * y0r - ai * y0i - y1r - bet[4]
y1i = ar * y0i + ai * y0r - y1i

z1r = ar * z0r - ai * z0i - z1r - twoT * 4 * bet[4]
z1i = ar * z0i + ai * z0r - z1i

y0r = ar * y1r - ai * y1i - y0r - bet[4-1]
y0i = ar * y1i + ai * y1r - y0i

z0r = ar * z1r - ai * z1i - z0r - twoT * (4 - 1) * bet[4-1]
z0i = ar * z1i + ai * z1r - z0i

# --- j = 2 ------
y1r = ar * y0r - ai * y0i - y1r - bet[2]
y1i = ar * y0i + ai * y0r - y1i

z1r = ar * z0r - ai * z0i - z1r - twoT * 2 * bet[2]
z1i = ar * z0i + ai * z0r - z1i

y0r = ar * y1r - ai * y1i - y0r - bet[2-1]
y0i = ar * y1i + ai * y1r - y0i

z0r = ar * z1r - ai * z1i - z0r - twoT * (2 - 1) * bet[2-1]
z0i = ar * z1i + ai * z1r - z0i

# ---
z1r = oneT - z1r + z0r * ar / twoT - z0i * ai / twoT
z1i= -z1i+ z0r* ai / twoT + z0i* ar / twoT

xip = xi + y0r * s0 * ch0 - y0i * c0 * sh0
etap = eta + y0r * c0 * sh0 + y0i * s0 * ch0

gam[i] = atan(z1i, z1r) / d2r
k[i] = b1 / sqrt(z1r^2 + z1i^2)

# ----
s = sinh(etap)
c = max(zeroT, cos(xip))
r = sqrt(s^2 + c^2)
lon[i] = atan(s, c) / d2r

sxip = sin(xip)
sr = sxip / r

gam[i] += atan(sxip * tanh(etap), c) / d2r

tau = sr / e2m

#numit = 5, iter = 1
tau1 = sqrt(tau^2 + oneT)
b = atanh(e * tau / tau1)
sig = sinh(e * b)
taupa = sqrt(sig^2 + oneT) * tau - sig * tau1

...
end
```