Commit f47fd3d8 authored by thomas dettbarn's avatar thomas dettbarn
Browse files

More deterministic "magic_shifts" in biDiagonalReductionLeft_64() and...

More deterministic "magic_shifts" in biDiagonalReductionLeft_64() and biDiagonalReductionRight_64().
parent f5737f6c
Loading
Loading
Loading
Loading
Loading
+19 −18
Original line number Diff line number Diff line
@@ -1004,8 +1004,9 @@ static void biDiagonalReductionLeft_64(
{

#define HEADROOM_LEFT_1 1 
#define HEADROOM_LEFT_2 ( HEADROOM_LEFT_1 + 1 )
#define HEADROOM_LEFT_2 2
#define HEADROOM_LEFT_3 2

    Word16 iCh, jCh;
    Word32 norm_x;
    Word16 norm_x_e;
@@ -1048,28 +1049,28 @@ static void biDiagonalReductionLeft_64(
            ( *g ) = L_negate( *g );
        }
        factor2 = W_extract_l( W_shr( singularVectors_Left_64[currChannel][currChannel], HEADROOM_LEFT_1 ) );
        tmp_e = shr( sub( norm_x_e, 1 ), 1 );
        tmp_e = sub( HEADROOM_LEFT_1, ( *g_e ) );
        tmpmul = W_mult0_32_32( ( *g ), factor2 );
        tmpmul = W_shr( tmpmul, tmp_e );
        r_64 = W_sub( tmpmul, norm_64 );
        r_e = W_norm( r_64 );
        r = W_extract_h( W_shl( r_64, r_e ) );


        invVal_e = sub( add( 1, ( HEADROOM_LEFT_1 + HEADROOM_LEFT_1 ) ), r_e );
        invVal = BASOP_Util_Divide3232_Scale_newton( MAXVAL_WORD32, r, &invVal_e );
        invVal_e = 0;
        move16();
        invVal = BASOP_Util_Divide3232_Scale_newton( MAXVAL_WORD32, maxWithSign_fx( r ), &invVal_e );


        tmp_e = add( 31, sub( HEADROOM_LEFT_1, *g_e ) );
        tmp_e = sub( 32, *g_e );
        singularVectors_Left_64[currChannel][currChannel] = W_sub( singularVectors_Left_64[currChannel][currChannel], W_shr( W_deposit32_h( *g ), tmp_e ) ); // here, the exponent goes up.

        FOR( iCh = add( currChannel, 1 ); iCh < nChannelsC; iCh++ )
        {
            Word16 magic_shift;
            Word32 factor1;
            Word32 factor2;
            Word32 f; // = norm / r
                      //            Word16 f_e;      // not really needed
            Word16 magic_shift;

            norm_64 = 0;
            for ( jCh = currChannel; jCh < nChannelsL; jCh++ )
@@ -1081,10 +1082,10 @@ static void biDiagonalReductionLeft_64(
            norm_x_e = W_norm( norm_64 );
            norm_x = W_extract_h( W_shl( norm_64, norm_x_e ) );
            f = Mpy_32_32( norm_x, invVal );
            magic_shift = 31 - HEADROOM_LEFT_3 + norm_x_e - r_e;
            magic_shift = norm_x_e-r_e+(28-HEADROOM_LEFT_3);	// works with headroom_left: 1,2,2
            FOR( jCh = currChannel; jCh < nChannelsL; jCh++ )
            {
                factor1 = W_extract_l( W_shr( singularVectors_Left_64[jCh][currChannel], HEADROOM_LEFT_3 ) );
                factor1 = W_extract_h( W_shl( singularVectors_Left_64[jCh][currChannel], 32-HEADROOM_LEFT_3 ) );
                singularVectors_Left_64[jCh][iCh] = W_add( singularVectors_Left_64[jCh][iCh], W_shr( W_mult0_32_32( f, factor1 ), magic_shift ) );
            }
        }
@@ -1112,8 +1113,8 @@ static void biDiagonalReductionRight_64(
    Word16 idx;

#define HEADROOM_RIGHT_1 2
#define HEADROOM_RIGHT_2 ( HEADROOM_RIGHT_1 + 1 )
#define HEADROOM_RIGHT_3 3
#define HEADROOM_RIGHT_2 2
#define HEADROOM_RIGHT_3 2


    ( *g ) = 0;
@@ -1159,9 +1160,9 @@ static void biDiagonalReductionRight_64(
            move32();
            move16();
            factor2 = W_extract_l( W_shr( singularVectors_Left_64[currChannel][idx], HEADROOM_RIGHT_1 ) );
            tmp_e = sub( *g_e, HEADROOM_RIGHT_1 );
            tmp_e = sub( HEADROOM_RIGHT_1, *g_e );
            tmpmul = W_mult0_32_32( *g, factor2 );
            tmpmul = W_shl( tmpmul, tmp_e );
            tmpmul = W_shr( tmpmul, tmp_e );
            r_64 = W_sub( tmpmul, norm_64 );
            r_e = W_norm( r_64 );
            r = W_extract_h( W_shl( r_64, r_e ) );
@@ -1170,8 +1171,8 @@ static void biDiagonalReductionRight_64(
            move16();
            invVal = BASOP_Util_Divide3232_Scale_newton( MAXVAL_WORD32, maxWithSign_fx( r ), &invVal_e );

            magic_shift = 32 - *g_e;
            singularVectors_Left_64[currChannel][idx] = W_sub( singularVectors_Left_64[currChannel][idx], W_shr( W_deposit32_h( *g ), magic_shift ) ); // here, the exponent goes up
            tmp_e = sub( 32, *g_e );
            singularVectors_Left_64[currChannel][idx] = W_sub( singularVectors_Left_64[currChannel][idx], W_shr( W_deposit32_h( *g ), tmp_e ) ); // here, the exponent goes up

            FOR( iCh = idx; iCh < nChannelsL; iCh++ )
            {
@@ -1188,11 +1189,11 @@ static void biDiagonalReductionRight_64(
                norm_x_e = W_norm( norm_64 );
                norm_x = W_extract_h( W_shl( norm_64, norm_x_e ) );
                f = Mpy_32_32( norm_x, invVal );
                magic_shift = 25 + norm_x_e - r_e; // headroom 3 FIXME: Why does this work?
                magic_shift = norm_x_e-r_e   +29-(HEADROOM_RIGHT_3); // works with headroom_right: 2,2,2

                FOR( jCh = idx; jCh < nChannelsC; jCh++ )
                {
                    factor2 = W_extract_l( W_shr( singularVectors_Left_64[currChannel][jCh], HEADROOM_RIGHT_3 ) );
                    factor2 = W_extract_h( W_shl( singularVectors_Left_64[currChannel][jCh], 32-HEADROOM_RIGHT_3 ) );
                    singularVectors_Left_64[iCh][jCh] = W_add( singularVectors_Left_64[iCh][jCh], W_shr( W_mult0_32_32( f, factor2 ), magic_shift ) );
                }
            }