ref: 02c88f6e0f1e29884c9d084c790e09e658d47890
parent: 442bfd1d2755876f902ba9edc7e51e516b1286fa
author: qwx <qwx@sciops.net>
date: Sun May 29 20:13:42 EDT 2022
slight reshuffling
--- a/bit.c
+++ /dev/null
@@ -1,85 +1,0 @@
-#include <u.h>
-#include <libc.h>
-#include "asif.h"
-
-/* 32bit round up to next power of 2 */
-u32int
-next32pow2(u32int v)
-{- v--;
- v |= v >> 1;
- v |= v >> 2;
- v |= v >> 4;
- v |= v >> 8;
- v |= v >> 16;
- return v + 1;
-}
-
-/* find least significant set bit in 64bit integers */
-int
-lsb64(uvlong v)
-{- int c;
-
- c = 0;
- if((v & 0xffffffff) == 0){- v >>= 32;
- c += 32;
- }
- if((v & 0xffff) == 0){- v >>= 16;
- c += 16;
- }
- if((v & 0xff) == 0){- v >>= 8;
- c += 8;
- }
- if((v & 0xf) == 0){- v >>= 4;
- c += 4;
- }
- if((v & 3) == 0){- v >>= 2;
- c += 2;
- }
- if((v & 1) == 0)
- c++;
- return c;
-}
-
-/* find first set bit in 64bit integers with lookup table */
-static void
-initffs(uchar *tab, int size)
-{- int i;
-
- tab[0] = 0;
- tab[1] = 0;
- for(i=2; i<size; i++)
- tab[i] = 1 + tab[i/2];
-}
-int
-msb64(uvlong v)
-{- static uchar ffstab[256];
- int n;
-
- if(ffstab[nelem(ffstab)-1] == 0)
- initffs(ffstab, nelem(ffstab));
- if(n = v >> 56)
- return 56 + ffstab[n];
- else if(n = v >> 48)
- return 48 + ffstab[n];
- else if(n = v >> 40)
- return 40 + ffstab[n];
- else if(n = v >> 32)
- return 32 + ffstab[n];
- else if(n = v >> 24)
- return 24 + ffstab[n];
- else if(n = v >> 16)
- return 16 + ffstab[n];
- else if(n = v >> 8)
- return 8 + ffstab[n];
- else
- return ffstab[v];
-}
--- /dev/null
+++ b/dsp/fft.c
@@ -1,0 +1,103 @@
+#include <u.h>
+#include <libc.h>
+#include "asif.h"
+
+/* numerical recipes 2e, pp. 510-514
+ * numerical recipes 3e, pp. 617-620
+ * argument n here is the number of complex values,
+ * not real array size
+ * isign is either 1 (forward) or -1 (inverse) */
+
+#define SWAP(a,b) {tempr=(a); (a)=(b); (b)=tempr;}+void
+four1(double *d, int n, int isign)
+{+ int nn, mmax, m, j, istep, i;
+ double wtemp, wr, wpr, wpi, wi, θ, tempr, tempi;
+
+ if(n < 2 || n & n - 1)
+ sysfatal("non power of two");+ nn = n << 1;
+ j = 1;
+ for(i=1; i<nn; i+=2){+ if(j > i){+ SWAP(d[j-1], d[i-1]);
+ SWAP(d[j], d[i]);
+ }
+ m = n;
+ while(m >= 2 && j > m){+ j -= m;
+ m >>= 1;
+ }
+ j += m;
+ }
+ mmax = 2;
+ while(nn > mmax){+ istep = mmax << 1;
+ θ = isign * (6.28318530717959 / mmax);
+ wtemp = sin(0.5 * θ);
+ wpr = -2.0 * wtemp * wtemp;
+ wpi = sin(θ);
+ wr = 1.0;
+ wi = 0.0;
+ for(m=1; m<mmax; m+=2){+ for(i=m; i<=nn; i+=istep){+ j = i + mmax;
+ tempr = wr * d[j-1] - wi * d[j];
+ tempi = wr * d[j] + wi * d[j-1];
+ d[j-1] = d[i-1] - tempr;
+ d[j] = d[i] - tempi;
+ d[i-1] += tempr;
+ d[i] += tempi;
+ }
+ wr = (wtemp = wr) * wpr - wi * wpi + wr;
+ wi = wi * wpr + wtemp * wpi + wi;
+ }
+ mmax = istep;
+ }
+}
+
+void
+realft(double *d, int n, int isign)
+{+ int i, i1, i2, i3, i4;
+ double θ, c1, c2, h1r, h1i, h2r, h2i, wr, wi, wpr, wpi, wtemp;
+
+ c1 = 0.5;
+ θ = PI / (double)(n >> 1);
+ if(isign == 1){+ c2 = -0.5;
+ four1(d, n>>1, 1);
+ }else{+ c2 = 0.5;
+ θ = -θ;
+ }
+ wtemp = sin(0.5 * θ);
+ wpr = -2.0 * wtemp * wtemp;
+ wpi = sin(θ);
+ wr = 1.0 + wpr;
+ wi = wpi;
+ h1r = 0;
+ for(i=1; i<n>>2; i++){+ i2 = 1 + (i1 = i+i);
+ i4 = 1 + (i3 = n - i1);
+ h1r = c1 * (d[i1] + d[i3]);
+ h1i = c1 * (d[i2] - d[i4]);
+ h2r = -c2 * (d[i2] + d[i4]);
+ h2i = c2 * (d[i1] - d[i3]);
+ d[i1] = h1r + wr * h2r - wi * h2i;
+ d[i2] = h1i + wr * h2i + wi * h2r;
+ d[i3] = h1r - wr * h2r + wi * h2i;
+ d[i4] = -h1i + wr * h2i + wi * h2r;
+ wr = (wtemp = wr) * wpr - wi * wpi + wr;
+ wi = wi * wpr + wtemp * wpi + wi;
+ }
+ if(isign == 1){+ d[0] = (h1r = d[0]) + d[1];
+ d[1] = h1r - d[1];
+ }else{+ d[0] = c1 * ((h1r == d[0]) + d[1]);
+ d[1] = c1 * (h1r - d[1]);
+ four1(d, n>>1, -1);
+ }
+}
--- a/fft.c
+++ /dev/null
@@ -1,103 +1,0 @@
-#include <u.h>
-#include <libc.h>
-#include "asif.h"
-
-/* numerical recipes 2e, pp. 510-514
- * numerical recipes 3e, pp. 617-620
- * argument n here is the number of complex values,
- * not real array size
- * isign is either 1 (forward) or -1 (inverse) */
-
-#define SWAP(a,b) {tempr=(a); (a)=(b); (b)=tempr;}-void
-four1(double *d, int n, int isign)
-{- int nn, mmax, m, j, istep, i;
- double wtemp, wr, wpr, wpi, wi, θ, tempr, tempi;
-
- if(n < 2 || n & n - 1)
- sysfatal("non power of two");- nn = n << 1;
- j = 1;
- for(i=1; i<nn; i+=2){- if(j > i){- SWAP(d[j-1], d[i-1]);
- SWAP(d[j], d[i]);
- }
- m = n;
- while(m >= 2 && j > m){- j -= m;
- m >>= 1;
- }
- j += m;
- }
- mmax = 2;
- while(nn > mmax){- istep = mmax << 1;
- θ = isign * (6.28318530717959 / mmax);
- wtemp = sin(0.5 * θ);
- wpr = -2.0 * wtemp * wtemp;
- wpi = sin(θ);
- wr = 1.0;
- wi = 0.0;
- for(m=1; m<mmax; m+=2){- for(i=m; i<=nn; i+=istep){- j = i + mmax;
- tempr = wr * d[j-1] - wi * d[j];
- tempi = wr * d[j] + wi * d[j-1];
- d[j-1] = d[i-1] - tempr;
- d[j] = d[i] - tempi;
- d[i-1] += tempr;
- d[i] += tempi;
- }
- wr = (wtemp = wr) * wpr - wi * wpi + wr;
- wi = wi * wpr + wtemp * wpi + wi;
- }
- mmax = istep;
- }
-}
-
-void
-realft(double *d, int n, int isign)
-{- int i, i1, i2, i3, i4;
- double θ, c1, c2, h1r, h1i, h2r, h2i, wr, wi, wpr, wpi, wtemp;
-
- c1 = 0.5;
- θ = PI / (double)(n >> 1);
- if(isign == 1){- c2 = -0.5;
- four1(d, n>>1, 1);
- }else{- c2 = 0.5;
- θ = -θ;
- }
- wtemp = sin(0.5 * θ);
- wpr = -2.0 * wtemp * wtemp;
- wpi = sin(θ);
- wr = 1.0 + wpr;
- wi = wpi;
- h1r = 0;
- for(i=1; i<n>>2; i++){- i2 = 1 + (i1 = i+i);
- i4 = 1 + (i3 = n - i1);
- h1r = c1 * (d[i1] + d[i3]);
- h1i = c1 * (d[i2] - d[i4]);
- h2r = -c2 * (d[i2] + d[i4]);
- h2i = c2 * (d[i1] - d[i3]);
- d[i1] = h1r + wr * h2r - wi * h2i;
- d[i2] = h1i + wr * h2i + wi * h2r;
- d[i3] = h1r - wr * h2r + wi * h2i;
- d[i4] = -h1i + wr * h2i + wi * h2r;
- wr = (wtemp = wr) * wpr - wi * wpi + wr;
- wi = wi * wpr + wtemp * wpi + wi;
- }
- if(isign == 1){- d[0] = (h1r = d[0]) + d[1];
- d[1] = h1r - d[1];
- }else{- d[0] = c1 * ((h1r == d[0]) + d[1]);
- d[1] = c1 * (h1r - d[1]);
- four1(d, n>>1, -1);
- }
-}
--- /dev/null
+++ b/math/bit.c
@@ -1,0 +1,85 @@
+#include <u.h>
+#include <libc.h>
+#include "asif.h"
+
+/* 32bit round up to next power of 2 */
+u32int
+next32pow2(u32int v)
+{+ v--;
+ v |= v >> 1;
+ v |= v >> 2;
+ v |= v >> 4;
+ v |= v >> 8;
+ v |= v >> 16;
+ return v + 1;
+}
+
+/* find least significant set bit in 64bit integers */
+int
+lsb64(uvlong v)
+{+ int c;
+
+ c = 0;
+ if((v & 0xffffffff) == 0){+ v >>= 32;
+ c += 32;
+ }
+ if((v & 0xffff) == 0){+ v >>= 16;
+ c += 16;
+ }
+ if((v & 0xff) == 0){+ v >>= 8;
+ c += 8;
+ }
+ if((v & 0xf) == 0){+ v >>= 4;
+ c += 4;
+ }
+ if((v & 3) == 0){+ v >>= 2;
+ c += 2;
+ }
+ if((v & 1) == 0)
+ c++;
+ return c;
+}
+
+/* find first set bit in 64bit integers with lookup table */
+static void
+initffs(uchar *tab, int size)
+{+ int i;
+
+ tab[0] = 0;
+ tab[1] = 0;
+ for(i=2; i<size; i++)
+ tab[i] = 1 + tab[i/2];
+}
+int
+msb64(uvlong v)
+{+ static uchar ffstab[256];
+ int n;
+
+ if(ffstab[nelem(ffstab)-1] == 0)
+ initffs(ffstab, nelem(ffstab));
+ if(n = v >> 56)
+ return 56 + ffstab[n];
+ else if(n = v >> 48)
+ return 48 + ffstab[n];
+ else if(n = v >> 40)
+ return 40 + ffstab[n];
+ else if(n = v >> 32)
+ return 32 + ffstab[n];
+ else if(n = v >> 24)
+ return 24 + ffstab[n];
+ else if(n = v >> 16)
+ return 16 + ffstab[n];
+ else if(n = v >> 8)
+ return 8 + ffstab[n];
+ else
+ return ffstab[v];
+}
--- /dev/null
+++ b/math/mod.c
@@ -1,0 +1,28 @@
+#include <u.h>
+#include <libc.h>
+#include "asif.h"
+
+/* modular exponentiation by repeated binary square and multiply,
+ * (addition chaining), right-to-left scan.
+ * assumptions: mod > 0, base >= 0, exp >= 0
+ * if base ≥ 0x10000, base² in the first iteration will overflow.
+ * if mod > 0x10000, base can be ≥ 0x10000 in subsequent iterations.
+ */
+int
+modpow(int base, int exp, int mod)
+{+ int r;
+
+ if(base == 0)
+ return 0;
+ assert(base < 0x10000);
+ assert(mod <= 0x10000);
+ r = 1;
+ while(exp > 0){+ if(exp & 1)
+ r = r * base % mod;
+ exp >>= 1;
+ base = base * base % mod;
+ }
+ return r;
+}
--- a/mod.c
+++ /dev/null
@@ -1,28 +1,0 @@
-#include <u.h>
-#include <libc.h>
-#include "asif.h"
-
-/* modular exponentiation by repeated binary square and multiply,
- * (addition chaining), right-to-left scan.
- * assumptions: mod > 0, base >= 0, exp >= 0
- * if base ≥ 0x10000, base² in the first iteration will overflow.
- * if mod > 0x10000, base can be ≥ 0x10000 in subsequent iterations.
- */
-int
-modpow(int base, int exp, int mod)
-{- int r;
-
- if(base == 0)
- return 0;
- assert(base < 0x10000);
- assert(mod <= 0x10000);
- r = 1;
- while(exp > 0){- if(exp & 1)
- r = r * base % mod;
- exp >>= 1;
- base = base * base % mod;
- }
- return r;
-}