X-Authentication-Warning: delorie.com: mail set sender to djgpp-bounces using -f
Message-Id: <3.0.1.16.20040302202549.3957c428@earthlink.net>
X-Sender: ethros AT earthlink DOT net
X-Mailer: Windows Eudora Light Version 3.0.1 (16)
Date: Tue, 02 Mar 2004 20:25:49 -0500
To: Hans-Bernhard Broeker <broeker AT physik DOT rwth-aachen DOT de>
From: Ethan Rosenberg <ethros AT earthlink DOT net>
Subject: Re: Cross Platform Incompatabilites? - code fragments
Cc: djgpp AT delorie DOT com
In-Reply-To: <200402261347.i1QDlDC31668@ac3b07.physik.rwth-aachen.de>
References: <3 DOT 0 DOT 1 DOT 16 DOT 20040216231142 DOT 38e7a7cc AT earthlink DOT net>
Mime-Version: 1.0
Content-Type: text/plain; charset="iso-8859-1"
Content-Transfer-Encoding: 8bit
X-MIME-Autoconverted: from quoted-printable to 8bit by delorie.com id i231QsoQ008570
Errors-To: nobody AT delorie DOT com
X-Mailing-List: djgpp AT delorie DOT com
X-Unsubscribes-To: listserv AT delorie DOT com
Precedence: bulk

Dear Hans and Listmembers -

>
>We haven't heard any feedback from you in ages... did you try any
>of the suggestions?  Did they work?
>[...]
>
>-- 
>Hans-Bernhard Broeker (broeker AT physik DOT rwth-aachen DOT de)
>Even if all the snow were burnt, ashes would remain.


Please pardon the delay.  I appreciate EVER SO MUCH your suggestions.  To
the best of my understanding, I have tried to implement them.  The answer
to the problem is still elusive...

I have compiled the program using double and long double.  The results follow:
[note HBR is my DOS computer; SARAH is my wife's Win2000 computer]
When the programs were run, on the windows OS, nothing else was running on
windows.  I have yet to run the programs on DOS 7.0 [the pure DOS shell of
Win 98]

[min is time, delta is the calculated value of interest, etot is the
integral under the FFT curve]

   HBR/float
             min            delta          etot
              0.45      141.79744  490616.065924
              3.15      -3.106526    1270.050473

    SARAH/float
              min            delta          etot
              0.45       5.886202      56.013115
              3.15      -3.106526    1270.050473

    HBR/double
               min            delta          etot
              0.45     102.283497  895621.961731
              3.15      -3.106527    1270.050546

    SARAH/double
             min            delta          etot
              0.45       0.399665     187.450401
              3.15      -3.106527    1270.050546

    HBR/long double
              min           delta          etot
              0.45     135.183805  650175.488059
              3.15      -3.106527    1270.050546

    SARAH/long double
                min             delta          etot
              0.45     135.183805  650175.488059
              3.15      -3.106527    1270.050546

The problem with the long double calculations, I think, is that the sqrt
function assumes double for input.

The FFT routine which I have not modified at all is:

/*
 *  Here is the real work-horse.
 *  It's a generic FFT, so nothing is lost or approximated.
 *  The samples in wave[] should be in order, and they
 *  will be decimated when fft() returns.
 */
void fft (complx wave[], int bits)
{
  register int	loop, loop1, loop2;
  unsigned	i1;		/* going to right shift this */
  int		i2, i3, i4, y;
  FLOAT		a1, a2, b1, b2, z1, z2;

  i1 = FFT_LEN / 2;
  i2 = 1;

  /* perform the butterfly's */

  for (loop = 0; loop < bits; loop++)
    {
      i3 = 0;
      i4 = i1;

      for (loop1 = 0; loop1 < i2; loop1++)
	{
	  y  = PERMUTE(i3 / (int)i1, bits);
	  z1 = CosineTable[y];
	  z2 = -SineTable[y];

	  for (loop2 = i3; loop2 < i4; loop2++)
	    {
	      a1 = REAL(loop2);
	      a2 = IMAG(loop2);

	      b1 = z1 * REAL(loop2+i1) - z2 * IMAG(loop2+i1);
	      b2 = z2 * REAL(loop2+i1) + z1 * IMAG(loop2+i1);

	      REAL(loop2) = a1 + b1;
	      IMAG(loop2) = a2 + b2;

	      REAL(loop2+i1) = a1 - b1;
	      IMAG(loop2+i1) = a2 - b2;
	    }

	  i3 += (i1 << 1);
	  i4 += (i1 << 1);
	}

      i1 >>= 1;
      i2 <<= 1;
    }
}


Recognizing the the problem may be from overflow of the float, certain
changes were made in the program:

The original amplitude calculating functions:

/*
 *  Calculate the hypotenuse of a right triangle.  The function in math.h
 *  assumes the the input and return values are double.  This uses FLOAT
 */
FLOAT hypot2(FLOAT a, FLOAT b)
{
  FLOAT c,asq,bsq;
  double aa,bb,cc;
  asq = a*a;
  bsq = b*b;
  aa = (double)asq;
  bb = (double)bsq;
  cc = sqrt(aa+bb);
  c = (FLOAT)sqrt((double)(a*a) + (double)(b*b));
  return c;
}


/*
 *  Calculate amplitude of component n in the decimated wave[] array.
 */

FLOAT amp (int n, complx wave[], int bits)
{
  n = PERMUTE (n, bits);
  return (hypot2 (REAL(n), IMAG(n)));

/*  return (REAL(n), IMAG(n)); */
}

Changed to:


/*
 *  Calculate the hypotenuse of a right triangle.  The function in math.h
 *  assumes the the input and return values are double.  This uses FLOAT
 */
FLOAT hypot2(FLOAT a, FLOAT b)
{
  FLOAT c,asq,bsq;
  double aa,bb,cc;
  asq = a*a;
  bsq = b*b;
  aa = (double)asq;
  bb = (double)bsq;
  cc = sqrt(aa+bb);
  c = (FLOAT)sqrt(((double)a*(double)a) + ((double)b*(double)b));
  return c;
}


/*
 *  Calculate amplitude of component n in the decimated wave[] array.
 */

FLOAT amp (int n, complx wave[], int bits)
{
  n = PERMUTE (n, bits);
/*  return (hypot2 (REAL(n), IMAG(n)));*/
  return (hypot (REAL(n), IMAG(n)));

/*  return (REAL(n), IMAG(n)); */
}

>>>NOTE  the hypot2 was ONLY used for float calls <<<<<
>>> If the problem is in the c*c = a*a + b*b, then I could use a routine
for that in Long Double

<<<<
The original integration routine:

/* Integrate the PSD function, using the trapaziodal rule */

FLOAT trap_integ(FLOAT temp2[],int start_pt,int end_pt,FLOAT
start_val,FLOAT end_val)
{
    FLOAT sum;
    int i;

    sum = 0.0;
    sum += temp2[start_pt];
    for(i = start_pt+1; i < end_pt; i++)
        sum += 2*temp2[i];
    sum += temp2[end_pt];
    sum = (end_val-start_val)/(2*((FLOAT)end_pt-(FLOAT)start_pt))*sum;
    return sum;
}


Changed to:

[to sort the numbers prior to the sum]


/* Integrate the PSD function, using the trapaziodal rule */

FLOAT trap_integ(FLOAT temp2[],int start_pt,int end_pt,FLOAT
start_val,FLOAT end_val)
{
    FLOAT sum, temp3[1450];
    int i, j, jj;

    sum = 0.0;
    for (i = 0; i < 1450; i++)
    	temp3[i] = 0.0;

    sum += temp2[start_pt];
/*    for(i = start_pt+1; i < end_pt; i++)
        sum += 2*temp2[i];*/
	jj = 0;
    for(i = start_pt+1, j = 0; i < end_pt; i++, j++)
    {
         temp3[j] = temp2[i];
         jj = j;
    }
    qsort(temp3, jj, sizeof(FLOAT), mycomp);
    for(i = 0; i < jj; i++)
        sum += 2*temp3[i];
    sum += temp2[end_pt];
    sum = (end_val-start_val)/(2*((FLOAT)end_pt-(FLOAT)start_pt))*sum;
    return sum;
}

=============

I read the paper "What every computer Scientist needs to know about
floating point..".  My background is NOT in numerical analysis, so, I do
not understand it totally.  I hope I can find a mathematician that can help.

MUCH THANKIS IN ADVANCE FOR ALL YOUR EXCELLENT HELP.

Ethan Rosenberg