blob: 54bf5875629fcbd15984e10d35dddae83fea18b0 [file] [log] [blame]
 // Copyright ©2015 The Gonum Authors. All rights reserved. // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. package gonum import "math" // Dlasq6 computes one dqd transform in ping-pong form with protection against // overflow and underflow. z has length at least 4*(n0+1) and holds the qd array. // i0 is the zero-based first index. // n0 is the zero-based last index. // // Dlasq6 is an internal routine. It is exported for testing purposes. func (impl Implementation) Dlasq6(i0, n0 int, z []float64, pp int) (dmin, dmin1, dmin2, dn, dnm1, dnm2 float64) { switch { case i0 < 0: panic(i0LT0) case n0 < 0: panic(n0LT0) case len(z) < 4*n0: panic(shortZ) case pp != 0 && pp != 1: panic(badPp) } if n0-i0-1 <= 0 { return dmin, dmin1, dmin2, dn, dnm1, dnm2 } safmin := dlamchS j4 := 4*(i0+1) + pp - 4 // -4 rather than -3 for zero indexing emin := z[j4+4] d := z[j4] dmin = d if pp == 0 { for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 { j4 := j4loop - 1 // Translate back to zero-indexed. z[j4-2] = d + z[j4-1] if z[j4-2] == 0 { z[j4] = 0 d = z[j4+1] dmin = d emin = 0 } else if safmin*z[j4+1] < z[j4-2] && safmin*z[j4-2] < z[j4+1] { tmp := z[j4+1] / z[j4-2] z[j4] = z[j4-1] * tmp d *= tmp } else { z[j4] = z[j4+1] * (z[j4-1] / z[j4-2]) d = z[j4+1] * (d / z[j4-2]) } dmin = math.Min(dmin, d) emin = math.Min(emin, z[j4]) } } else { for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 { j4 := j4loop - 1 z[j4-3] = d + z[j4] if z[j4-3] == 0 { z[j4-1] = 0 d = z[j4+2] dmin = d emin = 0 } else if safmin*z[j4+2] < z[j4-3] && safmin*z[j4-3] < z[j4+2] { tmp := z[j4+2] / z[j4-3] z[j4-1] = z[j4] * tmp d *= tmp } else { z[j4-1] = z[j4+2] * (z[j4] / z[j4-3]) d = z[j4+2] * (d / z[j4-3]) } dmin = math.Min(dmin, d) emin = math.Min(emin, z[j4-1]) } } // Unroll last two steps. dnm2 = d dmin2 = dmin j4 = 4*(n0-1) - pp - 1 j4p2 := j4 + 2*pp - 1 z[j4-2] = dnm2 + z[j4p2] if z[j4-2] == 0 { z[j4] = 0 dnm1 = z[j4p2+2] dmin = dnm1 emin = 0 } else if safmin*z[j4p2+2] < z[j4-2] && safmin*z[j4-2] < z[j4p2+2] { tmp := z[j4p2+2] / z[j4-2] z[j4] = z[j4p2] * tmp dnm1 = dnm2 * tmp } else { z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2]) dnm1 = z[j4p2+2] * (dnm2 / z[j4-2]) } dmin = math.Min(dmin, dnm1) dmin1 = dmin j4 += 4 j4p2 = j4 + 2*pp - 1 z[j4-2] = dnm1 + z[j4p2] if z[j4-2] == 0 { z[j4] = 0 dn = z[j4p2+2] dmin = dn emin = 0 } else if safmin*z[j4p2+2] < z[j4-2] && safmin*z[j4-2] < z[j4p2+2] { tmp := z[j4p2+2] / z[j4-2] z[j4] = z[j4p2] * tmp dn = dnm1 * tmp } else { z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2]) dn = z[j4p2+2] * (dnm1 / z[j4-2]) } dmin = math.Min(dmin, dn) z[j4+2] = dn z[4*(n0+1)-pp-1] = emin return dmin, dmin1, dmin2, dn, dnm1, dnm2 }