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

the calculation of the exponents is stable enough to detect overflows in the mantissas.

parent 1eb687a6
Loading
Loading
Loading
Loading
+63 −4
Original line number Diff line number Diff line
@@ -1011,32 +1011,49 @@ static void HouseholderReduction_fx(
    {
	    int i,j;
	    int magicnum;
	    static int replacecnt=0;
	    static int bettercnt=0;
	    static int totalcnt=0;

	    magicnum=singularVectors_Left_fx_e[0][0]+W_norm(singularVectors_Left_64[0][0])+1;
	    for (i=0;i<nChannelsL;i++)
	    {
			printf("EXPONENT%02d: IN%3d ",i,singularVectors_Left_e);
		    for (j=0;j<nChannelsC;j++)
		    {
			    Word16 n;
			    Word32 tmp;
			    unsigned int x,y;
				int mine,theirs;
			    int k;
			    int cnt;
			    n=norm_l(singularVectors_Left_fx[i][j]);
			    tmp=singularVectors_Left_fx[i][j]<<n;
                            x=(unsigned int)tmp;
				printf("[(%2d)",singularVectors_Left_fx_e[i][j]-n);
				theirs=singularVectors_Left_fx_e[i][j]-n;

			    n=W_norm(singularVectors_Left_64[i][j]);
			    tmp=W_extract_h(W_shl(singularVectors_Left_64[i][j],n));
			    y=(unsigned int)tmp;

				n=32+singularVectors_Left_e-n;
				mine=n;
				printf("MINE:%2d]",n);
				printf("%08X/%08X ",x,y);
			    cnt=0;
				if (!((x^y)&0xf0000000)) bettercnt++;
				if (!((x^y)&0xff000000)) bettercnt++;
			//	if (mine==theirs)
				{
					singularVectors_Left_fx[i][j]=x;
					singularVectors_Left_fx_e[i][j]=mine;
					replacecnt++;
				}
			    totalcnt++;
		    }
			printf("\n");
	    }
	    printf("\nbetter %d/%d\n",bettercnt,totalcnt);
	    printf("\nbetter %d replace:%d  /%d\n",bettercnt,replacecnt,totalcnt);
    }

    /* SingularVecotr Accumulation */
@@ -1064,25 +1081,44 @@ static void biDiagonalReductionLeft_64(
    const Word16 currChannel /* Q0 */
)
{

	static FILE *f_debug=NULL;


    Word16 iCh, jCh;
    Word32 norm_x, g;
    Word16 norm_x_e, g_e;
    Word64 norm_64;
    norm_x=0;
    move32();
	if (f_debug==NULL) f_debug=fopen("left64.bin","wb");
	{
		unsigned int magic;
		magic=0xd00faffe;
		fwrite(&magic,sizeof(int),1,f_debug);
	}
	fwrite(&nChannelsL,sizeof(short),1,f_debug);
	fwrite(&nChannelsC,sizeof(short),1,f_debug);
	fwrite(&currChannel,sizeof(short),1,f_debug);
    IF( LT_16( currChannel, nChannelsL ) ) /* i <= m */
    {

        Word32 tmp;
        norm_64=0;
        move64();
	printf("\nNORM%02d: \x1b[1;34m",currChannel);
        FOR( jCh = currChannel; jCh < nChannelsL; jCh++ ) /* nChannelsL */
        {
            tmp=W_extract_l(W_shr(singularVectors_Left_64[jCh][currChannel],bitwindow));
		fwrite(&tmp,sizeof(int),1,f_debug);
            norm_64=W_add(norm_64,W_mult0_32_32(tmp,tmp));
		fwrite(&norm_64,sizeof(long long),1,f_debug);
		printf("%016llX ",norm_64);
        }
        norm_x_e=W_norm(norm_64);
        norm_x=W_extract_h(W_shl(norm_64, norm_x_e ));
        norm_x_e = add(sub(shl(bitwindow, 1), norm_x_e), 1 );	
	printf(" --> %08X<%2d\x1b[0m\n",norm_x,norm_x_e);
    }
    IF ( norm_x )
    {
@@ -1316,6 +1352,7 @@ static void biDiagonalReductionLeft_fx(
    Word16 norm_x_e, f_e, r_e;
    Word32 L_temp;
    Word16 L_temp_e;
	static FILE *f_debug2=NULL;

    secDiag[currChannel] = Mpy_32_32( *sig_x, *g ); /* exp(sig_x_e) */
    move32();
@@ -1328,6 +1365,16 @@ static void biDiagonalReductionLeft_fx(
    ( *g ) = 0;
    move32();

	if(f_debug2==NULL) f_debug2=fopen("leftfx.bin","wb");
	{
		unsigned int magic;
		magic=0xd00faffe;
		fwrite(&magic,sizeof(int),1,f_debug2);
	}
	fwrite(&nChannelsL,sizeof(short),1,f_debug2);
	fwrite(&nChannelsC,sizeof(short),1,f_debug2);
	fwrite(&currChannel,sizeof(short),1,f_debug2);
	
    IF( LT_16( currChannel, nChannelsL ) ) /* i <= m */
    {
        idx = currChannel;
@@ -1337,25 +1384,37 @@ static void biDiagonalReductionLeft_fx(
        {
            ( *sig_x ) = BASOP_Util_Add_Mant32Exp( *sig_x, *sig_x_e, L_abs( singularVectors[jCh][currChannel] ), singularVectors2_e[jCh][currChannel], sig_x_e ); /* exp(sig_x_e) */
        }

	fwrite(sig_x,sizeof(int),1,f_debug2);
	fwrite(sig_x_e,sizeof(short),1,f_debug2);
        IF( ( *sig_x ) ) /*(fabsf(*sig_x) > EPSILON * fabsf(*sig_x)) { */
        {
            Word16 invVal_e;
            Word32 invVal;
            invVal = BASOP_Util_Divide3232_Scale_newton( MAXVAL_WORD32, maxWithSign_fx( *sig_x ), &invVal_e );
	fwrite(&invVal,sizeof(int),1,f_debug2);
	fwrite(&invVal_e,sizeof(short),1,f_debug2);
            norm_x = 0;
            move32();
            norm_x_e = 0;
            move16();
		printf("\nNORM%02d: ",currChannel, norm_x,norm_x_e,(signed long long)norm_x*((signed long long)(*sig_x)*(signed long long)(*sig_x)),(*sig_x_e));
            FOR( jCh = idx; jCh < nChannelsL; jCh++ ) /* nChannelsL */
            {
                Word16 temp_e = norm_l( singularVectors[jCh][currChannel] );
		Word32 tmp;
		tmp=L_shl( singularVectors[jCh][currChannel],temp_e);
		fwrite(&tmp,sizeof(int),1,f_debug2);
                singularVectors[jCh][currChannel] = Mpy_32_32( L_shl( singularVectors[jCh][currChannel], temp_e ), invVal ); /* exp(sing_exp + (singularVectors_e - sig_x_e) */
		fwrite(&singularVectors[jCh][currChannel],sizeof(int),1,f_debug2);
                move32();
                singularVectors2_e[jCh][currChannel] = sub( add( invVal_e, sub( singularVectors2_e[jCh][currChannel], *sig_x_e ) ), temp_e );
                move16();
                norm_x = BASOP_Util_Add_Mant32Exp( norm_x, norm_x_e, Mpy_32_32( singularVectors[jCh][currChannel], singularVectors[jCh][currChannel] ), shl( singularVectors2_e[jCh][currChannel], 1 ), &norm_x_e ); /* exp(norm_x_e) */
		printf("%08X<%2d ",norm_x,norm_x_e);
		fwrite(&norm_x,sizeof(int),1,f_debug2);
		fwrite(&norm_x_e,sizeof(short),1,f_debug2);
            }
		printf("--> %08X<%2d  %016llX<%2d\n", norm_x,norm_x_e,(signed long long)norm_x*((signed long long)(*sig_x)*(signed long long)(*sig_x)),(*sig_x_e));
            IF( GT_16( norm_x_e, 0 ) )
            {
                norm_x = MAX_32;