Skip to content

Commit f1a0dc9

Browse files
dodomorandidodomorandi
committed
Convert drfti1 to a safe version
I am still not satisfied with this, but at least it can be a decent starting point for improvements. My points are: 1. Variable names are still unintelligible, mainly because I do not have a clear vision of the algorithm. Maybe someone better than me can be able to grasp the meaning of each variable and to git them a decent name. 2. The first loop is still unidiomatic. Everything I thought involves an additional check to exit the outer loop, and at the current state (no benchmarks) I am not sure about the best approch. Nevertheless, cast checks that were not present in the original code should be helpful to find strange edge cases, and they should have a negligible cost due to the cpu branch predictor.
1 parent f7abf1e commit f1a0dc9

File tree

1 file changed

+62
-91
lines changed

1 file changed

+62
-91
lines changed

‎fft/src/smallft.rs‎

Lines changed: 62 additions & 91 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,10 @@
88
unused_mut
99
)]
1010

11-
use std::os::raw::{c_double, c_float, c_int};
11+
use std::{
12+
convert::TryInto,
13+
os::raw::{c_double, c_float, c_int},
14+
};
1215

1316
/* *******************************************************************
1417
* *
@@ -86,109 +89,77 @@ last mod: $Id: smallft.c,v 1.19 2003/10/08 05:12:37 jm Exp $
8689
* it follows R_0, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, I_n like the
8790
* FORTRAN version
8891
*/
89-
unsafe extern "C" fn drfti1(wa: &mut [f32], ifac: &mut [i32]) {
92+
extern "C" fn drfti1(wa: &mut [f32], ifac: &mut [i32]) {
93+
const NTRYH: [i32; 4] = [4, 2, 3, 5];
94+
const TPI: f32 = 6.283_185_5;
9095

91-
static mut ntryh: [c_int; 4] = [4 as c_int, 2 as c_int, 3 as c_int, 5 as c_int];
92-
static mut tpi: c_float = 6.28318530717958648f32;
96+
let n = wa.len().try_into().unwrap();
97+
let ifac: &mut [i32; 32] = ifac.try_into().unwrap();
9398

94-
let mut n = wa.len() as i32;
95-
let mut wa = wa.as_mut_ptr();
96-
let mut ifac = ifac.as_mut_ptr();
99+
let mut nf = 0;
100+
let mut nl = n;
97101

98-
let mut arg: c_float = 0.;
99-
let mut argh: c_float = 0.;
100-
let mut argld: c_float = 0.;
101-
let mut fi: c_float = 0.;
102-
let mut ntry: c_int = 0 as c_int;
103-
let mut i: c_int = 0;
104-
let mut j: c_int = -(1 as c_int);
105-
let mut k1: c_int = 0;
106-
let mut l1: c_int = 0;
107-
let mut l2: c_int = 0;
108-
let mut ib: c_int = 0;
109-
let mut ld: c_int = 0;
110-
let mut ii: c_int = 0;
111-
let mut ip: c_int = 0;
112-
let mut is: c_int = 0;
113-
let mut nq: c_int = 0;
114-
let mut nr: c_int = 0;
115-
let mut ido: c_int = 0;
116-
let mut ipm: c_int = 0;
117-
let mut nfm1: c_int = 0;
118-
let mut nl: c_int = n;
119-
let mut nf: c_int = 0 as c_int;
120-
'c_10244: loop {
121-
j += 1;
122-
if j < 4 as c_int {
123-
ntry = ntryh[j as usize]
124-
} else {
125-
ntry += 2 as c_int
126-
}
102+
'out: for ntry in NTRYH
103+
.iter()
104+
.copied()
105+
.chain(((NTRYH.last().unwrap() + 2)..).step_by(2))
106+
{
127107
loop {
128-
nq = nl / ntry;
129-
nr = nl - ntry * nq;
130-
if nr != 0 as c_int {
108+
let nq = nl / ntry;
109+
let nr = nl - ntry * nq;
110+
if nr != 0 {
131111
break;
132112
}
113+
133114
nf += 1;
134-
*ifac.offset((nf + 1 as c_int) as isize) = ntry;
135-
nl = nq;
136-
if !(ntry != 2 as c_int) {
137-
if !(nf == 1 as c_int) {
138-
i = 1 as c_int;
139-
while i < nf {
140-
ib = nf - i + 1 as c_int;
141-
*ifac.offset((ib + 1 as c_int) as isize) = *ifac.offset(ib as isize);
142-
i += 1
143-
}
144-
*ifac.offset(2 as c_int as isize) = 2 as c_int
145-
}
115+
ifac[nf + 1] = ntry;
116+
if ntry == 2 && nf != 1 {
117+
ifac[1..nf]
118+
.rchunks_exact_mut(2)
119+
.for_each(|chunk| chunk[1] = chunk[0]);
120+
ifac[2] = 2;
146121
}
147-
if !(nl != 1 as c_int) {
148-
break 'c_10244;
122+
123+
nl = nq;
124+
if nl == 1 {
125+
break 'out;
149126
}
150127
}
151128
}
152-
*ifac.offset(0 as c_int as isize) = n;
153-
*ifac.offset(1 as c_int as isize) = nf;
154-
argh = tpi / n as c_float;
155-
is = 0 as c_int;
156-
nfm1 = nf - 1 as c_int;
157-
l1 = 1 as c_int;
158-
if nfm1 == 0 as c_int {
129+
let nf = nf;
130+
131+
ifac[0] = n;
132+
ifac[1] = nf.try_into().unwrap();
133+
if nf == 0 {
159134
return;
160135
}
161-
k1 = 0 as c_int;
162-
while k1 < nfm1 {
163-
ip = *ifac.offset((k1 + 2 as c_int) as isize);
164-
ld = 0 as c_int;
165-
l2 = l1 * ip;
166-
ido = n / l2;
167-
ipm = ip - 1 as c_int;
168-
j = 0 as c_int;
169-
while j < ipm {
170-
ld += l1;
171-
i = is;
172-
argld = ld as c_float * argh;
173-
fi = 0.0f32;
174-
ii = 2 as c_int;
175-
while ii < ido {
176-
fi += 1.0f32;
177-
arg = fi * argld;
178-
let fresh0 = i;
179-
i = i + 1;
180-
*wa.offset(fresh0 as isize) = f64::cos(arg as c_double) as c_float;
181-
let fresh1 = i;
182-
i = i + 1;
183-
*wa.offset(fresh1 as isize) = f64::sin(arg as c_double) as c_float;
184-
ii += 2 as c_int
185-
}
186-
is += ido;
187-
j += 1
188-
}
189-
l1 = l2;
190-
k1 += 1
191-
}
136+
137+
let argh = TPI / wa.len() as f32;
138+
ifac[2..=nf]
139+
.iter()
140+
.map(|&ip| ip.try_into().unwrap())
141+
.fold((0, 1), |(is, l1), ip: usize| {
142+
let l2 = l1 * ip;
143+
let ido = wa.len() / l2;
144+
let ipm = ip - 1;
145+
146+
wa[is..]
147+
.chunks_exact_mut(ido)
148+
.zip((l1..).step_by(l1).map(|ld| ld as f32 * argh))
149+
.take(ipm)
150+
.for_each(|(chunk, argld)| {
151+
chunk[..(ido - 2).try_into().unwrap_or(0)]
152+
.chunks_exact_mut(2)
153+
.fold(1f32, |fi, chunk| {
154+
let arg = (fi * argld) as f64;
155+
chunk[0] = f64::cos(arg) as f32;
156+
chunk[1] = f64::sin(arg) as f32;
157+
fi + 1.
158+
});
159+
});
160+
161+
(is + ipm * ido, l2)
162+
});
192163
}
193164
unsafe extern "C" fn fdrffti(n: usize, wsave: &mut [f32], mut ifac: &mut [i32]) {
194165
if n == 1 {

0 commit comments

Comments
 (0)