X-Authentication-Warning: delorie.com: mail set sender to djgpp-bounces using -f
Message-Id: <3.0.1.16.20040216231142.38e7a7cc@earthlink.net>
X-Sender: ethros AT earthlink DOT net
X-Mailer: Windows Eudora Light Version 3.0.1 (16)
Date: Mon, 16 Feb 2004 23:11:42 -0500
To: djgpp AT delorie DOT com, Eli Zaretskii <eliz AT elta DOT co DOT il>,
   Charles Sandmann <sandmann AT clio DOT rice DOT edu>
From: Ethan Rosenberg <ethros AT earthlink DOT net>
Subject: Re: Cross Platform Incompatabilites? - code fragments
Mime-Version: 1.0
Content-Type: text/plain; charset="us-ascii"

Dear Listmembers -


MUCH THANKS IN ADVANCE FOR ANY AVICE WHICH WILL BE GIVEN.

>Does that mean that the modified program produces _exactly_ the same
>results as before, on each platform?  Or do you see any changes in the
>results with the _CRT0_FLAG_FILL_DEADBEEF bit set?

No change wantsoever!!!


The following are code FRAGMENTS from the program.  My feelings are that
the error is in the integration routine.  etot is defined as FLOAT.

I have extracted these fragments for a program of approx. 3000 lines of
code.  I hope I am not burdening you with too much.

************


Ethan Rosenberg
*********************
/* fft.h: Defs for fft routines
 * [Part of simple-fft-1.5.tar.Z (ver 1.5)]*/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include <conio.h>
#include <stdlib.h>
#include <pc.h>

/*
 * If you want double-precision floats used everywhere,
 * define this to be "double"
 */
#define FLOAT 		float

/*
 * Since every math.h seems to have their own idea of what the
 * elements in struct complex are called, we define our own.
 */
struct _complx
{
  FLOAT re;
  FLOAT im;
};
typedef struct _complx complx;


/*
 *  This PERMUTE could be defined as a table lookup (Sampson did this),
 *  but I found that to be significantly slower than just computing it
 *  inline.  Of course, this is w/ GCC, and not MS-DOS Turbo C.
 *  I'll leave it as a macro; you can tweak with it.
 */

#define PERMUTE(x, y)	reverse((x), (y))

#define REAL(x)		wave[(x)].re
#define IMAG(x)		wave[(x)].im
#define ALPHA		0.54	/* ala hanning */
#define FFT_LEN     16384


/*
 *  Bit reverser for unsigned ints
 *  Reverses 'bits' bits.
 */
inline const unsigned int reverse (unsigned int val, int bits)
{
  unsigned int retn = 0;

  while (bits--)
    {
      retn <<= 1;
      retn |= (val & 1);
      val >>= 1;
    }
  return (retn);
}


/*
 *  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)); */
}

/* 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;
}
...snip ....
for (l = 0; l < FFT_LEN; l++)
            fft_out[l] = amp(l, wave, 14);
        for(l = 0; l < FFT_LEN/2; l++)
            temp2[l] += fft_out[l]*fft_out[l]*2;
e_tot = trap_integ(temp2,103,1537,0.2,3.0);