precision for very large NSIDE parameters > 2^25

[ntessore]

When the NSIDE parameter is larger than 2^25 for pix2vec/vec2pix or 2^26 for pix2ang/ang2pix, respectively, the transforms suffer from a loss of precision.

There is a short script tools/nside_precision.py to check. It appears that only the first row of base pixels suffers, so this needs further investigation.

[ntessore]

I believe the loss of precision happens in hpd2loc:

static t_loc hpd2loc(int64_t nside, t_hpd hpd, double dx, double dy) {
    double z, s, phi;
    const double x = (hpd.a - hpd.b + dx)/nside;
    const double y = (hpd.a + hpd.b - nside + dy)/nside;
    const int32_t r = 1 - hpd.f/4;
    const double m = 1 - r*y;
    if (m < 1) {
        // polar cap
        const double tmp = m*m*(1./3.);
        z = r*(1. - tmp);
        s = sqrt(tmp*(2.-tmp));
        phi = (PI/4.)*(jpll[hpd.f] + x/m);
    } else {
        // equatorial region
        z = (y+r)*(2./3.);
        s = sqrt((1.+z)*(1.-z));
        phi = (PI/4.)*(jpll[hpd.f] + x);
    }
    return (t_loc){z, s, phi};
}

Near the poles fabs(y) is only about 1./nside from unity.

[mreineck]

Hi, just came across this by accident.

If I run tools/nside_precision.py, I only see failures in the RING variants,not the NEST ones. That would seem to exclude hpd2loc as the culprit. Do you see the failures in nested ordering too?

[...]
2^24      angXnest ....  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          angXring ....  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          vecXnest ....  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          vecXring ....  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          ringXnest ...  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          nestXring ...  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
2^25      angXnest ....  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          angXring ....  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          vecXnest ....  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          vecXring ....  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          ringXnest ...  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          nestXring ...  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
2^26      angXnest ....  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          angXring ....  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          vecXnest ....  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          vecXring ....  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          ringXnest ...  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          nestXring ...  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
2^27      angXnest ....  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          angXring ....  [0 0 0 0 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          vecXnest ....  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          vecXring ....  [0 0 0 0 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          ringXnest ...  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          nestXring ...  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
2^28      angXnest ....  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          angXring ....  [0 0 0 0 0 0 0 0 0 0 0 0]  [1 1 1 1 1 1 1 1 1 1 1 1]  [0 0 0 0 0 0 0 0 1 1 1 1]
          vecXnest ....  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          vecXring ....  [0 0 0 0 0 0 0 0 0 0 0 0]  [1 1 1 1 1 1 1 1 1 1 1 1]  [0 0 0 0 0 0 0 0 1 1 1 1]
          ringXnest ...  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          nestXring ...  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
2^29      angXnest ....  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          angXring ....  [0 0 0 0 0 0 0 0 0 0 0 0]  [1 1 1 1 1 1 1 1 1 1 1 1]  [0 0 0 0 0 0 0 0 1 1 1 1]
          vecXnest ....  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          vecXring ....  [0 0 0 0 0 0 0 0 0 0 0 0]  [1 1 1 1 1 1 1 1 1 1 1 1]  [0 0 0 0 0 0 0 0 1 1 1 1]
          ringXnest ...  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]
          nestXring ...  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1 1 1 1 1]

[mreineck]

Sorry, I was stupid ... for nested ordering the pixels with the extreme indices are typically not close to the poles, so it's natural that they won't fail.

[ntessore]

Hi @mreineck, I'm glad you are having a look. And you are raising a very good point, the problematic pixels (if the problem is indeed as above, which I still think is true) are still in there, just somewhere else. One could probably find them with just an ordering conversion?

[mreineck]

You should be able to find the problematic pixels in nested scheme just by running ring2nest on ipixv.

If the problem does not occur with the reference implementation (and I don't think it does), your hypothesis about hpd2loc is probably right.

[ntessore]

@mreineck You were absolutely right, this does not happen in the healpix_bare implementation. Part of the loss in precision was indeed in hpd2loc, and part was in an optimisation introduced in loc2hpd. All fixed now, thanks a lot for your insight!