changeset 922:7e14701aef54 trunk

[svn] - replace random number generator in dithering code with SIMD-oriented Fast Mersenne Twister (SFMT). it reduces CPU load on SSE2 or AltiVec capable platform.

--- a/ChangeLog	Sat Apr 07 11:06:36 2007 -0700
+++ b/ChangeLog	Sun Apr 08 21:30:22 2007 -0700
@@ -1,3 +1,15 @@
[phone]:06:36 +0000  Yoshiki Yazawa <[e-mail]>
+  revision [1970]
+  - now wav plugin can handle remaining buffered data at the end of playing.
+  - reduce cpu usage significantly.
+  - add cleanup code.
+  - remove decode_start_mutex and decode_start_cond.
+  
+  trunk/src/wav/wav-sndfile.c |   69 +++++++++++++++++++-------------------------
+  trunk/src/wav/wav-sndfile.h |    1 
+  2 files changed, 31 insertions(+), 39 deletions(-)
+
+
 2007-04-06 19:38:11 +0000  Giacomo Lozito <[e-mail]>
   revision [1968]
   - statusicon: fix preferences window name

--- a/configure.ac	Sat Apr 07 11:06:36 2007 -0700
+++ b/configure.ac	Sun Apr 08 21:30:22 2007 -0700
@@ -250,11 +250,31 @@
     [AC_DEFINE(HAVE_ALTIVEC, 1, [Define to 1 if your system has AltiVec.])
      AC_DEFINE(HAVE_ALTIVEC_H, 1, [Define to 1 if your system has an altivec.h file.])
      AC_DEFINE(ARCH_POWERPC, 1, [Define to 1 if your system is a PowerPC.])
-     DCT64=dct64_altivec.c],
+     DCT64=dct64_altivec.c
+     SIMD_CFLAGS=-maltivec
+     AC_SUBST(SIMD_CFLAGS)],
     [DCT64=dct64.c]
 )    
 AC_SUBST(DCT64)
 
+dnl *** SSE2
+
+AC_LANG(C)
+AC_LANG_CONFTEST([AC_LANG_SOURCE([[#include <xmmintrin.h>
+int main(){
+#if __GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ >= 2)
+return 0;
+#else
+#error no vector builtins
+#endif
+}]])])
+gcc -E -dD -msse2 -o - conftest.c
+if test $? = 0 ; then
+   AC_DEFINE(HAVE_SSE2, 1, [Define to 1 if your system has SSE2])
+   SIMD_CFLAGS=-msse2
+   AC_SUBST(SIMD_CFLAGS)
+fi
+
 dnl *** MP3
 
 AC_ARG_ENABLE(mp3,

--- a/mk/rules.mk.in	Sat Apr 07 11:06:36 2007 -0700
+++ b/mk/rules.mk.in	Sun Apr 08 21:30:22 2007 -0700
@@ -357,3 +357,4 @@
 MAD_LIBS ?= @MAD_LIBS@
 IMLIB2_CFLAGS ?= @IMLIB2_CFLAGS@
 IMLIB2_LIBS ?= @IMLIB2_LIBS@
+SIMD_CFLAGS ?= @SIMD_CFLAGS@

--- a/src/madplug/Makefile	Sat Apr 07 11:06:36 2007 -0700
+++ b/src/madplug/Makefile	Sun Apr 08 21:30:22 2007 -0700
@@ -17,7 +17,7 @@
 
 OBJECTS = ${SOURCES:.c=.o}
 
-CFLAGS   += $(PICFLAGS) $(GTK_CFLAGS) $(GLIB_CFLAGS) $(PANGO_CFLAGS) $(ARCH_DEFINES) -I../../intl -I../.. -Wall
+CFLAGS   += $(PICFLAGS) $(GTK_CFLAGS) $(GLIB_CFLAGS) $(PANGO_CFLAGS) $(ARCH_DEFINES) $(SIMD_CFLAGS) -I../../intl -I../.. -Wall
 
 LDFLAGS += -Wl,-rpath=$(plugindir)
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/madplug/SFMT-alti.c	Sun Apr 08 21:30:22 2007 -0700
@@ -0,0 +1,126 @@
+/**
+ * This function represents the recursion formula in AltiVec and BIG ENDIAN.
+ * @param a a 128-bit part of the interal state array
+ * @param b a 128-bit part of the interal state array
+ * @param c a 128-bit part of the interal state array
+ * @param d a 128-bit part of the interal state array
+ * @return output
+ */
+inline static vector unsigned int vec_recursion(vector unsigned int a,
+						vector unsigned int b,
+						vector unsigned int c,
+						vector unsigned int d) {
+
+    const vector unsigned int sl1 = (vector unsigned int)(SL1, SL1, SL1, SL1);
+    const vector unsigned int sr1 = (vector unsigned int)(SR1, SR1, SR1, SR1);
+#ifdef ONLY64
+    const vector unsigned int mask = (vector unsigned int)
+    (MSK2, MSK1, MSK4, MSK3);
+    const vector unsigned char perm_sl = ALTI_SL2_PERM64;
+    const vector unsigned char perm_sr = ALTI_SR2_PERM64;
+#else
+    const vector unsigned int mask = (vector unsigned int)
+    (MSK1, MSK2, MSK3, MSK4);
+    const vector unsigned char perm_sl = ALTI_SL2_PERM;
+    const vector unsigned char perm_sr = ALTI_SR2_PERM;
+#endif
+    vector unsigned int v, w, x, y, z;
+    x = vec_perm(a, (vector unsigned int)perm_sl, perm_sl);
+    v = a;
+    y = vec_sr(b, sr1);
+    z = vec_perm(c, (vector unsigned int)perm_sr, perm_sr);
+    w = vec_sl(d, sl1);
+    z = vec_xor(z, w);
+    y = vec_and(y, mask);
+    v = vec_xor(v, x);
+    z = vec_xor(z, y);
+    z = vec_xor(z, v);
+    return z;
+}
+
+/**
+ * This function fills the internal state array with psedorandom
+ * integers.
+ */
+inline static void gen_rand_all(void) {
+    int i;
+    vector unsigned int r, r1, r2;
+
+    r1 = sfmt[N - 2].s;
+    r2 = sfmt[N - 1].s;
+    for (i = 0; i < N - POS1; i++) {
+	r = vec_recursion(sfmt[i].s, sfmt[i + POS1].s, r1, r2);
+	sfmt[i].s = r;
+	r1 = r2;
+	r2 = r;
+    }
+    for (; i < N; i++) {
+	r = vec_recursion(sfmt[i].s, sfmt[i + POS1 - N].s, r1, r2);
+	sfmt[i].s = r;
+	r1 = r2;
+	r2 = r;
+    }
+}
+
+/**
+ * This function fills the user-specified array with psedorandom
+ * integers.
+ *
+ * @param array an 128-bit array to be filled by pseudorandom numbers.  
+ * @param size number of 128-bit pesudorandom numbers to be generated.
+ */
+inline static void gen_rand_array(w128_t array[], int size) {
+    int i, j;
+    vector unsigned int r, r1, r2;
+
+    r1 = sfmt[N - 2].s;
+    r2 = sfmt[N - 1].s;
+    for (i = 0; i < N - POS1; i++) {
+	r = vec_recursion(sfmt[i].s, sfmt[i + POS1].s, r1, r2);
+	array[i].s = r;
+	r1 = r2;
+	r2 = r;
+    }
+    for (; i < N; i++) {
+	r = vec_recursion(sfmt[i].s, array[i + POS1 - N].s, r1, r2);
+	array[i].s = r;
+	r1 = r2;
+	r2 = r;
+    }
+    /* main loop */
+    for (; i < size - N; i++) {
+	r = vec_recursion(array[i - N].s, array[i + POS1 - N].s, r1, r2);
+	array[i].s = r;
+	r1 = r2;
+	r2 = r;
+    }
+    for (j = 0; j < 2 * N - size; j++) {
+	sfmt[j].s = array[j + size - N].s;
+    }
+    for (; i < size; i++) {
+	r = vec_recursion(array[i - N].s, array[i + POS1 - N].s, r1, r2);
+	array[i].s = r;
+	sfmt[j++].s = r;
+	r1 = r2;
+	r2 = r;
+    }
+}
+
+#ifndef ONLY64
+/**
+ * This function swaps high and low 32-bit of 64-bit integers in user
+ * specified array.
+ *
+ * @param array an 128-bit array to be swaped.
+ * @param size size of 128-bit array.
+ */
+inline static void swap(w128_t array[], int size) {
+    int i;
+    const vector unsigned char perm = (vector unsigned char)
+	(4, 5, 6, 7, 0, 1, 2, 3, 12, 13, 14, 15, 8, 9, 10, 11);
+
+    for (i = 0; i < size; i++) {
+	array[i].s = vec_perm(array[i].s, (vector unsigned int)perm, perm);
+    }
+}
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/madplug/SFMT-alti.h	Sun Apr 08 21:30:22 2007 -0700
@@ -0,0 +1,40 @@
+/** 
+ * @file SFMT-alti.h 
+ *
+ * @brief SIMD oriented Fast Mersenne Twister(SFMT)
+ * pseudorandom number generator
+ *
+ * @author Mutsuo Saito (Hiroshima University)
+ * @author Makoto Matsumoto (Hiroshima University)
+ *
+ * Copyright (C) 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
+ * University. All rights reserved.
+ *
+ * The new BSD License is applied to this software.
+ * see LICENSE.txt
+ */
+
+#ifndef SFMT_ALTI_H
+#define SFMT_ALTI_H
+
+union W128_T {
+    vector unsigned int s;
+    uint32_t u[4];
+};
+
+typedef union W128_T w128_t;
+
+#ifdef __GNUC__
+inline static vector unsigned int vec_recursion(vector unsigned int a,
+						vector unsigned int b,
+						vector unsigned int c,
+						vector unsigned int d)
+    __attribute__((always_inline));
+#else
+inline static vector unsigned int vec_recursion(vector unsigned int a,
+						vector unsigned int b,
+						vector unsigned int c,
+						vector unsigned int d);
+#endif
+
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/madplug/SFMT-params.h	Sun Apr 08 21:30:22 2007 -0700
@@ -0,0 +1,90 @@
+#if !defined(MEXP)
+#ifdef __GNUC__
+  #warning "MEXP is not defined. I assume MEXP is 19937."
+#endif
+  #define MEXP 19937
+#endif
+/*-----------------
+  BASIC DEFINITIONS
+  -----------------*/
+/** Mersenne Exponent. The period of the sequence 
+ *  is a multiple of 2^MEXP-1.
+ * #define MEXP 19937 */
+/** SFMT generator has an internal state array of 128-bit integers,
+ * and N is its size. */
+#define N (MEXP / 128 + 1)
+/** N32 is the size of internal state array when regarded as an array
+ * of 32-bit integers.*/
+#define N32 (N * 4)
+/** N64 is the size of internal state array when regarded as an array
+ * of 64-bit integers.*/
+#define N64 (N * 2)
+
+/*----------------------
+  the parameters of SFMT
+  following definitions are in paramsXXXX.h file.
+  ----------------------*/
+/** the pick up position of the array.
+#define POS1 122 
+*/
+
+/** the parameter of shift left as four 32-bit registers.
+#define SL1 18
+ */
+
+/** the parameter of shift left as one 128-bit register. 
+ * The 128-bit integer is shifted by (SL2 * 8) bits. 
+#define SL2 1 
+*/
+
+/** the parameter of shift right as four 32-bit registers.
+#define SR1 11
+*/
+
+/** the parameter of shift right as one 128-bit register. 
+ * The 128-bit integer is shifted by (SL2 * 8) bits. 
+#define SR2 1 
+*/
+
+/** A bitmask, used in the recursion.  These parameters are introduced
+ * to break symmetry of SIMD.
+#define MSK1 0xdfffffefU
+#define MSK2 0xddfecb7fU
+#define MSK3 0xbffaffffU
+#define MSK4 0xbffffff6U 
+*/
+
+/** These definitions are part of a 128-bit period certification vector.
+#define PARITY1	0x00000001U
+#define PARITY2	0x00000000U
+#define PARITY3	0x00000000U
+#define PARITY4	0xc98e126aU
+*/
+
+#if MEXP == 607
+  #include "SFMT-params607.h"
+#elif MEXP == 1279
+  #include "SFMT-params1279.h"
+#elif MEXP == 2281
+  #include "SFMT-params2281.h"
+#elif MEXP == 4253
+  #include "SFMT-params4253.h"
+#elif MEXP == 11213
+  #include "SFMT-params11213.h"
+#elif MEXP == 19937
+  #include "SFMT-params19937.h"
+#elif MEXP == 44497
+  #include "SFMT-params44497.h"
+#elif MEXP == 86243
+  #include "SFMT-params86243.h"
+#elif MEXP == 132049
+  #include "SFMT-params132049.h"
+#else
+#ifdef __GNUC__
+  #error "MEXP is not valid."
+  #undef MEXP
+#else
+  #undef MEXP
+#endif
+
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/madplug/SFMT-params19937.h	Sun Apr 08 21:30:22 2007 -0700
@@ -0,0 +1,22 @@
+#define POS1	122
+#define SL1	18
+#define SL2	1
+#define SR1	11
+#define SR2	1
+#define MSK1	0xdfffffefU
+#define MSK2	0xddfecb7fU
+#define MSK3	0xbffaffffU
+#define MSK4	0xbffffff6U
+#define PARITY1	0x00000001U
+#define PARITY2	0x00000000U
+#define PARITY3	0x00000000U
+#define PARITY4	0x13c9e684U
+#define ALTI_SL2_PERM \
+(vector unsigned char)(1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8)
+#define ALTI_SL2_PERM64 \
+(vector unsigned char)(1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0)
+#define ALTI_SR2_PERM \
+(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)
+#define ALTI_SR2_PERM64 \
+(vector unsigned char)(15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14)
+#define IDSTR	"SFMT-19937:122-18-1-11-1:dfffffef-ddfecb7f-bffaffff-bffffff6"

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/madplug/SFMT-sse2.c	Sun Apr 08 21:30:22 2007 -0700
@@ -0,0 +1,113 @@
+/** 
+ * @file  SFMT-sse2.c
+ * @brief SIMD oriented Fast Mersenne Twister(SFMT) for intel SSE2
+ *
+ * @author Mutsuo Saito (Hiroshima University)
+ * @author Makoto Matsumoto (Hiroshima University)
+ *
+ * @note We assume LITTLE ENDIAN in this file
+ *
+ * Copyright (C) 2006, 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
+ * University. All rights reserved.
+ *
+ * The new BSD License is applied to this software, see LICENSE.txt
+ */
+
+/**
+ * This function represents the recursion formula.
+ * @param a a 128-bit part of the interal state array
+ * @param b a 128-bit part of the interal state array
+ * @param c a 128-bit part of the interal state array
+ * @param d a 128-bit part of the interal state array
+ * @param mask 128-bit mask
+ * @return output
+ */
+inline static __m128i mm_recursion(__m128i *a, __m128i *b, 
+				   __m128i c, __m128i d, __m128i mask) {
+    __m128i v, x, y, z;
+    
+    x = _mm_load_si128(a);
+    y = _mm_srli_epi32(*b, SR1);
+    z = _mm_srli_si128(c, SR2);
+    v = _mm_slli_epi32(d, SL1);
+    z = _mm_xor_si128(z, x);
+    z = _mm_xor_si128(z, v);
+    x = _mm_slli_si128(x, SL2);
+    y = _mm_and_si128(y, mask);
+    z = _mm_xor_si128(z, x);
+    z = _mm_xor_si128(z, y);
+    return z;
+}
+
+/**
+ * This function fills the internal state array with psedorandom
+ * integers.
+ */
+inline void gen_rand_all(void) {
+    int i;
+    __m128i r, r1, r2, mask;
+    mask = _mm_set_epi32(MSK4, MSK3, MSK2, MSK1);
+
+    r1 = _mm_load_si128(&sfmt[N - 2].si);
+    r2 = _mm_load_si128(&sfmt[N - 1].si);
+    for (i = 0; i < N - POS1; i++) {
+	r = mm_recursion(&sfmt[i].si, &sfmt[i + POS1].si, r1, r2, mask);
+	_mm_store_si128(&sfmt[i].si, r);
+	r1 = r2;
+	r2 = r;
+    }
+    for (; i < N; i++) {
+	r = mm_recursion(&sfmt[i].si, &sfmt[i + POS1 - N].si, r1, r2, mask);
+	_mm_store_si128(&sfmt[i].si, r);
+	r1 = r2;
+	r2 = r;
+    }
+}
+
+/**
+ * This function fills the user-specified array with psedorandom
+ * integers.
+ *
+ * @param array an 128-bit array to be filled by pseudorandom numbers.  
+ * @param size number of 128-bit pesudorandom numbers to be generated.
+ */
+inline static void gen_rand_array(w128_t array[], int size) {
+    int i, j;
+    __m128i r, r1, r2, mask;
+    mask = _mm_set_epi32(MSK4, MSK3, MSK2, MSK1);
+
+    r1 = _mm_load_si128(&sfmt[N - 2].si);
+    r2 = _mm_load_si128(&sfmt[N - 1].si);
+    for (i = 0; i < N - POS1; i++) {
+	r = mm_recursion(&sfmt[i].si, &sfmt[i + POS1].si, r1, r2, mask);
+	_mm_store_si128(&array[i].si, r);
+	r1 = r2;
+	r2 = r;
+    }
+    for (; i < N; i++) {
+	r = mm_recursion(&sfmt[i].si, &array[i + POS1 - N].si, r1, r2, mask);
+	_mm_store_si128(&array[i].si, r);
+	r1 = r2;
+	r2 = r;
+    }
+    /* main loop */
+    for (; i < size - N; i++) {
+	r = mm_recursion(&array[i - N].si, &array[i + POS1 - N].si, r1, r2,
+			 mask);
+	_mm_store_si128(&array[i].si, r);
+	r1 = r2;
+	r2 = r;
+    }
+    for (j = 0; j < 2 * N - size; j++) {
+	r = _mm_load_si128(&array[j + size - N].si);
+	_mm_store_si128(&sfmt[j].si, r);
+    }
+    for (; i < size; i++) {
+	r = mm_recursion(&array[i - N].si, &array[i + POS1 - N].si, r1, r2,
+			 mask);
+	_mm_store_si128(&array[i].si, r);
+	_mm_store_si128(&sfmt[j++].si, r);
+	r1 = r2;
+	r2 = r;
+    }
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/madplug/SFMT-sse2.h	Sun Apr 08 21:30:22 2007 -0700
@@ -0,0 +1,36 @@
+/** 
+ * @file SFMT-sse2.h 
+ *
+ * @brief SIMD oriented Fast Mersenne Twister(SFMT)
+ * pseudorandom number generator
+ *
+ * @author Mutsuo Saito (Hiroshima University)
+ * @author Makoto Matsumoto (Hiroshima University)
+ *
+ * Copyright (C) 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
+ * University. All rights reserved.
+ *
+ * The new BSD License is applied to this software.
+ * see LICENSE.txt
+ */
+
+#ifndef SFMT_SSE2_H
+#define SFMT_SSE2_H
+#include <emmintrin.h>
+
+union W128_T {
+    __m128i si;
+    uint32_t u[4];
+};
+
+typedef union W128_T w128_t;
+
+#ifdef __GNUC__
+inline static __m128i mm_recursion(__m128i *a, __m128i *b, __m128i c,
+				   __m128i d, __m128i mask)
+    __attribute__((always_inline));
+#else
+inline static __m128i mm_recursion(__m128i *a, __m128i *b,
+				   __m128i c, __m128i d, __m128i mask);
+#endif
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/madplug/SFMT.c	Sun Apr 08 21:30:22 2007 -0700
@@ -0,0 +1,589 @@
+/** 
+ * @file  SFMT.c
+ * @brief SIMD oriented Fast Mersenne Twister(SFMT)
+ *
+ * @author Mutsuo Saito (Hiroshima University)
+ * @author Makoto Matsumoto (Hiroshima University)
+ *
+ * Copyright (C) 2006,2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
+ * University. All rights reserved.
+ *
+ * The new BSD License is applied to this software, see LICENSE.txt
+ */
+#include <string.h>
+#include <assert.h>
+#include "SFMT.h"
+#include "SFMT-params.h"
+
+#if defined(ALTIVEC)
+  #include "SFMT-alti.h"
+#elif defined(SSE2)
+  #include "SFMT-sse2.h"
+#else
+/*------------------------------------------
+  128-bit SIMD like data type for standard C
+  ------------------------------------------*/
+/** 128-bit data structure */
+struct W128_T {
+    uint32_t u[4];
+};
+
+/** 128-bit data type */
+typedef struct W128_T w128_t;
+
+#endif
+
+/*--------------------------------------
+  FILE GLOBAL VARIABLES
+  internal state, index counter and flag 
+  --------------------------------------*/
+/** the 128-bit internal state array */
+static w128_t sfmt[N];
+/** the 32bit integer pointer to the 128-bit internal state array */
+static uint32_t *psfmt32 = &sfmt[0].u[0];
+#if !defined(BIG_ENDIAN64) || defined(ONLY64)
+/** the 64bit integer pointer to the 128-bit internal state array */
+static uint64_t *psfmt64 = (uint64_t *)&sfmt[0].u[0];
+#endif
+/** index counter to the 32-bit internal state array */
+static int idx;
+/** a flag: it is 0 if and only if the internal state is not yet
+ * initialized. */
+static int initialized = 0;
+/** a parity check vector which certificate the period of 2^{MEXP} */
+static uint32_t parity[4] = {PARITY1, PARITY2, PARITY3, PARITY4};
+
+/*----------------
+  STATIC FUNCTIONS
+  ----------------*/
+inline static int idxof(int i);
+inline static void rshift128(w128_t *out,  w128_t const *in, int shift);
+inline static void lshift128(w128_t *out,  w128_t const *in, int shift);
+inline static void gen_rand_all(void);
+inline static void gen_rand_array(w128_t array[], int size);
+inline static uint32_t func1(uint32_t x);
+inline static uint32_t func2(uint32_t x);
+static void period_certification(void);
+#if defined(BIG_ENDIAN64) && !defined(ONLY64)
+inline static void swap(w128_t array[], int size);
+#endif
+
+#if defined(ALTIVEC)
+  #include "SFMT-alti.c"
+#elif defined(SSE2)
+  #include "SFMT-sse2.c"
+#endif
+
+/**
+ * This function simulate a 64-bit index of LITTLE ENDIAN 
+ * in BIG ENDIAN machine.
+ */
+#ifdef ONLY64
+inline static int idxof(int i) {
+    return i ^ 1;
+}
+#else
+inline static int idxof(int i) {
+    return i;
+}
+#endif
+/**
+ * This function simulates SIMD 128-bit right shift by the standard C.
+ * The 128-bit integer given in in is shifted by (shift * 8) bits.
+ * This function simulates the LITTLE ENDIAN SIMD.
+ * @param out the output of this function
+ * @param in the 128-bit data to be shifted
+ * @param shift the shift value
+ */
+#ifdef ONLY64
+inline static void rshift128(w128_t *out, w128_t const *in, int shift) {
+    uint64_t th, tl, oh, ol;
+
+    th = ((uint64_t)in->u[2] << 32) | ((uint64_t)in->u[3]);
+    tl = ((uint64_t)in->u[0] << 32) | ((uint64_t)in->u[1]);
+
+    oh = th >> (shift * 8);
+    ol = tl >> (shift * 8);
+    ol |= th << (64 - shift * 8);
+    out->u[0] = (uint32_t)(ol >> 32);
+    out->u[1] = (uint32_t)ol;
+    out->u[2] = (uint32_t)(oh >> 32);
+    out->u[3] = (uint32_t)oh;
+}
+#else
+inline static void rshift128(w128_t *out, w128_t const *in, int shift) {
+    uint64_t th, tl, oh, ol;
+
+    th = ((uint64_t)in->u[3] << 32) | ((uint64_t)in->u[2]);
+    tl = ((uint64_t)in->u[1] << 32) | ((uint64_t)in->u[0]);
+
+    oh = th >> (shift * 8);
+    ol = tl >> (shift * 8);
+    ol |= th << (64 - shift * 8);
+    out->u[1] = (uint32_t)(ol >> 32);
+    out->u[0] = (uint32_t)ol;
+    out->u[3] = (uint32_t)(oh >> 32);
+    out->u[2] = (uint32_t)oh;
+}
+#endif
+/**
+ * This function simulates SIMD 128-bit left shift by the standard C.
+ * The 128-bit integer given in in is shifted by (shift * 8) bits.
+ * This function simulates the LITTLE ENDIAN SIMD.
+ * @param out the output of this function
+ * @param in the 128-bit data to be shifted
+ * @param shift the shift value
+ */
+#ifdef ONLY64
+inline static void lshift128(w128_t *out, w128_t const *in, int shift) {
+    uint64_t th, tl, oh, ol;
+
+    th = ((uint64_t)in->u[2] << 32) | ((uint64_t)in->u[3]);
+    tl = ((uint64_t)in->u[0] << 32) | ((uint64_t)in->u[1]);
+
+    oh = th << (shift * 8);
+    ol = tl << (shift * 8);
+    oh |= tl >> (64 - shift * 8);
+    out->u[0] = (uint32_t)(ol >> 32);
+    out->u[1] = (uint32_t)ol;
+    out->u[2] = (uint32_t)(oh >> 32);
+    out->u[3] = (uint32_t)oh;
+}
+#else
+inline static void lshift128(w128_t *out, w128_t const *in, int shift) {
+    uint64_t th, tl, oh, ol;
+
+    th = ((uint64_t)in->u[3] << 32) | ((uint64_t)in->u[2]);
+    tl = ((uint64_t)in->u[1] << 32) | ((uint64_t)in->u[0]);
+
+    oh = th << (shift * 8);
+    ol = tl << (shift * 8);
+    oh |= tl >> (64 - shift * 8);
+    out->u[1] = (uint32_t)(ol >> 32);
+    out->u[0] = (uint32_t)ol;
+    out->u[3] = (uint32_t)(oh >> 32);
+    out->u[2] = (uint32_t)oh;
+}
+#endif
+
+/**
+ * This function represents the recursion formula.
+ * @param r output
+ * @param a a 128-bit part of the internal state array
+ * @param b a 128-bit part of the internal state array
+ * @param c a 128-bit part of the internal state array
+ * @param d a 128-bit part of the internal state array
+ */
+#ifdef ONLY64
+inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *c,
+				w128_t *d) {
+    w128_t x;
+    w128_t y;
+
+    lshift128(&x, a, SL2);
+    rshift128(&y, c, SR2);
+    r->u[0] = a->u[0] ^ x.u[0] ^ ((b->u[0] >> SR1) & MSK2) ^ y.u[0] 
+	^ (d->u[0] << SL1);
+    r->u[1] = a->u[1] ^ x.u[1] ^ ((b->u[1] >> SR1) & MSK1) ^ y.u[1] 
+	^ (d->u[1] << SL1);
+    r->u[2] = a->u[2] ^ x.u[2] ^ ((b->u[2] >> SR1) & MSK4) ^ y.u[2] 
+	^ (d->u[2] << SL1);
+    r->u[3] = a->u[3] ^ x.u[3] ^ ((b->u[3] >> SR1) & MSK3) ^ y.u[3] 
+	^ (d->u[3] << SL1);
+}
+#else
+inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *c,
+				w128_t *d) {
+    w128_t x;
+    w128_t y;
+
+    lshift128(&x, a, SL2);
+    rshift128(&y, c, SR2);
+    r->u[0] = a->u[0] ^ x.u[0] ^ ((b->u[0] >> SR1) & MSK1) ^ y.u[0] 
+	^ (d->u[0] << SL1);
+    r->u[1] = a->u[1] ^ x.u[1] ^ ((b->u[1] >> SR1) & MSK2) ^ y.u[1] 
+	^ (d->u[1] << SL1);
+    r->u[2] = a->u[2] ^ x.u[2] ^ ((b->u[2] >> SR1) & MSK3) ^ y.u[2] 
+	^ (d->u[2] << SL1);
+    r->u[3] = a->u[3] ^ x.u[3] ^ ((b->u[3] >> SR1) & MSK4) ^ y.u[3] 
+	^ (d->u[3] << SL1);
+}
+#endif
+
+#if (!defined(ALTIVEC)) && (!defined(SSE2))
+/**
+ * This function fills the internal state array with pseudorandom
+ * integers.
+ */
+inline static void gen_rand_all(void) {
+    int i;
+    w128_t *r1, *r2;
+
+    r1 = &sfmt[N - 2];
+    r2 = &sfmt[N - 1];
+    for (i = 0; i < N - POS1; i++) {
+	do_recursion(&sfmt[i], &sfmt[i], &sfmt[i + POS1], r1, r2);
+	r1 = r2;
+	r2 = &sfmt[i];
+    }
+    for (; i < N; i++) {
+	do_recursion(&sfmt[i], &sfmt[i], &sfmt[i + POS1 - N], r1, r2);
+	r1 = r2;
+	r2 = &sfmt[i];
+    }
+}
+
+/**
+ * This function fills the user-specified array with pseudorandom
+ * integers.
+ *
+ * @param array an 128-bit array to be filled by pseudorandom numbers.  
+ * @param size number of 128-bit pseudorandom numbers to be generated.
+ */
+inline static void gen_rand_array(w128_t array[], int size) {
+    int i, j;
+    w128_t *r1, *r2;
+
+    r1 = &sfmt[N - 2];
+    r2 = &sfmt[N - 1];
+    for (i = 0; i < N - POS1; i++) {
+	do_recursion(&array[i], &sfmt[i], &sfmt[i + POS1], r1, r2);
+	r1 = r2;
+	r2 = &array[i];
+    }
+    for (; i < N; i++) {
+	do_recursion(&array[i], &sfmt[i], &array[i + POS1 - N], r1, r2);
+	r1 = r2;
+	r2 = &array[i];
+    }
+    for (; i < size - N; i++) {
+	do_recursion(&array[i], &array[i - N], &array[i + POS1 - N], r1, r2);
+	r1 = r2;
+	r2 = &array[i];
+    }
+    for (j = 0; j < 2 * N - size; j++) {
+	sfmt[j] = array[j + size - N];
+    }
+    for (; i < size; i++, j++) {
+	do_recursion(&array[i], &array[i - N], &array[i + POS1 - N], r1, r2);
+	r1 = r2;
+	r2 = &array[i];
+	sfmt[j] = array[i];
+    }
+}
+#endif
+
+#if defined(BIG_ENDIAN64) && !defined(ONLY64) && !defined(ALTIVEC)
+inline static void swap(w128_t array[], int size) {
+    int i;
+    uint32_t x, y;
+
+    for (i = 0; i < size; i++) {
+	x = array[i].u[0];
+	y = array[i].u[2];
+	array[i].u[0] = array[i].u[1];
+	array[i].u[2] = array[i].u[3];
+	array[i].u[1] = x;
+	array[i].u[3] = y;
+    }
+}
+#endif
+/**
+ * This function represents a function used in the initialization
+ * by init_by_array
+ * @param x 32-bit integer
+ * @return 32-bit integer
+ */
+static uint32_t func1(uint32_t x) {
+    return (x ^ (x >> 27)) * (uint32_t)1664525UL;
+}
+
+/**
+ * This function represents a function used in the initialization
+ * by init_by_array
+ * @param x 32-bit integer
+ * @return 32-bit integer
+ */
+static uint32_t func2(uint32_t x) {
+    return (x ^ (x >> 27)) * (uint32_t)1566083941UL;
+}
+
+/**
+ * This function certificate the period of 2^{MEXP}
+ */
+static void period_certification(void) {
+    int inner = 0;
+    int i, j;
+    uint32_t work;
+
+    for (i = 0; i < 4; i++) {
+	work = psfmt32[idxof(i)] & parity[i];
+	for (j = 0; j < 32; j++) {
+	    inner ^= work & 1;
+	    work = work >> 1;
+	}
+    }
+    /* check OK */
+    if (inner == 1) {
+	return;
+    }
+    /* check NG, and modification */
+    for (i = 0; i < 4; i++) {
+	work = 1;
+	for (j = 0; j < 32; j++) {
+	    if ((work & parity[i]) != 0) {
+		psfmt32[idxof(i)] ^= work;
+		return;
+	    }
+	    work = work << 1;
+	}
+    }
+}
+
+/*----------------
+  PUBLIC FUNCTIONS
+  ----------------*/
+/**
+ * This function returns the identification string.
+ * The string shows the word size, the Mersenne exponent,
+ * and all parameters of this generator.
+ */
+char *get_idstring(void) {
+    return IDSTR;
+}
+
+/**
+ * This function returns the minimum size of array used for \b
+ * fill_array32() function.
+ * @return minimum size of array used for fill_array32() function.
+ */
+int get_min_array_size32(void) {
+    return N32;
+}
+
+/**
+ * This function returns the minimum size of array used for \b
+ * fill_array64() function.
+ * @return minimum size of array used for fill_array64() function.
+ */
+int get_min_array_size64(void) {
+    return N64;
+}
+
+#ifndef ONLY64
+/**
+ * This function generates and returns 32-bit pseudorandom number.
+ * init_gen_rand or init_by_array must be called before this function.
+ * @return 32-bit pseudorandom number
+ */
+inline uint32_t gen_rand32(void) {
+    uint32_t r;
+
+    assert(initialized);
+    if (idx >= N32) {
+	gen_rand_all();
+	idx = 0;
+    }
+    r = psfmt32[idx++];
+    return r;
+}
+#endif
+/**
+ * This function generates and returns 64-bit pseudorandom number.
+ * init_gen_rand or init_by_array must be called before this function.
+ * The function gen_rand64 should not be called after gen_rand32,
+ * unless an initialization is again executed. 
+ * @return 64-bit pseudorandom number
+ */
+inline uint64_t gen_rand64(void) {
+#if defined(BIG_ENDIAN64) && !defined(ONLY64)
+    uint32_t r1, r2;
+#else
+    uint64_t r;
+#endif
+
+    assert(initialized);
+    assert(idx % 2 == 0);
+
+    if (idx >= N32) {
+	gen_rand_all();
+	idx = 0;
+    }
+#if defined(BIG_ENDIAN64) && !defined(ONLY64)
+    r1 = psfmt32[idx];
+    r2 = psfmt32[idx + 1];
+    idx += 2;
+    return ((uint64_t)r2 << 32) | r1;
+#else
+    r = psfmt64[idx / 2];
+    idx += 2;
+    return r;
+#endif
+}
+
+#ifndef ONLY64
+/**
+ * This function generates pseudorandom 32-bit integers in the
+ * specified array[] by one call. The number of pseudorandom integers
+ * is specified by the argument size, which must be at least 624 and a
+ * multiple of four.  The generation by this function is much faster
+ * than the following gen_rand function.
+ *
+ * For initialization, init_gen_rand or init_by_array must be called
+ * before the first call of this function. This function can not be
+ * used after calling gen_rand function, without initialization.
+ *
+ * @param array an array where pseudorandom 32-bit integers are filled
+ * by this function.  The pointer to the array must be \b "aligned"
+ * (namely, must be a multiple of 16) in the SIMD version, since it
+ * refers to the address of a 128-bit integer.  In the standard C
+ * version, the pointer is arbitrary.
+ *
+ * @param size the number of 32-bit pseudorandom integers to be
+ * generated.  size must be a multiple of 4, and greater than or equal
+ * to (MEXP / 128 + 1) * 4.
+ *
+ * @note \b memalign or \b posix_memalign is available to get aligned
+ * memory. Mac OSX doesn't have these functions, but \b malloc of OSX
+ * returns the pointer to the aligned memory block.
+ */
+inline void fill_array32(uint32_t array[], int size) {
+    assert(initialized);
+    assert(idx == N32);
+    assert(size % 4 == 0);
+    assert(size >= N32);
+
+    gen_rand_array((w128_t *)array, size / 4);
+    idx = N32;
+}
+#endif
+
+/**
+ * This function generates pseudorandom 64-bit integers in the
+ * specified array[] by one call. The number of pseudorandom integers
+ * is specified by the argument size, which must be at least 312 and a
+ * multiple of two.  The generation by this function is much faster
+ * than the following gen_rand function.
+ *
+ * For initialization, init_gen_rand or init_by_array must be called
+ * before the first call of this function. This function can not be
+ * used after calling gen_rand function, without initialization.
+ *
+ * @param array an array where pseudorandom 64-bit integers are filled
+ * by this function.  The pointer to the array must be "aligned"
+ * (namely, must be a multiple of 16) in the SIMD version, since it
+ * refers to the address of a 128-bit integer.  In the standard C
+ * version, the pointer is arbitrary.
+ *
+ * @param size the number of 64-bit pseudorandom integers to be
+ * generated.  size must be a multiple of 2, and greater than or equal
+ * to (MEXP / 128 + 1) * 2
+ *
+ * @note \b memalign or \b posix_memalign is available to get aligned
+ * memory. Mac OSX doesn't have these functions, but \b malloc of OSX
+ * returns the pointer to the aligned memory block.
+ */
+inline void fill_array64(uint64_t array[], int size) {
+    assert(initialized);
+    assert(idx == N32);
+    assert(size % 2 == 0);
+    assert(size >= N64);
+
+    gen_rand_array((w128_t *)array, size / 2);
+    idx = N32;
+
+#if defined(BIG_ENDIAN64) && !defined(ONLY64)
+    swap((w128_t *)array, size /2);
+#endif
+}
+
+/**
+ * This function initializes the internal state array with a 32-bit
+ * integer seed.
+ *
+ * @param seed a 32-bit integer used as the seed.
+ */
+void init_gen_rand(uint32_t seed) {
+    int i;
+
+    psfmt32[idxof(0)] = seed;
+    for (i = 1; i < N32; i++) {
+	psfmt32[idxof(i)] = 1812433253UL * (psfmt32[idxof(i - 1)] 
+					    ^ (psfmt32[idxof(i - 1)] >> 30))
+	    + i;
+    }
+    idx = N32;
+    period_certification();
+    initialized = 1;
+}
+
+/**
+ * This function initializes the internal state array,
+ * with an array of 32-bit integers used as the seeds
+ * @param init_key the array of 32-bit integers, used as a seed.
+ * @param key_length the length of init_key.
+ */
+void init_by_array(uint32_t init_key[], int key_length) {
+    int i, j, count;
+    uint32_t r;
+    int lag;
+    int mid;
+    int size = N * 4;
+
+    if (size >= 623) {
+	lag = 11;
+    } else if (size >= 68) {
+	lag = 7;
+    } else if (size >= 39) {
+	lag = 5;
+    } else {
+	lag = 3;
+    }
+    mid = (size - lag) / 2;
+
+    memset(sfmt, 0x8b, sizeof(sfmt));
+    if (key_length + 1 > N32) {
+	count = key_length + 1;
+    } else {
+	count = N32;
+    }
+    r = func1(psfmt32[idxof(0)] ^ psfmt32[idxof(mid)] 
+	      ^ psfmt32[idxof(N32 - 1)]);
+    psfmt32[idxof(mid)] += r;
+    r += key_length;
+    psfmt32[idxof(mid + lag)] += r;
+    psfmt32[idxof(0)] = r;
+    i = 1;
+    count--;
+    for (i = 1, j = 0; (j < count) && (j < key_length); j++) {
+	r = func1(psfmt32[idxof(i)] ^ psfmt32[idxof((i + mid) % N32)] 
+		  ^ psfmt32[idxof((i + N32 - 1) % N32)]);
+	psfmt32[idxof((i + mid) % N32)] += r;
+	r += init_key[j] + i;
+	psfmt32[idxof((i + mid + lag) % N32)] += r;
+	psfmt32[idxof(i)] = r;
+	i = (i + 1) % N32;
+    }
+    for (; j < count; j++) {
+	r = func1(psfmt32[idxof(i)] ^ psfmt32[idxof((i + mid) % N32)] 
+		  ^ psfmt32[idxof((i + N32 - 1) % N32)]);
+	psfmt32[idxof((i + mid) % N32)] += r;
+	r += i;
+	psfmt32[idxof((i + mid + lag) % N32)] += r;
+	psfmt32[idxof(i)] = r;
+	i = (i + 1) % N32;
+    }
+    for (j = 0; j < N32; j++) {
+	r = func2(psfmt32[idxof(i)] + psfmt32[idxof((i + mid) % N32)] 
+		  + psfmt32[idxof((i + N32 - 1) % N32)]);
+	psfmt32[idxof((i + mid) % N32)] ^= r;
+	r -= i;
+	psfmt32[idxof((i + mid + lag) % N32)] ^= r;
+	psfmt32[idxof(i)] = r;
+	i = (i + 1) % N32;
+    }
+
+    idx = N32;
+    period_certification();
+    initialized = 1;
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/madplug/SFMT.h	Sun Apr 08 21:30:22 2007 -0700
@@ -0,0 +1,121 @@
+/** 
+ * @file SFMT.h 
+ *
+ * @brief SIMD oriented Fast Mersenne Twister(SFMT) pseudorandom
+ * number generator
+ *
+ * @author Mutsuo Saito (Hiroshima University)
+ * @author Makoto Matsumoto (Hiroshima University)
+ *
+ * Copyright (C) 2006, 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
+ * University. All rights reserved.
+ *
+ * The new BSD License is applied to this software.
+ * see LICENSE.txt
+ *
+ * @note We assume that your system has inttypes.h.  If your system
+ * doesn't have inttypes.h, you have to typedef uint32_t and uint64_t,
+ * and you have to define PRIu64 and PRIx64 in this file as follows:
+ * @verbatim
+ typedef unsigned int uint32_t
+ typedef unsigned long long uint64_t  
+ #define PRIu64 "llu"
+ #define PRIx64 "llx"
+@endverbatim
+ * uint32_t must be exactly 32-bit unsigned integer type (no more, no
+ * less), and uint64_t must be exactly 64-bit unsigned integer type.
+ * PRIu64 and PRIx64 are used for printf function to print 64-bit
+ * unsigned int and 64-bit unsigned int in hexadecimal format.
+ */
+
+#ifndef SFMT_H
+#define SFMT_H
+
+#include <stdio.h>
+
+#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L)
+  #include <inttypes.h>
+#elif defined(_MSC_VER)
+  typedef unsigned int uint32_t;
+  typedef unsigned long long uint64_t;
+  #define inline
+#else
+  #include <inttypes.h>
+  #if defined(__GNUC__)
+    #define inline __inline__
+  #endif
+#endif
+
+#ifndef PRIu64
+  #if defined(_MSC_VER)
+    #define PRIu64 "I64u"
+    #define PRIx64 "I64x"
+  #else
+    #define PRIu64 "llu"
+    #define PRIx64 "llx"
+  #endif
+#endif
+
+inline uint32_t gen_rand32(void);
+inline uint64_t gen_rand64(void);
+inline void fill_array32(uint32_t array[], int size);
+inline void fill_array64(uint64_t array[], int size);
+void init_gen_rand(uint32_t seed);
+void init_by_array(uint32_t init_key[], int key_length);
+char *get_idstring(void);
+int get_min_array_size32(void);
+int get_min_array_size64(void);
+
+/* These real versions are due to Isaku Wada */
+/** generates a random number on [0,1]-real-interval */
+inline static double to_real1(uint32_t v)
+{
+    return v * (1.0/4294967295.0); 
+    /* divided by 2^32-1 */ 
+}
+
+/** generates a random number on [0,1]-real-interval */
+inline static double genrand_real1(void)
+{
+    return to_real1(gen_rand32());
+}
+
+/** generates a random number on [0,1)-real-interval */
+inline static double to_real2(uint32_t v)
+{
+    return v * (1.0/4294967296.0); 
+    /* divided by 2^32 */
+}
+
+/** generates a random number on [0,1)-real-interval */
+inline static double genrand_real2(void)
+{
+    return to_real2(gen_rand32());
+}
+
+/** generates a random number on (0,1)-real-interval */
+inline static double to_real3(uint32_t v)
+{
+    return (((double)v) + 0.5)*(1.0/4294967296.0); 
+    /* divided by 2^32 */
+}
+
+/** generates a random number on (0,1)-real-interval */
+inline static double genrand_real3(void)
+{
+    return to_real3(gen_rand32());
+}
+/** These real versions are due to Isaku Wada */
+
+/** generates a random number on [0,1) with 53-bit resolution*/
+inline static double to_res53(uint64_t v) 
+{ 
+    return v * (1.0/18446744073709551616.0L);
+}
+
+/** generates a random number on [0,1) with 53-bit resolution*/
+inline static double genrand_res53(void) 
+{ 
+    return to_res53(gen_rand64());
+} 
+#endif

--- a/src/madplug/dither.c	Sat Apr 07 11:06:36 2007 -0700
+++ b/src/madplug/dither.c	Sun Apr 08 21:30:22 2007 -0700
@@ -1,100 +1,18 @@
-/* A C-program for MT19937: Integer     version                   */
-/*  genrand() generates one pseudorandom unsigned integer (32bit) */
-/* which is uniformly distributed among 0 to 2^32-1  for each     */
-/* call. sgenrand(seed) set initial values to the working area    */
-/* of 624 words. Before genrand(), sgenrand(seed) must be         */
-/* called once. (seed is any 32-bit integer except for 0).        */
-/*   Coded by Takuji Nishimura, considering the suggestions by    */
-/* Topher Cooper and Marc Rieffel in July-Aug. 1997.              */
-
-/* This library is free software; you can redistribute it and/or   */
-/* modify it under the terms of the GNU Library General Public     */
-/* License as published by the Free Software Foundation; either    */
-/* version 2 of the License, or (at your option) any later         */
-/* version.                                                        */
-/* This library is distributed in the hope that it will be useful, */
-/* but WITHOUT ANY WARRANTY; without even the implied warranty of  */
-/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.            */
-/* See the GNU Library General Public License for more details.    */
-/* You should have received a copy of the GNU Library General      */
-/* Public License along with this library; if not, write to the    */
-/* Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA   */
-/* 02111-1307  USA                                                 */
-
-/* Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.       */
-/* Any feedback is very welcome. For any question, comments,       */
-/* see http://www.math.keio.ac.jp/matumoto/emt.html or email       */
-/* [e-mail]                                        */
-
 #include <stdio.h>
 #include <assert.h>
-
-
-/* Period parameters */
-#define N 624
-#define M 397
-#define MATRIX_A 0x9908b0df     /* constant vector a */
-#define UPPER_MASK 0x80000000   /* most significant w-r bits */
-#define LOWER_MASK 0x7fffffff   /* least significant r bits */
+#include "../../config.h"
 
-/* Tempering parameters */
-#define TEMPERING_MASK_B 0x9d2c5680
-#define TEMPERING_MASK_C 0xefc60000
-#define TEMPERING_SHIFT_U(y)  (y >> 11)
-#define TEMPERING_SHIFT_S(y)  (y << 7)
-#define TEMPERING_SHIFT_T(y)  (y << 15)
-#define TEMPERING_SHIFT_L(y)  (y >> 18)
-
-static unsigned long mt[N];     /* the array for the state vector  */
-static int mti = N + 1;         /* mti==N+1 means mt[N] is not initialized */
-
-/* initializing the array with a NONZERO seed */
-void sgenrand(seed)
-unsigned long seed;
-{
-    /* setting initial seeds to mt[N] using         */
-    /* the generator Line 25 of Table 1 in          */
-    /* [KNUTH 1981, The Art of Computer Programming */
-    /*    Vol. 2 (2nd Ed.), pp102]                  */
-    mt[0] = seed & 0xffffffff;
-    for (mti = 1; mti < N; mti++)
-        mt[mti] = (69069 * mt[mti - 1]) & 0xffffffff;
-}
+#define MEXP 19937
 
-unsigned long genrand()
-{
-    unsigned long y;
-    static unsigned long mag01[2] = { 0x0, MATRIX_A };
-    /* mag01[x] = x * MATRIX_A  for x=0,1 */
-
-    if (mti >= N) {             /* generate N words at one time */
-        int kk;
-
-        if (mti == N + 1)       /* if sgenrand() has not been called, */
-            sgenrand(4357);     /* a default initial seed is used   */
+#ifdef HAVE_SSE2
+#define SSE2 1
+#endif
+#ifdef HAVE_ALTIVEC
+#define ALTIVEC 1
+#endif
 
-        for (kk = 0; kk < N - M; kk++) {
-            y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
-            mt[kk] = mt[kk + M] ^ (y >> 1) ^ mag01[y & 0x1];
-        }
-        for (; kk < N - 1; kk++) {
-            y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
-            mt[kk] = mt[kk + (M - N)] ^ (y >> 1) ^ mag01[y & 0x1];
-        }
-        y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
-        mt[N - 1] = mt[M - 1] ^ (y >> 1) ^ mag01[y & 0x1];
-
-        mti = 0;
-    }
-
-    y = mt[mti++];
-    y ^= TEMPERING_SHIFT_U(y);
-    y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
-    y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
-    y ^= TEMPERING_SHIFT_L(y);
-
-    return y;
-}
+#include "SFMT.h"
+#include "SFMT.c"
 
 long triangular_dither_noise(int nbits)
 {
@@ -105,7 +23,7 @@
     // see The Theory of Dithered Quantization by Robert Alexander Wannamaker
     // for complete proof of why that's optimal
 
-    long v = (genrand() / 2 - genrand() / 2);   // in ]-2^31, 2^31[
+    long v = (gen_rand32() / 2 - gen_rand32() / 2);   // in ]-2^31, 2^31[
     //int signe = (v>0) ? 1 : -1;
     long P = 1 << (32 - nbits); // the power of 2
     v /= P;

--- a/src/madplug/plugin.c	Sat Apr 07 11:06:36 2007 -0700
+++ b/src/madplug/plugin.c	Sun Apr 08 21:30:22 2007 -0700
@@ -32,6 +32,7 @@
 #include <fcntl.h>
 #include <audacious/vfs.h>
 #include <sys/stat.h>
+#include "SFMT.h"
 
 /*
  * Global variables
@@ -174,6 +175,9 @@
 
     if (!audmad_config.id3_format)
         audmad_config.id3_format = g_strdup("");
+
+    init_gen_rand(4357);
+
 }
 
 static void audmad_cleanup()