On 1/24/07, Greg Chicares <chicares@cox.net> wrote:
> On 2007-1-24 7:52 UTC, Aron Rubin wrote:
> > I have been pulling my hair out for several weeks here trying to
> > figure out why there are strange artifacts in my software. Tonight I
> > finally figured it out. The math library is broken (in my install).
>
> But MinGW has no real 'libm'--no math library per se; it gets
> C math functions from the msvc runtime library provided as
> part of the operating system. A few functions are overridden
> by 'libmingwex', but by no means all.

I realize this but lets say for example that the actual functions in msvc were single precision float but the prototypes in math.h were double it would cause behavior like this.

>
> > All math functions return bogus numbers.
>
> It would help to have a testcase so we can see which math
> functions you're using, the options you use for compiling,
> and a copy of the results you see.

<pre>/* ------------------------------
   Function: geo_geod_to_utm_hifi
     Transform a geographic (<latitude> <longitude>) coordinate to <UTM>.

   Parameters:
     ellipsoid_id - [in]  the identifier for the ellipsoid to use found by
                    <geo_ellipsoid_by_name>
     lat          - [in]  <latitude> (degrees)
     lon          - [in]  <longitude> (degrees)
     utm_zone     - [out] <UTM> zone name consists of a numerical longitude zone and a
                    letter identifying the latitude zone. The zone number
                    actually specifies the reference meridian for <easting>.
     easting      - [out] <easting> (meters)
     northing     - [out] <northing> (meters)

   Notes:
     - Equations from USGS Bulletin 1532
     - East Longitudes are positive, West longitudes are negative.
     - North latitudes are positive, South latitudes are negative.
     - Lat and lon are in decimal degrees.
     - Originally written by Chuck Gantz <chuck.gantz@globalstar.com > with no
       stated license.
   ------------------------------ */
void geo_geod_to_utm_hifi( int ellipsoid_id, double lat, double lon,
               char *utm_zone, double *easting, double *northing ) {
  double a = geo_ellipsoid_table[ellipsoid_id].equatorial_radius;
  double eccSquared = geo_ellipsoid_table[ellipsoid_id].eccentricity_squared;
  double k0 = 0.9996;

  double lon_origin;
  double eccPrimeSquared;
  double N, T, C, A, M;
    
  //Make sure the longitude is between -180.00 .. 179.9
  double lontmp = (lon + 180.0) - floor( (lon + 180.0)*(1.0/360.0) )*360.0 - 180.0;

  double lat_rad = lat*DEGTORAD;
  double lon_rad = lontmp*DEGTORAD;
  double lon_origin_rad;
  int    zone_num;

  printf( "floor val = %g\n", floor( (lon + 180.0)*(1.0/360.0) ) );

  zone_num = (int)((lontmp + 180.0)*(1.0 /6.0) + 1.0);
  
  if( lat >= 56.0 && lat < 64.0 && lontmp >= 3.0 && lontmp < 12.0 )
    zone_num = 32;

  // Special zones for Svalbard
  if( lat >= 72.0 && lat < 84.0 ) {
    if(      lontmp >= 0.0  && lontmp <  9.0 ) zone_num = 31;
    else if( lontmp >= 9.0  && lontmp < 21.0 ) zone_num = 33;
    else if( lontmp >= 21.0 && lontmp < 33.0 ) zone_num = 35;
    else if( lontmp >= 33.0 && lontmp < 42.0 ) zone_num = 37;
  }
  lon_origin = (zone_num - 1)*6 - 180 + 3;  //+3 puts origin in middle of zone
  lon_origin_rad = lon_origin*DEGTORAD;

  // set the UTM Zone from the latitude and longitude
  sprintf( utm_zone, "%d%c", zone_num, geo_utm_lat_zone_by_lat( lat ) );

  eccPrimeSquared = (eccSquared)/(1-eccSquared);

  N = a/sqrt(1 - eccSquared*sin(lat_rad)*sin(lat_rad));
  T = tan(lat_rad)*tan(lat_rad);
  C = eccPrimeSquared*cos(lat_rad)*cos(lat_rad);
  A = cos(lat_rad)*(lon_rad - lon_origin_rad);

  M = a*((1 - eccSquared/4 - 3*eccSquared*eccSquared/64    - 5*eccSquared*eccSquared*eccSquared/256)*lat_rad
     - (3*eccSquared/8 + 3*eccSquared*eccSquared/32    + 45*eccSquared*eccSquared*eccSquared/1024)*sin(2*lat_rad)
     +                (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(4*lat_rad)
     -                                               (35*eccSquared*eccSquared*eccSquared/3072)*sin(6*lat_rad));

  *easting = (k0*N*(A+(1-T+C)*A*A*A/6 + (5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120) +
          500000.0);

  *northing = (k0*(M+N*tan(lat_rad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
                     + (61-58*T+T*T+600*C-330*eccPrimeSquared)*A*A*A*A*A*A/720)));
  if( lat < 0 )
    *northing += 10000000.0; //10,000,000 meter offset for southern hemisphere
}

</pre>
>
> > Just to get it out of the
> > way, yes I am including math.h. I have the latest mingw-runtime
> > installed (11/18/06).
> >
> > For the moment I have made macro versions of some functions and
> > -ffast-math seems to take care of the rest (I hope). What could be
> > causing this?
>
> It depends on the sort of bogosity you're observing.
>
> -------------------------------------------------------------------------
> Take Surveys. Earn Cash. Influence the Future of IT
> Join SourceForge.net's Techsay panel and you'll get the chance to share your
> opinions on IT & business topics through brief surveys - and earn cash
> http://www.techsay.com/default.php?page=join.php&p=sourceforge&CID=DEVDEV
> _______________________________________________
> Mingw-msys mailing list
> Mingw-msys@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/mingw-msys
>