@@ -76,7 +76,6 @@ floatmin2(::Type{T}) where {T} = (twopar = 2one(T); twopar^trunc(Integer,log(flo
76
76
# NAG Ltd.
77
77
function givensAlgorithm (f:: T , g:: T ) where T<: AbstractFloat
78
78
onepar = one (T)
79
- twopar = 2 one (T)
80
79
T0 = typeof (onepar) # dimensionless
81
80
zeropar = T0 (zero (T)) # must be dimensionless
82
81
@@ -105,7 +104,7 @@ function givensAlgorithm(f::T, g::T) where T<:AbstractFloat
105
104
f1 *= safmn2
106
105
g1 *= safmn2
107
106
scalepar = max (abs (f1), abs (g1))
108
- if scalepar < safmx2u break end
107
+ if scalepar < safmx2u || count >= 20 break end
109
108
end
110
109
r = sqrt (f1* f1 + g1* g1)
111
110
cs = f1/ r
149
148
# Univ. of Colorado Denver
150
149
# NAG Ltd.
151
150
function givensAlgorithm (f:: Complex{T} , g:: Complex{T} ) where T<: AbstractFloat
152
- twopar, onepar = 2 one (T), one (T)
151
+ onepar = one (T)
153
152
T0 = typeof (onepar) # dimensionless
154
153
zeropar = T0 (zero (T)) # must be dimensionless
155
154
czero = complex (zeropar)
@@ -170,7 +169,7 @@ function givensAlgorithm(f::Complex{T}, g::Complex{T}) where T<:AbstractFloat
170
169
fs *= safmn2
171
170
gs *= safmn2
172
171
scalepar *= safmn2
173
- if scalepar < safmx2u break end
172
+ if scalepar < safmx2u || count >= 20 break end
174
173
end
175
174
elseif scalepar <= safmn2u
176
175
if g == 0
@@ -193,13 +192,13 @@ function givensAlgorithm(f::Complex{T}, g::Complex{T}) where T<:AbstractFloat
193
192
# This is a rare case: F is very small.
194
193
if f == 0
195
194
cs = zero (T)
196
- r = complex (hypot ( real (g), imag (g) ))
195
+ r = complex (abs (g ))
197
196
# do complex/real division explicitly with two real divisions
198
- d = hypot ( real (gs), imag (gs) )
197
+ d = abs (gs )
199
198
sn = complex (real (gs)/ d, - imag (gs)/ d)
200
199
return cs, sn, r
201
200
end
202
- f2s = hypot ( real (fs), imag (fs) )
201
+ f2s = abs (fs )
203
202
# g2 and g2s are accurate
204
203
# g2 is at least safmin, and g2s is at least safmn2
205
204
g2s = sqrt (g2)
@@ -214,7 +213,7 @@ function givensAlgorithm(f::Complex{T}, g::Complex{T}) where T<:AbstractFloat
214
213
# make sure abs(ff) = 1
215
214
# do complex/real division explicitly with 2 real divisions
216
215
if abs1 (f) > 1
217
- d = hypot ( real (f), imag (f) )
216
+ d = abs (f )
218
217
ff = complex (real (f)/ d, imag (f)/ d)
219
218
else
220
219
dr = safmx2* real (f)
0 commit comments