fourier: remove refactoring support types and make arrays zero-based
MODIFY
fourier/array.go
fourier/array.go
@@ -8,29 +8,10 @@ import "fmt"
// The types in array.go implement Fortran-like arrays for bootstrapping
// the implementation of the FFT functions translated from FFTPACK; they
-// are column-major one-based.
-// TODO(kortschak): Make types zero-based and remove linear array types.
+// are column-major.
// TODO(kortschak): Make a bounds-checking and non-bounds-checking version
// for testing/production.
-type oneArray []float64
-
-func (a oneArray) at(i int) float64 {
- i--
- if i < 0 || len(a) <= i {
- panic(fmt.Sprintf("out of bounds at(%d): bounds i=%d", i, len(a)))
- }
- return a[i]
-}
-
-func (a oneArray) set(i int, v float64) {
- i--
- if i < 0 || len(a) <= i {
- panic(fmt.Sprintf("out of bounds set(%d): bounds i=%d", i, len(a)))
- }
- a[i] = v
-}
-
type twoArray struct {
i, j int
jStride int
@@ -50,8 +31,6 @@ func newTwoArray(i, j int, data []float64) twoArray {
}
func (a twoArray) at(i, j int) float64 {
- i--
- j--
if i < 0 || a.i <= i || j < 0 || a.j <= j {
panic(fmt.Sprintf("out of bounds at(%d, %d): bounds i=%d, j=%d", i, j, a.i, a.j))
}
@@ -59,8 +38,6 @@ func (a twoArray) at(i, j int) float64 {
}
func (a twoArray) set(i, j int, v float64) {
- i--
- j--
if i < 0 || a.i <= i || j < 0 || a.j <= j {
panic(fmt.Sprintf("out of bounds set(%d, %d): bounds i=%d, j=%d", i, j, a.i, a.j))
}
@@ -88,9 +65,6 @@ func newThreeArray(i, j, k int, data []float64) threeArray {
}
func (a threeArray) at(i, j, k int) float64 {
- i--
- j--
- k--
if i < 0 || a.i <= i || j < 0 || a.j <= j || k < 0 || a.k <= k {
panic(fmt.Sprintf("out of bounds at(%d, %d, %d): bounds i=%d, j=%d, k=%d", i, j, k, a.i, a.j, a.k))
}
@@ -98,9 +72,6 @@ func (a threeArray) at(i, j, k int) float64 {
}
func (a threeArray) set(i, j, k int, v float64) {
- i--
- j--
- k--
if i < 0 || a.i <= i || j < 0 || a.j <= j || k < 0 || a.k <= k {
panic(fmt.Sprintf("out of bounds set(%d, %d, %d): bounds i=%d, j=%d, k=%d", i, j, k, a.i, a.j, a.k))
}

MODIFY
fourier/cfft.go
fourier/cfft.go
@@ -135,9 +135,9 @@ outer:
//
// output parameters
//
-// c for j=1,...,n
-// c(j)=the sum from k=1,...,n of
-// c(k)*exp(-i*(j-1)*(k-1)*2*pi/n)
+// c for j=0, ..., n-1
+// c[j]=the sum from k=0, ..., n-1 of
+// c[k]*exp(-i*j*k*2*pi/n)
//
// where i=sqrt(-1)
//
@@ -228,8 +228,8 @@ func cfftf1(n int, c, ch, wa []float64, ifac []int) {
}
func passf2(ido, l1 int, cc, ch, wa1 []float64) {
- cc3 := newThreeArrayZ(ido, 2, l1, cc)
- ch3 := newThreeArrayZ(ido, l1, 2, ch)
+ cc3 := newThreeArray(ido, 2, l1, cc)
+ ch3 := newThreeArray(ido, l1, 2, ch)
if ido <= 2 {
for k := 0; k < l1; k++ {
@@ -258,8 +258,8 @@ func passf3(ido, l1 int, cc, ch, wa1, wa2 []float64) {
taui = -0.866025403784439 // -sqrt(3)/2
)
- cc3 := newThreeArrayZ(ido, 3, l1, cc)
- ch3 := newThreeArrayZ(ido, l1, 3, ch)
+ cc3 := newThreeArray(ido, 3, l1, cc)
+ ch3 := newThreeArray(ido, l1, 3, ch)
if ido == 2 {
for k := 0; k < l1; k++ {
@@ -301,8 +301,8 @@ func passf3(ido, l1 int, cc, ch, wa1, wa2 []float64) {
}
func passf4(ido, l1 int, cc, ch, wa1, wa2, wa3 []float64) {
- cc3 := newThreeArrayZ(ido, 4, l1, cc)
- ch3 := newThreeArrayZ(ido, l1, 4, ch)
+ cc3 := newThreeArray(ido, 4, l1, cc)
+ ch3 := newThreeArray(ido, l1, 4, ch)
if ido == 2 {
for k := 0; k < l1; k++ {
@@ -361,8 +361,8 @@ func passf5(ido, l1 int, cc, ch, wa1, wa2, wa3, wa4 []float64) {
ti12 = -0.587785252292473
)
- cc3 := newThreeArrayZ(ido, 5, l1, cc)
- ch3 := newThreeArrayZ(ido, l1, 5, ch)
+ cc3 := newThreeArray(ido, 5, l1, cc)
+ ch3 := newThreeArray(ido, l1, 5, ch)
if ido == 2 {
for k := 0; k < l1; k++ {
@@ -436,11 +436,11 @@ func passf5(ido, l1 int, cc, ch, wa1, wa2, wa3, wa4 []float64) {
}
func passf(ido, ip, l1, idl1 int, cc, c1, c2, ch, ch2, wa []float64) (nac bool) {
- cc3 := newThreeArrayZ(ido, ip, l1, cc)
- c13 := newThreeArrayZ(ido, l1, ip, c1)
- ch3 := newThreeArrayZ(ido, l1, ip, ch)
- c2m := newTwoArrayZ(idl1, ip, c2)
- ch2m := newTwoArrayZ(idl1, ip, ch2)
+ cc3 := newThreeArray(ido, ip, l1, cc)
+ c13 := newThreeArray(ido, l1, ip, c1)
+ ch3 := newThreeArray(ido, l1, ip, ch)
+ c2m := newTwoArray(idl1, ip, c2)
+ ch2m := newTwoArray(idl1, ip, ch2)
idot := ido / 2
ipp2 := ip + 1
@@ -599,9 +599,9 @@ func passf(ido, ip, l1, idl1 int, cc, c1, c2, ch, ch2, wa []float64) (nac bool)
//
// output parameters
//
-// c for j=1,...,n
-// c(j)=the sum from k=1,...,n of
-// c(k)*exp(i*(j-1)*(k-1)*2*pi/n)
+// c for j=0, ..., n-1
+// c[j]=the sum from k=0, ..., n-1 of
+// c[k]*exp(i*j*k*2*pi/n)
//
// where i=sqrt(-1)
//
@@ -692,8 +692,8 @@ func cfftb1(n int, c, ch, wa []float64, ifac []int) {
}
func passb2(ido, l1 int, cc, ch, wa1 []float64) {
- cc3 := newThreeArrayZ(ido, 2, l1, cc)
- ch3 := newThreeArrayZ(ido, l1, 2, ch)
+ cc3 := newThreeArray(ido, 2, l1, cc)
+ ch3 := newThreeArray(ido, l1, 2, ch)
if ido <= 2 {
for k := 0; k < l1; k++ {
@@ -722,8 +722,8 @@ func passb3(ido, l1 int, cc, ch, wa1, wa2 []float64) {
taui = 0.866025403784439 // sqrt(3)/2
)
- cc3 := newThreeArrayZ(ido, 3, l1, cc)
- ch3 := newThreeArrayZ(ido, l1, 3, ch)
+ cc3 := newThreeArray(ido, 3, l1, cc)
+ ch3 := newThreeArray(ido, l1, 3, ch)
if ido == 2 {
for k := 0; k < l1; k++ {
@@ -765,8 +765,8 @@ func passb3(ido, l1 int, cc, ch, wa1, wa2 []float64) {
}
func passb4(ido, l1 int, cc, ch, wa1, wa2, wa3 []float64) {
- cc3 := newThreeArrayZ(ido, 4, l1, cc)
- ch3 := newThreeArrayZ(ido, l1, 4, ch)
+ cc3 := newThreeArray(ido, 4, l1, cc)
+ ch3 := newThreeArray(ido, l1, 4, ch)
if ido == 2 {
for k := 0; k < l1; k++ {
@@ -825,8 +825,8 @@ func passb5(ido, l1 int, cc, ch, wa1, wa2, wa3, wa4 []float64) {
ti12 = 0.587785252292473
)
- cc3 := newThreeArrayZ(ido, 5, l1, cc)
- ch3 := newThreeArrayZ(ido, l1, 5, ch)
+ cc3 := newThreeArray(ido, 5, l1, cc)
+ ch3 := newThreeArray(ido, l1, 5, ch)
if ido == 2 {
for k := 0; k < l1; k++ {
@@ -900,11 +900,11 @@ func passb5(ido, l1 int, cc, ch, wa1, wa2, wa3, wa4 []float64) {
}
func passb(ido, ip, l1, idl1 int, cc, c1, c2, ch, ch2, wa []float64) (nac bool) {
- cc3 := newThreeArrayZ(ido, ip, l1, cc)
- c13 := newThreeArrayZ(ido, l1, ip, c1)
- ch3 := newThreeArrayZ(ido, l1, ip, ch)
- c2m := newTwoArrayZ(idl1, ip, c2)
- ch2m := newTwoArrayZ(idl1, ip, ch2)
+ cc3 := newThreeArray(ido, ip, l1, cc)
+ c13 := newThreeArray(ido, l1, ip, c1)
+ ch3 := newThreeArray(ido, l1, ip, ch)
+ c2m := newTwoArray(idl1, ip, c2)
+ ch2m := newTwoArray(idl1, ip, ch2)
idot := ido / 2
ipp2 := ip + 1

MODIFY
fourier/rfft.go
fourier/rfft.go
@@ -131,18 +131,18 @@ outer:
//
// output parameters
//
-// r r(1) = the sum from i=1 to i=n of r(i)
+// r r[0] = the sum from i=0 to i=n-1 of r[i]
//
// if n is even set l=n/2, if n is odd set l = (n+1)/2
-// then for k = 2,...,l
-// r(2*k-2) = the sum from i = 1 to i = n of
-// r(i)*cos((k-1)*(i-1)*2*pi/n)
-// r(2*k-1) = the sum from i = 1 to i = n of
-// -r(i)*sin((k-1)*(i-1)*2*pi/n)
+// then for k = 1, ..., l-1
+// r[2*k-1] = the sum from i = 0 to i = n-1 of
+// r[i]*cos(k*i*2*pi/n)
+// r[2*k] = the sum from i = 0 to i = n-1 of
+// -r[i]*sin(k*i*2*pi/n)
//
// if n is even
-// r(n) = the sum from i = 1 to i = n of
-// (-1)**(i-1)*r(i)
+// r[n-1] = the sum from i = 0 to i = n-1 of
+// (-1)^i*r[i]
//
// This transform is unnormalized since a call of rfftf
// followed by a call of rfftb will multiply the input
@@ -230,8 +230,8 @@ func rfftf1(n int, c, ch, wa []float64, ifac []int) {
}
func radf2(ido, l1 int, cc, ch, wa1 []float64) {
- cc3 := newThreeArrayZ(ido, l1, 2, cc)
- ch3 := newThreeArrayZ(ido, 2, l1, ch)
+ cc3 := newThreeArray(ido, l1, 2, cc)
+ ch3 := newThreeArray(ido, 2, l1, ch)
for k := 0; k < l1; k++ {
ch3.set(0, 0, k, cc3.at(0, k, 0)+cc3.at(0, k, 1))
@@ -269,8 +269,8 @@ func radf3(ido, l1 int, cc, ch, wa1, wa2 []float64) {
taui = 0.866025403784439 // sqrt(3)/2
)
- cc3 := newThreeArrayZ(ido, l1, 3, cc)
- ch3 := newThreeArrayZ(ido, 3, l1, ch)
+ cc3 := newThreeArray(ido, l1, 3, cc)
+ ch3 := newThreeArray(ido, 3, l1, ch)
for k := 0; k < l1; k++ {
cr2 := cc3.at(0, k, 1) + cc3.at(0, k, 2)
@@ -308,8 +308,8 @@ func radf3(ido, l1 int, cc, ch, wa1, wa2 []float64) {
func radf4(ido, l1 int, cc, ch, wa1, wa2, wa3 []float64) {
const hsqt2 = math.Sqrt2 / 2
- cc3 := newThreeArrayZ(ido, l1, 4, cc)
- ch3 := newThreeArrayZ(ido, 4, l1, ch)
+ cc3 := newThreeArray(ido, l1, 4, cc)
+ ch3 := newThreeArray(ido, 4, l1, ch)
for k := 0; k < l1; k++ {
tr1 := cc3.at(0, k, 1) + cc3.at(0, k, 3)
@@ -374,8 +374,8 @@ func radf5(ido, l1 int, cc, ch, wa1, wa2, wa3, wa4 []float64) {
ti12 = 0.587785252292473
)
- cc3 := newThreeArrayZ(ido, l1, 5, cc)
- ch3 := newThreeArrayZ(ido, 5, l1, ch)
+ cc3 := newThreeArray(ido, l1, 5, cc)
+ ch3 := newThreeArray(ido, 5, l1, ch)
for k := 0; k < l1; k++ {
cr2 := cc3.at(0, k, 4) + cc3.at(0, k, 1)
@@ -638,16 +638,16 @@ func radfg(ido, ip, l1, idl1 int, cc, c1, c2, ch, ch2, wa []float64) {
//
// output parameters
//
-// r for n even and for i = 1,...,n
-// r(i) = r(1)+(-1)**(i-1)*r(n)
-// plus the sum from k=2 to k=n/2 of
-// 2.*r(2*k-2)*cos((k-1)*(i-1)*2*pi/n)
-// -2.*r(2*k-1)*sin((k-1)*(i-1)*2*pi/n)
+// r for n even and for i = 0, ..., n
+// r[i] = r[0]+(-1)^i*r[n-1]
+// plus the sum from k=1 to k=n/2-1 of
+// 2*r(2*k-1)*cos(k*i*2*pi/n)
+// -2*r(2*k)*sin(k*i*2*pi/n)
//
-// for n odd and for i = 1,...,n
-// r(i) = r(1) plus the sum from k=2 to k=(n+1)/2 of
-// 2.*r(2*k-2)*cos((k-1)*(i-1)*2*pi/n)
-// -2.*r(2*k-1)*sin((k-1)*(i-1)*2*pi/n)
+// for n odd and for i = 0, ..., n-1
+// r[i] = r[0] plus the sum from k=1 to k=(n-1)/2 of
+// 2*r(2*k-1)*cos(k*i*2*pi/n)
+// -2*r(2*k)*sin(k*i*2*pi/n)
//
// This transform is unnormalized since a call of rfftf
// followed by a call of rfftb will multiply the input
@@ -734,8 +734,8 @@ func rfftb1(n int, c, ch, wa []float64, ifac []int) {
}
func radb2(ido, l1 int, cc, ch, wa1 []float64) {
- cc3 := newThreeArrayZ(ido, 2, l1, cc)
- ch3 := newThreeArrayZ(ido, l1, 2, ch)
+ cc3 := newThreeArray(ido, 2, l1, cc)
+ ch3 := newThreeArray(ido, l1, 2, ch)
for k := 0; k < l1; k++ {
ch3.set(0, k, 0, cc3.at(0, 0, k)+cc3.at(ido-1, 1, k))
@@ -775,8 +775,8 @@ func radb3(ido, l1 int, cc, ch, wa1, wa2 []float64) {
taui = 0.866025403784439 // sqrt(3)/2
)
- cc3 := newThreeArrayZ(ido, 3, l1, cc)
- ch3 := newThreeArrayZ(ido, l1, 3, ch)
+ cc3 := newThreeArray(ido, 3, l1, cc)
+ ch3 := newThreeArray(ido, l1, 3, ch)
for k := 0; k < l1; k++ {
tr2 := cc3.at(ido-1, 1, k) + cc3.at(ido-1, 1, k)
@@ -816,8 +816,8 @@ func radb3(ido, l1 int, cc, ch, wa1, wa2 []float64) {
}
func radb4(ido, l1 int, cc, ch, wa1, wa2, wa3 []float64) {
- cc3 := newThreeArrayZ(ido, 4, l1, cc)
- ch3 := newThreeArrayZ(ido, l1, 4, ch)
+ cc3 := newThreeArray(ido, 4, l1, cc)
+ ch3 := newThreeArray(ido, l1, 4, ch)
for k := 0; k < l1; k++ {
tr1 := cc3.at(0, 0, k) - cc3.at(ido-1, 3, k)
@@ -887,8 +887,8 @@ func radb5(ido, l1 int, cc, ch, wa1, wa2, wa3, wa4 []float64) {
ti12 = 0.587785252292473
)
- cc3 := newThreeArrayZ(ido, 5, l1, cc)
- ch3 := newThreeArrayZ(ido, l1, 5, ch)
+ cc3 := newThreeArray(ido, 5, l1, cc)
+ ch3 := newThreeArray(ido, l1, 5, ch)
for k := 0; k < l1; k++ {
ti5 := cc3.at(0, 2, k) + cc3.at(0, 2, k)
@@ -953,11 +953,11 @@ func radb5(ido, l1 int, cc, ch, wa1, wa2, wa3, wa4 []float64) {
}
func radbg(ido, ip, l1, idl1 int, cc, c1, c2, ch, ch2, wa []float64) {
- cc3 := newThreeArrayZ(ido, ip, l1, cc)
- c13 := newThreeArrayZ(ido, l1, ip, c1)
- ch3 := newThreeArrayZ(ido, l1, ip, ch)
- c2m := newTwoArrayZ(idl1, ip, c2)
- ch2m := newTwoArrayZ(idl1, ip, ch2)
+ cc3 := newThreeArray(ido, ip, l1, cc)
+ c13 := newThreeArray(ido, l1, ip, c1)
+ ch3 := newThreeArray(ido, l1, ip, ch)
+ c2m := newTwoArray(idl1, ip, c2)
+ ch2m := newTwoArray(idl1, ip, ch2)
arg := 2 * math.Pi / float64(ip)
dcp := math.Cos(arg)

fourier: remove refactoring support types and make arrays zero-based