[Proj] Polar stereographic, different values on different platforms?

[Glynn Clements]

Frank Warmerdam wrote:

> > However, I see two areas for improvement here:
> > 
> > 1)  a rigorous test script to test the output of proj against known
> > correct values (within some tolerance).  It should test forward and
> > reverse transformations, as many of the projection/datum combinations
> > as is practical, using points from all eight octants of the globe.  If
> > this existed (and covered this particular projection), the bug would
> > never have been released.  I've offered to make a start at this if
> > this doesn't exist already.
> 
> We have a modest test suite in the proj/nad directory (do "make check")
> but unfortunately it is not convenient to run on win32/msvc builds
> since it depends on a unix style shell.
> 
> It is my hope that the MetaCRS project will publish "test data files"
> at some point which each of the software packages (proj.4, csmap,
> proj4js,geotools) can use as a base for a substantial regression test.

Even if you lack external test data, testing different versions and
different builds of the program against each other can often detect
problems; it would certainly have caught this one.

Comparing the results of debug and optimised builds is always a good
idea, particularly for numerical software. If you look at the errata
when a new version of a compiler is released, it's invariably the case
that most of bugs only applied when optimisation was enabled.

FWIW, i'm really curious about how the lack of that option came to
result in an error of 32km.

The description of the /Op option suggests that it's the same as
-ffloat-store on gcc. But gcc gives the same answer (to 2 decimal
places, i.e. centimetre precision) with or without that option. It
also gives the same answer with -funsafe-math-optimizations, which can
introduce more substantial errors (e.g. it allows replacing x/y with
x*(1/y)).

That suggests that it isn't just a case of minor rounding errors
exploding due to algorithmic instability.

I also note that the documentation for /Op says:

	http://msdn.microsoft.com/en-us/library/aa984742(VS.71).aspx

	Note The /Op option disables inline generation of
	floating-point functions. The standard run-time library
	routines are used instead. For more information, see the /Oi
	option.

The documentation for /Oi says:

	http://msdn.microsoft.com/en-us/library/f99tchzc(VS.71).aspx

	The floating-point functions listed below do not have true
	intrinsic forms. If you use the Generate Intrinsic Functions
	option, the listed functions are replaced with versions that
	pass arguments directly to the floating-point chip rather than
	pushing them onto the program stack:

	acos 	cosh 	pow 	tanh
	asin 	fmod 	sinh 	 

	The floating-point functions listed below have true intrinsic
	forms when you specify both /Oi and /Og (or any option that
	includes /Og: /Ox, /O1, and /O2):

	atan 	exp 	sin
	atan2 	log 	sqrt
	cos 	log10 	tan

Switching to an inline implementation would seem to provide more scope
for substantial errors. It would also mean that the issue would never
be observed with gcc, as the problem would be in code generated by
MSVC rather than in the run-time.

-- 
Glynn Clements <glynn at gclements.plus.com>