#include <stdio.h>
#include <math.h>

float X [18] = {
    18.000000f,         // 0.
    11.981016f,
     8.965634f,
     7.172353f,
     5.962646f,
     5.109376f,
     4.470957f,
     3.966050f,
     3.569583f,
     3.241481f,
     2.970905f,
     2.740498f,
     2.547215f,
     2.374246f,
     2.224069f,
     2.092965f,
     1.978325f,
     1.873230f,         // 1.
} ;

#define M_PI    3.14159265358979

static const  unsigned char    Parity [256] = {  // parity
    0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,
    1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,
    1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,
    0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,
    1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,
    0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,
    0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,
    1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0
};

static unsigned int  __r1 = 1;
static unsigned int  __r2 = 1;

static void
set_seed ( unsigned int  seed )
{
    __r1 = seed != 0 ? seed : 1;                // 0 is a forbidden value which locks the generator
    __r2 = 1;                                   // the are several (8389119) forbidden codes. 1 is not one of them
}


unsigned int
random_int ( void )
{
    unsigned int  t1, t2, t3, t4;

    t3   = t1 = __r1;   t4   = t2 = __r2;       // Parity calculation is done via table lookup, this is also available
    t1  &= 0xF5;        t2 >>= 25;              // on CPUs without parity, can be implemented in C and avoid unpredictable
    t1   = Parity [t1]; t2  &= 0x63;            // jumps and slow rotate through the carry flag operations.
    t1 <<= 31;          t2   = Parity [t2];

    return (__r1 = (t3 >> 1) | t1 ) ^ (__r2 = (t4 + t4) | t2 );
}



int  usePlevel;
int  MinPNSSubband;

static float  A [36] [36];
static float  B [1801];

static void
Init ( void )
{
        int     i;
        int     j;
        double  w;

        for ( i = 0; i < 36; i += 2 ) {
                for ( j = 0; j < 36; j++ ) {
                        w = 2 * M_PI * (i+1) / 36 * j;
                        A [i+0] [j] = cos (w);
                        A [i+1] [j] = sin (w);
                }
        }
        for ( i = 0; i <= 1800; i++ )
                B [i] = 1.0 ;
        for ( i = 0; i < 18000; i++ ) {
                j     =
                B [j] = 0.01
        }
}


        static double T [4];
        static short k;


static float
Penalties ( float x, int n )                                            // Value between 1 and 18
{                                                                                                       //             =>1    =>0
// Values between 1.0 and 2.0 are noise       -> return 0.0 to 0.1
// Values between 6.0 und 18.0 are not noise  -> return 0.9 to 1.0

//      static double T [4];
        //static short k;


        T [n+0] += x;
        T [n+2] += 1;

        if ( k++ )
                return 0.;
        printf ( "%8.6f%c", T[n]/T[n+2], n ? '\n' : '\t'  );
        return B [ (int) (100. * x) ] ;
}


int
TestForPNSUsage ( int          Subband,
                                  float        SMR,                                     // 1 = 0 dB, 4 = 6 dB
                                  const float  Samples[36] )            // Samples are freed of SCF and decolored upon prediction
{
        float  x;
        float  xx;
        float  tmp;
        int    i;
        int    j;
        float  Result [36];
        float  Ptime;
        float  Pfreq;
        float  Plevel;

        if ( Subband < MinPNSSubband )
                return 0;

        x = xx = 0. ;
        for ( i = 0; i < 36; i += 2 ) {
                tmp = Samples[i] * Samples[i] + Samples[i+1] * Samples[i+1];
                x  += tmp;
                xx += tmp * tmp;
        }
        if ( x == 0. )
                return 0;

        Plevel = fabs (x / 18.) - 1. ;
        Ptime  = Penalties ( 18.*xx / (x*x), 0 );       // possible values between 1 and 18

        return 0;
        for ( i = 0; i < 36; i++ ) {
                x = 0. ;
                for ( j = 0; j < 36; j++ )
                        x += Samples[j] * A[i][j];
                Result [i] = x;
        }

        x = xx = 0. ;
        for ( i = 0; i < 36; i += 2 ) {
                tmp = Result[i] * Result[i] + Result[i+1] * Result[i+1];
                x  += tmp;
                xx += tmp * tmp;
        }

        Pfreq  = Penalties ( 18.*xx / (x*x), 1 );       // possible values between 1 and 18

        if ( usePlevel ) {
                tmp = (1. - Plevel) * Ptime * Pfreq ;
        }
        else {
                tmp = 1. * Ptime * Pfreq ;
        }

        if ( tmp < 1 - 1. / SMR )
                return 0;

        return 1;
}


int
main ( void )
{
        int           i;
        int           j;
        static float  X [36];
        int l;

        Init ();
        for ( l = 2; l <= 36; l+=2 ) {
                printf ("%2u: ", l );
                memset ( T, 0, sizeof(T) );
        for ( i = 0; i < 0x60000; i++ ) {
                for ( j = 0; j < l; j++ ) {
                        X [j] = ((double) (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int () + (int) random_int ()) * 1.e-10 ;
                        //X[j] = j;
                }
                TestForPNSUsage ( 0, 1, X );
        }
        printf ("\n");
        }

        return 0 ;
}