Changeset 109 for mppenc/trunk/src


Ignore:
Timestamp:
11/10/06 22:04:46 (18 years ago)
Author:
zorg
Message:

Remove math aliasing issues which were breaking gcc

Location:
mppenc/trunk/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • mppenc/trunk/src/cvd.c

    r97 r109  
    219219logfast ( float x )
    220220{
    221     double  y = x * x;
    222     y *= y;
    223     y *= y;
    224     return (((int*)(&y))[1] + (45127.5 - 1072693248.)) * ( M_LN2 / (1L<<23) );
     221    union { double d; Int32_t n[2]; } tmp;
     222    tmp.d  = x * x;
     223    tmp.d *= tmp.d;
     224    tmp.d *= tmp.d;
     225#if ENDIAN == HAVE_LITTLE_ENDIAN
     226    return (tmp.n[1] + (45127.5 - 1072693248.)) * ( M_LN2 / (1L<<23) );
     227#else
     228    return (tmp.n[0] + (45127.5 - 1072693248.)) * ( M_LN2 / (1L<<23) );
     229#endif
    225230}
    226231

  • mppenc/trunk/src/fastmath.c

    r97 r109  
    1919
    2020#include "mppenc.h"
     21#include "fastmath.h"
    2122
    2223#ifdef FAST_MATH
    … …  
    3031void   Init_FastMath ( void )
    3132{
    32     int     i;
    33     float   X;
    34     float   Y;
    35     double  xm;
    36     double  x0;
    37     double  xp;
    38     double  x;
    39     double  y;
    40     float*  p;
     33    int i; mpc_floatint X, Y; double xm, x0, xp, x, y; float* p;
    4134
    4235    p = (float*) tabatan2;
    … …  
    6457    p = (float*) tabsqrt_ex;
    6558    for ( i = 0; i < 255; i++ ) {
    66         *(int*)&X = (i << 23);
    67         *(int*)&Y = (i << 23) + (1<<23) - 1;
    68         *p++ = sqrt(X);
     59        X.n = (i << 23);
     60        Y.n = (i << 23) + (1<<23) - 1;
     61        *p++ = sqrt(X.f);
    6962    }
    70     *(int*)&X = (255 << 23) - 1;
    71     *p++ = sqrt(X);
     63    X.n = (255 << 23) - 1;
     64    *p++ = sqrt(X.f);
    7265
    7366    p = (float*) tabsqrt_m;

  • mppenc/trunk/src/fastmath.h

    r97 r109  
    1818 */
    1919
    20 #if 1
    21 # define ROUND32(x)   ( floattmp = (x) + (int)0x00FD8000L, *(int*)(&floattmp) - (int)0x4B7D8000L )
    22 #else
    23 # define ROUND32(x)   ( (int) floor ((x) + 0.5) )
    24 #endif
     20typedef union mpc_floatint
     21{
     22    float   f;
     23    Int32_t n;
     24} mpc_floatint;
     25
     26static __inline Int32_t mpc_nearbyintf(float fVal)
     27{
     28    mpc_floatint tmp;
     29    tmp.f = fVal  + 0x00FF8000;
     30    tmp.n = tmp.n - 0x4B7F8000;
     31    return tmp.n;
     32}
     33#define mpc_lrintf mpc_nearbyintf
     34static __inline Int32_t mpc_round32(float fVal)
     35{
     36    mpc_floatint tmp;
     37    tmp.f = fVal  + 0x00FD8000;
     38    tmp.n = tmp.n - 0x4B7D8000;
     39    return tmp.n;
     40}
    2541
    2642#ifdef FAST_MATH
    … …  
    2945my_atan2 ( float x, float y )
    3046{
    31     float  t;
    32     int    i;
    33     float  ret;
    34     float  floattmp;
     47    float t, ret; int i; mpc_floatint mx, my;
    3548
    36     if ( (*(int*)&x & 0x7FFFFFFF) < (*(int*)&y & 0x7FFFFFFF) ) {
    37         i   = ROUND32 (t = TABSTEP * (x / y));
     49    mx.f = x;
     50    my.f = y;
     51    if ( (mx.n & 0x7FFFFFFF) < (my.n & 0x7FFFFFFF) ) {
     52        i   = mpc_round32 (t = TABSTEP * (mx.f / my.f));
    3853        ret = tabatan2 [1*TABSTEP+i][0] + tabatan2 [1*TABSTEP+i][1] * (t-i);
    39         if ( *(int*)&y < 0 )
     54        if ( my.n < 0 )
    4055           ret = (float)(ret - M_PI);
    4156    }
    42     else if ( *(int*)&x < 0) {
    43         i   = ROUND32 (t = TABSTEP * (y / x));
     57    else if ( mx.n < 0 ) {
     58        i   = mpc_round32 (t = TABSTEP * (my.f / mx.f));
    4459        ret = - M_PI/2 - tabatan2 [1*TABSTEP+i][0] + tabatan2 [1*TABSTEP+i][1] * (i-t);
    4560    }
    46     else if ( *(int*)&x > 0) {
    47         i   = ROUND32 (t = TABSTEP * (y / x));
     61    else if ( mx.n > 0 ) {
     62        i   = mpc_round32 (t = TABSTEP * (my.f / mx.f));
    4863        ret = + M_PI/2 - tabatan2 [1*TABSTEP+i][0] + tabatan2 [1*TABSTEP+i][1] * (i-t);
    4964    }
    … …  
    5873my_cos ( float x )
    5974{
    60     float  t;
    61     int    i;
    62     float  ret;
    63     float  floattmp;
    64 
    65     i   = ROUND32 (t = TABSTEP * x);
     75    float t, ret; int i;
     76    i   = mpc_round32 (t = TABSTEP * x);
    6677    ret = tabcos [13*TABSTEP+i][0] + tabcos [13*TABSTEP+i][1] * (t-i);
    6778    return ret;
    … …  
    7283my_ifloor ( float x )
    7384{
    74     x = x + (0x0C00000L + 0.500000001);
    75     return *(int*)&x - 1262485505;
     85    mpc_floatint mx;
     86    mx.f = (float) (x + (0x0C00000L + 0.500000001));
     87    return mx.n - 1262485505;
    7688}
    7789
    … …  
    8092my_sqrt ( float x )
    8193{
    82     float  ret;
    83     int    i;
    84     int    ex = *(int*)&x >> 23;                                // get the exponent
    85     float  floattmp;
    86 
    87     *(int*)&x = (*(int*)&x & 0x7FFFFF) | 0x42800000;            // delete the exponent
    88     i    = ROUND32 (x);                                         // Integer-part of the mantissa  (round ????????????)
    89     ret  = tabsqrt_m [i-TABSTEP][0] + tabsqrt_m [i-TABSTEP][1] * (x-i); // calculate value
     94    float  ret; int i, ex; mpc_floatint mx;
     95    mx.f = x;
     96    ex   = mx.n >> 23;                     // get the exponent
     97    mx.n = (mx.n & 0x7FFFFF) | 0x42800000; // delete the exponent
     98    i    = mpc_round32 (mx.f);             // Integer-part of the mantissa  (round ????????????)
     99    ret  = tabsqrt_m [i-TABSTEP][0] + tabsqrt_m [i-TABSTEP][1] * (mx.f-i); // calculate value
    90100    ret *= tabsqrt_ex [ex];
    91101    return ret;

  • mppenc/trunk/src/quant.c

    r97 r109  
    1919
    2020#include "mppenc.h"
     21#include "fastmath.h"
    2122
    2223/* V A R I A B L E S */
    … …  
    149150ISNR_Schaetzer ( const float* input, const float SNRcomp, const int res )
    150151{
    151     int    k;
    152     float  fac    = A [res];
    153     float  invfac = C [res];
    154     float  Signal = 1.e-30f;
    155     float  Fehler = 1.e-30f;
    156     float  tmp ;
    157     float  tmp2;
    158     float  tmp3;
     152    int k; float signal, error, sig, err;
     153    const float fac       = A [res];
     154    const float invfac    = C [res];
     155    const float noiseComp = NoiseInjectionCompensation1D [res];
    159156
    160157    // Summation of the absolute power and the quadratic error
    161     for ( k = 0; k < 36; k++ ) {
    162         tmp2    = input[k] * NoiseInjectionCompensation1D [res];
    163         // q = ftol(in), correct rounding
    164         tmp  = tmp2 * fac + 0xFF8000;
    165         tmp3 = (*(int*) & tmp - 0x4B7F8000) * invfac;
    166         tmp  = tmp3 - tmp2;
    167 
    168         Fehler += tmp * tmp;
    169         Signal += tmp2 * tmp2;
     158    signal = error = 1.e-30f;
     159    for ( k = 0; k < 36; k++ )
     160    {
     161        sig     = input[k] * noiseComp;
     162        err     = mpc_nearbyintf(sig * fac) * invfac - sig;
     163        signal += sig * sig;
     164        error  += err * err;
    170165    }
    171166
    172167    // Utilization of SNRcomp only if SNR > 1 !!!
    173     return Signal > Fehler  ?  Fehler / (SNRcomp * Signal)  :  Fehler / Signal;
     168    return signal > error ? error / (signal * SNRcomp) : error / signal;
    174169}
    175170
    … …  
    178173ISNR_Schaetzer_Trans ( const float* input, const float SNRcomp, const int res )
    179174{
    180     int    k;
    181     float  fac    = A [res];
    182     float  invfac = C [res];
    183     float  Signal;
    184     float  Fehler;
    185     float  ret ;
    186     float  tmp ;
    187     float  tmp2;
    188     float  tmp3;
    189 
     175    int k; float  signal, error, result, sig, err;
     176    const float fac       = A [res];
     177    const float invfac    = C [res];
     178    const float noiseComp = NoiseInjectionCompensation1D [res];
     179   
    190180    // Summation of the absolute power and the quadratic error
    191     k = 0;
    192     Signal = Fehler = 1.e-30f;
    193     for ( ; k < 12; k++ ) {
    194         tmp2    = input[k] * NoiseInjectionCompensation1D [res];
    195         // q = ftol(in), correct rounding
    196         tmp  = tmp2 * fac + 0xFF8000;
    197         tmp3 = (*(int*) & tmp - 0x4B7F8000) * invfac;
    198         tmp  = tmp3 - tmp2;
    199 
    200         Fehler += tmp * tmp;
    201         Signal += tmp2 * tmp2;
    202     }
    203     tmp = Signal > Fehler  ?  Fehler / (SNRcomp * Signal)  :  Fehler / Signal;
    204     ret = tmp;
    205     Signal = Fehler = 1.e-30f;
    206     for ( ; k < 24; k++ ) {
    207         tmp2    = input[k] * NoiseInjectionCompensation1D [res];
    208         // q = ftol(in), correct rounding
    209         tmp  = tmp2 * fac + 0xFF8000;
    210         tmp3 = (*(int*) & tmp - 0x4B7F8000) * invfac;
    211         tmp  = tmp3 - tmp2;
    212 
    213         Fehler += tmp * tmp;
    214         Signal += tmp2 * tmp2;
    215     }
    216     tmp = Signal > Fehler  ?  Fehler / (SNRcomp * Signal)  :  Fehler / Signal;
    217     if ( tmp > ret ) ret = tmp;
    218     //ret += tmp;
    219     Signal = Fehler = 1.e-30f;
    220     for ( ; k < 36; k++ ) {
    221         tmp2    = input[k] * NoiseInjectionCompensation1D [res];
    222         // q = ftol(in), correct rounding
    223         tmp  = tmp2 * fac + 0xFF8000;
    224         tmp3 = (*(int*) & tmp - 0x4B7F8000) * invfac;
    225         tmp  = tmp3 - tmp2;
    226 
    227         Fehler += tmp * tmp;
    228         Signal += tmp2 * tmp2;
    229     }
    230     tmp = Signal > Fehler  ?  Fehler / (SNRcomp * Signal)  :  Fehler / Signal;
    231     if ( tmp > ret ) ret = tmp;
    232     //ret += tmp;
    233     //ret *= 0.33333333333f;
    234 
    235     return ret;
     181    signal = error = 1.e-30f;
     182    for ( k = 0 ; k < 12; k++ )
     183    {
     184        sig     = input[k] * noiseComp;
     185        err     = mpc_nearbyintf(sig * fac) * invfac - sig;
     186        signal += sig * sig;
     187        error  += err * err;
     188    }
     189    result = signal > error ? error / (signal * SNRcomp) : error / signal;
     190
     191    signal = error = 1.e-30f;
     192    for ( ; k < 24; k++ )
     193    {
     194        sig     = input[k] * noiseComp;
     195        err     = mpc_nearbyintf(sig * fac) * invfac - sig;
     196        signal += sig * sig;
     197        error  += err * err;
     198    }
     199    sig = signal > error ? error / (signal * SNRcomp) : error / signal;
     200    if ( sig > result ) result = sig;
     201
     202    signal = error = 1.e-30f;
     203    for ( ; k < 36; k++ )
     204    {
     205        sig     = input[k] * noiseComp;
     206        err     = mpc_nearbyintf(sig * fac) * invfac - sig;
     207        signal += sig * sig;
     208        error  += err * err;
     209    }
     210    sig = signal > error ? error / (signal * SNRcomp) : error / signal;
     211    if ( sig > result ) result = sig;
     212
     213    return result;
    236214}
    237215
    … …  
    241219QuantizeSubband ( unsigned int* qu_output, const float* input, const int res, float* errors )
    242220{
    243     int    n;
    244     int    offset  = D [res];
    245     float  mult    = A [res] * NoiseInjectionCompensation1D [res];
    246     float  invmult = C [res];
    247     float  tmp;
    248     int    quant;
    249     float  signal;
    250 
    251     for ( n = 0; n < 36 - MAX_NS_ORDER; n++, input++, qu_output++ ) {
    252         // q = ftol(in), correct rounding
    253         tmp   = *input * mult + 0xFF8000;
    254         quant = (unsigned int)(*(int*) & tmp - 0x4B7F8000 + offset);
    255 
     221    int n, quant; float signal;
     222    const int   offset    = D [res];
     223    const float noiseComp = NoiseInjectionCompensation1D [res];
     224    const float mult      = A [res] * noiseComp;
     225    const float invmult   = C [res];
     226
     227    for ( n = 0; n < 36 - MAX_NS_ORDER; n++)
     228    {
     229        quant = (unsigned int) ( mpc_lrintf( input[n] * mult ) + offset );
    256230        // limitation to 0...2D
    257         if ((unsigned int)quant > (unsigned int)2*offset ) {
    258             quant = mini ( quant, 2*offset );
    259             quant = maxi ( quant,        0 );
     231        if ( quant > offset * 2 )
     232        {
     233            quant = mini ( quant, offset * 2 );
     234            quant = maxi ( quant, 0 );
    260235        }
    261         *qu_output  = quant;
    262     }
    263 
    264     for ( ; n < 36; n++, input++, qu_output++ ) {
    265         // q = ftol(in), correct rounding
    266         signal = *input * mult;
    267         tmp   =  signal + 0xFF8000;
    268         quant = (unsigned int)(*(int*) & tmp - 0x4B7F8000 + offset);
     236        qu_output[n] = quant;
     237    }
     238
     239    for ( ; n < 36; n++ )
     240    {
     241        signal = input[n] * mult;
     242        quant  = (unsigned int) ( mpc_lrintf( signal ) + offset );
    269243
    270244        // calculate the current error and save it for error refeeding
    271         errors [n + 6] = invmult * (quant - offset) - signal * NoiseInjectionCompensation1D [res];
     245        errors [n + 6] = invmult * (quant - offset) - signal * noiseComp;
    272246
    273247        // limitation to 0...2D
    274         if ((unsigned int)quant > (unsigned int)2*offset ) {
    275             quant = mini ( quant, 2*offset );
    276             quant = maxi ( quant,        0 );
     248        if ( quant > offset * 2 )
     249        {
     250            quant = mini ( quant, offset * 2);
     251            quant = maxi ( quant, 0 );
    277252        }
    278         *qu_output = quant;
     253        qu_output[n] = quant;
    279254    }
    280255}
    … …  
    285260QuantizeSubbandWithNoiseShaping ( unsigned int* qu_output, const float* input, const int res, float* errors, const float* FIR )
    286261{
    287 #define E(x) *((int*)errors+(x))
    288 
    289     float  signal;
    290     float  tmp;
    291     float  mult    = A [res];
    292     float  invmult = C [res];
    293     int    offset  = D [res];
    294     int    n;
    295     int    quant;
    296 
    297     E(0) = E(1) = E(2) = E(3) = E(4) = E(5) = 0;       // arghh, it produces pops on each frame boundary!
    298 
    299     for ( n = 0; n < 36; n++, input++, qu_output++ ) {
    300         signal = *input * NoiseInjectionCompensation1D [res] - (FIR[5]*errors[n+0] + FIR[4]*errors[n+1] + FIR[3]*errors[n+2] + FIR[2]*errors[n+3] + FIR[1]*errors[n+4] + FIR[0]*errors[n+5]);
    301 
    302         // quant = ftol(signal), correct rounding
    303         tmp   = signal * mult + 0xFF8000;
    304         quant = *(int*) & tmp - 0x4B7F8000;
     262    int n, quant; float signal;
     263    const float mult      = A [res];
     264    const float invmult   = C [res];
     265    const int   offset    = D [res];
     266    const float noiseComp = NoiseInjectionCompensation1D [res];
     267
     268    memset(errors, 0, 6 * sizeof *errors); // arghh, it produces pops on each frame boundary!
     269
     270    for ( n = 0; n < 36; n++ )
     271    {
     272        signal = input[n] * noiseComp -
     273            (FIR[5]*errors[n+0] + FIR[4]*errors[n+1] + FIR[3]*errors[n+2] +
     274             FIR[2]*errors[n+3] + FIR[1]*errors[n+4] + FIR[0]*errors[n+5]);
     275        quant = mpc_lrintf(signal * mult);
    305276
    306277        // calculate the current error and save it for error refeeding
    307         errors [n + 6] = invmult * quant - signal * NoiseInjectionCompensation1D [res];
     278        errors [n + 6] = invmult * quant - signal * noiseComp;
    308279
    309280        // limitation to +/-D
    … …  
    311282        quant = maxf ( quant, -offset );
    312283
    313         *qu_output = (unsigned int)(quant + offset);
     284        qu_output[n] = (unsigned int) ( quant + offset );
    314285    }
    315286}
Note: See TracChangeset for help on using the changeset viewer.