[loris] Add C++ LorisPeakDetection.

4 files changed, 295 insertions, 3631 deletions

diff --git a/src/loris/Analyzer.C b/src/loris/Analyzer.C
index 5b4996c..deaf178 100644
--- a/src/loris/Analyzer.C
+++ b/src/loris/Analyzer.C
@@ -89,78 +89,6 @@ static double compare2nd( const Pair & p1, const Pair & p2 )
     return p1.second < p2.second;
 }
 
-// ---------------------------------------------------------------------------
-//  LinearEnvelopeBuilder
-// ---------------------------------------------------------------------------
-//  Base class for envelope builders that add a point (possibly) at each
-//  analysis frame. 
-//
-//  TODO: make a dictionary of these things and allow clients to add their
-//  own envelope builders and builder functions, and retrieve them after
-//  analysis.
-class LinearEnvelopeBuilder
-{
-public:
-    virtual ~LinearEnvelopeBuilder( void ) {}
-    virtual LinearEnvelopeBuilder * clone( void ) const = 0;
-    virtual void build( const Peaks & peaks, double frameTime ) = 0;
-    
-    const LinearEnvelope & envelope( void ) const { return mEnvelope; }
-    
-    //  reset (clear) envelope, override if necesssary:
-    virtual void reset( void ) { mEnvelope.clear(); }
-    
-protected:
-
-    LinearEnvelope mEnvelope;   //  build this
-};
-
-// ---------------------------------------------------------------------------
-//  FundamentalBuilder - for constructing an F0 envelope during analysis
-// ---------------------------------------------------------------------------
-class FundamentalBuilder : public LinearEnvelopeBuilder
-{
-    std::auto_ptr< Envelope > mFminEnv;
-    std::auto_ptr< Envelope > mFmaxEnv; 
-    
-    double mAmpThresh, mFreqThresh;
-    
-    std::vector< double > amplitudes, frequencies;
-    
-    const double mMinConfidence;    // 0.9, this could be made a parameter, 
-                                    // or raised to make estimates smoother
-    
-public:
-    FundamentalBuilder( double fmin, double fmax, double threshDb = -60, double threshHz = 8000 ) :
-        mFminEnv( new LinearEnvelope( fmin ) ), 
-        mFmaxEnv( new LinearEnvelope( fmax ) ), 
-        mAmpThresh( std::pow( 10., 0.05*(threshDb) ) ),
-        mFreqThresh( threshHz ),
-        mMinConfidence( 0.9 )
-        {}
-
-    FundamentalBuilder( const Envelope & fmin, const Envelope & fmax, 
-    					double threshDb = -60, double threshHz = 8000 ) :
-        mFminEnv( fmin.clone() ), 
-        mFmaxEnv( fmax.clone() ), 
-        mAmpThresh( std::pow( 10., 0.05*(threshDb) ) ),
-        mFreqThresh( threshHz ),
-        mMinConfidence( 0.9 )
-        {}
-        	
-    FundamentalBuilder( const FundamentalBuilder & rhs ) :
-        mFminEnv( rhs.mFminEnv->clone() ), 
-        mFmaxEnv( rhs.mFmaxEnv->clone() ), 
-        mAmpThresh( rhs.mAmpThresh ),
-        mFreqThresh( rhs.mFreqThresh ),
-        mMinConfidence( rhs.mMinConfidence )
-        {}
-    
-        
-	FundamentalBuilder * clone( void ) const { return new FundamentalBuilder(*this); }
-	
-    void build( const Peaks & peaks, double frameTime );
-};
 
 // ---------------------------------------------------------------------------
 //  FundamentalBuilder::build
@@ -198,19 +126,6 @@ void FundamentalBuilder::build( const Peaks & peaks, double frameTime )
     
 }
 
-// ---------------------------------------------------------------------------
-//  AmpEnvBuilder - for constructing an amplitude envelope during analysis
-// ---------------------------------------------------------------------------
-class AmpEnvBuilder : public LinearEnvelopeBuilder
-{
-public:
-    AmpEnvBuilder( void ) {}
-		
-	AmpEnvBuilder * clone( void ) const { return new AmpEnvBuilder(*this); }
-	
-    void build( const Peaks & peaks, double frameTime );
-
-};
 
 // ---------------------------------------------------------------------------
 //  AmpEnvBuilder::build
@@ -732,11 +647,6 @@ Analyzer::analyze( const double * bufBegin, const double * bufEnd, double srate,
     }
 }
 
-void 
-Analyzer::analyze_peaks( const std::vector<double> & vec, double srate )
-{
-    printf("analyze_peaks\n");
-}
 // -- parameter access --
 
 // ---------------------------------------------------------------------------
diff --git a/src/loris/Analyzer.h b/src/loris/Analyzer.h
index 58ece87..9effdc5 100644
--- a/src/loris/Analyzer.h
+++ b/src/loris/Analyzer.h
@@ -1,8 +1,8 @@
 #ifndef INCLUDE_ANALYZER_H
 #define INCLUDE_ANALYZER_H
 /*
- * This is the Loris C++ Class Library, implementing analysis, 
- * manipulation, and synthesis of digitized sounds using the Reassigned 
+ * This is the Loris C++ Class Library, implementing analysis,
+ * manipulation, and synthesis of digitized sounds using the Reassigned
  * Bandwidth-Enhanced Additive Sound Model.
  *
  * Loris is Copyright (c) 1999-2010 by Kelly Fitz and Lippold Haken
@@ -34,6 +34,7 @@
  */
 #include <memory>
 #include <vector>
+#include <cmath>
 #include "LinearEnvelope.h"
 #include "Partial.h"
 #include "PartialList.h"
@@ -43,42 +44,127 @@
 namespace Loris {
 
 class Envelope;
-class LinearEnvelopeBuilder;
-//  oooo, this is nasty, need to fix it!
-// class SpectralPeak;
-// typedef std::vector< SpectralPeak > Peaks;
+
+// ---------------------------------------------------------------------------
+//  LinearEnvelopeBuilder
+// ---------------------------------------------------------------------------
+//  Base class for envelope builders that add a point (possibly) at each
+//  analysis frame.
+//
+//  TODO: make a dictionary of these things and allow clients to add their
+//  own envelope builders and builder functions, and retrieve them after
+//  analysis.
+class LinearEnvelopeBuilder
+{
+public:
+    virtual ~LinearEnvelopeBuilder( void ) {}
+    virtual LinearEnvelopeBuilder * clone( void ) const = 0;
+    virtual void build( const Peaks & peaks, double frameTime ) = 0;
+
+    const LinearEnvelope & envelope( void ) const { return mEnvelope; }
+
+    //  reset (clear) envelope, override if necesssary:
+    virtual void reset( void ) { mEnvelope.clear(); }
+
+protected:
+
+    LinearEnvelope mEnvelope;   //  build this
+};
+
+
+// ---------------------------------------------------------------------------
+//  FundamentalBuilder - for constructing an F0 envelope during analysis
+// ---------------------------------------------------------------------------
+class FundamentalBuilder : public LinearEnvelopeBuilder
+{
+    std::auto_ptr< Envelope > mFminEnv;
+    std::auto_ptr< Envelope > mFmaxEnv;
+
+    double mAmpThresh, mFreqThresh;
+
+    std::vector< double > amplitudes, frequencies;
+
+    const double mMinConfidence;    // 0.9, this could be made a parameter,
+                                    // or raised to make estimates smoother
+
+public:
+    FundamentalBuilder( double fmin, double fmax, double threshDb = -60, double threshHz = 8000 ) :
+        mFminEnv( new LinearEnvelope( fmin ) ),
+        mFmaxEnv( new LinearEnvelope( fmax ) ),
+        mAmpThresh( std::pow( 10., 0.05*(threshDb) ) ),
+        mFreqThresh( threshHz ),
+        mMinConfidence( 0.9 )
+        {}
+
+    FundamentalBuilder( const Envelope & fmin, const Envelope & fmax,
+    					double threshDb = -60, double threshHz = 8000 ) :
+        mFminEnv( fmin.clone() ),
+        mFmaxEnv( fmax.clone() ),
+        mAmpThresh( std::pow( 10., 0.05*(threshDb) ) ),
+        mFreqThresh( threshHz ),
+        mMinConfidence( 0.9 )
+        {}
+
+    FundamentalBuilder( const FundamentalBuilder & rhs ) :
+        mFminEnv( rhs.mFminEnv->clone() ),
+        mFmaxEnv( rhs.mFmaxEnv->clone() ),
+        mAmpThresh( rhs.mAmpThresh ),
+        mFreqThresh( rhs.mFreqThresh ),
+        mMinConfidence( rhs.mMinConfidence )
+        {}
+
+
+	FundamentalBuilder * clone( void ) const { return new FundamentalBuilder(*this); }
+
+    void build( const Peaks & peaks, double frameTime );
+};
+
+
+// ---------------------------------------------------------------------------
+//  AmpEnvBuilder - for constructing an amplitude envelope during analysis
+// ---------------------------------------------------------------------------
+class AmpEnvBuilder : public LinearEnvelopeBuilder
+{
+public:
+    AmpEnvBuilder( void ) {}
+
+	AmpEnvBuilder * clone( void ) const { return new AmpEnvBuilder(*this); }
+
+    void build( const Peaks & peaks, double frameTime );
+
+};
 
 // ---------------------------------------------------------------------------
 //  class Analyzer
 //
 //! Class Analyzer represents a configuration of parameters for
 //! performing Reassigned Bandwidth-Enhanced Additive Analysis
-//! of sampled sounds. The analysis process yields a collection 
+//! of sampled sounds. The analysis process yields a collection
 //! of Partials, each having a trio of synchronous, non-uniformly-
-//! sampled breakpoint envelopes representing the time-varying 
+//! sampled breakpoint envelopes representing the time-varying
 //! frequency, amplitude, and noisiness of a single bandwidth-
 //! enhanced sinusoid. These Partials are accumulated in the
 //! Analyzer.
 //!
-//! The core analysis parameter is the frequency resolution, the 
-//! minimum instantaneous frequency spacing between partials. Most 
-//! other parameters are initially configured according to this 
+//! The core analysis parameter is the frequency resolution, the
+//! minimum instantaneous frequency spacing between partials. Most
+//! other parameters are initially configured according to this
 //! parameter (and the analysis window width, if specified).
 //! Subsequent parameter mutations are independent.
-//! 
+//!
 //! Bandwidth enhancement:
 //! Two different strategies are available for computing bandwidth
-//! (or noisiness) envelope: 
-//! 
+//! (or noisiness) envelope:
+//!
 //! One strategy is to construct bandwidth envelopes during analysis
-//! by associating residual energy in the spectrum (after peak 
-//! extraction) with the selected spectral peaks that are used 
+//! by associating residual energy in the spectrum (after peak
+//! extraction) with the selected spectral peaks that are used
 //! to construct Partials. This is the original bandwidth enhancement
 //! algorithm, and bandwidth envelopes constructed in this way may
-//! be suitable for use in bandwidth-enhanced synthesis. 
-//! 
-//! Another stategy is to construct bandwidth envelopes during 
-//! analysis by storing the mixed derivative of short-time phase, 
+//! be suitable for use in bandwidth-enhanced synthesis.
+//!
+//! Another stategy is to construct bandwidth envelopes during
+//! analysis by storing the mixed derivative of short-time phase,
 //! scaled and shifted so that a value of 0 corresponds
 //! to a pure sinusoid, and a value of 1 corresponds to a
 //! bandwidth-enhanced sinusoid with maximal energy spread
@@ -86,17 +172,17 @@ class LinearEnvelopeBuilder;
 //! are not suitable for bandwidth-enhanced synthesis, be sure
 //! to set the bandwidth to 0, or to disable bandwidth enhancement
 //! before rendering.
-//! 
+//!
 //! The Analyzer may be configured to use either of these two
 //! strategies for bandwidth-enhanced analysis, or to construct
 //! no bandwidth envelopes at all. If unspecified, the default
 //! Analyzer configuration uses spectral residue to construct
 //! bandwidth envelopes.
-//! 
+//!
 //! \sa storeResidueBandwidth, storeConvergenceBandwidth, storeNoBandwidth
-//! 
-//! For more information about Reassigned Bandwidth-Enhanced 
-//! Analysis and the Reassigned Bandwidth-Enhanced Additive Sound 
+//!
+//! For more information about Reassigned Bandwidth-Enhanced
+//! Analysis and the Reassigned Bandwidth-Enhanced Additive Sound
 //! Model, refer to the Loris website: www.cerlsoundgroup.org/Loris/.
 //
 class Analyzer
@@ -106,74 +192,74 @@ public:
 
 //  -- construction --
 
-    //! Construct a new Analyzer configured with the given  
-    //! frequency resolution (minimum instantaneous frequency   
-    //! difference between Partials). All other Analyzer parameters     
-    //! are computed from the specified frequency resolution.   
-    //! 
+    //! Construct a new Analyzer configured with the given
+    //! frequency resolution (minimum instantaneous frequency
+    //! difference between Partials). All other Analyzer parameters
+    //! are computed from the specified frequency resolution.
+    //!
     //! \param resolutionHz is the frequency resolution in Hz.
     explicit Analyzer( double resolutionHz );
-    
-    //! Construct a new Analyzer configured with the given  
-    //! frequency resolution (minimum instantaneous frequency   
+
+    //! Construct a new Analyzer configured with the given
+    //! frequency resolution (minimum instantaneous frequency
     //! difference between Partials) and analysis window width
-    //! (main lobe, zero-to-zero). All other Analyzer parameters    
-    //! are computed from the specified resolution and window width.    
-    //! 
+    //! (main lobe, zero-to-zero). All other Analyzer parameters
+    //! are computed from the specified resolution and window width.
+    //!
     //! \param resolutionHz is the frequency resolution in Hz.
     //! \param windowWidthHz is the main lobe width of the Kaiser
     //! analysis window in Hz.
     Analyzer( double resolutionHz, double windowWidthHz );
 
     //! Construct a new Analyzer configured with the given time-varying
-    //! frequency resolution (minimum instantaneous frequency   
+    //! frequency resolution (minimum instantaneous frequency
     //! difference between Partials) and analysis window width
-    //! (main lobe, zero-to-zero). All other Analyzer parameters    
-    //! are computed from the specified resolution and window width.    
-    //! 
+    //! (main lobe, zero-to-zero). All other Analyzer parameters
+    //! are computed from the specified resolution and window width.
+    //!
     //! \param resolutionHz is the frequency resolution in Hz.
     //! \param windowWidthHz is the main lobe width of the Kaiser
     //! analysis window in Hz.
     Analyzer( const Envelope & resolutionEnv, double windowWidthHz );
 
     //! Construct  a new Analyzer having identical
-    //! parameter configuration to another Analyzer. 
-    //! The list of collected Partials is not copied.       
-    //! 
-    //! \param other is the Analyzer to copy.   
+    //! parameter configuration to another Analyzer.
+    //! The list of collected Partials is not copied.
+    //!
+    //! \param other is the Analyzer to copy.
     Analyzer( const Analyzer & other );
 
     //! Destroy this Analyzer.
     ~Analyzer( void );
 
     //! Construct  a new Analyzer having identical
-    //! parameter configuration to another Analyzer. 
-    //! The list of collected Partials is not copied.       
-    //! 
-    //! \param rhs is the Analyzer to copy. 
+    //! parameter configuration to another Analyzer.
+    //! The list of collected Partials is not copied.
+    //!
+    //! \param rhs is the Analyzer to copy.
     Analyzer & operator=( const Analyzer & rhs );
 
 //  -- configuration --
 
-    //! Configure this Analyzer with the given frequency resolution 
-    //! (minimum instantaneous frequency difference between Partials, 
-    //! in Hz). All other Analyzer parameters are (re-)computed from the 
+    //! Configure this Analyzer with the given frequency resolution
+    //! (minimum instantaneous frequency difference between Partials,
+    //! in Hz). All other Analyzer parameters are (re-)computed from the
     //! frequency resolution, including the window width, which is
-    //! twice the resolution.      
-    //! 
+    //! twice the resolution.
+    //!
     //! \param resolutionHz is the frequency resolution in Hz.
     void configure( double resolutionHz );
 
-    //! Configure this Analyzer with the given frequency resolution 
+    //! Configure this Analyzer with the given frequency resolution
     //! (minimum instantaneous frequency difference between Partials)
-    //! and analysis window width (main lobe, zero-to-zero, in Hz). 
-    //! All other Analyzer parameters are (re-)computed from the 
-    //! frequency resolution and window width.      
-    //! 
+    //! and analysis window width (main lobe, zero-to-zero, in Hz).
+    //! All other Analyzer parameters are (re-)computed from the
+    //! frequency resolution and window width.
+    //!
     //! \param resolutionHz is the frequency resolution in Hz.
     //! \param windowWidthHz is the main lobe width of the Kaiser
     //! analysis window in Hz.
-    //!     
+    //!
     //! There are three categories of analysis parameters:
     //! - the resolution, and params that are usually related to (or
     //! identical to) the resolution (frequency floor and drift)
@@ -182,17 +268,17 @@ public:
     //! - independent parameters (bw region width and amp floor)
     void configure( double resolutionHz, double windowWidthHz );
 
-	//! Configure this Analyzer with the given time-varying frequency resolution 
+	//! Configure this Analyzer with the given time-varying frequency resolution
 	//! (minimum instantaneous frequency difference between Partials)
-	//! and analysis window width (main lobe, zero-to-zero, in Hz). 
-	//! All other Analyzer parameters are (re-)computed from the 
-	//! frequency resolution and window width.      
-	//! 
-	//! \param resolutionEnv is the time-varying frequency resolution 
+	//! and analysis window width (main lobe, zero-to-zero, in Hz).
+	//! All other Analyzer parameters are (re-)computed from the
+	//! frequency resolution and window width.
+	//!
+	//! \param resolutionEnv is the time-varying frequency resolution
 	//!	in Hz.
 	//! \param windowWidthHz is the main lobe width of the Kaiser
 	//! analysis window in Hz.
-	//!     
+	//!
 	//! There are three categories of analysis parameters:
 	//! - the resolution, and params that are usually related to (or
 	//! identical to) the resolution (frequency floor and drift)
@@ -200,31 +286,31 @@ public:
 	//! identical to) the window width (hop and crop times)
 	//! - independent parameters (bw region width and amp floor)
 	//
-	void configure( const Envelope & resolutionEnv, double windowWidthHz );    
-    
+	void configure( const Envelope & resolutionEnv, double windowWidthHz );
+
 //  -- analysis --
 
-    //! Analyze a vector of (mono) samples at the given sample rate         
-    //! (in Hz) and store the extracted Partials in the Analyzer's 
-    //! PartialList (std::list of Partials).    
-    //! 
+    //! Analyze a vector of (mono) samples at the given sample rate
+    //! (in Hz) and store the extracted Partials in the Analyzer's
+    //! PartialList (std::list of Partials).
+    //!
     //! \param  vec is a vector of floating point samples
-    //! \param  srate is the sample rate of the samples in the vector 
+    //! \param  srate is the sample rate of the samples in the vector
     void analyze( const std::vector<double> & vec, double srate );
-    
-    //! Analyze a range of (mono) samples at the given sample rate      
+
+    //! Analyze a range of (mono) samples at the given sample rate
     //! (in Hz) and store the extracted Partials in the Analyzer's
-    //! PartialList (std::list of Partials).    
-    //! 
+    //! PartialList (std::list of Partials).
+    //!
     //! \param  bufBegin is a pointer to a buffer of floating point samples
-    //! \param  bufEnd is (one-past) the end of a buffer of floating point 
+    //! \param  bufEnd is (one-past) the end of a buffer of floating point
     //!         samples
     //! \param  srate is the sample rate of the samples in the buffer
     void analyze( const double * bufBegin, const double * bufEnd, double srate );
-    
+
 //  -- tracking analysis --
 
-    //! Analyze a vector of (mono) samples at the given sample rate         
+    //! Analyze a vector of (mono) samples at the given sample rate
     //! (in Hz) and store the extracted Partials in the Analyzer's
     //! PartialList (std::list of Partials). Use the specified envelope
     //! as a frequency reference for Partial tracking.
@@ -233,29 +319,27 @@ public:
     //! \param  srate is the sample rate of the samples in the vector
     //! \param  reference is an Envelope having the approximate
     //!         frequency contour expected of the resulting Partials.
-    void analyze( const std::vector<double> & vec, double srate, 
+    void analyze( const std::vector<double> & vec, double srate,
                   const Envelope & reference );
-    
-    //! Analyze a range of (mono) samples at the given sample rate      
+
+    //! Analyze a range of (mono) samples at the given sample rate
     //! (in Hz) and store the extracted Partials in the Analyzer's
     //! PartialList (std::list of Partials). Use the specified envelope
     //! as a frequency reference for Partial tracking.
-    //! 
+    //!
     //! \param  bufBegin is a pointer to a buffer of floating point samples
-    //! \param  bufEnd is (one-past) the end of a buffer of floating point 
+    //! \param  bufEnd is (one-past) the end of a buffer of floating point
     //!         samples
     //! \param  srate is the sample rate of the samples in the buffer
     //! \param  reference is an Envelope having the approximate
     //!         frequency contour expected of the resulting Partials.
-    void analyze( const double * bufBegin, const double * bufEnd, double srate,
-                  const Envelope & reference );
+    virtual void analyze( const double * bufBegin, const double * bufEnd, double srate,
+                          const Envelope & reference );
 
-    void analyze_peaks( const std::vector<double> & vec, double srate );
-    
 //  -- parameter access --
 
-    //! Return the amplitude floor (lowest detected spectral amplitude),            
-    //! in (negative) dB, for this Analyzer.                
+    //! Return the amplitude floor (lowest detected spectral amplitude),
+    //! in (negative) dB, for this Analyzer.
     double ampFloor( void ) const;
 
     //! Return the crop time (maximum temporal displacement of a time-
@@ -264,40 +348,40 @@ public:
     //! for this Analyzer.
     double cropTime( void ) const;
 
-    //! Return the maximum allowable frequency difference between                   
-    //! consecutive Breakpoints in a Partial envelope for this Analyzer.                
+    //! Return the maximum allowable frequency difference between
+    //! consecutive Breakpoints in a Partial envelope for this Analyzer.
     double freqDrift( void ) const;
 
-    //! Return the frequency floor (minimum instantaneous Partial               
-    //! frequency), in Hz, for this Analyzer.               
+    //! Return the frequency floor (minimum instantaneous Partial
+    //! frequency), in Hz, for this Analyzer.
     double freqFloor( void ) const;
-	
-	//! Return the frequency resolution (minimum instantaneous frequency        
+
+	//! Return the frequency resolution (minimum instantaneous frequency
 	//! difference between Partials) for this Analyzer at the specified
 	//! time in seconds. If no time is specified, then the initial resolution
 	//!	(at 0 seconds) is returned.
-	//! 
-	//! \param time is the time in seconds at which to evaluate the 
+	//!
+	//! \param time is the time in seconds at which to evaluate the
 	//!		   frequency resolution
-	double freqResolution( double time = 0.0 ) const; 
+	double freqResolution( double time = 0.0 ) const;
 
-    //! Return the hop time (which corresponds approximately to the 
-    //! average density of Partial envelope Breakpoint data) for this 
+    //! Return the hop time (which corresponds approximately to the
+    //! average density of Partial envelope Breakpoint data) for this
     //! Analyzer.
     double hopTime( void ) const;
 
     //! Return the sidelobe attenutation level for the Kaiser analysis window in
-    //! positive dB. Larger numbers (e.g. 90) give very good sidelobe 
-    //! rejection but cause the window to be longer in time. Smaller numbers 
+    //! positive dB. Larger numbers (e.g. 90) give very good sidelobe
+    //! rejection but cause the window to be longer in time. Smaller numbers
     //! (like 60) raise the level of the sidelobes, increasing the likelihood
     //! of frequency-domain interference, but allow the window to be shorter
     //! in time.
     double sidelobeLevel( void ) const;
 
-    //! Return the frequency-domain main lobe width (measured between 
-    //! zero-crossings) of the analysis window used by this Analyzer.               
+    //! Return the frequency-domain main lobe width (measured between
+    //! zero-crossings) of the analysis window used by this Analyzer.
     double windowWidth( void ) const;
-     
+
     //! Return true if the phases and frequencies of the constructed
     //! partials should be modified to be consistent at the end of the
     //! analysis, and false otherwise. (Default is true.)
@@ -306,66 +390,66 @@ public:
 
 //  -- parameter mutation --
 
-    //! Set the amplitude floor (lowest detected spectral amplitude), in            
-    //! (negative) dB, for this Analyzer. 
-    //! 
-    //! \param x is the new value of this parameter.            
+    //! Set the amplitude floor (lowest detected spectral amplitude), in
+    //! (negative) dB, for this Analyzer.
+    //!
+    //! \param x is the new value of this parameter.
     void setAmpFloor( double x );
 
     //! Set the crop time (maximum temporal displacement of a time-
     //! frequency data point from the time-domain center of the analysis
     //! window, beyond which data points are considered "unreliable")
     //! for this Analyzer.
-    //! 
-    //! \param x is the new value of this parameter.            
+    //!
+    //! \param x is the new value of this parameter.
     void setCropTime( double x );
 
-    //! Set the maximum allowable frequency difference between                  
-    //! consecutive Breakpoints in a Partial envelope for this Analyzer.                
-    //! 
-    //! \param x is the new value of this parameter.            
+    //! Set the maximum allowable frequency difference between
+    //! consecutive Breakpoints in a Partial envelope for this Analyzer.
+    //!
+    //! \param x is the new value of this parameter.
     void setFreqDrift( double x );
 
-    //! Set the frequency floor (minimum instantaneous Partial                  
+    //! Set the frequency floor (minimum instantaneous Partial
     //! frequency), in Hz, for this Analyzer.
-    //! 
-    //! \param x is the new value of this parameter.            
+    //!
+    //! \param x is the new value of this parameter.
     void setFreqFloor( double x );
 
-    //! Set the frequency resolution (minimum instantaneous frequency       
-    //! difference between Partials) for this Analyzer. (Does not cause     
-    //! other parameters to be recomputed.)                                     
-    //! 
-    //! \param x is the new value of this parameter.            
+    //! Set the frequency resolution (minimum instantaneous frequency
+    //! difference between Partials) for this Analyzer. (Does not cause
+    //! other parameters to be recomputed.)
+    //!
+    //! \param x is the new value of this parameter.
     void setFreqResolution( double x );
-    
-	//! Set the time-varying frequency resolution (minimum instantaneous frequency       
-	//! difference between Partials) for this Analyzer. (Does not cause     
-	//! other parameters to be recomputed.)                                     
-	//! 
-	//! \param e is the envelope to copy for this parameter.                                        
-	void setFreqResolution( const Envelope & e );    
+
+	//! Set the time-varying frequency resolution (minimum instantaneous frequency
+	//! difference between Partials) for this Analyzer. (Does not cause
+	//! other parameters to be recomputed.)
+	//!
+	//! \param e is the envelope to copy for this parameter.
+	void setFreqResolution( const Envelope & e );
 
     //! Set the hop time (which corresponds approximately to the average
     //! density of Partial envelope Breakpoint data) for this Analyzer.
-    //! 
-    //! \param x is the new value of this parameter.            
+    //!
+    //! \param x is the new value of this parameter.
     void setHopTime( double x );
 
     //! Set the sidelobe attenutation level for the Kaiser analysis window in
-    //! positive dB. More negative numbers (e.g. -90) give very good sidelobe 
-    //! rejection but cause the window to be longer in time. Less negative 
+    //! positive dB. More negative numbers (e.g. -90) give very good sidelobe
+    //! rejection but cause the window to be longer in time. Less negative
     //! numbers raise the level of the sidelobes, increasing the likelihood
     //! of frequency-domain interference, but allow the window to be shorter
     //! in time.
-    //! 
-    //! \param x is the new value of this parameter.            
+    //!
+    //! \param x is the new value of this parameter.
     void setSidelobeLevel( double x );
 
-    //! Set the frequency-domain main lobe width (measured between 
-    //! zero-crossings) of the analysis window used by this Analyzer.   
-    //! 
-    //! \param x is the new value of this parameter.            
+    //! Set the frequency-domain main lobe width (measured between
+    //! zero-crossings) of the analysis window used by this Analyzer.
+    //!
+    //! \param x is the new value of this parameter.
     void setWindowWidth( double x );
 
     //! Indicate whether the phases and frequencies of the constructed
@@ -375,54 +459,54 @@ public:
     //! \param  TF is a flag indicating whether or not to construct
     //!         phase-corrected Partials
     void setPhaseCorrect( bool TF = true );
-    
-    
+
+
 //  -- bandwidth envelope specification --
 
     enum { Default_ResidueBandwidth_RegionWidth = 2000,
            Default_ConvergenceBandwidth_TolerancePct = 10 };
-           
+
     //! Construct Partial bandwidth envelopes during analysis
     //! by associating residual energy in the spectrum (after
     //! peak extraction) with the selected spectral peaks that
-    //! are used to construct Partials. 
+    //! are used to construct Partials.
     //!
     //!	This is the default bandwidth-enhancement strategy.
-    //! 
-    //! \param regionWidth is the width (in Hz) of the bandwidth 
+    //!
+    //! \param regionWidth is the width (in Hz) of the bandwidth
     //! association regions used by this process, must be positive.
     //! If unspecified, a default value is used.
     void storeResidueBandwidth( double regionWidth = Default_ResidueBandwidth_RegionWidth );
-    
+
 	//!	Construct Partial bandwidth envelopes during analysis
-	//!	by storing the mixed derivative of short-time phase, 
+	//!	by storing the mixed derivative of short-time phase,
 	//!	scaled and shifted so that a value of 0 corresponds
 	//!	to a pure sinusoid, and a value of 1 corresponds to a
 	//! bandwidth-enhanced sinusoid with maximal energy spread
 	//! (minimum sinusoidal convergence).
 	//!
-	//!	\param tolerance is the amount of range over which the 
-	//!	mixed derivative indicator should be allowed to drift away 
+	//!	\param tolerance is the amount of range over which the
+	//!	mixed derivative indicator should be allowed to drift away
 	//!	from a pure sinusoid before saturating. This range is mapped
-	//!	to bandwidth values on the range [0,1]. Must be positive and 
+	//!	to bandwidth values on the range [0,1]. Must be positive and
 	//!	not greater than 1. If unspecified, a default value is used.
-    void storeConvergenceBandwidth( double tolerancePct = 
+    void storeConvergenceBandwidth( double tolerancePct =
     		0.01 * (double)Default_ConvergenceBandwidth_TolerancePct );
-    
-    //! Disable bandwidth envelope construction. Bandwidth 
+
+    //! Disable bandwidth envelope construction. Bandwidth
     //! will be zero for all Breakpoints in all Partials.
     void storeNoBandwidth( void );
-    
+
     //! Return true if this Analyzer is configured to compute
     //! bandwidth envelopes using the spectral residue after
     //! peaks have been identified, and false otherwise.
     bool bandwidthIsResidue( void ) const;
-    
+
     //! Return true if this Analyzer is configured to compute
     //! bandwidth envelopes using the mixed derivative convergence
     //! indicator, and false otherwise.
     bool bandwidthIsConvergence( void ) const;
-    
+
     //! Return the width (in Hz) of the Bandwidth Association regions
     //! used by this Analyzer, only if the spectral residue method is
     //! used to compute bandwidth envelopes. Return zero if the mixed
@@ -436,45 +520,45 @@ public:
     double bwConvergenceTolerance( void ) const;
 
     //! Return true if bandwidth envelopes are to be constructed
-    //!	by any means, that is, if either bandwidthIsResidue() or 
+    //!	by any means, that is, if either bandwidthIsResidue() or
     //!	bandwidthIsConvergence() are true. Otherwise, return
     //! false.
     bool associateBandwidth( void ) const
         { return bandwidthIsResidue() || bandwidthIsConvergence(); }
 
-    //! Deprecated, use storeResidueBandwidth or storeNoBandwidth instead.            
-    void setBwRegionWidth( double x ) 
-        { 
+    //! Deprecated, use storeResidueBandwidth or storeNoBandwidth instead.
+    void setBwRegionWidth( double x )
+        {
             if ( x != 0 )
             {
-                storeResidueBandwidth( x ); 
+                storeResidueBandwidth( x );
             }
             else
             {
                 storeNoBandwidth();
             }
         }
-           
+
 
 //  -- PartialList access --
 
-    //! Return a mutable reference to this Analyzer's list of 
-    //! analyzed Partials. 
+    //! Return a mutable reference to this Analyzer's list of
+    //! analyzed Partials.
     PartialList & partials( void );
 
-    //! Return an immutable (const) reference to this Analyzer's 
-    //! list of analyzed Partials. 
+    //! Return an immutable (const) reference to this Analyzer's
+    //! list of analyzed Partials.
     const PartialList & partials( void ) const;
 
 //  -- envelope access --
 
-    enum { Default_FundamentalEnv_ThreshDb = -60, 
+    enum { Default_FundamentalEnv_ThreshDb = -60,
            Default_FundamentalEnv_ThreshHz = 8000 };
 
-    //! Specify parameters for constructing a fundamental frequency 
-    //! envelope for the analyzed sound during analysis. The fundamental 
-    //! frequency estimate can be accessed by fundamentalEnv() after the 
-    //! analysis is complete. 
+    //! Specify parameters for constructing a fundamental frequency
+    //! envelope for the analyzed sound during analysis. The fundamental
+    //! frequency estimate can be accessed by fundamentalEnv() after the
+    //! analysis is complete.
     //!
     //! By default, a fundamental envelope is estimated during analysis
     //! between the frequency resolution and 1.5 times the resolution.
@@ -488,9 +572,9 @@ public:
     //! \param  threshHz is the upper bound on the frequency of a spectral
     //!         peak that will constribute to the fundamental frequency estimate.
     //!         Default is 8 kHz.
-    void buildFundamentalEnv( double fmin, double fmax, 
-                              double threshDb = Default_FundamentalEnv_ThreshDb, 
-                              double threshHz = Default_FundamentalEnv_ThreshHz );                                     
+    void buildFundamentalEnv( double fmin, double fmax,
+                              double threshDb = Default_FundamentalEnv_ThreshDb,
+                              double threshHz = Default_FundamentalEnv_ThreshHz );
 
 
     //! Return the fundamental frequency estimate envelope constructed
@@ -503,10 +587,10 @@ public:
     //! Return the overall amplitude estimate envelope constructed
     //! during the most recent analysis performed by this Analyzer.
     const LinearEnvelope & ampEnv( void ) const;
-    
-    
+
+
 //  -- legacy support --
-    
+
     //  Fundamental and amplitude envelopes are always constructed during
     //  analysis, these members do nothing, and are retained for backwards
     //  compatibility.
@@ -515,51 +599,51 @@ public:
 
 //  -- private member variables --
 
-private:
+protected:
 
-    std::auto_ptr< Envelope > m_freqResolutionEnv;    
+    std::auto_ptr< Envelope > m_freqResolutionEnv;
     							//!  in Hz, minimum instantaneous frequency distance;
                                 //!  this is the core parameter, others are, by default,
                                 //!  computed from this one
-    
+
     double m_ampFloor;          //!  dB, relative to full amplitude sine wave, absolute
                                 //!  amplitude threshold (negative)
-    
+
     double m_windowWidth;       //!  in Hz, width of main lobe; this might be more
                                 //!  conveniently presented as window length, but
                                 //!  the main lobe width more explicitly highlights
                                 //!  the critical interaction with resolution
-    
-    // std::auto_ptr< Envelope > m_freqFloorEnv; 
+
+    // std::auto_ptr< Envelope > m_freqFloorEnv;
     double m_freqFloor;         //!  lowest frequency (Hz) component extracted
                                 //!  in spectral analysis
-    
-    double m_freqDrift;         //!  the maximum frequency (Hz) difference between two 
+
+    double m_freqDrift;         //!  the maximum frequency (Hz) difference between two
                                 //!  consecutive Breakpoints that will be linked to
                                 //!  form a Partial
-    
+
     double m_hopTime;           //!  in seconds, time between analysis windows in
                                 //!  successive spectral analyses
-    
+
     double m_cropTime;          //!  in seconds, maximum time correction for a spectral
                                 //!  component to be considered reliable, and to be eligible
                                 //!  for extraction and for Breakpoint formation
-    
-    double m_bwAssocParam;      //!  formerly, width in Hz of overlapping bandwidth 
+
+    double m_bwAssocParam;      //!  formerly, width in Hz of overlapping bandwidth
                                 //!  association regions, or zero if bandwidth association
                                 //!  is disabled, now a catch-all bandwidth association
                                 //!  parameter that, if negative, indicates the tolerance (%)
                                 //!  level used to construct bandwidth envelopes from the
                                 //!  mixed phase derivative indicator
-                                                        
-    double m_sidelobeLevel;     //!  sidelobe attenutation level for the Kaiser analysis 
+
+    double m_sidelobeLevel;     //!  sidelobe attenutation level for the Kaiser analysis
                                 //!  window, in positive dB
-                                
+
     bool m_phaseCorrect;        //!  flag indicating that phases/frequencies should be
                                 //!  made consistent at the end of the analysis
-                            
+
     PartialList m_partials;     //!  collect Partials here
-        
+
     //! builder object for constructing a fundamental frequency
     //! estimate during analysis
     std::auto_ptr< LinearEnvelopeBuilder > m_f0Builder;
@@ -573,22 +657,22 @@ private:
 /*
 
 	//	These members make up the sequence of operations in an
-	//	analysis. If analysis were ever to be made into a 
+	//	analysis. If analysis were ever to be made into a
 	//	template method, these would be the operations that
 	//	derived classes could override. Or each of these could
 	//	be represented by a strategy class.
 
 	//!	Compute the spectrum of the next sequence of samples.
 	void computeSpectrum( void );
-	
+
 	//!	Identify and select the spectral components that will be
-	//!	used to form Partials. 
+	//!	used to form Partials.
 	void selectPeaks( void );
-	
+
 	//!	Compute the bandwidth coefficients for the Breakpoints
-	//!	that are going to be used to form Partials. 
+	//!	that are going to be used to form Partials.
 	void associateBandwidth( void );
-	
+
 	//!	Construct Partials from extracted spectral components.
 	//!	Partials are built up frame by frame by appending
 	//!	Breakpoints to Partials under construction, and giving
@@ -598,21 +682,21 @@ private:
     //  Reject peaks that are too close in frequency to a louder peak that is
     //  being retained, and peaks that are too quiet. Peaks that are retained,
     //  but are quiet enough to be in the specified fadeRange should be faded.
-    //  
+    //
     //  Rejected peaks are placed at the end of the peak collection.
     //  Return the first position in the collection containing a rejected peak,
     //  or the end of the collection if no peaks are rejected.
     Peaks::iterator thinPeaks( Peaks & peaks, double frameTime  );
-                
-    //  Fix the bandwidth value stored in the specified Peaks. 
+
+    //  Fix the bandwidth value stored in the specified Peaks.
     //  This function is invoked if the spectral residue method is
     //  not used to compute bandwidth (that method overwrites the
-    //  bandwidth already). If the convergence method is used to 
+    //  bandwidth already). If the convergence method is used to
     //  compute bandwidth, the appropriate scaling is applied
     //  to the stored mixed phase derivative. Otherwise, the
     //  Peak bandwidth is set to zero.
     void fixBandwidth( Peaks & peaks );
-                    
+
 };  //  end of class Analyzer
 
 }   //  end of namespace Loris
diff --git a/src/loris/FourierTransform.C b/src/loris/FourierTransform.C
index 2a29238..b67e6de 100644
--- a/src/loris/FourierTransform.C
+++ b/src/loris/FourierTransform.C
@@ -328,7 +328,7 @@ public:
 #define SORRY_NO_FFTW  1
 
 //  function prototype, definition in fftsg.c
-extern "C" void cdft(int, int, double *, int *, double *);
+/* extern "C" void cdft(int, int, double *, int *, double *); */
 
 //  function prototype, definition below
 static void slowDFT( double * in, double * out, int N );
@@ -444,7 +444,8 @@ public:
     {        
         if ( mIsPO2 )
         {
-            cdft( 2*N, -1, mTxInOut, mWorkspace, mTwiddle );
+            /* cdft( 2*N, -1, mTxInOut, mWorkspace, mTwiddle ); */
+            slowDFT( mTxInOut, mTwiddle, N );
         }
         else
         {
diff --git a/src/loris/fftsg.c b/src/loris/fftsg.c
deleted file mode 100644
index dc3e237..0000000
--- a/src/loris/fftsg.c
+++ /dev/null
@@ -1,3331 +0,0 @@
-/*  
-    double-precision floating point transform functions
-    adapted from the tranform library of Ooura.
-    http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html
-
-    Identical to fftsg.c, wrapped in extern "C"
-    to prevent name-mangling by a C++ compiler.
-
-    Kelly Fitz 21 July 2006
-    kfitz@cerlsoundgroup.org
-*/
-/*
-Fast Fourier/Cosine/Sine Transform
-    dimension   :one
-    data length :power of 2
-    decimation  :frequency
-    radix       :split-radix
-    data        :inplace
-    table       :use
-functions
-    cdft: Complex Discrete Fourier Transform
-    rdft: Real Discrete Fourier Transform
-    ddct: Discrete Cosine Transform
-    ddst: Discrete Sine Transform
-    dfct: Cosine Transform of RDFT (Real Symmetric DFT)
-    dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
-function prototypes
-    void cdft(int, int, double *, int *, double *);
-    void rdft(int, int, double *, int *, double *);
-    void ddct(int, int, double *, int *, double *);
-    void ddst(int, int, double *, int *, double *);
-    void dfct(int, double *, double *, int *, double *);
-    void dfst(int, double *, double *, int *, double *);
-macro definitions
-    USE_CDFT_PTHREADS : default=not defined
-        CDFT_THREADS_BEGIN_N  : must be >= 512, default=8192
-        CDFT_4THREADS_BEGIN_N : must be >= 512, default=65536
-    USE_CDFT_WINTHREADS : default=not defined
-        CDFT_THREADS_BEGIN_N  : must be >= 512, default=32768
-        CDFT_4THREADS_BEGIN_N : must be >= 512, default=524288
-
-
--------- Complex DFT (Discrete Fourier Transform) --------
-    [definition]
-        <case1>
-            X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
-        <case2>
-            X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
-        (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
-    [usage]
-        <case1>
-            ip[0] = 0; // first time only
-            cdft(2*n, 1, a, ip, w);
-        <case2>
-            ip[0] = 0; // first time only
-            cdft(2*n, -1, a, ip, w);
-    [parameters]
-        2*n            :data length (int)
-                        n >= 1, n = power of 2
-        a[0...2*n-1]   :input/output data (double *)
-                        input data
-                            a[2*j] = Re(x[j]), 
-                            a[2*j+1] = Im(x[j]), 0<=j<n
-                        output data
-                            a[2*k] = Re(X[k]), 
-                            a[2*k+1] = Im(X[k]), 0<=k<n
-        ip[0...*]      :work area for bit reversal (int *)
-                        length of ip >= 2+sqrt(n)
-                        strictly, 
-                        length of ip >= 
-                            2+(1<<(int)(log(n+0.5)/log(2))/2).
-                        ip[0],ip[1] are pointers of the cos/sin table.
-        w[0...n/2-1]   :cos/sin table (double *)
-                        w[],ip[] are initialized if ip[0] == 0.
-    [remark]
-        Inverse of 
-            cdft(2*n, -1, a, ip, w);
-        is 
-            cdft(2*n, 1, a, ip, w);
-            for (j = 0; j <= 2 * n - 1; j++) {
-                a[j] *= 1.0 / n;
-            }
-        .
-
-
--------- Real DFT / Inverse of Real DFT --------
-    [definition]
-        <case1> RDFT
-            R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
-            I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
-        <case2> IRDFT (excluding scale)
-            a[k] = (R[0] + R[n/2]*cos(pi*k))/2 + 
-                   sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) + 
-                   sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
-    [usage]
-        <case1>
-            ip[0] = 0; // first time only
-            rdft(n, 1, a, ip, w);
-        <case2>
-            ip[0] = 0; // first time only
-            rdft(n, -1, a, ip, w);
-    [parameters]
-        n              :data length (int)
-                        n >= 2, n = power of 2
-        a[0...n-1]     :input/output data (double *)
-                        <case1>
-                            output data
-                                a[2*k] = R[k], 0<=k<n/2
-                                a[2*k+1] = I[k], 0<k<n/2
-                                a[1] = R[n/2]
-                        <case2>
-                            input data
-                                a[2*j] = R[j], 0<=j<n/2
-                                a[2*j+1] = I[j], 0<j<n/2
-                                a[1] = R[n/2]
-        ip[0...*]      :work area for bit reversal (int *)
-                        length of ip >= 2+sqrt(n/2)
-                        strictly, 
-                        length of ip >= 
-                            2+(1<<(int)(log(n/2+0.5)/log(2))/2).
-                        ip[0],ip[1] are pointers of the cos/sin table.
-        w[0...n/2-1]   :cos/sin table (double *)
-                        w[],ip[] are initialized if ip[0] == 0.
-    [remark]
-        Inverse of 
-            rdft(n, 1, a, ip, w);
-        is 
-            rdft(n, -1, a, ip, w);
-            for (j = 0; j <= n - 1; j++) {
-                a[j] *= 2.0 / n;
-            }
-        .
-
-
--------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
-    [definition]
-        <case1> IDCT (excluding scale)
-            C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
-        <case2> DCT
-            C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
-    [usage]
-        <case1>
-            ip[0] = 0; // first time only
-            ddct(n, 1, a, ip, w);
-        <case2>
-            ip[0] = 0; // first time only
-            ddct(n, -1, a, ip, w);
-    [parameters]
-        n              :data length (int)
-                        n >= 2, n = power of 2
-        a[0...n-1]     :input/output data (double *)
-                        output data
-                            a[k] = C[k], 0<=k<n
-        ip[0...*]      :work area for bit reversal (int *)
-                        length of ip >= 2+sqrt(n/2)
-                        strictly, 
-                        length of ip >= 
-                            2+(1<<(int)(log(n/2+0.5)/log(2))/2).
-                        ip[0],ip[1] are pointers of the cos/sin table.
-        w[0...n*5/4-1] :cos/sin table (double *)
-                        w[],ip[] are initialized if ip[0] == 0.
-    [remark]
-        Inverse of 
-            ddct(n, -1, a, ip, w);
-        is 
-            a[0] *= 0.5;
-            ddct(n, 1, a, ip, w);
-            for (j = 0; j <= n - 1; j++) {
-                a[j] *= 2.0 / n;
-            }
-        .
-
-
--------- DST (Discrete Sine Transform) / Inverse of DST --------
-    [definition]
-        <case1> IDST (excluding scale)
-            S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
-        <case2> DST
-            S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
-    [usage]
-        <case1>
-            ip[0] = 0; // first time only
-            ddst(n, 1, a, ip, w);
-        <case2>
-            ip[0] = 0; // first time only
-            ddst(n, -1, a, ip, w);
-    [parameters]
-        n              :data length (int)
-                        n >= 2, n = power of 2
-        a[0...n-1]     :input/output data (double *)
-                        <case1>
-                            input data
-                                a[j] = A[j], 0<j<n
-                                a[0] = A[n]
-                            output data
-                                a[k] = S[k], 0<=k<n
-                        <case2>
-                            output data
-                                a[k] = S[k], 0<k<n
-                                a[0] = S[n]
-        ip[0...*]      :work area for bit reversal (int *)
-                        length of ip >= 2+sqrt(n/2)
-                        strictly, 
-                        length of ip >= 
-                            2+(1<<(int)(log(n/2+0.5)/log(2))/2).
-                        ip[0],ip[1] are pointers of the cos/sin table.
-        w[0...n*5/4-1] :cos/sin table (double *)
-                        w[],ip[] are initialized if ip[0] == 0.
-    [remark]
-        Inverse of 
-            ddst(n, -1, a, ip, w);
-        is 
-            a[0] *= 0.5;
-            ddst(n, 1, a, ip, w);
-            for (j = 0; j <= n - 1; j++) {
-                a[j] *= 2.0 / n;
-            }
-        .
-
-
--------- Cosine Transform of RDFT (Real Symmetric DFT) --------
-    [definition]
-        C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
-    [usage]
-        ip[0] = 0; // first time only
-        dfct(n, a, t, ip, w);
-    [parameters]
-        n              :data length - 1 (int)
-                        n >= 2, n = power of 2
-        a[0...n]       :input/output data (double *)
-                        output data
-                            a[k] = C[k], 0<=k<=n
-        t[0...n/2]     :work area (double *)
-        ip[0...*]      :work area for bit reversal (int *)
-                        length of ip >= 2+sqrt(n/4)
-                        strictly, 
-                        length of ip >= 
-                            2+(1<<(int)(log(n/4+0.5)/log(2))/2).
-                        ip[0],ip[1] are pointers of the cos/sin table.
-        w[0...n*5/8-1] :cos/sin table (double *)
-                        w[],ip[] are initialized if ip[0] == 0.
-    [remark]
-        Inverse of 
-            a[0] *= 0.5;
-            a[n] *= 0.5;
-            dfct(n, a, t, ip, w);
-        is 
-            a[0] *= 0.5;
-            a[n] *= 0.5;
-            dfct(n, a, t, ip, w);
-            for (j = 0; j <= n; j++) {
-                a[j] *= 2.0 / n;
-            }
-        .
-
-
--------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
-    [definition]
-        S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
-    [usage]
-        ip[0] = 0; // first time only
-        dfst(n, a, t, ip, w);
-    [parameters]
-        n              :data length + 1 (int)
-                        n >= 2, n = power of 2
-        a[0...n-1]     :input/output data (double *)
-                        output data
-                            a[k] = S[k], 0<k<n
-                        (a[0] is used for work area)
-        t[0...n/2-1]   :work area (double *)
-        ip[0...*]      :work area for bit reversal (int *)
-                        length of ip >= 2+sqrt(n/4)
-                        strictly, 
-                        length of ip >= 
-                            2+(1<<(int)(log(n/4+0.5)/log(2))/2).
-                        ip[0],ip[1] are pointers of the cos/sin table.
-        w[0...n*5/8-1] :cos/sin table (double *)
-                        w[],ip[] are initialized if ip[0] == 0.
-    [remark]
-        Inverse of 
-            dfst(n, a, t, ip, w);
-        is 
-            dfst(n, a, t, ip, w);
-            for (j = 1; j <= n - 1; j++) {
-                a[j] *= 2.0 / n;
-            }
-        .
-
-
-Appendix :
-    The cos/sin table is recalculated when the larger table required.
-    w[] and ip[] are compatible with all routines.
-*/
-
-#if defined(__cplusplus)
-extern "C" {
-#endif
-
-void cdft(int n, int isgn, double *a, int *ip, double *w)
-{
-    void makewt(int nw, int *ip, double *w);
-    void cftfsub(int n, double *a, int *ip, int nw, double *w);
-    void cftbsub(int n, double *a, int *ip, int nw, double *w);
-    int nw;
-    
-    nw = ip[0];
-    if (n > (nw << 2)) {
-        nw = n >> 2;
-        makewt(nw, ip, w);
-    }
-    if (isgn >= 0) {
-        cftfsub(n, a, ip, nw, w);
-    } else {
-        cftbsub(n, a, ip, nw, w);
-    }
-}
-
-
-void rdft(int n, int isgn, double *a, int *ip, double *w)
-{
-    void makewt(int nw, int *ip, double *w);
-    void makect(int nc, int *ip, double *c);
-    void cftfsub(int n, double *a, int *ip, int nw, double *w);
-    void cftbsub(int n, double *a, int *ip, int nw, double *w);
-    void rftfsub(int n, double *a, int nc, double *c);
-    void rftbsub(int n, double *a, int nc, double *c);
-    int nw, nc;
-    double xi;
-    
-    nw = ip[0];
-    if (n > (nw << 2)) {
-        nw = n >> 2;
-        makewt(nw, ip, w);
-    }
-    nc = ip[1];
-    if (n > (nc << 2)) {
-        nc = n >> 2;
-        makect(nc, ip, w + nw);
-    }
-    if (isgn >= 0) {
-        if (n > 4) {
-            cftfsub(n, a, ip, nw, w);
-            rftfsub(n, a, nc, w + nw);
-        } else if (n == 4) {
-            cftfsub(n, a, ip, nw, w);
-        }
-        xi = a[0] - a[1];
-        a[0] += a[1];
-        a[1] = xi;
-    } else {
-        a[1] = 0.5 * (a[0] - a[1]);
-        a[0] -= a[1];
-        if (n > 4) {
-            rftbsub(n, a, nc, w + nw);
-            cftbsub(n, a, ip, nw, w);
-        } else if (n == 4) {
-            cftbsub(n, a, ip, nw, w);
-        }
-    }
-}
-
-
-void ddct(int n, int isgn, double *a, int *ip, double *w)
-{
-    void makewt(int nw, int *ip, double *w);
-    void makect(int nc, int *ip, double *c);
-    void cftfsub(int n, double *a, int *ip, int nw, double *w);
-    void cftbsub(int n, double *a, int *ip, int nw, double *w);
-    void rftfsub(int n, double *a, int nc, double *c);
-    void rftbsub(int n, double *a, int nc, double *c);
-    void dctsub(int n, double *a, int nc, double *c);
-    int j, nw, nc;
-    double xr;
-    
-    nw = ip[0];
-    if (n > (nw << 2)) {
-        nw = n >> 2;
-        makewt(nw, ip, w);
-    }
-    nc = ip[1];
-    if (n > nc) {
-        nc = n;
-        makect(nc, ip, w + nw);
-    }
-    if (isgn < 0) {
-        xr = a[n - 1];
-        for (j = n - 2; j >= 2; j -= 2) {
-            a[j + 1] = a[j] - a[j - 1];
-            a[j] += a[j - 1];
-        }
-        a[1] = a[0] - xr;
-        a[0] += xr;
-        if (n > 4) {
-            rftbsub(n, a, nc, w + nw);
-            cftbsub(n, a, ip, nw, w);
-        } else if (n == 4) {
-            cftbsub(n, a, ip, nw, w);
-        }
-    }
-    dctsub(n, a, nc, w + nw);
-    if (isgn >= 0) {
-        if (n > 4) {
-            cftfsub(n, a, ip, nw, w);
-            rftfsub(n, a, nc, w + nw);
-        } else if (n == 4) {
-            cftfsub(n, a, ip, nw, w);
-        }
-        xr = a[0] - a[1];
-        a[0] += a[1];
-        for (j = 2; j < n; j += 2) {
-            a[j - 1] = a[j] - a[j + 1];
-            a[j] += a[j + 1];
-        }
-        a[n - 1] = xr;
-    }
-}
-
-
-void ddst(int n, int isgn, double *a, int *ip, double *w)
-{
-    void makewt(int nw, int *ip, double *w);
-    void makect(int nc, int *ip, double *c);
-    void cftfsub(int n, double *a, int *ip, int nw, double *w);
-    void cftbsub(int n, double *a, int *ip, int nw, double *w);
-    void rftfsub(int n, double *a, int nc, double *c);
-    void rftbsub(int n, double *a, int nc, double *c);
-    void dstsub(int n, double *a, int nc, double *c);
-    int j, nw, nc;
-    double xr;
-    
-    nw = ip[0];
-    if (n > (nw << 2)) {
-        nw = n >> 2;
-        makewt(nw, ip, w);
-    }
-    nc = ip[1];
-    if (n > nc) {
-        nc = n;
-        makect(nc, ip, w + nw);
-    }
-    if (isgn < 0) {
-        xr = a[n - 1];
-        for (j = n - 2; j >= 2; j -= 2) {
-            a[j + 1] = -a[j] - a[j - 1];
-            a[j] -= a[j - 1];
-        }
-        a[1] = a[0] + xr;
-        a[0] -= xr;
-        if (n > 4) {
-            rftbsub(n, a, nc, w + nw);
-            cftbsub(n, a, ip, nw, w);
-        } else if (n == 4) {
-            cftbsub(n, a, ip, nw, w);
-        }
-    }
-    dstsub(n, a, nc, w + nw);
-    if (isgn >= 0) {
-        if (n > 4) {
-            cftfsub(n, a, ip, nw, w);
-            rftfsub(n, a, nc, w + nw);
-        } else if (n == 4) {
-            cftfsub(n, a, ip, nw, w);
-        }
-        xr = a[0] - a[1];
-        a[0] += a[1];
-        for (j = 2; j < n; j += 2) {
-            a[j - 1] = -a[j] - a[j + 1];
-            a[j] -= a[j + 1];
-        }
-        a[n - 1] = -xr;
-    }
-}
-
-
-void dfct(int n, double *a, double *t, int *ip, double *w)
-{
-    void makewt(int nw, int *ip, double *w);
-    void makect(int nc, int *ip, double *c);
-    void cftfsub(int n, double *a, int *ip, int nw, double *w);
-    void rftfsub(int n, double *a, int nc, double *c);
-    void dctsub(int n, double *a, int nc, double *c);
-    int j, k, l, m, mh, nw, nc;
-    double xr, xi, yr, yi;
-    
-    nw = ip[0];
-    if (n > (nw << 3)) {
-        nw = n >> 3;
-        makewt(nw, ip, w);
-    }
-    nc = ip[1];
-    if (n > (nc << 1)) {
-        nc = n >> 1;
-        makect(nc, ip, w + nw);
-    }
-    m = n >> 1;
-    yi = a[m];
-    xi = a[0] + a[n];
-    a[0] -= a[n];
-    t[0] = xi - yi;
-    t[m] = xi + yi;
-    if (n > 2) {
-        mh = m >> 1;
-        for (j = 1; j < mh; j++) {
-            k = m - j;
-            xr = a[j] - a[n - j];
-            xi = a[j] + a[n - j];
-            yr = a[k] - a[n - k];
-            yi = a[k] + a[n - k];
-            a[j] = xr;
-            a[k] = yr;
-            t[j] = xi - yi;
-            t[k] = xi + yi;
-        }
-        t[mh] = a[mh] + a[n - mh];
-        a[mh] -= a[n - mh];
-        dctsub(m, a, nc, w + nw);
-        if (m > 4) {
-            cftfsub(m, a, ip, nw, w);
-            rftfsub(m, a, nc, w + nw);
-        } else if (m == 4) {
-            cftfsub(m, a, ip, nw, w);
-        }
-        a[n - 1] = a[0] - a[1];
-        a[1] = a[0] + a[1];
-        for (j = m - 2; j >= 2; j -= 2) {
-            a[2 * j + 1] = a[j] + a[j + 1];
-            a[2 * j - 1] = a[j] - a[j + 1];
-        }
-        l = 2;
-        m = mh;
-        while (m >= 2) {
-            dctsub(m, t, nc, w + nw);
-            if (m > 4) {
-                cftfsub(m, t, ip, nw, w);
-                rftfsub(m, t, nc, w + nw);
-            } else if (m == 4) {
-                cftfsub(m, t, ip, nw, w);
-            }
-            a[n - l] = t[0] - t[1];
-            a[l] = t[0] + t[1];
-            k = 0;
-            for (j = 2; j < m; j += 2) {
-                k += l << 2;
-                a[k - l] = t[j] - t[j + 1];
-                a[k + l] = t[j] + t[j + 1];
-            }
-            l <<= 1;
-            mh = m >> 1;
-            for (j = 0; j < mh; j++) {
-                k = m - j;
-                t[j] = t[m + k] - t[m + j];
-                t[k] = t[m + k] + t[m + j];
-            }
-            t[mh] = t[m + mh];
-            m = mh;
-        }
-        a[l] = t[0];
-        a[n] = t[2] - t[1];
-        a[0] = t[2] + t[1];
-    } else {
-        a[1] = a[0];
-        a[2] = t[0];
-        a[0] = t[1];
-    }
-}
-
-
-void dfst(int n, double *a, double *t, int *ip, double *w)
-{
-    void makewt(int nw, int *ip, double *w);
-    void makect(int nc, int *ip, double *c);
-    void cftfsub(int n, double *a, int *ip, int nw, double *w);
-    void rftfsub(int n, double *a, int nc, double *c);
-    void dstsub(int n, double *a, int nc, double *c);
-    int j, k, l, m, mh, nw, nc;
-    double xr, xi, yr, yi;
-    
-    nw = ip[0];
-    if (n > (nw << 3)) {
-        nw = n >> 3;
-        makewt(nw, ip, w);
-    }
-    nc = ip[1];
-    if (n > (nc << 1)) {
-        nc = n >> 1;
-        makect(nc, ip, w + nw);
-    }
-    if (n > 2) {
-        m = n >> 1;
-        mh = m >> 1;
-        for (j = 1; j < mh; j++) {
-            k = m - j;
-            xr = a[j] + a[n - j];
-            xi = a[j] - a[n - j];
-            yr = a[k] + a[n - k];
-            yi = a[k] - a[n - k];
-            a[j] = xr;
-            a[k] = yr;
-            t[j] = xi + yi;
-            t[k] = xi - yi;
-        }
-        t[0] = a[mh] - a[n - mh];
-        a[mh] += a[n - mh];
-        a[0] = a[m];
-        dstsub(m, a, nc, w + nw);
-        if (m > 4) {
-            cftfsub(m, a, ip, nw, w);
-            rftfsub(m, a, nc, w + nw);
-        } else if (m == 4) {
-            cftfsub(m, a, ip, nw, w);
-        }
-        a[n - 1] = a[1] - a[0];
-        a[1] = a[0] + a[1];
-        for (j = m - 2; j >= 2; j -= 2) {
-            a[2 * j + 1] = a[j] - a[j + 1];
-            a[2 * j - 1] = -a[j] - a[j + 1];
-        }
-        l = 2;
-        m = mh;
-        while (m >= 2) {
-            dstsub(m, t, nc, w + nw);
-            if (m > 4) {
-                cftfsub(m, t, ip, nw, w);
-                rftfsub(m, t, nc, w + nw);
-            } else if (m == 4) {
-                cftfsub(m, t, ip, nw, w);
-            }
-            a[n - l] = t[1] - t[0];
-            a[l] = t[0] + t[1];
-            k = 0;
-            for (j = 2; j < m; j += 2) {
-                k += l << 2;
-                a[k - l] = -t[j] - t[j + 1];
-                a[k + l] = t[j] - t[j + 1];
-            }
-            l <<= 1;
-            mh = m >> 1;
-            for (j = 1; j < mh; j++) {
-                k = m - j;
-                t[j] = t[m + k] + t[m + j];
-                t[k] = t[m + k] - t[m + j];
-            }
-            t[0] = t[m + mh];
-            m = mh;
-        }
-        a[l] = t[0];
-    }
-    a[0] = 0;
-}
-
-
-/* -------- initializing routines -------- */
-
-
-#include <math.h>
-
-void makewt(int nw, int *ip, double *w)
-{
-    void makeipt(int nw, int *ip);
-    int j, nwh, nw0, nw1;
-    double delta, wn4r, wk1r, wk1i, wk3r, wk3i;
-    
-    ip[0] = nw;
-    ip[1] = 1;
-    if (nw > 2) {
-        nwh = nw >> 1;
-        delta = atan(1.0) / nwh;
-        wn4r = cos(delta * nwh);
-        w[0] = 1;
-        w[1] = wn4r;
-        if (nwh == 4) {
-            w[2] = cos(delta * 2);
-            w[3] = sin(delta * 2);
-        } else if (nwh > 4) {
-            makeipt(nw, ip);
-            w[2] = 0.5 / cos(delta * 2);
-            w[3] = 0.5 / cos(delta * 6);
-            for (j = 4; j < nwh; j += 4) {
-                w[j] = cos(delta * j);
-                w[j + 1] = sin(delta * j);
-                w[j + 2] = cos(3 * delta * j);
-                w[j + 3] = -sin(3 * delta * j);
-            }
-        }
-        nw0 = 0;
-        while (nwh > 2) {
-            nw1 = nw0 + nwh;
-            nwh >>= 1;
-            w[nw1] = 1;
-            w[nw1 + 1] = wn4r;
-            if (nwh == 4) {
-                wk1r = w[nw0 + 4];
-                wk1i = w[nw0 + 5];
-                w[nw1 + 2] = wk1r;
-                w[nw1 + 3] = wk1i;
-            } else if (nwh > 4) {
-                wk1r = w[nw0 + 4];
-                wk3r = w[nw0 + 6];
-                w[nw1 + 2] = 0.5 / wk1r;
-                w[nw1 + 3] = 0.5 / wk3r;
-                for (j = 4; j < nwh; j += 4) {
-                    wk1r = w[nw0 + 2 * j];
-                    wk1i = w[nw0 + 2 * j + 1];
-                    wk3r = w[nw0 + 2 * j + 2];
-                    wk3i = w[nw0 + 2 * j + 3];
-                    w[nw1 + j] = wk1r;
-                    w[nw1 + j + 1] = wk1i;
-                    w[nw1 + j + 2] = wk3r;
-                    w[nw1 + j + 3] = wk3i;
-                }
-            }
-            nw0 = nw1;
-        }
-    }
-}
-
-
-void makeipt(int nw, int *ip)
-{
-    int j, l, m, m2, p, q;
-    
-    ip[2] = 0;
-    ip[3] = 16;
-    m = 2;
-    for (l = nw; l > 32; l >>= 2) {
-        m2 = m << 1;
-        q = m2 << 3;
-        for (j = m; j < m2; j++) {
-            p = ip[j] << 2;
-            ip[m + j] = p;
-            ip[m2 + j] = p + q;
-        }
-        m = m2;
-    }
-}
-
-
-void makect(int nc, int *ip, double *c)
-{
-    int j, nch;
-    double delta;
-    
-    ip[1] = nc;
-    if (nc > 1) {
-        nch = nc >> 1;
-        delta = atan(1.0) / nch;
-        c[0] = cos(delta * nch);
-        c[nch] = 0.5 * c[0];
-        for (j = 1; j < nch; j++) {
-            c[j] = 0.5 * cos(delta * j);
-            c[nc - j] = 0.5 * sin(delta * j);
-        }
-    }
-}
-
-
-/* -------- child routines -------- */
-
-
-#ifdef USE_CDFT_PTHREADS
-#define USE_CDFT_THREADS
-#ifndef CDFT_THREADS_BEGIN_N
-#define CDFT_THREADS_BEGIN_N 8192
-#endif
-#ifndef CDFT_4THREADS_BEGIN_N
-#define CDFT_4THREADS_BEGIN_N 65536
-#endif
-#include <pthread.h>
-#include <stdio.h>
-#include <stdlib.h>
-#define cdft_thread_t pthread_t
-#define cdft_thread_create(thp,func,argp) { \
-    if (pthread_create(thp, NULL, func, (void *) argp) != 0) { \
-        fprintf(stderr, "cdft thread error\n"); \
-        exit(1); \
-    } \
-}
-#define cdft_thread_wait(th) { \
-    if (pthread_join(th, NULL) != 0) { \
-        fprintf(stderr, "cdft thread error\n"); \
-        exit(1); \
-    } \
-}
-#endif /* USE_CDFT_PTHREADS */
-
-
-#ifdef USE_CDFT_WINTHREADS
-#define USE_CDFT_THREADS
-#ifndef CDFT_THREADS_BEGIN_N
-#define CDFT_THREADS_BEGIN_N 32768
-#endif
-#ifndef CDFT_4THREADS_BEGIN_N
-#define CDFT_4THREADS_BEGIN_N 524288
-#endif
-#include <windows.h>
-#include <stdio.h>
-#include <stdlib.h>
-#define cdft_thread_t HANDLE
-#define cdft_thread_create(thp,func,argp) { \
-    DWORD thid; \
-    *(thp) = CreateThread(NULL, 0, (LPTHREAD_START_ROUTINE) func, (LPVOID) argp, 0, &thid); \
-    if (*(thp) == 0) { \
-        fprintf(stderr, "cdft thread error\n"); \
-        exit(1); \
-    } \
-}
-#define cdft_thread_wait(th) { \
-    WaitForSingleObject(th, INFINITE); \
-    CloseHandle(th); \
-}
-#endif /* USE_CDFT_WINTHREADS */
-
-
-void cftfsub(int n, double *a, int *ip, int nw, double *w)
-{
-    void bitrv2(int n, int *ip, double *a);
-    void bitrv216(double *a);
-    void bitrv208(double *a);
-    void cftf1st(int n, double *a, double *w);
-    void cftrec4(int n, double *a, int nw, double *w);
-    void cftleaf(int n, int isplt, double *a, int nw, double *w);
-    void cftfx41(int n, double *a, int nw, double *w);
-    void cftf161(double *a, double *w);
-    void cftf081(double *a, double *w);
-    void cftf040(double *a);
-    void cftx020(double *a);
-#ifdef USE_CDFT_THREADS
-    void cftrec4_th(int n, double *a, int nw, double *w);
-#endif /* USE_CDFT_THREADS */
-    
-    if (n > 8) {
-        if (n > 32) {
-            cftf1st(n, a, &w[nw - (n >> 2)]);
-#ifdef USE_CDFT_THREADS
-            if (n > CDFT_THREADS_BEGIN_N) {
-                cftrec4_th(n, a, nw, w);
-            } else 
-#endif /* USE_CDFT_THREADS */
-            if (n > 512) {
-                cftrec4(n, a, nw, w);
-            } else if (n > 128) {
-                cftleaf(n, 1, a, nw, w);
-            } else {
-                cftfx41(n, a, nw, w);
-            }
-            bitrv2(n, ip, a);
-        } else if (n == 32) {
-            cftf161(a, &w[nw - 8]);
-            bitrv216(a);
-        } else {
-            cftf081(a, w);
-            bitrv208(a);
-        }
-    } else if (n == 8) {
-        cftf040(a);
-    } else if (n == 4) {
-        cftx020(a);
-    }
-}
-
-
-void cftbsub(int n, double *a, int *ip, int nw, double *w)
-{
-    void bitrv2conj(int n, int *ip, double *a);
-    void bitrv216neg(double *a);
-    void bitrv208neg(double *a);
-    void cftb1st(int n, double *a, double *w);
-    void cftrec4(int n, double *a, int nw, double *w);
-    void cftleaf(int n, int isplt, double *a, int nw, double *w);
-    void cftfx41(int n, double *a, int nw, double *w);
-    void cftf161(double *a, double *w);
-    void cftf081(double *a, double *w);
-    void cftb040(double *a);
-    void cftx020(double *a);
-#ifdef USE_CDFT_THREADS
-    void cftrec4_th(int n, double *a, int nw, double *w);
-#endif /* USE_CDFT_THREADS */
-    
-    if (n > 8) {
-        if (n > 32) {
-            cftb1st(n, a, &w[nw - (n >> 2)]);
-#ifdef USE_CDFT_THREADS
-            if (n > CDFT_THREADS_BEGIN_N) {
-                cftrec4_th(n, a, nw, w);
-            } else 
-#endif /* USE_CDFT_THREADS */
-            if (n > 512) {
-                cftrec4(n, a, nw, w);
-            } else if (n > 128) {
-                cftleaf(n, 1, a, nw, w);
-            } else {
-                cftfx41(n, a, nw, w);
-            }
-            bitrv2conj(n, ip, a);
-        } else if (n == 32) {
-            cftf161(a, &w[nw - 8]);
-            bitrv216neg(a);
-        } else {
-            cftf081(a, w);
-            bitrv208neg(a);
-        }
-    } else if (n == 8) {
-        cftb040(a);
-    } else if (n == 4) {
-        cftx020(a);
-    }
-}
-
-
-void bitrv2(int n, int *ip, double *a)
-{
-    int j, j1, k, k1, l, m, nh, nm;
-    double xr, xi, yr, yi;
-    
-    m = 1;
-    for (l = n >> 2; l > 8; l >>= 2) {
-        m <<= 1;
-    }
-    nh = n >> 1;
-    nm = 4 * m;
-    if (l == 8) {
-        for (k = 0; k < m; k++) {
-            for (j = 0; j < k; j++) {
-                j1 = 4 * j + 2 * ip[m + k];
-                k1 = 4 * k + 2 * ip[m + j];
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nm;
-                k1 += 2 * nm;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nm;
-                k1 -= nm;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nm;
-                k1 += 2 * nm;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nh;
-                k1 += 2;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nm;
-                k1 -= 2 * nm;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nm;
-                k1 += nm;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nm;
-                k1 -= 2 * nm;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += 2;
-                k1 += nh;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nm;
-                k1 += 2 * nm;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nm;
-                k1 -= nm;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nm;
-                k1 += 2 * nm;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nh;
-                k1 -= 2;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nm;
-                k1 -= 2 * nm;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nm;
-                k1 += nm;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nm;
-                k1 -= 2 * nm;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-            }
-            k1 = 4 * k + 2 * ip[m + k];
-            j1 = k1 + 2;
-            k1 += nh;
-            xr = a[j1];
-            xi = a[j1 + 1];
-            yr = a[k1];
-            yi = a[k1 + 1];
-            a[j1] = yr;
-            a[j1 + 1] = yi;
-            a[k1] = xr;
-            a[k1 + 1] = xi;
-            j1 += nm;
-            k1 += 2 * nm;
-            xr = a[j1];
-            xi = a[j1 + 1];
-            yr = a[k1];
-            yi = a[k1 + 1];
-            a[j1] = yr;
-            a[j1 + 1] = yi;
-            a[k1] = xr;
-            a[k1 + 1] = xi;
-            j1 += nm;
-            k1 -= nm;
-            xr = a[j1];
-            xi = a[j1 + 1];
-            yr = a[k1];
-            yi = a[k1 + 1];
-            a[j1] = yr;
-            a[j1 + 1] = yi;
-            a[k1] = xr;
-            a[k1 + 1] = xi;
-            j1 -= 2;
-            k1 -= nh;
-            xr = a[j1];
-            xi = a[j1 + 1];
-            yr = a[k1];
-            yi = a[k1 + 1];
-            a[j1] = yr;
-            a[j1 + 1] = yi;
-            a[k1] = xr;
-            a[k1 + 1] = xi;
-            j1 += nh + 2;
-            k1 += nh + 2;
-            xr = a[j1];
-            xi = a[j1 + 1];
-            yr = a[k1];
-            yi = a[k1 + 1];
-            a[j1] = yr;
-            a[j1 + 1] = yi;
-            a[k1] = xr;
-            a[k1 + 1] = xi;
-            j1 -= nh - nm;
-            k1 += 2 * nm - 2;
-            xr = a[j1];
-            xi = a[j1 + 1];
-            yr = a[k1];
-            yi = a[k1 + 1];
-            a[j1] = yr;
-            a[j1 + 1] = yi;
-            a[k1] = xr;
-            a[k1 + 1] = xi;
-        }
-    } else {
-        for (k = 0; k < m; k++) {
-            for (j = 0; j < k; j++) {
-                j1 = 4 * j + ip[m + k];
-                k1 = 4 * k + ip[m + j];
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nm;
-                k1 += nm;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nh;
-                k1 += 2;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nm;
-                k1 -= nm;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += 2;
-                k1 += nh;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nm;
-                k1 += nm;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nh;
-                k1 -= 2;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nm;
-                k1 -= nm;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-            }
-            k1 = 4 * k + ip[m + k];
-            j1 = k1 + 2;
-            k1 += nh;
-            xr = a[j1];
-            xi = a[j1 + 1];
-            yr = a[k1];
-            yi = a[k1 + 1];
-            a[j1] = yr;
-            a[j1 + 1] = yi;
-            a[k1] = xr;
-            a[k1 + 1] = xi;
-            j1 += nm;
-            k1 += nm;
-            xr = a[j1];
-            xi = a[j1 + 1];
-            yr = a[k1];
-            yi = a[k1 + 1];
-            a[j1] = yr;
-            a[j1 + 1] = yi;
-            a[k1] = xr;
-            a[k1 + 1] = xi;
-        }
-    }
-}
-
-
-void bitrv2conj(int n, int *ip, double *a)
-{
-    int j, j1, k, k1, l, m, nh, nm;
-    double xr, xi, yr, yi;
-    
-    m = 1;
-    for (l = n >> 2; l > 8; l >>= 2) {
-        m <<= 1;
-    }
-    nh = n >> 1;
-    nm = 4 * m;
-    if (l == 8) {
-        for (k = 0; k < m; k++) {
-            for (j = 0; j < k; j++) {
-                j1 = 4 * j + 2 * ip[m + k];
-                k1 = 4 * k + 2 * ip[m + j];
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nm;
-                k1 += 2 * nm;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nm;
-                k1 -= nm;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nm;
-                k1 += 2 * nm;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nh;
-                k1 += 2;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nm;
-                k1 -= 2 * nm;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nm;
-                k1 += nm;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nm;
-                k1 -= 2 * nm;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += 2;
-                k1 += nh;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nm;
-                k1 += 2 * nm;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nm;
-                k1 -= nm;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nm;
-                k1 += 2 * nm;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nh;
-                k1 -= 2;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nm;
-                k1 -= 2 * nm;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nm;
-                k1 += nm;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nm;
-                k1 -= 2 * nm;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-            }
-            k1 = 4 * k + 2 * ip[m + k];
-            j1 = k1 + 2;
-            k1 += nh;
-            a[j1 - 1] = -a[j1 - 1];
-            xr = a[j1];
-            xi = -a[j1 + 1];
-            yr = a[k1];
-            yi = -a[k1 + 1];
-            a[j1] = yr;
-            a[j1 + 1] = yi;
-            a[k1] = xr;
-            a[k1 + 1] = xi;
-            a[k1 + 3] = -a[k1 + 3];
-            j1 += nm;
-            k1 += 2 * nm;
-            xr = a[j1];
-            xi = -a[j1 + 1];
-            yr = a[k1];
-            yi = -a[k1 + 1];
-            a[j1] = yr;
-            a[j1 + 1] = yi;
-            a[k1] = xr;
-            a[k1 + 1] = xi;
-            j1 += nm;
-            k1 -= nm;
-            xr = a[j1];
-            xi = -a[j1 + 1];
-            yr = a[k1];
-            yi = -a[k1 + 1];
-            a[j1] = yr;
-            a[j1 + 1] = yi;
-            a[k1] = xr;
-            a[k1 + 1] = xi;
-            j1 -= 2;
-            k1 -= nh;
-            xr = a[j1];
-            xi = -a[j1 + 1];
-            yr = a[k1];
-            yi = -a[k1 + 1];
-            a[j1] = yr;
-            a[j1 + 1] = yi;
-            a[k1] = xr;
-            a[k1 + 1] = xi;
-            j1 += nh + 2;
-            k1 += nh + 2;
-            xr = a[j1];
-            xi = -a[j1 + 1];
-            yr = a[k1];
-            yi = -a[k1 + 1];
-            a[j1] = yr;
-            a[j1 + 1] = yi;
-            a[k1] = xr;
-            a[k1 + 1] = xi;
-            j1 -= nh - nm;
-            k1 += 2 * nm - 2;
-            a[j1 - 1] = -a[j1 - 1];
-            xr = a[j1];
-            xi = -a[j1 + 1];
-            yr = a[k1];
-            yi = -a[k1 + 1];
-            a[j1] = yr;
-            a[j1 + 1] = yi;
-            a[k1] = xr;
-            a[k1 + 1] = xi;
-            a[k1 + 3] = -a[k1 + 3];
-        }
-    } else {
-        for (k = 0; k < m; k++) {
-            for (j = 0; j < k; j++) {
-                j1 = 4 * j + ip[m + k];
-                k1 = 4 * k + ip[m + j];
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nm;
-                k1 += nm;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nh;
-                k1 += 2;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nm;
-                k1 -= nm;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += 2;
-                k1 += nh;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += nm;
-                k1 += nm;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nh;
-                k1 -= 2;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 -= nm;
-                k1 -= nm;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-            }
-            k1 = 4 * k + ip[m + k];
-            j1 = k1 + 2;
-            k1 += nh;
-            a[j1 - 1] = -a[j1 - 1];
-            xr = a[j1];
-            xi = -a[j1 + 1];
-            yr = a[k1];
-            yi = -a[k1 + 1];
-            a[j1] = yr;
-            a[j1 + 1] = yi;
-            a[k1] = xr;
-            a[k1 + 1] = xi;
-            a[k1 + 3] = -a[k1 + 3];
-            j1 += nm;
-            k1 += nm;
-            a[j1 - 1] = -a[j1 - 1];
-            xr = a[j1];
-            xi = -a[j1 + 1];
-            yr = a[k1];
-            yi = -a[k1 + 1];
-            a[j1] = yr;
-            a[j1 + 1] = yi;
-            a[k1] = xr;
-            a[k1 + 1] = xi;
-            a[k1 + 3] = -a[k1 + 3];
-        }
-    }
-}
-
-
-void bitrv216(double *a)
-{
-    double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
-        x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i, 
-        x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i;
-    
-    x1r = a[2];
-    x1i = a[3];
-    x2r = a[4];
-    x2i = a[5];
-    x3r = a[6];
-    x3i = a[7];
-    x4r = a[8];
-    x4i = a[9];
-    x5r = a[10];
-    x5i = a[11];
-    x7r = a[14];
-    x7i = a[15];
-    x8r = a[16];
-    x8i = a[17];
-    x10r = a[20];
-    x10i = a[21];
-    x11r = a[22];
-    x11i = a[23];
-    x12r = a[24];
-    x12i = a[25];
-    x13r = a[26];
-    x13i = a[27];
-    x14r = a[28];
-    x14i = a[29];
-    a[2] = x8r;
-    a[3] = x8i;
-    a[4] = x4r;
-    a[5] = x4i;
-    a[6] = x12r;
-    a[7] = x12i;
-    a[8] = x2r;
-    a[9] = x2i;
-    a[10] = x10r;
-    a[11] = x10i;
-    a[14] = x14r;
-    a[15] = x14i;
-    a[16] = x1r;
-    a[17] = x1i;
-    a[20] = x5r;
-    a[21] = x5i;
-    a[22] = x13r;
-    a[23] = x13i;
-    a[24] = x3r;
-    a[25] = x3i;
-    a[26] = x11r;
-    a[27] = x11i;
-    a[28] = x7r;
-    a[29] = x7i;
-}
-
-
-void bitrv216neg(double *a)
-{
-    double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
-        x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i, 
-        x9r, x9i, x10r, x10i, x11r, x11i, x12r, x12i, 
-        x13r, x13i, x14r, x14i, x15r, x15i;
-    
-    x1r = a[2];
-    x1i = a[3];
-    x2r = a[4];
-    x2i = a[5];
-    x3r = a[6];
-    x3i = a[7];
-    x4r = a[8];
-    x4i = a[9];
-    x5r = a[10];
-    x5i = a[11];
-    x6r = a[12];
-    x6i = a[13];
-    x7r = a[14];
-    x7i = a[15];
-    x8r = a[16];
-    x8i = a[17];
-    x9r = a[18];
-    x9i = a[19];
-    x10r = a[20];
-    x10i = a[21];
-    x11r = a[22];
-    x11i = a[23];
-    x12r = a[24];
-    x12i = a[25];
-    x13r = a[26];
-    x13i = a[27];
-    x14r = a[28];
-    x14i = a[29];
-    x15r = a[30];
-    x15i = a[31];
-    a[2] = x15r;
-    a[3] = x15i;
-    a[4] = x7r;
-    a[5] = x7i;
-    a[6] = x11r;
-    a[7] = x11i;
-    a[8] = x3r;
-    a[9] = x3i;
-    a[10] = x13r;
-    a[11] = x13i;
-    a[12] = x5r;
-    a[13] = x5i;
-    a[14] = x9r;
-    a[15] = x9i;
-    a[16] = x1r;
-    a[17] = x1i;
-    a[18] = x14r;
-    a[19] = x14i;
-    a[20] = x6r;
-    a[21] = x6i;
-    a[22] = x10r;
-    a[23] = x10i;
-    a[24] = x2r;
-    a[25] = x2i;
-    a[26] = x12r;
-    a[27] = x12i;
-    a[28] = x4r;
-    a[29] = x4i;
-    a[30] = x8r;
-    a[31] = x8i;
-}
-
-
-void bitrv208(double *a)
-{
-    double x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i;
-    
-    x1r = a[2];
-    x1i = a[3];
-    x3r = a[6];
-    x3i = a[7];
-    x4r = a[8];
-    x4i = a[9];
-    x6r = a[12];
-    x6i = a[13];
-    a[2] = x4r;
-    a[3] = x4i;
-    a[6] = x6r;
-    a[7] = x6i;
-    a[8] = x1r;
-    a[9] = x1i;
-    a[12] = x3r;
-    a[13] = x3i;
-}
-
-
-void bitrv208neg(double *a)
-{
-    double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
-        x5r, x5i, x6r, x6i, x7r, x7i;
-    
-    x1r = a[2];
-    x1i = a[3];
-    x2r = a[4];
-    x2i = a[5];
-    x3r = a[6];
-    x3i = a[7];
-    x4r = a[8];
-    x4i = a[9];
-    x5r = a[10];
-    x5i = a[11];
-    x6r = a[12];
-    x6i = a[13];
-    x7r = a[14];
-    x7i = a[15];
-    a[2] = x7r;
-    a[3] = x7i;
-    a[4] = x3r;
-    a[5] = x3i;
-    a[6] = x5r;
-    a[7] = x5i;
-    a[8] = x1r;
-    a[9] = x1i;
-    a[10] = x6r;
-    a[11] = x6i;
-    a[12] = x2r;
-    a[13] = x2i;
-    a[14] = x4r;
-    a[15] = x4i;
-}
-
-
-void cftf1st(int n, double *a, double *w)
-{
-    int j, j0, j1, j2, j3, k, m, mh;
-    double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i, 
-        wd1r, wd1i, wd3r, wd3i;
-    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
-        y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
-    
-    mh = n >> 3;
-    m = 2 * mh;
-    j1 = m;
-    j2 = j1 + m;
-    j3 = j2 + m;
-    x0r = a[0] + a[j2];
-    x0i = a[1] + a[j2 + 1];
-    x1r = a[0] - a[j2];
-    x1i = a[1] - a[j2 + 1];
-    x2r = a[j1] + a[j3];
-    x2i = a[j1 + 1] + a[j3 + 1];
-    x3r = a[j1] - a[j3];
-    x3i = a[j1 + 1] - a[j3 + 1];
-    a[0] = x0r + x2r;
-    a[1] = x0i + x2i;
-    a[j1] = x0r - x2r;
-    a[j1 + 1] = x0i - x2i;
-    a[j2] = x1r - x3i;
-    a[j2 + 1] = x1i + x3r;
-    a[j3] = x1r + x3i;
-    a[j3 + 1] = x1i - x3r;
-    wn4r = w[1];
-    csc1 = w[2];
-    csc3 = w[3];
-    wd1r = 1;
-    wd1i = 0;
-    wd3r = 1;
-    wd3i = 0;
-    k = 0;
-    for (j = 2; j < mh - 2; j += 4) {
-        k += 4;
-        wk1r = csc1 * (wd1r + w[k]);
-        wk1i = csc1 * (wd1i + w[k + 1]);
-        wk3r = csc3 * (wd3r + w[k + 2]);
-        wk3i = csc3 * (wd3i + w[k + 3]);
-        wd1r = w[k];
-        wd1i = w[k + 1];
-        wd3r = w[k + 2];
-        wd3i = w[k + 3];
-        j1 = j + m;
-        j2 = j1 + m;
-        j3 = j2 + m;
-        x0r = a[j] + a[j2];
-        x0i = a[j + 1] + a[j2 + 1];
-        x1r = a[j] - a[j2];
-        x1i = a[j + 1] - a[j2 + 1];
-        y0r = a[j + 2] + a[j2 + 2];
-        y0i = a[j + 3] + a[j2 + 3];
-        y1r = a[j + 2] - a[j2 + 2];
-        y1i = a[j + 3] - a[j2 + 3];
-        x2r = a[j1] + a[j3];
-        x2i = a[j1 + 1] + a[j3 + 1];
-        x3r = a[j1] - a[j3];
-        x3i = a[j1 + 1] - a[j3 + 1];
-        y2r = a[j1 + 2] + a[j3 + 2];
-        y2i = a[j1 + 3] + a[j3 + 3];
-        y3r = a[j1 + 2] - a[j3 + 2];
-        y3i = a[j1 + 3] - a[j3 + 3];
-        a[j] = x0r + x2r;
-        a[j + 1] = x0i + x2i;
-        a[j + 2] = y0r + y2r;
-        a[j + 3] = y0i + y2i;
-        a[j1] = x0r - x2r;
-        a[j1 + 1] = x0i - x2i;
-        a[j1 + 2] = y0r - y2r;
-        a[j1 + 3] = y0i - y2i;
-        x0r = x1r - x3i;
-        x0i = x1i + x3r;
-        a[j2] = wk1r * x0r - wk1i * x0i;
-        a[j2 + 1] = wk1r * x0i + wk1i * x0r;
-        x0r = y1r - y3i;
-        x0i = y1i + y3r;
-        a[j2 + 2] = wd1r * x0r - wd1i * x0i;
-        a[j2 + 3] = wd1r * x0i + wd1i * x0r;
-        x0r = x1r + x3i;
-        x0i = x1i - x3r;
-        a[j3] = wk3r * x0r + wk3i * x0i;
-        a[j3 + 1] = wk3r * x0i - wk3i * x0r;
-        x0r = y1r + y3i;
-        x0i = y1i - y3r;
-        a[j3 + 2] = wd3r * x0r + wd3i * x0i;
-        a[j3 + 3] = wd3r * x0i - wd3i * x0r;
-        j0 = m - j;
-        j1 = j0 + m;
-        j2 = j1 + m;
-        j3 = j2 + m;
-        x0r = a[j0] + a[j2];
-        x0i = a[j0 + 1] + a[j2 + 1];
-        x1r = a[j0] - a[j2];
-        x1i = a[j0 + 1] - a[j2 + 1];
-        y0r = a[j0 - 2] + a[j2 - 2];
-        y0i = a[j0 - 1] + a[j2 - 1];
-        y1r = a[j0 - 2] - a[j2 - 2];
-        y1i = a[j0 - 1] - a[j2 - 1];
-        x2r = a[j1] + a[j3];
-        x2i = a[j1 + 1] + a[j3 + 1];
-        x3r = a[j1] - a[j3];
-        x3i = a[j1 + 1] - a[j3 + 1];
-        y2r = a[j1 - 2] + a[j3 - 2];
-        y2i = a[j1 - 1] + a[j3 - 1];
-        y3r = a[j1 - 2] - a[j3 - 2];
-        y3i = a[j1 - 1] - a[j3 - 1];
-        a[j0] = x0r + x2r;
-        a[j0 + 1] = x0i + x2i;
-        a[j0 - 2] = y0r + y2r;
-        a[j0 - 1] = y0i + y2i;
-        a[j1] = x0r - x2r;
-        a[j1 + 1] = x0i - x2i;
-        a[j1 - 2] = y0r - y2r;
-        a[j1 - 1] = y0i - y2i;
-        x0r = x1r - x3i;
-        x0i = x1i + x3r;
-        a[j2] = wk1i * x0r - wk1r * x0i;
-        a[j2 + 1] = wk1i * x0i + wk1r * x0r;
-        x0r = y1r - y3i;
-        x0i = y1i + y3r;
-        a[j2 - 2] = wd1i * x0r - wd1r * x0i;
-        a[j2 - 1] = wd1i * x0i + wd1r * x0r;
-        x0r = x1r + x3i;
-        x0i = x1i - x3r;
-        a[j3] = wk3i * x0r + wk3r * x0i;
-        a[j3 + 1] = wk3i * x0i - wk3r * x0r;
-        x0r = y1r + y3i;
-        x0i = y1i - y3r;
-        a[j3 - 2] = wd3i * x0r + wd3r * x0i;
-        a[j3 - 1] = wd3i * x0i - wd3r * x0r;
-    }
-    wk1r = csc1 * (wd1r + wn4r);
-    wk1i = csc1 * (wd1i + wn4r);
-    wk3r = csc3 * (wd3r - wn4r);
-    wk3i = csc3 * (wd3i - wn4r);
-    j0 = mh;
-    j1 = j0 + m;
-    j2 = j1 + m;
-    j3 = j2 + m;
-    x0r = a[j0 - 2] + a[j2 - 2];
-    x0i = a[j0 - 1] + a[j2 - 1];
-    x1r = a[j0 - 2] - a[j2 - 2];
-    x1i = a[j0 - 1] - a[j2 - 1];
-    x2r = a[j1 - 2] + a[j3 - 2];
-    x2i = a[j1 - 1] + a[j3 - 1];
-    x3r = a[j1 - 2] - a[j3 - 2];
-    x3i = a[j1 - 1] - a[j3 - 1];
-    a[j0 - 2] = x0r + x2r;
-    a[j0 - 1] = x0i + x2i;
-    a[j1 - 2] = x0r - x2r;
-    a[j1 - 1] = x0i - x2i;
-    x0r = x1r - x3i;
-    x0i = x1i + x3r;
-    a[j2 - 2] = wk1r * x0r - wk1i * x0i;
-    a[j2 - 1] = wk1r * x0i + wk1i * x0r;
-    x0r = x1r + x3i;
-    x0i = x1i - x3r;
-    a[j3 - 2] = wk3r * x0r + wk3i * x0i;
-    a[j3 - 1] = wk3r * x0i - wk3i * x0r;
-    x0r = a[j0] + a[j2];
-    x0i = a[j0 + 1] + a[j2 + 1];
-    x1r = a[j0] - a[j2];
-    x1i = a[j0 + 1] - a[j2 + 1];
-    x2r = a[j1] + a[j3];
-    x2i = a[j1 + 1] + a[j3 + 1];
-    x3r = a[j1] - a[j3];
-    x3i = a[j1 + 1] - a[j3 + 1];
-    a[j0] = x0r + x2r;
-    a[j0 + 1] = x0i + x2i;
-    a[j1] = x0r - x2r;
-    a[j1 + 1] = x0i - x2i;
-    x0r = x1r - x3i;
-    x0i = x1i + x3r;
-    a[j2] = wn4r * (x0r - x0i);
-    a[j2 + 1] = wn4r * (x0i + x0r);
-    x0r = x1r + x3i;
-    x0i = x1i - x3r;
-    a[j3] = -wn4r * (x0r + x0i);
-    a[j3 + 1] = -wn4r * (x0i - x0r);
-    x0r = a[j0 + 2] + a[j2 + 2];
-    x0i = a[j0 + 3] + a[j2 + 3];
-    x1r = a[j0 + 2] - a[j2 + 2];
-    x1i = a[j0 + 3] - a[j2 + 3];
-    x2r = a[j1 + 2] + a[j3 + 2];
-    x2i = a[j1 + 3] + a[j3 + 3];
-    x3r = a[j1 + 2] - a[j3 + 2];
-    x3i = a[j1 + 3] - a[j3 + 3];
-    a[j0 + 2] = x0r + x2r;
-    a[j0 + 3] = x0i + x2i;
-    a[j1 + 2] = x0r - x2r;
-    a[j1 + 3] = x0i - x2i;
-    x0r = x1r - x3i;
-    x0i = x1i + x3r;
-    a[j2 + 2] = wk1i * x0r - wk1r * x0i;
-    a[j2 + 3] = wk1i * x0i + wk1r * x0r;
-    x0r = x1r + x3i;
-    x0i = x1i - x3r;
-    a[j3 + 2] = wk3i * x0r + wk3r * x0i;
-    a[j3 + 3] = wk3i * x0i - wk3r * x0r;
-}
-
-
-void cftb1st(int n, double *a, double *w)
-{
-    int j, j0, j1, j2, j3, k, m, mh;
-    double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i, 
-        wd1r, wd1i, wd3r, wd3i;
-    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
-        y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
-    
-    mh = n >> 3;
-    m = 2 * mh;
-    j1 = m;
-    j2 = j1 + m;
-    j3 = j2 + m;
-    x0r = a[0] + a[j2];
-    x0i = -a[1] - a[j2 + 1];
-    x1r = a[0] - a[j2];
-    x1i = -a[1] + a[j2 + 1];
-    x2r = a[j1] + a[j3];
-    x2i = a[j1 + 1] + a[j3 + 1];
-    x3r = a[j1] - a[j3];
-    x3i = a[j1 + 1] - a[j3 + 1];
-    a[0] = x0r + x2r;
-    a[1] = x0i - x2i;
-    a[j1] = x0r - x2r;
-    a[j1 + 1] = x0i + x2i;
-    a[j2] = x1r + x3i;
-    a[j2 + 1] = x1i + x3r;
-    a[j3] = x1r - x3i;
-    a[j3 + 1] = x1i - x3r;
-    wn4r = w[1];
-    csc1 = w[2];
-    csc3 = w[3];
-    wd1r = 1;
-    wd1i = 0;
-    wd3r = 1;
-    wd3i = 0;
-    k = 0;
-    for (j = 2; j < mh - 2; j += 4) {
-        k += 4;
-        wk1r = csc1 * (wd1r + w[k]);
-        wk1i = csc1 * (wd1i + w[k + 1]);
-        wk3r = csc3 * (wd3r + w[k + 2]);
-        wk3i = csc3 * (wd3i + w[k + 3]);
-        wd1r = w[k];
-        wd1i = w[k + 1];
-        wd3r = w[k + 2];
-        wd3i = w[k + 3];
-        j1 = j + m;
-        j2 = j1 + m;
-        j3 = j2 + m;
-        x0r = a[j] + a[j2];
-        x0i = -a[j + 1] - a[j2 + 1];
-        x1r = a[j] - a[j2];
-        x1i = -a[j + 1] + a[j2 + 1];
-        y0r = a[j + 2] + a[j2 + 2];
-        y0i = -a[j + 3] - a[j2 + 3];
-        y1r = a[j + 2] - a[j2 + 2];
-        y1i = -a[j + 3] + a[j2 + 3];
-        x2r = a[j1] + a[j3];
-        x2i = a[j1 + 1] + a[j3 + 1];
-        x3r = a[j1] - a[j3];
-        x3i = a[j1 + 1] - a[j3 + 1];
-        y2r = a[j1 + 2] + a[j3 + 2];
-        y2i = a[j1 + 3] + a[j3 + 3];
-        y3r = a[j1 + 2] - a[j3 + 2];
-        y3i = a[j1 + 3] - a[j3 + 3];
-        a[j] = x0r + x2r;
-        a[j + 1] = x0i - x2i;
-        a[j + 2] = y0r + y2r;
-        a[j + 3] = y0i - y2i;
-        a[j1] = x0r - x2r;
-        a[j1 + 1] = x0i + x2i;
-        a[j1 + 2] = y0r - y2r;
-        a[j1 + 3] = y0i + y2i;
-        x0r = x1r + x3i;
-        x0i = x1i + x3r;
-        a[j2] = wk1r * x0r - wk1i * x0i;
-        a[j2 + 1] = wk1r * x0i + wk1i * x0r;
-        x0r = y1r + y3i;
-        x0i = y1i + y3r;
-        a[j2 + 2] = wd1r * x0r - wd1i * x0i;
-        a[j2 + 3] = wd1r * x0i + wd1i * x0r;
-        x0r = x1r - x3i;
-        x0i = x1i - x3r;
-        a[j3] = wk3r * x0r + wk3i * x0i;
-        a[j3 + 1] = wk3r * x0i - wk3i * x0r;
-        x0r = y1r - y3i;
-        x0i = y1i - y3r;
-        a[j3 + 2] = wd3r * x0r + wd3i * x0i;
-        a[j3 + 3] = wd3r * x0i - wd3i * x0r;
-        j0 = m - j;
-        j1 = j0 + m;
-        j2 = j1 + m;
-        j3 = j2 + m;
-        x0r = a[j0] + a[j2];
-        x0i = -a[j0 + 1] - a[j2 + 1];
-        x1r = a[j0] - a[j2];
-        x1i = -a[j0 + 1] + a[j2 + 1];
-        y0r = a[j0 - 2] + a[j2 - 2];
-        y0i = -a[j0 - 1] - a[j2 - 1];
-        y1r = a[j0 - 2] - a[j2 - 2];
-        y1i = -a[j0 - 1] + a[j2 - 1];
-        x2r = a[j1] + a[j3];
-        x2i = a[j1 + 1] + a[j3 + 1];
-        x3r = a[j1] - a[j3];
-        x3i = a[j1 + 1] - a[j3 + 1];
-        y2r = a[j1 - 2] + a[j3 - 2];
-        y2i = a[j1 - 1] + a[j3 - 1];
-        y3r = a[j1 - 2] - a[j3 - 2];
-        y3i = a[j1 - 1] - a[j3 - 1];
-        a[j0] = x0r + x2r;
-        a[j0 + 1] = x0i - x2i;
-        a[j0 - 2] = y0r + y2r;
-        a[j0 - 1] = y0i - y2i;
-        a[j1] = x0r - x2r;
-        a[j1 + 1] = x0i + x2i;
-        a[j1 - 2] = y0r - y2r;
-        a[j1 - 1] = y0i + y2i;
-        x0r = x1r + x3i;
-        x0i = x1i + x3r;
-        a[j2] = wk1i * x0r - wk1r * x0i;
-        a[j2 + 1] = wk1i * x0i + wk1r * x0r;
-        x0r = y1r + y3i;
-        x0i = y1i + y3r;
-        a[j2 - 2] = wd1i * x0r - wd1r * x0i;
-        a[j2 - 1] = wd1i * x0i + wd1r * x0r;
-        x0r = x1r - x3i;
-        x0i = x1i - x3r;
-        a[j3] = wk3i * x0r + wk3r * x0i;
-        a[j3 + 1] = wk3i * x0i - wk3r * x0r;
-        x0r = y1r - y3i;
-        x0i = y1i - y3r;
-        a[j3 - 2] = wd3i * x0r + wd3r * x0i;
-        a[j3 - 1] = wd3i * x0i - wd3r * x0r;
-    }
-    wk1r = csc1 * (wd1r + wn4r);
-    wk1i = csc1 * (wd1i + wn4r);
-    wk3r = csc3 * (wd3r - wn4r);
-    wk3i = csc3 * (wd3i - wn4r);
-    j0 = mh;
-    j1 = j0 + m;
-    j2 = j1 + m;
-    j3 = j2 + m;
-    x0r = a[j0 - 2] + a[j2 - 2];
-    x0i = -a[j0 - 1] - a[j2 - 1];
-    x1r = a[j0 - 2] - a[j2 - 2];
-    x1i = -a[j0 - 1] + a[j2 - 1];
-    x2r = a[j1 - 2] + a[j3 - 2];
-    x2i = a[j1 - 1] + a[j3 - 1];
-    x3r = a[j1 - 2] - a[j3 - 2];
-    x3i = a[j1 - 1] - a[j3 - 1];
-    a[j0 - 2] = x0r + x2r;
-    a[j0 - 1] = x0i - x2i;
-    a[j1 - 2] = x0r - x2r;
-    a[j1 - 1] = x0i + x2i;
-    x0r = x1r + x3i;
-    x0i = x1i + x3r;
-    a[j2 - 2] = wk1r * x0r - wk1i * x0i;
-    a[j2 - 1] = wk1r * x0i + wk1i * x0r;
-    x0r = x1r - x3i;
-    x0i = x1i - x3r;
-    a[j3 - 2] = wk3r * x0r + wk3i * x0i;
-    a[j3 - 1] = wk3r * x0i - wk3i * x0r;
-    x0r = a[j0] + a[j2];
-    x0i = -a[j0 + 1] - a[j2 + 1];
-    x1r = a[j0] - a[j2];
-    x1i = -a[j0 + 1] + a[j2 + 1];
-    x2r = a[j1] + a[j3];
-    x2i = a[j1 + 1] + a[j3 + 1];
-    x3r = a[j1] - a[j3];
-    x3i = a[j1 + 1] - a[j3 + 1];
-    a[j0] = x0r + x2r;
-    a[j0 + 1] = x0i - x2i;
-    a[j1] = x0r - x2r;
-    a[j1 + 1] = x0i + x2i;
-    x0r = x1r + x3i;
-    x0i = x1i + x3r;
-    a[j2] = wn4r * (x0r - x0i);
-    a[j2 + 1] = wn4r * (x0i + x0r);
-    x0r = x1r - x3i;
-    x0i = x1i - x3r;
-    a[j3] = -wn4r * (x0r + x0i);
-    a[j3 + 1] = -wn4r * (x0i - x0r);
-    x0r = a[j0 + 2] + a[j2 + 2];
-    x0i = -a[j0 + 3] - a[j2 + 3];
-    x1r = a[j0 + 2] - a[j2 + 2];
-    x1i = -a[j0 + 3] + a[j2 + 3];
-    x2r = a[j1 + 2] + a[j3 + 2];
-    x2i = a[j1 + 3] + a[j3 + 3];
-    x3r = a[j1 + 2] - a[j3 + 2];
-    x3i = a[j1 + 3] - a[j3 + 3];
-    a[j0 + 2] = x0r + x2r;
-    a[j0 + 3] = x0i - x2i;
-    a[j1 + 2] = x0r - x2r;
-    a[j1 + 3] = x0i + x2i;
-    x0r = x1r + x3i;
-    x0i = x1i + x3r;
-    a[j2 + 2] = wk1i * x0r - wk1r * x0i;
-    a[j2 + 3] = wk1i * x0i + wk1r * x0r;
-    x0r = x1r - x3i;
-    x0i = x1i - x3r;
-    a[j3 + 2] = wk3i * x0r + wk3r * x0i;
-    a[j3 + 3] = wk3i * x0i - wk3r * x0r;
-}
-
-
-#ifdef USE_CDFT_THREADS
-struct cdft_arg_st {
-    int n0;
-    int n;
-    double *a;
-    int nw;
-    double *w;
-};
-typedef struct cdft_arg_st cdft_arg_t;
-
-
-void cftrec4_th(int n, double *a, int nw, double *w)
-{
-    void *cftrec1_th(void *p);
-    void *cftrec2_th(void *p);
-    int i, idiv4, m, nthread;
-    cdft_thread_t th[4];
-    cdft_arg_t ag[4];
-    
-    nthread = 2;
-    idiv4 = 0;
-    m = n >> 1;
-    if (n > CDFT_4THREADS_BEGIN_N) {
-        nthread = 4;
-        idiv4 = 1;
-        m >>= 1;
-    }
-    for (i = 0; i < nthread; i++) {
-        ag[i].n0 = n;
-        ag[i].n = m;
-        ag[i].a = &a[i * m];
-        ag[i].nw = nw;
-        ag[i].w = w;
-        if (i != idiv4) {
-            cdft_thread_create(&th[i], cftrec1_th, &ag[i]);
-        } else {
-            cdft_thread_create(&th[i], cftrec2_th, &ag[i]);
-        }
-    }
-    for (i = 0; i < nthread; i++) {
-        cdft_thread_wait(th[i]);
-    }
-}
-
-
-void *cftrec1_th(void *p)
-{
-    int cfttree(int n, int j, int k, double *a, int nw, double *w);
-    void cftleaf(int n, int isplt, double *a, int nw, double *w);
-    void cftmdl1(int n, double *a, double *w);
-    int isplt, j, k, m, n, n0, nw;
-    double *a, *w;
-    
-    n0 = ((cdft_arg_t *) p)->n0;
-    n = ((cdft_arg_t *) p)->n;
-    a = ((cdft_arg_t *) p)->a;
-    nw = ((cdft_arg_t *) p)->nw;
-    w = ((cdft_arg_t *) p)->w;
-    m = n0;
-    while (m > 512) {
-        m >>= 2;
-        cftmdl1(m, &a[n - m], &w[nw - (m >> 1)]);
-    }
-    cftleaf(m, 1, &a[n - m], nw, w);
-    k = 0;
-    for (j = n - m; j > 0; j -= m) {
-        k++;
-        isplt = cfttree(m, j, k, a, nw, w);
-        cftleaf(m, isplt, &a[j - m], nw, w);
-    }
-    return (void *) 0;
-}
-
-
-void *cftrec2_th(void *p)
-{
-    int cfttree(int n, int j, int k, double *a, int nw, double *w);
-    void cftleaf(int n, int isplt, double *a, int nw, double *w);
-    void cftmdl2(int n, double *a, double *w);
-    int isplt, j, k, m, n, n0, nw;
-    double *a, *w;
-    
-    n0 = ((cdft_arg_t *) p)->n0;
-    n = ((cdft_arg_t *) p)->n;
-    a = ((cdft_arg_t *) p)->a;
-    nw = ((cdft_arg_t *) p)->nw;
-    w = ((cdft_arg_t *) p)->w;
-    k = 1;
-    m = n0;
-    while (m > 512) {
-        m >>= 2;
-        k <<= 2;
-        cftmdl2(m, &a[n - m], &w[nw - m]);
-    }
-    cftleaf(m, 0, &a[n - m], nw, w);
-    k >>= 1;
-    for (j = n - m; j > 0; j -= m) {
-        k++;
-        isplt = cfttree(m, j, k, a, nw, w);
-        cftleaf(m, isplt, &a[j - m], nw, w);
-    }
-    return (void *) 0;
-}
-#endif /* USE_CDFT_THREADS */
-
-
-void cftrec4(int n, double *a, int nw, double *w)
-{
-    int cfttree(int n, int j, int k, double *a, int nw, double *w);
-    void cftleaf(int n, int isplt, double *a, int nw, double *w);
-    void cftmdl1(int n, double *a, double *w);
-    int isplt, j, k, m;
-    
-    m = n;
-    while (m > 512) {
-        m >>= 2;
-        cftmdl1(m, &a[n - m], &w[nw - (m >> 1)]);
-    }
-    cftleaf(m, 1, &a[n - m], nw, w);
-    k = 0;
-    for (j = n - m; j > 0; j -= m) {
-        k++;
-        isplt = cfttree(m, j, k, a, nw, w);
-        cftleaf(m, isplt, &a[j - m], nw, w);
-    }
-}
-
-
-int cfttree(int n, int j, int k, double *a, int nw, double *w)
-{
-    void cftmdl1(int n, double *a, double *w);
-    void cftmdl2(int n, double *a, double *w);
-    int i, isplt, m;
-    
-    if ((k & 3) != 0) {
-        isplt = k & 1;
-        if (isplt != 0) {
-            cftmdl1(n, &a[j - n], &w[nw - (n >> 1)]);
-        } else {
-            cftmdl2(n, &a[j - n], &w[nw - n]);
-        }
-    } else {
-        m = n;
-        for (i = k; (i & 3) == 0; i >>= 2) {
-            m <<= 2;
-        }
-        isplt = i & 1;
-        if (isplt != 0) {
-            while (m > 128) {
-                cftmdl1(m, &a[j - m], &w[nw - (m >> 1)]);
-                m >>= 2;
-            }
-        } else {
-            while (m > 128) {
-                cftmdl2(m, &a[j - m], &w[nw - m]);
-                m >>= 2;
-            }
-        }
-    }
-    return isplt;
-}
-
-
-void cftleaf(int n, int isplt, double *a, int nw, double *w)
-{
-    void cftmdl1(int n, double *a, double *w);
-    void cftmdl2(int n, double *a, double *w);
-    void cftf161(double *a, double *w);
-    void cftf162(double *a, double *w);
-    void cftf081(double *a, double *w);
-    void cftf082(double *a, double *w);
-    
-    if (n == 512) {
-        cftmdl1(128, a, &w[nw - 64]);
-        cftf161(a, &w[nw - 8]);
-        cftf162(&a[32], &w[nw - 32]);
-        cftf161(&a[64], &w[nw - 8]);
-        cftf161(&a[96], &w[nw - 8]);
-        cftmdl2(128, &a[128], &w[nw - 128]);
-        cftf161(&a[128], &w[nw - 8]);
-        cftf162(&a[160], &w[nw - 32]);
-        cftf161(&a[192], &w[nw - 8]);
-        cftf162(&a[224], &w[nw - 32]);
-        cftmdl1(128, &a[256], &w[nw - 64]);
-        cftf161(&a[256], &w[nw - 8]);
-        cftf162(&a[288], &w[nw - 32]);
-        cftf161(&a[320], &w[nw - 8]);
-        cftf161(&a[352], &w[nw - 8]);
-        if (isplt != 0) {
-            cftmdl1(128, &a[384], &w[nw - 64]);
-            cftf161(&a[480], &w[nw - 8]);
-        } else {
-            cftmdl2(128, &a[384], &w[nw - 128]);
-            cftf162(&a[480], &w[nw - 32]);
-        }
-        cftf161(&a[384], &w[nw - 8]);
-        cftf162(&a[416], &w[nw - 32]);
-        cftf161(&a[448], &w[nw - 8]);
-    } else {
-        cftmdl1(64, a, &w[nw - 32]);
-        cftf081(a, &w[nw - 8]);
-        cftf082(&a[16], &w[nw - 8]);
-        cftf081(&a[32], &w[nw - 8]);
-        cftf081(&a[48], &w[nw - 8]);
-        cftmdl2(64, &a[64], &w[nw - 64]);
-        cftf081(&a[64], &w[nw - 8]);
-        cftf082(&a[80], &w[nw - 8]);
-        cftf081(&a[96], &w[nw - 8]);
-        cftf082(&a[112], &w[nw - 8]);
-        cftmdl1(64, &a[128], &w[nw - 32]);
-        cftf081(&a[128], &w[nw - 8]);
-        cftf082(&a[144], &w[nw - 8]);
-        cftf081(&a[160], &w[nw - 8]);
-        cftf081(&a[176], &w[nw - 8]);
-        if (isplt != 0) {
-            cftmdl1(64, &a[192], &w[nw - 32]);
-            cftf081(&a[240], &w[nw - 8]);
-        } else {
-            cftmdl2(64, &a[192], &w[nw - 64]);
-            cftf082(&a[240], &w[nw - 8]);
-        }
-        cftf081(&a[192], &w[nw - 8]);
-        cftf082(&a[208], &w[nw - 8]);
-        cftf081(&a[224], &w[nw - 8]);
-    }
-}
-
-
-void cftmdl1(int n, double *a, double *w)
-{
-    int j, j0, j1, j2, j3, k, m, mh;
-    double wn4r, wk1r, wk1i, wk3r, wk3i;
-    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
-    
-    mh = n >> 3;
-    m = 2 * mh;
-    j1 = m;
-    j2 = j1 + m;
-    j3 = j2 + m;
-    x0r = a[0] + a[j2];
-    x0i = a[1] + a[j2 + 1];
-    x1r = a[0] - a[j2];
-    x1i = a[1] - a[j2 + 1];
-    x2r = a[j1] + a[j3];
-    x2i = a[j1 + 1] + a[j3 + 1];
-    x3r = a[j1] - a[j3];
-    x3i = a[j1 + 1] - a[j3 + 1];
-    a[0] = x0r + x2r;
-    a[1] = x0i + x2i;
-    a[j1] = x0r - x2r;
-    a[j1 + 1] = x0i - x2i;
-    a[j2] = x1r - x3i;
-    a[j2 + 1] = x1i + x3r;
-    a[j3] = x1r + x3i;
-    a[j3 + 1] = x1i - x3r;
-    wn4r = w[1];
-    k = 0;
-    for (j = 2; j < mh; j += 2) {
-        k += 4;
-        wk1r = w[k];
-        wk1i = w[k + 1];
-        wk3r = w[k + 2];
-        wk3i = w[k + 3];
-        j1 = j + m;
-        j2 = j1 + m;
-        j3 = j2 + m;
-        x0r = a[j] + a[j2];
-        x0i = a[j + 1] + a[j2 + 1];
-        x1r = a[j] - a[j2];
-        x1i = a[j + 1] - a[j2 + 1];
-        x2r = a[j1] + a[j3];
-        x2i = a[j1 + 1] + a[j3 + 1];
-        x3r = a[j1] - a[j3];
-        x3i = a[j1 + 1] - a[j3 + 1];
-        a[j] = x0r + x2r;
-        a[j + 1] = x0i + x2i;
-        a[j1] = x0r - x2r;
-        a[j1 + 1] = x0i - x2i;
-        x0r = x1r - x3i;
-        x0i = x1i + x3r;
-        a[j2] = wk1r * x0r - wk1i * x0i;
-        a[j2 + 1] = wk1r * x0i + wk1i * x0r;
-        x0r = x1r + x3i;
-        x0i = x1i - x3r;
-        a[j3] = wk3r * x0r + wk3i * x0i;
-        a[j3 + 1] = wk3r * x0i - wk3i * x0r;
-        j0 = m - j;
-        j1 = j0 + m;
-        j2 = j1 + m;
-        j3 = j2 + m;
-        x0r = a[j0] + a[j2];
-        x0i = a[j0 + 1] + a[j2 + 1];
-        x1r = a[j0] - a[j2];
-        x1i = a[j0 + 1] - a[j2 + 1];
-        x2r = a[j1] + a[j3];
-        x2i = a[j1 + 1] + a[j3 + 1];
-        x3r = a[j1] - a[j3];
-        x3i = a[j1 + 1] - a[j3 + 1];
-        a[j0] = x0r + x2r;
-        a[j0 + 1] = x0i + x2i;
-        a[j1] = x0r - x2r;
-        a[j1 + 1] = x0i - x2i;
-        x0r = x1r - x3i;
-        x0i = x1i + x3r;
-        a[j2] = wk1i * x0r - wk1r * x0i;
-        a[j2 + 1] = wk1i * x0i + wk1r * x0r;
-        x0r = x1r + x3i;
-        x0i = x1i - x3r;
-        a[j3] = wk3i * x0r + wk3r * x0i;
-        a[j3 + 1] = wk3i * x0i - wk3r * x0r;
-    }
-    j0 = mh;
-    j1 = j0 + m;
-    j2 = j1 + m;
-    j3 = j2 + m;
-    x0r = a[j0] + a[j2];
-    x0i = a[j0 + 1] + a[j2 + 1];
-    x1r = a[j0] - a[j2];
-    x1i = a[j0 + 1] - a[j2 + 1];
-    x2r = a[j1] + a[j3];
-    x2i = a[j1 + 1] + a[j3 + 1];
-    x3r = a[j1] - a[j3];
-    x3i = a[j1 + 1] - a[j3 + 1];
-    a[j0] = x0r + x2r;
-    a[j0 + 1] = x0i + x2i;
-    a[j1] = x0r - x2r;
-    a[j1 + 1] = x0i - x2i;
-    x0r = x1r - x3i;
-    x0i = x1i + x3r;
-    a[j2] = wn4r * (x0r - x0i);
-    a[j2 + 1] = wn4r * (x0i + x0r);
-    x0r = x1r + x3i;
-    x0i = x1i - x3r;
-    a[j3] = -wn4r * (x0r + x0i);
-    a[j3 + 1] = -wn4r * (x0i - x0r);
-}
-
-
-void cftmdl2(int n, double *a, double *w)
-{
-    int j, j0, j1, j2, j3, k, kr, m, mh;
-    double wn4r, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i;
-    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y2r, y2i;
-    
-    mh = n >> 3;
-    m = 2 * mh;
-    wn4r = w[1];
-    j1 = m;
-    j2 = j1 + m;
-    j3 = j2 + m;
-    x0r = a[0] - a[j2 + 1];
-    x0i = a[1] + a[j2];
-    x1r = a[0] + a[j2 + 1];
-    x1i = a[1] - a[j2];
-    x2r = a[j1] - a[j3 + 1];
-    x2i = a[j1 + 1] + a[j3];
-    x3r = a[j1] + a[j3 + 1];
-    x3i = a[j1 + 1] - a[j3];
-    y0r = wn4r * (x2r - x2i);
-    y0i = wn4r * (x2i + x2r);
-    a[0] = x0r + y0r;
-    a[1] = x0i + y0i;
-    a[j1] = x0r - y0r;
-    a[j1 + 1] = x0i - y0i;
-    y0r = wn4r * (x3r - x3i);
-    y0i = wn4r * (x3i + x3r);
-    a[j2] = x1r - y0i;
-    a[j2 + 1] = x1i + y0r;
-    a[j3] = x1r + y0i;
-    a[j3 + 1] = x1i - y0r;
-    k = 0;
-    kr = 2 * m;
-    for (j = 2; j < mh; j += 2) {
-        k += 4;
-        wk1r = w[k];
-        wk1i = w[k + 1];
-        wk3r = w[k + 2];
-        wk3i = w[k + 3];
-        kr -= 4;
-        wd1i = w[kr];
-        wd1r = w[kr + 1];
-        wd3i = w[kr + 2];
-        wd3r = w[kr + 3];
-        j1 = j + m;
-        j2 = j1 + m;
-        j3 = j2 + m;
-        x0r = a[j] - a[j2 + 1];
-        x0i = a[j + 1] + a[j2];
-        x1r = a[j] + a[j2 + 1];
-        x1i = a[j + 1] - a[j2];
-        x2r = a[j1] - a[j3 + 1];
-        x2i = a[j1 + 1] + a[j3];
-        x3r = a[j1] + a[j3 + 1];
-        x3i = a[j1 + 1] - a[j3];
-        y0r = wk1r * x0r - wk1i * x0i;
-        y0i = wk1r * x0i + wk1i * x0r;
-        y2r = wd1r * x2r - wd1i * x2i;
-        y2i = wd1r * x2i + wd1i * x2r;
-        a[j] = y0r + y2r;
-        a[j + 1] = y0i + y2i;
-        a[j1] = y0r - y2r;
-        a[j1 + 1] = y0i - y2i;
-        y0r = wk3r * x1r + wk3i * x1i;
-        y0i = wk3r * x1i - wk3i * x1r;
-        y2r = wd3r * x3r + wd3i * x3i;
-        y2i = wd3r * x3i - wd3i * x3r;
-        a[j2] = y0r + y2r;
-        a[j2 + 1] = y0i + y2i;
-        a[j3] = y0r - y2r;
-        a[j3 + 1] = y0i - y2i;
-        j0 = m - j;
-        j1 = j0 + m;
-        j2 = j1 + m;
-        j3 = j2 + m;
-        x0r = a[j0] - a[j2 + 1];
-        x0i = a[j0 + 1] + a[j2];
-        x1r = a[j0] + a[j2 + 1];
-        x1i = a[j0 + 1] - a[j2];
-        x2r = a[j1] - a[j3 + 1];
-        x2i = a[j1 + 1] + a[j3];
-        x3r = a[j1] + a[j3 + 1];
-        x3i = a[j1 + 1] - a[j3];
-        y0r = wd1i * x0r - wd1r * x0i;
-        y0i = wd1i * x0i + wd1r * x0r;
-        y2r = wk1i * x2r - wk1r * x2i;
-        y2i = wk1i * x2i + wk1r * x2r;
-        a[j0] = y0r + y2r;
-        a[j0 + 1] = y0i + y2i;
-        a[j1] = y0r - y2r;
-        a[j1 + 1] = y0i - y2i;
-        y0r = wd3i * x1r + wd3r * x1i;
-        y0i = wd3i * x1i - wd3r * x1r;
-        y2r = wk3i * x3r + wk3r * x3i;
-        y2i = wk3i * x3i - wk3r * x3r;
-        a[j2] = y0r + y2r;
-        a[j2 + 1] = y0i + y2i;
-        a[j3] = y0r - y2r;
-        a[j3 + 1] = y0i - y2i;
-    }
-    wk1r = w[m];
-    wk1i = w[m + 1];
-    j0 = mh;
-    j1 = j0 + m;
-    j2 = j1 + m;
-    j3 = j2 + m;
-    x0r = a[j0] - a[j2 + 1];
-    x0i = a[j0 + 1] + a[j2];
-    x1r = a[j0] + a[j2 + 1];
-    x1i = a[j0 + 1] - a[j2];
-    x2r = a[j1] - a[j3 + 1];
-    x2i = a[j1 + 1] + a[j3];
-    x3r = a[j1] + a[j3 + 1];
-    x3i = a[j1 + 1] - a[j3];
-    y0r = wk1r * x0r - wk1i * x0i;
-    y0i = wk1r * x0i + wk1i * x0r;
-    y2r = wk1i * x2r - wk1r * x2i;
-    y2i = wk1i * x2i + wk1r * x2r;
-    a[j0] = y0r + y2r;
-    a[j0 + 1] = y0i + y2i;
-    a[j1] = y0r - y2r;
-    a[j1 + 1] = y0i - y2i;
-    y0r = wk1i * x1r - wk1r * x1i;
-    y0i = wk1i * x1i + wk1r * x1r;
-    y2r = wk1r * x3r - wk1i * x3i;
-    y2i = wk1r * x3i + wk1i * x3r;
-    a[j2] = y0r - y2r;
-    a[j2 + 1] = y0i - y2i;
-    a[j3] = y0r + y2r;
-    a[j3 + 1] = y0i + y2i;
-}
-
-
-void cftfx41(int n, double *a, int nw, double *w)
-{
-    void cftf161(double *a, double *w);
-    void cftf162(double *a, double *w);
-    void cftf081(double *a, double *w);
-    void cftf082(double *a, double *w);
-    
-    if (n == 128) {
-        cftf161(a, &w[nw - 8]);
-        cftf162(&a[32], &w[nw - 32]);
-        cftf161(&a[64], &w[nw - 8]);
-        cftf161(&a[96], &w[nw - 8]);
-    } else {
-        cftf081(a, &w[nw - 8]);
-        cftf082(&a[16], &w[nw - 8]);
-        cftf081(&a[32], &w[nw - 8]);
-        cftf081(&a[48], &w[nw - 8]);
-    }
-}
-
-
-void cftf161(double *a, double *w)
-{
-    double wn4r, wk1r, wk1i, 
-        x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
-        y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
-        y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i, 
-        y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i, 
-        y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
-    
-    wn4r = w[1];
-    wk1r = w[2];
-    wk1i = w[3];
-    x0r = a[0] + a[16];
-    x0i = a[1] + a[17];
-    x1r = a[0] - a[16];
-    x1i = a[1] - a[17];
-    x2r = a[8] + a[24];
-    x2i = a[9] + a[25];
-    x3r = a[8] - a[24];
-    x3i = a[9] - a[25];
-    y0r = x0r + x2r;
-    y0i = x0i + x2i;
-    y4r = x0r - x2r;
-    y4i = x0i - x2i;
-    y8r = x1r - x3i;
-    y8i = x1i + x3r;
-    y12r = x1r + x3i;
-    y12i = x1i - x3r;
-    x0r = a[2] + a[18];
-    x0i = a[3] + a[19];
-    x1r = a[2] - a[18];
-    x1i = a[3] - a[19];
-    x2r = a[10] + a[26];
-    x2i = a[11] + a[27];
-    x3r = a[10] - a[26];
-    x3i = a[11] - a[27];
-    y1r = x0r + x2r;
-    y1i = x0i + x2i;
-    y5r = x0r - x2r;
-    y5i = x0i - x2i;
-    x0r = x1r - x3i;
-    x0i = x1i + x3r;
-    y9r = wk1r * x0r - wk1i * x0i;
-    y9i = wk1r * x0i + wk1i * x0r;
-    x0r = x1r + x3i;
-    x0i = x1i - x3r;
-    y13r = wk1i * x0r - wk1r * x0i;
-    y13i = wk1i * x0i + wk1r * x0r;
-    x0r = a[4] + a[20];
-    x0i = a[5] + a[21];
-    x1r = a[4] - a[20];
-    x1i = a[5] - a[21];
-    x2r = a[12] + a[28];
-    x2i = a[13] + a[29];
-    x3r = a[12] - a[28];
-    x3i = a[13] - a[29];
-    y2r = x0r + x2r;
-    y2i = x0i + x2i;
-    y6r = x0r - x2r;
-    y6i = x0i - x2i;
-    x0r = x1r - x3i;
-    x0i = x1i + x3r;
-    y10r = wn4r * (x0r - x0i);
-    y10i = wn4r * (x0i + x0r);
-    x0r = x1r + x3i;
-    x0i = x1i - x3r;
-    y14r = wn4r * (x0r + x0i);
-    y14i = wn4r * (x0i - x0r);
-    x0r = a[6] + a[22];
-    x0i = a[7] + a[23];
-    x1r = a[6] - a[22];
-    x1i = a[7] - a[23];
-    x2r = a[14] + a[30];
-    x2i = a[15] + a[31];
-    x3r = a[14] - a[30];
-    x3i = a[15] - a[31];
-    y3r = x0r + x2r;
-    y3i = x0i + x2i;
-    y7r = x0r - x2r;
-    y7i = x0i - x2i;
-    x0r = x1r - x3i;
-    x0i = x1i + x3r;
-    y11r = wk1i * x0r - wk1r * x0i;
-    y11i = wk1i * x0i + wk1r * x0r;
-    x0r = x1r + x3i;
-    x0i = x1i - x3r;
-    y15r = wk1r * x0r - wk1i * x0i;
-    y15i = wk1r * x0i + wk1i * x0r;
-    x0r = y12r - y14r;
-    x0i = y12i - y14i;
-    x1r = y12r + y14r;
-    x1i = y12i + y14i;
-    x2r = y13r - y15r;
-    x2i = y13i - y15i;
-    x3r = y13r + y15r;
-    x3i = y13i + y15i;
-    a[24] = x0r + x2r;
-    a[25] = x0i + x2i;
-    a[26] = x0r - x2r;
-    a[27] = x0i - x2i;
-    a[28] = x1r - x3i;
-    a[29] = x1i + x3r;
-    a[30] = x1r + x3i;
-    a[31] = x1i - x3r;
-    x0r = y8r + y10r;
-    x0i = y8i + y10i;
-    x1r = y8r - y10r;
-    x1i = y8i - y10i;
-    x2r = y9r + y11r;
-    x2i = y9i + y11i;
-    x3r = y9r - y11r;
-    x3i = y9i - y11i;
-    a[16] = x0r + x2r;
-    a[17] = x0i + x2i;
-    a[18] = x0r - x2r;
-    a[19] = x0i - x2i;
-    a[20] = x1r - x3i;
-    a[21] = x1i + x3r;
-    a[22] = x1r + x3i;
-    a[23] = x1i - x3r;
-    x0r = y5r - y7i;
-    x0i = y5i + y7r;
-    x2r = wn4r * (x0r - x0i);
-    x2i = wn4r * (x0i + x0r);
-    x0r = y5r + y7i;
-    x0i = y5i - y7r;
-    x3r = wn4r * (x0r - x0i);
-    x3i = wn4r * (x0i + x0r);
-    x0r = y4r - y6i;
-    x0i = y4i + y6r;
-    x1r = y4r + y6i;
-    x1i = y4i - y6r;
-    a[8] = x0r + x2r;
-    a[9] = x0i + x2i;
-    a[10] = x0r - x2r;
-    a[11] = x0i - x2i;
-    a[12] = x1r - x3i;
-    a[13] = x1i + x3r;
-    a[14] = x1r + x3i;
-    a[15] = x1i - x3r;
-    x0r = y0r + y2r;
-    x0i = y0i + y2i;
-    x1r = y0r - y2r;
-    x1i = y0i - y2i;
-    x2r = y1r + y3r;
-    x2i = y1i + y3i;
-    x3r = y1r - y3r;
-    x3i = y1i - y3i;
-    a[0] = x0r + x2r;
-    a[1] = x0i + x2i;
-    a[2] = x0r - x2r;
-    a[3] = x0i - x2i;
-    a[4] = x1r - x3i;
-    a[5] = x1i + x3r;
-    a[6] = x1r + x3i;
-    a[7] = x1i - x3r;
-}
-
-
-void cftf162(double *a, double *w)
-{
-    double wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i, 
-        x0r, x0i, x1r, x1i, x2r, x2i, 
-        y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
-        y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i, 
-        y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i, 
-        y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
-    
-    wn4r = w[1];
-    wk1r = w[4];
-    wk1i = w[5];
-    wk3r = w[6];
-    wk3i = -w[7];
-    wk2r = w[8];
-    wk2i = w[9];
-    x1r = a[0] - a[17];
-    x1i = a[1] + a[16];
-    x0r = a[8] - a[25];
-    x0i = a[9] + a[24];
-    x2r = wn4r * (x0r - x0i);
-    x2i = wn4r * (x0i + x0r);
-    y0r = x1r + x2r;
-    y0i = x1i + x2i;
-    y4r = x1r - x2r;
-    y4i = x1i - x2i;
-    x1r = a[0] + a[17];
-    x1i = a[1] - a[16];
-    x0r = a[8] + a[25];
-    x0i = a[9] - a[24];
-    x2r = wn4r * (x0r - x0i);
-    x2i = wn4r * (x0i + x0r);
-    y8r = x1r - x2i;
-    y8i = x1i + x2r;
-    y12r = x1r + x2i;
-    y12i = x1i - x2r;
-    x0r = a[2] - a[19];
-    x0i = a[3] + a[18];
-    x1r = wk1r * x0r - wk1i * x0i;
-    x1i = wk1r * x0i + wk1i * x0r;
-    x0r = a[10] - a[27];
-    x0i = a[11] + a[26];
-    x2r = wk3i * x0r - wk3r * x0i;
-    x2i = wk3i * x0i + wk3r * x0r;
-    y1r = x1r + x2r;
-    y1i = x1i + x2i;
-    y5r = x1r - x2r;
-    y5i = x1i - x2i;
-    x0r = a[2] + a[19];
-    x0i = a[3] - a[18];
-    x1r = wk3r * x0r - wk3i * x0i;
-    x1i = wk3r * x0i + wk3i * x0r;
-    x0r = a[10] + a[27];
-    x0i = a[11] - a[26];
-    x2r = wk1r * x0r + wk1i * x0i;
-    x2i = wk1r * x0i - wk1i * x0r;
-    y9r = x1r - x2r;
-    y9i = x1i - x2i;
-    y13r = x1r + x2r;
-    y13i = x1i + x2i;
-    x0r = a[4] - a[21];
-    x0i = a[5] + a[20];
-    x1r = wk2r * x0r - wk2i * x0i;
-    x1i = wk2r * x0i + wk2i * x0r;
-    x0r = a[12] - a[29];
-    x0i = a[13] + a[28];
-    x2r = wk2i * x0r - wk2r * x0i;
-    x2i = wk2i * x0i + wk2r * x0r;
-    y2r = x1r + x2r;
-    y2i = x1i + x2i;
-    y6r = x1r - x2r;
-    y6i = x1i - x2i;
-    x0r = a[4] + a[21];
-    x0i = a[5] - a[20];
-    x1r = wk2i * x0r - wk2r * x0i;
-    x1i = wk2i * x0i + wk2r * x0r;
-    x0r = a[12] + a[29];
-    x0i = a[13] - a[28];
-    x2r = wk2r * x0r - wk2i * x0i;
-    x2i = wk2r * x0i + wk2i * x0r;
-    y10r = x1r - x2r;
-    y10i = x1i - x2i;
-    y14r = x1r + x2r;
-    y14i = x1i + x2i;
-    x0r = a[6] - a[23];
-    x0i = a[7] + a[22];
-    x1r = wk3r * x0r - wk3i * x0i;
-    x1i = wk3r * x0i + wk3i * x0r;
-    x0r = a[14] - a[31];
-    x0i = a[15] + a[30];
-    x2r = wk1i * x0r - wk1r * x0i;
-    x2i = wk1i * x0i + wk1r * x0r;
-    y3r = x1r + x2r;
-    y3i = x1i + x2i;
-    y7r = x1r - x2r;
-    y7i = x1i - x2i;
-    x0r = a[6] + a[23];
-    x0i = a[7] - a[22];
-    x1r = wk1i * x0r + wk1r * x0i;
-    x1i = wk1i * x0i - wk1r * x0r;
-    x0r = a[14] + a[31];
-    x0i = a[15] - a[30];
-    x2r = wk3i * x0r - wk3r * x0i;
-    x2i = wk3i * x0i + wk3r * x0r;
-    y11r = x1r + x2r;
-    y11i = x1i + x2i;
-    y15r = x1r - x2r;
-    y15i = x1i - x2i;
-    x1r = y0r + y2r;
-    x1i = y0i + y2i;
-    x2r = y1r + y3r;
-    x2i = y1i + y3i;
-    a[0] = x1r + x2r;
-    a[1] = x1i + x2i;
-    a[2] = x1r - x2r;
-    a[3] = x1i - x2i;
-    x1r = y0r - y2r;
-    x1i = y0i - y2i;
-    x2r = y1r - y3r;
-    x2i = y1i - y3i;
-    a[4] = x1r - x2i;
-    a[5] = x1i + x2r;
-    a[6] = x1r + x2i;
-    a[7] = x1i - x2r;
-    x1r = y4r - y6i;
-    x1i = y4i + y6r;
-    x0r = y5r - y7i;
-    x0i = y5i + y7r;
-    x2r = wn4r * (x0r - x0i);
-    x2i = wn4r * (x0i + x0r);
-    a[8] = x1r + x2r;
-    a[9] = x1i + x2i;
-    a[10] = x1r - x2r;
-    a[11] = x1i - x2i;
-    x1r = y4r + y6i;
-    x1i = y4i - y6r;
-    x0r = y5r + y7i;
-    x0i = y5i - y7r;
-    x2r = wn4r * (x0r - x0i);
-    x2i = wn4r * (x0i + x0r);
-    a[12] = x1r - x2i;
-    a[13] = x1i + x2r;
-    a[14] = x1r + x2i;
-    a[15] = x1i - x2r;
-    x1r = y8r + y10r;
-    x1i = y8i + y10i;
-    x2r = y9r - y11r;
-    x2i = y9i - y11i;
-    a[16] = x1r + x2r;
-    a[17] = x1i + x2i;
-    a[18] = x1r - x2r;
-    a[19] = x1i - x2i;
-    x1r = y8r - y10r;
-    x1i = y8i - y10i;
-    x2r = y9r + y11r;
-    x2i = y9i + y11i;
-    a[20] = x1r - x2i;
-    a[21] = x1i + x2r;
-    a[22] = x1r + x2i;
-    a[23] = x1i - x2r;
-    x1r = y12r - y14i;
-    x1i = y12i + y14r;
-    x0r = y13r + y15i;
-    x0i = y13i - y15r;
-    x2r = wn4r * (x0r - x0i);
-    x2i = wn4r * (x0i + x0r);
-    a[24] = x1r + x2r;
-    a[25] = x1i + x2i;
-    a[26] = x1r - x2r;
-    a[27] = x1i - x2i;
-    x1r = y12r + y14i;
-    x1i = y12i - y14r;
-    x0r = y13r - y15i;
-    x0i = y13i + y15r;
-    x2r = wn4r * (x0r - x0i);
-    x2i = wn4r * (x0i + x0r);
-    a[28] = x1r - x2i;
-    a[29] = x1i + x2r;
-    a[30] = x1r + x2i;
-    a[31] = x1i - x2r;
-}
-
-
-void cftf081(double *a, double *w)
-{
-    double wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
-        y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
-        y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
-    
-    wn4r = w[1];
-    x0r = a[0] + a[8];
-    x0i = a[1] + a[9];
-    x1r = a[0] - a[8];
-    x1i = a[1] - a[9];
-    x2r = a[4] + a[12];
-    x2i = a[5] + a[13];
-    x3r = a[4] - a[12];
-    x3i = a[5] - a[13];
-    y0r = x0r + x2r;
-    y0i = x0i + x2i;
-    y2r = x0r - x2r;
-    y2i = x0i - x2i;
-    y1r = x1r - x3i;
-    y1i = x1i + x3r;
-    y3r = x1r + x3i;
-    y3i = x1i - x3r;
-    x0r = a[2] + a[10];
-    x0i = a[3] + a[11];
-    x1r = a[2] - a[10];
-    x1i = a[3] - a[11];
-    x2r = a[6] + a[14];
-    x2i = a[7] + a[15];
-    x3r = a[6] - a[14];
-    x3i = a[7] - a[15];
-    y4r = x0r + x2r;
-    y4i = x0i + x2i;
-    y6r = x0r - x2r;
-    y6i = x0i - x2i;
-    x0r = x1r - x3i;
-    x0i = x1i + x3r;
-    x2r = x1r + x3i;
-    x2i = x1i - x3r;
-    y5r = wn4r * (x0r - x0i);
-    y5i = wn4r * (x0r + x0i);
-    y7r = wn4r * (x2r - x2i);
-    y7i = wn4r * (x2r + x2i);
-    a[8] = y1r + y5r;
-    a[9] = y1i + y5i;
-    a[10] = y1r - y5r;
-    a[11] = y1i - y5i;
-    a[12] = y3r - y7i;
-    a[13] = y3i + y7r;
-    a[14] = y3r + y7i;
-    a[15] = y3i - y7r;
-    a[0] = y0r + y4r;
-    a[1] = y0i + y4i;
-    a[2] = y0r - y4r;
-    a[3] = y0i - y4i;
-    a[4] = y2r - y6i;
-    a[5] = y2i + y6r;
-    a[6] = y2r + y6i;
-    a[7] = y2i - y6r;
-}
-
-
-void cftf082(double *a, double *w)
-{
-    double wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i, 
-        y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
-        y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
-    
-    wn4r = w[1];
-    wk1r = w[2];
-    wk1i = w[3];
-    y0r = a[0] - a[9];
-    y0i = a[1] + a[8];
-    y1r = a[0] + a[9];
-    y1i = a[1] - a[8];
-    x0r = a[4] - a[13];
-    x0i = a[5] + a[12];
-    y2r = wn4r * (x0r - x0i);
-    y2i = wn4r * (x0i + x0r);
-    x0r = a[4] + a[13];
-    x0i = a[5] - a[12];
-    y3r = wn4r * (x0r - x0i);
-    y3i = wn4r * (x0i + x0r);
-    x0r = a[2] - a[11];
-    x0i = a[3] + a[10];
-    y4r = wk1r * x0r - wk1i * x0i;
-    y4i = wk1r * x0i + wk1i * x0r;
-    x0r = a[2] + a[11];
-    x0i = a[3] - a[10];
-    y5r = wk1i * x0r - wk1r * x0i;
-    y5i = wk1i * x0i + wk1r * x0r;
-    x0r = a[6] - a[15];
-    x0i = a[7] + a[14];
-    y6r = wk1i * x0r - wk1r * x0i;
-    y6i = wk1i * x0i + wk1r * x0r;
-    x0r = a[6] + a[15];
-    x0i = a[7] - a[14];
-    y7r = wk1r * x0r - wk1i * x0i;
-    y7i = wk1r * x0i + wk1i * x0r;
-    x0r = y0r + y2r;
-    x0i = y0i + y2i;
-    x1r = y4r + y6r;
-    x1i = y4i + y6i;
-    a[0] = x0r + x1r;
-    a[1] = x0i + x1i;
-    a[2] = x0r - x1r;
-    a[3] = x0i - x1i;
-    x0r = y0r - y2r;
-    x0i = y0i - y2i;
-    x1r = y4r - y6r;
-    x1i = y4i - y6i;
-    a[4] = x0r - x1i;
-    a[5] = x0i + x1r;
-    a[6] = x0r + x1i;
-    a[7] = x0i - x1r;
-    x0r = y1r - y3i;
-    x0i = y1i + y3r;
-    x1r = y5r - y7r;
-    x1i = y5i - y7i;
-    a[8] = x0r + x1r;
-    a[9] = x0i + x1i;
-    a[10] = x0r - x1r;
-    a[11] = x0i - x1i;
-    x0r = y1r + y3i;
-    x0i = y1i - y3r;
-    x1r = y5r + y7r;
-    x1i = y5i + y7i;
-    a[12] = x0r - x1i;
-    a[13] = x0i + x1r;
-    a[14] = x0r + x1i;
-    a[15] = x0i - x1r;
-}
-
-
-void cftf040(double *a)
-{
-    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
-    
-    x0r = a[0] + a[4];
-    x0i = a[1] + a[5];
-    x1r = a[0] - a[4];
-    x1i = a[1] - a[5];
-    x2r = a[2] + a[6];
-    x2i = a[3] + a[7];
-    x3r = a[2] - a[6];
-    x3i = a[3] - a[7];
-    a[0] = x0r + x2r;
-    a[1] = x0i + x2i;
-    a[2] = x1r - x3i;
-    a[3] = x1i + x3r;
-    a[4] = x0r - x2r;
-    a[5] = x0i - x2i;
-    a[6] = x1r + x3i;
-    a[7] = x1i - x3r;
-}
-
-
-void cftb040(double *a)
-{
-    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
-    
-    x0r = a[0] + a[4];
-    x0i = a[1] + a[5];
-    x1r = a[0] - a[4];
-    x1i = a[1] - a[5];
-    x2r = a[2] + a[6];
-    x2i = a[3] + a[7];
-    x3r = a[2] - a[6];
-    x3i = a[3] - a[7];
-    a[0] = x0r + x2r;
-    a[1] = x0i + x2i;
-    a[2] = x1r + x3i;
-    a[3] = x1i - x3r;
-    a[4] = x0r - x2r;
-    a[5] = x0i - x2i;
-    a[6] = x1r - x3i;
-    a[7] = x1i + x3r;
-}
-
-
-void cftx020(double *a)
-{
-    double x0r, x0i;
-    
-    x0r = a[0] - a[2];
-    x0i = a[1] - a[3];
-    a[0] += a[2];
-    a[1] += a[3];
-    a[2] = x0r;
-    a[3] = x0i;
-}
-
-
-void rftfsub(int n, double *a, int nc, double *c)
-{
-    int j, k, kk, ks, m;
-    double wkr, wki, xr, xi, yr, yi;
-    
-    m = n >> 1;
-    ks = 2 * nc / m;
-    kk = 0;
-    for (j = 2; j < m; j += 2) {
-        k = n - j;
-        kk += ks;
-        wkr = 0.5 - c[nc - kk];
-        wki = c[kk];
-        xr = a[j] - a[k];
-        xi = a[j + 1] + a[k + 1];
-        yr = wkr * xr - wki * xi;
-        yi = wkr * xi + wki * xr;
-        a[j] -= yr;
-        a[j + 1] -= yi;
-        a[k] += yr;
-        a[k + 1] -= yi;
-    }
-}
-
-
-void rftbsub(int n, double *a, int nc, double *c)
-{
-    int j, k, kk, ks, m;
-    double wkr, wki, xr, xi, yr, yi;
-    
-    m = n >> 1;
-    ks = 2 * nc / m;
-    kk = 0;
-    for (j = 2; j < m; j += 2) {
-        k = n - j;
-        kk += ks;
-        wkr = 0.5 - c[nc - kk];
-        wki = c[kk];
-        xr = a[j] - a[k];
-        xi = a[j + 1] + a[k + 1];
-        yr = wkr * xr + wki * xi;
-        yi = wkr * xi - wki * xr;
-        a[j] -= yr;
-        a[j + 1] -= yi;
-        a[k] += yr;
-        a[k + 1] -= yi;
-    }
-}
-
-
-void dctsub(int n, double *a, int nc, double *c)
-{
-    int j, k, kk, ks, m;
-    double wkr, wki, xr;
-    
-    m = n >> 1;
-    ks = nc / n;
-    kk = 0;
-    for (j = 1; j < m; j++) {
-        k = n - j;
-        kk += ks;
-        wkr = c[kk] - c[nc - kk];
-        wki = c[kk] + c[nc - kk];
-        xr = wki * a[j] - wkr * a[k];
-        a[j] = wkr * a[j] + wki * a[k];
-        a[k] = xr;
-    }
-    a[m] *= c[0];
-}
-
-
-void dstsub(int n, double *a, int nc, double *c)
-{
-    int j, k, kk, ks, m;
-    double wkr, wki, xr;
-    
-    m = n >> 1;
-    ks = nc / n;
-    kk = 0;
-    for (j = 1; j < m; j++) {
-        k = n - j;
-        kk += ks;
-        wkr = c[kk] - c[nc - kk];
-        wki = c[kk] + c[nc - kk];
-        xr = wki * a[k] - wkr * a[j];
-        a[k] = wkr * a[k] + wki * a[j];
-        a[j] = xr;
-    }
-    a[m] *= c[0];
-}
-
-#if defined(__cplusplus)
-} // end of extern "C"
-#endif