if(i && !(i % 4))
		fprintf(f, "\n");
	    fprintf(f, "%.12g, ", nbrs[i + j * len2]);
	}
	fprintf(f, "\n},\n");
    }
    fprintf(f, "\n};\n\n");

    fclose(f);
}

#endif


#ifdef USE_DATA
#include "pow2table1.h"
#include "pow2table2.h"
#include "powerscales.h"
#include "imdct_win.h"
#include "globgaintable.h"
#else
static mpfloat Granule_pow2table_m2[32]; /* (2^(-2))^i */
static mpfloat Granule_pow2table_m05[32]; /* (2^(-0.5))^i */
static mpfloat Granule_powerscaling[200]; /* i^(4/3) */
mpfloat Granule_imdct_win[4][36];
static mpfloat Granule_globgaintable[256];
#endif

static float Granule_alias_c[8] = { -0.6f, -0.535f, -0.33f, -0.185f,
				    -0.095f, -0.041f, -0.0142f, -0.0037f };

mpfloat Granule_alias_cs[8], Granule_alias_ca[8];

#ifdef LEE_IMDCT
#ifdef USE_DATA
#include "idct_9x9.h"
#else
mpfloat Granule_9x9_idct[72];
#endif
#else
static mpfloat Granule_imdct_bigCOS[18*18];
static mpfloat Granule_imdct_bigCOS2[12*6];
#endif

mpfloat Granule_imdct_previous[2][576];  /* used for overlapping */

int Granule_sbsynth_Vptr[2] = { 64, 64 };
#include "mp3dec_D.h"

void
Granule_init()
{
    int i, j, k, m, p, x1, x2;
	int odd_i, two_odd_i, four_odd_i, eight_odd_i;
    float f;
    static int inited = 0;

    if(inited)
	return;

    inited = 1;

    /* initialize some powertables used in Granule_requantize */
    /* ## alternative is to make a .h file with these calculated
       directly in the source */
#ifndef USE_DATA
    for(i = 0; i < 32; i++) {
	Granule_pow2table_m2[i] = (mpfloat)pow(2.0, -2.0 * i);
	Granule_pow2table_m05[i] = (mpfloat)pow(2.0, -0.5 * i);
    }
#endif

#ifdef MAKE_DATA
    make_data_file("pow2table1.h", "Granule_pow2table_m2", 
		   Granule_pow2table_m2, 32);
    make_data_file("pow2table2.h", "Granule_pow2table_m05", 
		   Granule_pow2table_m05, 32);
#endif

#ifndef USE_DATA
    for(i = 0; i < 200; i++)
	Granule_powerscaling[i] = (mpfloat)pow( i, (4.0 / 3.0) );
#endif

#ifdef MAKE_DATA
    make_data_file("powerscales.h", "Granule_powerscaling", 
		   Granule_powerscaling, 200);
#endif

    /* calculate the anti-aliasing butterfly coefficients */

    for(i = 0; i < 8; i++) {
	f = sqrt(1.0 + Granule_alias_c[i] * Granule_alias_c[i]);
	Granule_alias_cs[i] = (mpfloat)(1.0f / f);
	Granule_alias_ca[i] = (mpfloat)(Granule_alias_c[i] / f);
    }

#ifndef USE_DATA
    /* calculate some window shapes and cos tables for the IMDCT */

    /* block_type 0 (normal window) */

    for(i = 0; i < 36; i++)
	Granule_imdct_win[0][i] = (mpfloat)sin(PI/36 * (i + 0.5));
    
    /* block_type 1 (start block) */

    for(i = 0; i < 18; i++)
	Granule_imdct_win[1][i] = (mpfloat)sin(PI/36 * (i + 0.5));

    for(i = 18; i < 24; i++)
	Granule_imdct_win[1][i] = (mpfloat)1.0f;

    for(i = 24; i < 30; i++)
	Granule_imdct_win[1][i] = (mpfloat)sin(PI/12 * (i - 18 + 0.5));

    for(i = 30; i < 36; i++)
	Granule_imdct_win[1][i] = (mpfloat)0.0f;
    
    /* block_type 3 (stop block) */

    for(i = 0; i < 6; i++)
	Granule_imdct_win[3][i] = (mpfloat)0.0f;

    for(i = 6; i < 12; i++)
	Granule_imdct_win[3][i] = (mpfloat)sin(PI/12 * (i - 6 + 0.5));

    for(i = 12; i < 18; i++)
	Granule_imdct_win[3][i] = (mpfloat)1.0f;

    for(i = 18; i < 36; i++)
	Granule_imdct_win[3][i] = (mpfloat)sin(PI/36 * (i + 0.5));
    
    /* block_type 2 (short block) */

    for(i = 0; i < 12; i++)
	Granule_imdct_win[2][i] = (mpfloat)sin(PI/12 * (i + 0.5));

    for(i = 12; i < 36; i++)
	Granule_imdct_win[2][i] = (mpfloat)0.0 ;

#endif

#ifdef MAKE_DATA
    make_data_file_2d("imdct_win.h", "Granule_imdct_win", 
		   &Granule_imdct_win[0][0], 4, 36);
#endif
          
#ifdef LEE_IMDCT
#ifndef USE_DATA
    j = 0;
    for(i = 0; i < 9; i++) {
	odd_i = (i << 1) + 1;
	two_odd_i = odd_i << 1;
	four_odd_i = odd_i << 2;
	Granule_9x9_idct[j++] = (mpfloat)cos(PI/18 * odd_i);
	Granule_9x9_idct[j++] = (mpfloat)cos(PI/18 * two_odd_i);
	eight_odd_i = two_odd_i << 2;
	Granule_9x9_idct[j++] = (mpfloat)cos(PI/18 * (four_odd_i - odd_i));
	Granule_9x9_idct[j++] = (mpfloat)cos(PI/18 * four_odd_i);
	Granule_9x9_idct[j++] = (mpfloat)cos(PI/18 * (four_odd_i + odd_i));
	Granule_9x9_idct[j++] = (mpfloat)cos(PI/18 * (four_odd_i + two_odd_i));
	Granule_9x9_idct[j++] = (mpfloat)cos(PI/18 * (eight_odd_i - odd_i));
	Granule_9x9_idct[j++] = (mpfloat)cos(PI/18 * eight_odd_i);
    }
#endif
#ifdef MAKE_DATA
    make_data_file("idct_9x9.h", "Granule_9x9_idct", 
		   Granule_9x9_idct, 72);
#endif

#else
    /* make the costable for the short blocks */
    for(i = 0, p = 0; i < 6; i++)
	for(j = 0; j < 6; j++)
	    Granule_imdct_bigCOS2[p++] = (mpfloat)-cos((2*(i-3)+7)*
						       (2*j+1)*PI/24);

    /* make a linear cos factor table, so the IMDCT can be calculated
       from it without messing with any modulos, weird index additions
       etc */

    for(i = 0, p = 0; i < 18; i++) {
	for(j = 0; j < 18; j++)
	    Granule_imdct_bigCOS[p++] = (mpfloat)-cos((2*(i-9)+19)*
						      (2*j+1)*PI/72); 
    }
#endif

    /* clear the previous data table used for overlapping in IMDCT */
    
    for(p = 0; p < 2; p++)
	for(i = 0; i < 576; i++)
	    Granule_imdct_previous[p][i] = (mpfloat)0.0f;

#ifndef USE_DATA
    fast_idct_init();
#endif

#ifndef INT_MATH
    /* rescale the D table to include the 32768 factor */
    for(i = 0; i < 512; i++)
#ifdef DSP_LOFI
	Granule_sbsynth_D[i] *= 8192;
#else
	Granule_sbsynth_D[i] *= 32768;
#endif
#endif
#if 0
    for(i = 0; i < 512; i++)
	printf("D[%d] = 0x%8x, ", i, Granule_sbsynth_D[i].rint());
#endif

    windowing_init();

#ifndef USE_DATA
    for(i = 0; i < 256; i++)
	Granule_globgaintable[i] = (mpfloat)pow(2.0 , 
						(0.25 * (i -
							 210.0)));
#endif
#ifdef MAKE_DATA
    make_data_file("globgaintable.h", "Granule_globgaintable", 
		   Granule_globgaintable, 256);
#endif

}

void
Granule_decode_info(Granule *gr, Bitstream *bs)
{
    int region, window;

    gr->part2_3_length = Bitstream32_get(bs, 12);
    gr->big_values = Bitstream32_get(bs, 9);
    gr->global_gain = Bitstream32_get(bs, 8);
    gr->scalefac_compress = Bitstream32_get(bs, 4);
      
    if((gr->window_switching_flag = Bitstream32_get1(bs))) {
	
	/* the block_type indicates the window type for the granule */
	
	gr->block_type = Bitstream32_get(bs, 2);

	gr->mixed_block_flag = Bitstream32_get1(bs);

	/* table selectors select different Huffman tables for each
	   frequency region */

	for(region = 0; region < 2; region++)
	    gr->table_select[region] = Bitstream32_get(bs, 5);

	/* subblock_gain multiplied by 4 indicates the gain offset for
	   each window, compared to the global_gain of the granule */
	
	for(window = 0; window < 3; window++)
	  gr->subblock_gain[window] = Bitstream32_get(bs, 3);

	/* region0count and region1count are set to default values */

	if(gr->block_type == BLOCKTYPE_3WIN &&
	   !gr->mixed_block_flag)
	    gr->region0_count = 8;
	else
	    gr->region0_count = 7;
	
	gr->region1_count = 20 - gr->region0_count;

    } else {
	
	gr->block_type = 0;
	
	for(region = 0; region < 3; region++)
	    gr->table_select[region] = Bitstream32_get(bs, 5);
	
	gr->region0_count = Bitstream32_get(bs, 4);
	gr->region1_count = Bitstream32_get(bs, 3);
    }

    gr->preflag = Bitstream32_get1(bs);
    gr->scalefac_scale = Bitstream32_get1(bs);
    gr->count1table_select = Bitstream32_get1(bs);   
}

void
Granule_decode_scalefactors(Granule *gr, Bitstream *bs, 
			    int channel, int grnum, Frame *f)
{
    int slen1, slen2, i, sfb, window;

    slen1 = Frame_slen[0][gr->scalefac_compress];
    slen2 = Frame_slen[1][gr->scalefac_compress];

    /* read the scalefactors - these are read in two different ways
       depending on the block type */

    if(gr->window_switching_flag &&
       gr->block_type == BLOCKTYPE_3WIN) {

	if(gr->mixed_block_flag) {
	    
	    /* mixed blocks - a mixed block uses
	       both long scalefactors (0-7) and short scalefactors (3-11)
	       which together makes up the whole spectrum 
	       */
	    
	    /* block_type 2 and mixed_block_flag 1, 
	       slen1: length of scalefactors for bands 0 to 7 of the
	       long sf and length of sf's 3 to 5
	       of the short scalefactor band 
	       slen2: length of scalefactors for bands 6 to 11 of the
	       short scalefactor band
	       */
	    
	    for(sfb = 0; sfb < 8; sfb++) 
		gr->scalefac_l[sfb] =
		    Bitstream_get(bs, slen1);
	    
	    for(sfb = 3; sfb < 6; sfb++)
		for(window = 0; window < 3; window++)
		    gr->scalefac_s[sfb][window] =
			Bitstream_get(bs, slen1);
	    
	    for(sfb = 6; sfb < 12; sfb++)
		for(window = 0; window < 3; window++)
		    gr->scalefac_s[sfb][window] =
			Bitstream_get(bs, slen2);
	    
	    for(window = 0; window < 3; window++)
		    gr->scalefac_s[sfb][window] = 0;
		
	} else {
	    
	    /* short block using only short scalefactors
	       slen1: length of scalefactors for short sf bands 0 to 5
	       slen2: length of scalefactors for short sf bands 6 to 11
	       */
	    
	    for(sfb = 0; sfb < 6; sfb++)
		for(window = 0; window < 3; window++)
		    gr->scalefac_s[sfb][window] =
			Bitstream_get(bs, slen1);
	    
	    for(sfb = 6; sfb < 12; sfb++)
		for(window = 0; window < 3; window++)
		    gr->scalefac_s[sfb][window] =
			Bitstream_get(bs, slen2);

	    for(window = 0; window < 3; window++)
		    gr->scalefac_s[sfb][window] = 0;
	    
	}  
	
    } else {
	
	/* long block
	   slen1: length of scalefactors for long sf bands 0 to 10
	   slen2: length of scalefactors for long sf bands 11 to 20
	   */
	
	/* the sf bands are divided in 4 parts, 0-5, 6-10, 11-15, 16-20
	   the standard seems to disagree with itself here, but apparently
	   it works this way:
	   */
	
	static int sfb_bound[5] = { 0, 6, 11, 16, 21 };
	
	for(i = 0; i < 4; i++)
	    if(f->scfsi[channel][i] == 0 || grnum == 0)
		for(sfb = sfb_bound[i]; sfb < sfb_bound[i + 1]; sfb++) 
		    gr->scalefac_l[sfb] =
			Bitstream_get(bs, sfb < 11 ? slen1 : slen2);
	
	gr->scalefac_l[21] = 0;
	gr->scalefac_l[22] = 0;

    }    
}


/* The so-called Requantizer takes the de-huffmanized energies of each
   frequency and applies the scalefactors to them to rescale them to 
   their former glory. 
 */

void
Granule_requantize(Granule *gr, Granule_intfreqs is,
		   Granule_floatfreqs xr, Frame *f)
{
    /* 
       We have to go through all frequency lines in all bands and apply
       a large formula to the is[freq] value. All bands have their own
       scalefactors that are incorporated in this formula, so we have to
       extract these per band before requantizing a band.
     
       The formula for short blocks is:

       xr[i] = sign(is[i]) * |is[i]|^(4/3) * 2 ^ 0.25 * (global_gain - 210 -
               8 * subblock_gain[window]) * 2 ^ -(scalefac_multiplier *
	       scalefac_s[band][window])

       and for long blocks:

       xr[i] = sign(is[i]) * | is[i] |^(4/3) * 2 ^ 0.25 * (global_gain - 210) *
               2 ^ -(scalefac_multiplier * (scalefac_l[band] + preflag *
	       pretab[band]))

       The preflag when set applies a predefined pre-emphasis value from
       the pretab table. 

       The 2 ^ operation is very costly, but fortunately, the range of the
       possible exponents is fairly small so we can do a table lookup and
       avoid all 2 ^ operations.

    */

    static int pretab[22] = { 0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,2,2,3,3,3,2,0 };
    int i, band, bandend, window, windowlength = 0, j = 0, t, sign, needc;
    mpfloat myxr = (mpfloat)0.0f;

    /* get the correct band-boundaries table */

    struct Frame_band *bands = &Frame_bands[f->sampling_frequency];

    /* get the global scaling factor - ## see if this can be table-ized */
/*     globscale = (mpfloat)pow(2.0 , 
		       	      (0.25 * (gr->global_gain -
			      		       210.0)));
					       */

    mpfloat globscale = Granule_globgaintable[gr->global_gain];
	
    /* In mixed blocks, the lowest two polyphase subbands are transformed
       using normal block bands, and the rest according to block_type.
       
       Remember that band boundaries in short blocks should be multiplied
       by 3 (since there are 3 windows in each band).

    */

    if(gr->window_switching_flag && gr->block_type == BLOCKTYPE_3WIN &&
       !gr->mixed_block_flag) {
	bandend = bands->s[1] * 3;  /* short block */
	j = windowlength = bands->s[1];
    } else
	bandend = bands->l[1];      /* normal block or mixed short block */

    band = window = 0;
    needc = 1;
    
    for(i = 0; i < 576; i++) {
	
	if(i == bandend) {

	    needc = 1;
	    /* advance to the next band */
	    
	    if(gr->window_switching_flag && gr->block_type == BLOCKTYPE_3WIN) {

		window = 0;   /* reset window number */

		if(gr->mixed_block_flag) {

		    /* a short, mixed block is coded as a long block for
		     the first two polyphase subbands, so we have 3 
		     situations here - either we are in the long part, we are
		     on the boundary, or we are in the short part 

		     bands->l[8] == bands->s[3] == 36 for all frequencies

		     */

		    if(i == bands->l[8]) {
			/* on the boundary */
			bandend = bands->s[4] * 3;
			band = 3;
			windowlength = 4;  /* s[4]-s[3] is same for all freq */
		    } else {
			if(i < bands->l[8]) 
			    /* in the long part */
			    bandend = bands->l[++band + 1];
			else {
			    /* in the short part */
			    bandend = bands->s[++band + 1] * 3;
			    windowlength = bands->s[band + 1] - bands->s[band];
			}
		    }
		} else {
		    /* pure short block */
		    bandend = bands->s[++band + 1] * 3;
		    windowlength = bands->s[band + 1] - bands->s[band];
		}

		j = windowlength;

		/* For short blocks, i.e blocks with 3 windows in each band,
		   the windowlength gives how many values are in each band
		   in each window. We put this value in j, and decrease j
		   for each value. When j hits 0, we restore j and go to the
		   next window. 
		 */

	    } else
		/* long block */
		bandend = bands->l[++band + 1];
	}
	
	if(i == 36)
	    needc = 1;

/* only recalculate the xr formula when we pass between bands and
	windows
*/
	if(needc) {
	    needc = 0;
	    
	    /* start with the global part (same for all values in the block) */

	    myxr = globscale;

	    /* apply the formula according to if its a short or long block */
	    /* remember that mixed blocks are long in the lowest 2 subbands */

	    if (gr->window_switching_flag && gr->block_type == BLOCKTYPE_3WIN &&
		(!gr->mixed_block_flag || (gr->mixed_block_flag &&
					   i >= 2 * NUM_DCTBANDS))) {

		/* its a short block, use this formula:

		   xr[i] = sign(is[i])*|is[i]|^(4/3)*2^0.25*(global_gain - 210 -
		   8 * subblock_gain[window]) * 2^-(scalefac_multiplier *
		   scalefac_s[band][window])

		   the global gain part is already in *xr
	
		   we are clever enough to use tables for the powers-of-2

		   scalefac_multiplier is 0.5 for scalefac_scale == 0 and 1
		   for scalefac_scale == 1
	       
		   */
		
		t = gr->subblock_gain[window];

		if(t >= 32) { /* when dividing by 2^32, you don't have much left */
		    myxr = (mpfloat)0.0f;
		}
		else
		    myxr *= Granule_pow2table_m2[t];
		ISCALE(myxr);
   		
		t = (1 + gr->scalefac_scale) * gr->scalefac_s[band][window];
	    
		if(t >= 32) { /* Only 32 values in the table */
		    myxr = (mpfloat)0.0f;
		}
		else
		    myxr *= Granule_pow2table_m05[t];
		ISCALE(myxr);