Precision issues near poles for large values of nside

[astrofrog]

At the moment, our longitude/latitude to pixel conversion has a lower precision than that of the official HEALPix library close to the poles. This is visible when making a map of pixel values for 2**28 for the following range of coordinates:

lat_min = 89.999996
lat_max = 89.999997
lon_min = -5
lon_max = +5

Original HEALPix library (via healpy)

This package

I think I've tracked down the issue to the xyztohp function, and I think the issue could be with these lines:

        // solve eqn 20: k = Ns - xx (in the northern hemi)
		root = (1.0 - vz*zfactor) * 3.0 * mysquare(Nside * (2.0 * phi_t - pi) / pi);
		kx = (root <= 0.0) ? 0.0 : sqrt(root);

        // solve eqn 19 for k = Ns - yy
		root = (1.0 - vz*zfactor) * 3.0 * mysquare(Nside * 2.0 * phi_t / pi);
		ky = (root <= 0.0) ? 0.0 : sqrt(root);

		if (north) {
			xx = Nside - kx;
			yy = Nside - ky;
		} else {
			xx = ky;
			yy = kx;
		}

The issue is that doubles actually have less precision information than int64 (only 52 bits are allocated to significant precision), so by going through double and doing operations such as sqrt etc which further degrade precision, I think we are losing precision there. So this may not be a 32-bit vs 64-bit int issue but rather an issue with floating point precision. We should investigate whether we can improve the precision here, for example by using long double or finding other ways to avoid precision loss. For now, I will update the documentation to mention the resolution to which the results can be trusted.

@dstndstn - if you have any ideas on improving precision here, please feel free to comment!

[dstndstn]

Hi, If this is at Dec = 89.99999, then I'm guessing you're seeing cancellation error in (1. - vz), where (I'm guessing) vz is the unit-sphere z vector -- which will be approaching 1 as Dec approaches 90 degrees. Maybe you could use a different trig identity here to avoid the cancellation error -- sin(Dec) = cos(pi/2 - Dec) -- for vz > 0.9 or 0.5 or whatever. cheers, --dustin

…

On Mon, Sep 25, 2017 at 9:30 AM, Thomas Robitaille ***@***.*** > wrote: At the moment, our longitude/latitude to pixel conversion has a lower precision than that of the official HEALPix library. This is visible when making a map of pixel values for 2**28: *Original HEALPix library (via healpy)* [image: healpy_true] <https://user-images.githubusercontent.com/314716/30810655-aeabe51e-a1fd-11e7-82d0-e956b59b1029.png> *This package* [image: healpy_false] <https://user-images.githubusercontent.com/314716/30810664-b4bbafac-a1fd-11e7-8cd9-bad400e11f53.png> I think I've tracked down the issue to the xyztohp function, and specifically these lines: // solve eqn 20: k = Ns - xx (in the northern hemi) root = (1.0 - vz*zfactor) * 3.0 * mysquare(Nside * (2.0 * phi_t - pi) / pi); kx = (root <= 0.0) ? 0.0 : sqrt(root); // solve eqn 19 for k = Ns - yy root = (1.0 - vz*zfactor) * 3.0 * mysquare(Nside * 2.0 * phi_t / pi); ky = (root <= 0.0) ? 0.0 : sqrt(root); if (north) { xx = Nside - kx; yy = Nside - ky; } else { xx = ky; yy = kx; } and in particular the issue is with the x value (y seems fine I think). The issue is that doubles actually have less precision information than int64 (only 52 bits are allocated to significant precision), so by going through double and doing operations such as sqrt etc which further degrade precision, I *think* we are losing precision there. So this may not be a 32-bit vs 64-bit int issue but rather an issue with floating point precision. We should investigate whether we can improve the precision here, for example by using long double or finding other ways to avoid precision loss. For now, I will update the documentation to mention the resolution to which the results can be trusted. @dstndstn <https://github.com/dstndstn> - if you have any ideas on improving precision here, please feel free to comment! — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <https://github.com/cdeil/healpix/issues/34>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABBD_Z5p3TkcPkeib_Fg2WSXTGlLaSPJks5sl6sDgaJpZM4PivH7> .

[dstndstn]

I had a further look at this, and indeed the precision loss is in converting to unit-sphere coordinates. In order to fix this, I think one would have to pass Dec, in addition to z, into xyztohp -- or maybe pass in (1-z) and z -- or event sqrt(1-z) and z if you want to get really fancy -- and then use trig identities.

[dstndstn]

PS, it's things like this that make me hate the healpix paper.

[dstndstn]

In e1df57b I used "long double" trig functions and the result seems to fix this resolution problem.

I only did that for the one function used by the library, but presumably we could use this approach more widely. I assume there is some computational cost for using the long double trig functions, so I don't know whether to use this everywhere or not.

[cdeil]

@dstndstn - Do you understand this? Is "double" 64 bit and "long double" more? Or are both 64 bit and just the C functions called are higher precision somehow?

Maybe open a PR with the commit you have, and wait for @astrofrog to comment whether to merge as-is or apply more widely?

[dstndstn]

"long double" is longer than 64 bits. How long depends on implementation;
https://en.wikipedia.org/wiki/Long_double

[dstndstn]

PR #93

[dstndstn]

PR #94 instead?

[fxpineau]

I don't know if it helps but in the case of the CDS Java HEALPix lib I developed, I replaced \sqrt{3(1-\sin\delta)} by the equivalent \sqrt{6}\cos(\delta/2 + \pi/4) for the Collignon projection (=> we also save a sqrt computation, which is a slow operation).
Basically, the costly operations are: a sine for the Cylindrical projection, a cosine for the Collignon projection.

[fxpineau]

(In your lib, see more efficient algo to compute the Morton code for not small orders here: https://graphics.stanford.edu/~seander/bithacks.html#InterleaveTableLookup)