Repository: CMB-S4/spt3g_software (language: python)

Issue:
Integer overflow in `maps.HealpixSkyMap.angle_to_pixel` for `nside` = 16384
```
m = maps.HealpixSkyMap(16384)
m.angle_to_pixel(0, -pi / 2)
>>> 18446744073709551615
```



Discussion / hints:
This issue stems from using int32 for all of the math.  There are int64 versions of all the chealpix functions, but `HealpixSkyMapInfo` class probably needs some TLC to make sure that all of the internal functions there are using int64 as well.

Fix patch:
diff --git a/maps/include/maps/FlatSkyMap.h b/maps/include/maps/FlatSkyMap.h
index b526b987..a2279d59 100644
--- a/maps/include/maps/FlatSkyMap.h
+++ b/maps/include/maps/FlatSkyMap.h
@@ -125,12 +125,12 @@ class FlatSkyMap : public G3FrameObject, public G3SkyMap {
 
 	std::vector<double> PixelToAngleGrad(size_t pixel, double h=0.001) const;
 
-	void GetRebinAngles(long pixel, size_t scale,
+	void GetRebinAngles(size_t pixel, size_t scale,
 	    std::vector<double> & alphas, std::vector<double> & deltas) const override;
 	void GetInterpPixelsWeights(double alpha, double delta,
-	    std::vector<long> & pixels, std::vector<double> & weights) const override;
+	    std::vector<size_t> & pixels, std::vector<double> & weights) const override;
 
-	std::vector<long> QueryDisc(double alpha, double delta,
+	std::vector<size_t> QueryDisc(double alpha, double delta,
 	    double radius) const override;
 
 	G3SkyMapPtr Rebin(size_t scale, bool norm = true) const override;
diff --git a/maps/include/maps/FlatSkyProjection.h b/maps/include/maps/FlatSkyProjection.h
index 6e9a58d9..e9e046aa 100644
--- a/maps/include/maps/FlatSkyProjection.h
+++ b/maps/include/maps/FlatSkyProjection.h
@@ -85,25 +85,27 @@ class FlatSkyProjection : public G3FrameObject {
 	double yres() const { return y_res_; };
 	double res() const { return y_res_; };
 
-	long XYToPixel(double x, double y) const;
-	std::vector<double> PixelToXY(long pixel) const;
+	size_t XYToPixel(double x, double y) const;
+	std::vector<double> PixelToXY(size_t pixel) const;
 	std::vector<double> XYToAngle(double x, double y, bool wrap_alpha=false) const;
 	std::vector<double> AngleToXY(double alpha, double delta) const;
-	std::vector<double> PixelToAngle(long pixel, bool wrap_alpha=false) const;
-	long AngleToPixel(double alpha, double delta) const;
+	std::vector<double> PixelToAngle(size_t pixel, bool wrap_alpha=false) const;
+	size_t AngleToPixel(double alpha, double delta) const;
 	std::vector<double> QuatToXY(quat q, bool local) const;
 	quat XYToQuat(double x, double y, bool local) const;
 
 	std::vector<double> XYToAngleGrad(double x, double y, double h=0.001) const;
-	std::vector<double> PixelToAngleGrad(long pixel, double h=0.001) const;
+	std::vector<double> PixelToAngleGrad(size_t pixel, double h=0.001) const;
 
-	void GetRebinAngles(long pixel, size_t scale,
+	size_t RebinPixel(size_t pixel, size_t scale) const;
+
+	void GetRebinAngles(size_t pixel, size_t scale,
 	    std::vector<double> & alphas, std::vector<double> & deltas,
 	    bool wrap_alpha=false) const;
 	void GetInterpPixelsWeights(double alpha, double delta,
-	    std::vector<long> & pixels, std::vector<double> & weights) const;
+	    std::vector<size_t> & pixels, std::vector<double> & weights) const;
 
-	std::vector<long> QueryDisc(double alpha, double delta, double radius,
+	std::vector<size_t> QueryDisc(double alpha, double delta, double radius,
 	    bool local) const;
 
 	FlatSkyProjection Rebin(size_t scale, double x_center = 0.0 / 0.0,
diff --git a/maps/include/maps/G3SkyMap.h b/maps/include/maps/G3SkyMap.h
index 5a7dc306..8e63347b 100644
--- a/maps/include/maps/G3SkyMap.h
+++ b/maps/include/maps/G3SkyMap.h
@@ -168,14 +168,14 @@ class G3SkyMap {
 	size_t QuatToPixel(quat q) const;
 
 	// Rebinning and interpolation
-	virtual void GetRebinAngles(long pixel, size_t scale,
+	virtual void GetRebinAngles(size_t pixel, size_t scale,
 	    std::vector<double> & alphas, std::vector<double> & deltas) const = 0;
-	G3VectorQuat GetRebinQuats(long pixel, size_t scale) const;
+	G3VectorQuat GetRebinQuats(size_t pixel, size_t scale) const;
 	virtual void GetInterpPixelsWeights(double alpha, double delta,
-	    std::vector<long> & pixels, std::vector<double> & weights) const = 0;
-	void GetInterpPixelsWeights(quat q, std::vector<long> & pixels,
+	    std::vector<size_t> & pixels, std::vector<double> & weights) const = 0;
+	void GetInterpPixelsWeights(quat q, std::vector<size_t> & pixels,
 	    std::vector<double> & weights) const;
-	double GetInterpPrecalc(const std::vector<long> &pixels,
+	double GetInterpPrecalc(const std::vector<size_t> &pixels,
 	    const std::vector<double> &weights) const;
 	double GetInterpValue(double alpha, double delta) const;
 	double GetInterpValue(quat q) const;
@@ -186,10 +186,10 @@ class G3SkyMap {
 	virtual boost::shared_ptr<G3SkyMap> Rebin(size_t scale, bool norm = true) const = 0;
 
 	/* Analogue to healpy.query_disc, returns list of pixels within a disc */
-	virtual std::vector<long> QueryDisc(double alpha, double delta, double radius) const {
+	virtual std::vector<size_t> QueryDisc(double alpha, double delta, double radius) const {
 		throw std::runtime_error("QueryDisc not implemented");
 	};
-	std::vector<long> QueryAlphaEllipse(double alpha, double delta, double a, double b) const;
+	std::vector<size_t> QueryAlphaEllipse(double alpha, double delta, double a, double b) const;
 
 	virtual bool IsDense() const {
 		throw std::runtime_error("Checking array density not implemented");
diff --git a/maps/include/maps/HealpixSkyMap.h b/maps/include/maps/HealpixSkyMap.h
index e306dea6..a28194d9 100644
--- a/maps/include/maps/HealpixSkyMap.h
+++ b/maps/include/maps/HealpixSkyMap.h
@@ -81,12 +81,12 @@ class HealpixSkyMap : public G3FrameObject, public G3SkyMap {
 	size_t AngleToPixel(double alpha, double delta) const override;
 	std::vector<double> PixelToAngle(size_t pixel) const override;
 
-	void GetRebinAngles(long pixel, size_t scale,
+	void GetRebinAngles(size_t pixel, size_t scale,
 	    std::vector<double> & alphas, std::vector<double> & deltas) const override;
 	void GetInterpPixelsWeights(double alpha, double delta,
-	    std::vector<long> & pixels, std::vector<double> & weights) const override;
+	    std::vector<size_t> & pixels, std::vector<double> & weights) const override;
 
-	std::vector<long> QueryDisc(double alpha, double delta,
+	std::vector<size_t> QueryDisc(double alpha, double delta,
 	    double radius) const override;
 
 	G3SkyMapPtr Rebin(size_t scale, bool norm = true) const override;
diff --git a/maps/include/maps/HealpixSkyMapInfo.h b/maps/include/maps/HealpixSkyMapInfo.h
index 4751ef67..b2461385 100644
--- a/maps/include/maps/HealpixSkyMapInfo.h
+++ b/maps/include/maps/HealpixSkyMapInfo.h
@@ -14,7 +14,8 @@ class HealpixSkyMapInfo : public G3FrameObject {
 	HealpixSkyMapInfo();
 	HealpixSkyMapInfo(const HealpixSkyMapInfo& info);
 
-	void initialize(size_t nside, bool nested=false, bool shifted=false);
+	void initialize(size_t nside_or_npix, bool nested=false,
+	    bool shifted=false, bool is_npix=false);
 
 	template <class A> void load(A &ar, unsigned v);
 	template <class A> void save(A &ar, unsigned v) const;
@@ -22,6 +23,7 @@ class HealpixSkyMapInfo : public G3FrameObject {
 	std::string Description() const override;
 
 	void SetNSide(size_t nside);
+	void SetNPix(size_t npix);
 	void SetNested(bool nested);
 	void SetShifted(bool shifted);
 
@@ -38,13 +40,15 @@ class HealpixSkyMapInfo : public G3FrameObject {
 	std::vector<double> PixelToAngle(size_t pixel) const;
 	size_t AngleToPixel(double alpha, double delta) const;
 
-	void GetRebinAngles(long pixel, size_t scale,
+	size_t RebinPixel(size_t pixel, size_t scale) const;
+
+	void GetRebinAngles(size_t pixel, size_t scale,
 	    std::vector<double> & alphas, std::vector<double> & deltas) const;
 
 	void GetInterpPixelsWeights(double alpha, double delta,
-	    std::vector<long> & pixels, std::vector<double> & weights) const;
+	    std::vector<size_t> & pixels, std::vector<double> & weights) const;
 
-	std::vector<long> QueryDisc(double alpha, double delta,
+	std::vector<size_t> QueryDisc(double alpha, double delta,
 	    double radius) const;
 
 private:
diff --git a/maps/src/FlatSkyMap.cxx b/maps/src/FlatSkyMap.cxx
index 7378b514..3ba0f5a8 100644
--- a/maps/src/FlatSkyMap.cxx
+++ b/maps/src/FlatSkyMap.cxx
@@ -738,19 +738,19 @@ std::vector<double> FlatSkyMap::PixelToAngleGrad(size_t pixel, double h) const {
 	return proj_info.PixelToAngleGrad(pixel, h);
 }
 
-void FlatSkyMap::GetRebinAngles(long pixel, size_t scale,
+void FlatSkyMap::GetRebinAngles(size_t pixel, size_t scale,
     std::vector<double> & alphas, std::vector<double> & deltas) const
 {
 	proj_info.GetRebinAngles(pixel, scale, alphas, deltas, false);
 }
 
 void FlatSkyMap::GetInterpPixelsWeights(double alpha, double delta,
-    std::vector<long> & pixels, std::vector<double> & weights) const
+    std::vector<size_t> & pixels, std::vector<double> & weights) const
 {
 	proj_info.GetInterpPixelsWeights(alpha, delta, pixels, weights);
 }
 
-std::vector<long>
+std::vector<size_t>
 FlatSkyMap::QueryDisc(double alpha, double delta, double radius) const
 {
 	return proj_info.QueryDisc(alpha, delta, radius, coord_ref == Local);
@@ -780,9 +780,7 @@ G3SkyMapPtr FlatSkyMap::Rebin(size_t scale, bool norm) const
 	for (auto i : *this) {
 		if (i.second == 0)
 			continue;
-		size_t x = (i.first % xpix_) / scale;
-		size_t y = ((size_t)(i.first / xpix_)) / scale;
-		size_t ip = x + y * out->xpix_;
+		size_t ip = proj_info.RebinPixel(i.first, scale);
 		(*out)[ip] += i.second / sqscal;
 	}
 
diff --git a/maps/src/FlatSkyProjection.cxx b/maps/src/FlatSkyProjection.cxx
index a250df49..744d505a 100644
--- a/maps/src/FlatSkyProjection.cxx
+++ b/maps/src/FlatSkyProjection.cxx
@@ -247,7 +247,7 @@ void FlatSkyProjection::SetRes(double res, double x_res)
 	SetXRes(x_res);
 }
 
-long
+size_t
 FlatSkyProjection::XYToPixel(double x, double y) const
 {
 	// Truncate X/Y coordinates to integer pixels and wrap to 1D.
@@ -255,20 +255,20 @@ FlatSkyProjection::XYToPixel(double x, double y) const
 	// The floor() is important to properly truncate negative numbers,
 	// which otherwise pile up at zero from both directions.
 
-	long ix = (long)floor(x + 0.5);
-	long iy = (long)floor(y + 0.5);
+	ssize_t ix = (ssize_t)floor(x + 0.5);
+	ssize_t iy = (ssize_t)floor(y + 0.5);
 	return (ix < 0 || iy < 0 || ix >= xpix_ || iy >= ypix_) ?
 	    -1 : ix + iy * xpix_;
 }
 
 std::vector<double>
-FlatSkyProjection::PixelToXY(long pixel) const
+FlatSkyProjection::PixelToXY(size_t pixel) const
 {
 	std::vector<double> out(2, -1);
 
 	if (pixel >= 0 && pixel < xpix_ * ypix_) {
-		out[0] = (int)(pixel % xpix_);
-		out[1] = (int)(pixel / xpix_);
+		out[0] = (ssize_t)(pixel % xpix_);
+		out[1] = (ssize_t)(pixel / xpix_);
 	}
 
 	return out;
@@ -544,7 +544,7 @@ FlatSkyProjection::XYToQuat(double x, double y, bool local) const
 }
 
 std::vector<double>
-FlatSkyProjection::PixelToAngle(long pixel, bool wrap_alpha) const
+FlatSkyProjection::PixelToAngle(size_t pixel, bool wrap_alpha) const
 {
 	if (pixel < 0 || pixel >= xpix_ * ypix_)
 		return {0., 0.};
@@ -552,14 +552,14 @@ FlatSkyProjection::PixelToAngle(long pixel, bool wrap_alpha) const
 	return XYToAngle(xy[0], xy[1], wrap_alpha);
 }
 
-long
+size_t
 FlatSkyProjection::AngleToPixel(double alpha, double delta) const
 {
 	std::vector<double> xy = AngleToXY(alpha, delta);
 	return XYToPixel(xy[0], xy[1]);
 }
 
-void FlatSkyProjection::GetRebinAngles(long pixel, size_t scale,
+void FlatSkyProjection::GetRebinAngles(size_t pixel, size_t scale,
     std::vector<double> & alphas, std::vector<double> & deltas,
     bool wrap_alpha) const
 {
@@ -620,7 +620,7 @@ FlatSkyProjection::XYToAngleGrad(double x, double y, double h) const
 }
 
 std::vector<double>
-FlatSkyProjection::PixelToAngleGrad(long pixel, double h) const
+FlatSkyProjection::PixelToAngleGrad(size_t pixel, double h) const
 {
 	if (pixel < 0 || pixel >= xpix_ * ypix_)
 		return {0., 0., 0., 0.};
@@ -629,19 +629,19 @@ FlatSkyProjection::PixelToAngleGrad(long pixel, double h) const
 }
 
 void FlatSkyProjection::GetInterpPixelsWeights(double alpha, double delta,
-    std::vector<long> & pixels, std::vector<double> & weights) const
+    std::vector<size_t> & pixels, std::vector<double> & weights) const
 {
 	std::vector<double> xy = AngleToXY(alpha, delta);
 	double x = xy[0];
 	double y = xy[1];
 
-	pixels = std::vector<long>(4, -1);
+	pixels = std::vector<size_t>(4, (size_t) -1);
 	weights = std::vector<double>(4, 0);
 
-	int x_1 = (int)floorf(x);
-	int x_2 = x_1 + 1;
-	int y_1 = (int)floorf(y);
-	int y_2 = y_1 + 1;
+	ssize_t x_1 = (ssize_t)floorf(x);
+	ssize_t x_2 = x_1 + 1;
+	ssize_t y_1 = (ssize_t)floorf(y);
+	ssize_t y_2 = y_1 + 1;
 	if (x_1 < 0 || y_1 < 0 || x_2 >= xpix_ || y_2 >= ypix_){
 		log_debug("Point lies outside of pixel grid\n");
 		return;
@@ -653,7 +653,7 @@ void FlatSkyProjection::GetInterpPixelsWeights(double alpha, double delta,
 	pixels[3] = x_2 + y_2 * xpix_;  weights[3] = (x - x_1) * (y - y_1);
 }
 
-std::vector<long>
+std::vector<size_t>
 FlatSkyProjection::QueryDisc(double alpha, double delta, double radius, bool local) const
 {
 	static const size_t npts = 72;
@@ -671,10 +671,10 @@ FlatSkyProjection::QueryDisc(double alpha, double delta, double radius, bool loc
 	double pvb = pv.R_component_3();
 	double pvc = pv.R_component_4();
 
-	long xmin = xpix_;
-	long xmax = 0;
-	long ymin = ypix_;
-	long ymax = 0;
+	ssize_t xmin = xpix_;
+	ssize_t xmax = 0;
+	ssize_t ymin = ypix_;
+	ssize_t ymax = 0;
 
 	double a, d;
 	double phi = 0;
@@ -687,10 +687,10 @@ FlatSkyProjection::QueryDisc(double alpha, double delta, double radius, bool loc
 		double s = sin(phi / 2.0);
 		quat q = quat(c, pva * s, pvb * s, pvc * s);
 		auto xy = QuatToXY(q * p / q, local);
-		long fx = std::floor(xy[0]);
-		long cx = std::ceil(xy[0]);
-		long fy = std::floor(xy[1]);
-		long cy = std::ceil(xy[1]);
+		ssize_t fx = std::floor(xy[0]);
+		ssize_t cx = std::ceil(xy[0]);
+		ssize_t fy = std::floor(xy[1]);
+		ssize_t cy = std::ceil(xy[1]);
 		if (fx < xmin)
 			xmin = fx < 0 ? 0 : fx;
 		if (cx > xmax)
@@ -701,10 +701,10 @@ FlatSkyProjection::QueryDisc(double alpha, double delta, double radius, bool loc
 			ymax = cy > ypix_ ? ypix_ : cy;
 	}
 
-	std::vector<long> pixels;
+	std::vector<size_t> pixels;
 	for (size_t x = xmin; x < xmax; x++) {
 		for (size_t y = ymin; y < ymax; y++) {
-			long pixel = y * xpix_ + x;
+			size_t pixel = y * xpix_ + x;
 			auto ang = PixelToAngle(pixel);
 			if (angular_distance(alpha, delta, ang[0], ang[1]) < radius)
 				pixels.push_back(pixel);
@@ -738,6 +738,13 @@ FlatSkyProjection FlatSkyProjection::Rebin(size_t scale, double x_center, double
 	return fp;
 }
 
+size_t FlatSkyProjection::RebinPixel(size_t pixel, size_t scale) const
+{
+	size_t x = (pixel % xpix_) / scale;
+	size_t y = ((size_t)(pixel / xpix_)) / scale;
+	return x + y * (xpix_ / scale);
+}
+
 FlatSkyProjection FlatSkyProjection::OverlayPatch(double x0, double y0,
     size_t width, size_t height) const
 {
diff --git a/maps/src/G3SkyMap.cxx b/maps/src/G3SkyMap.cxx
index 8309d401..4c105eaa 100644
--- a/maps/src/G3SkyMap.cxx
+++ b/maps/src/G3SkyMap.cxx
@@ -334,7 +334,7 @@ G3SkyMap &G3SkyMap::operator/=(double rhs)
 
 
 void
-G3SkyMap::GetInterpPixelsWeights(quat q, std::vector<long> & pixels,
+G3SkyMap::GetInterpPixelsWeights(quat q, std::vector<size_t> & pixels,
     std::vector<double> & weights) const
 {
 	double alpha, delta;
@@ -346,7 +346,7 @@ G3SkyMap::GetInterpPixelsWeights(quat q, std::vector<long> & pixels,
 }
 
 G3VectorQuat
-G3SkyMap::GetRebinQuats(long pixel, size_t scale) const
+G3SkyMap::GetRebinQuats(size_t pixel, size_t scale) const
 {
 	std::vector<double> alphas, deltas;
 
@@ -363,7 +363,7 @@ G3SkyMap::GetRebinQuats(long pixel, size_t scale) const
 }
 
 double
-G3SkyMap::GetInterpPrecalc(const std::vector<long> & pix,
+G3SkyMap::GetInterpPrecalc(const std::vector<size_t> & pix,
     const std::vector<double> & weight) const
 {
 	double outval = 0;
@@ -376,7 +376,7 @@ G3SkyMap::GetInterpPrecalc(const std::vector<long> & pix,
 double
 G3SkyMap::GetInterpValue(double alpha, double delta) const
 {
-	std::vector<long> pix;
+	std::vector<size_t> pix;
 	std::vector<double> weight;
 	GetInterpPixelsWeights(alpha, delta, pix, weight);
 	return GetInterpPrecalc(pix, weight);
@@ -417,7 +417,7 @@ G3SkyMap::GetInterpValues(const G3VectorQuat & quats) const
 	return outvals;
 }
 
-std::vector<long>
+std::vector<size_t>
 G3SkyMap::QueryAlphaEllipse(double alpha, double delta, double a, double b) const
 {
 	double rmaj = a > b ? a : b;
@@ -428,7 +428,7 @@ G3SkyMap::QueryAlphaEllipse(double alpha, double delta, double a, double b) cons
 
 	auto disc = QueryDisc(alpha, delta, rmaj);
 
-	std::vector<long> pixels;
+	std::vector<size_t> pixels;
 	for (auto i: disc) {
 		auto ang = PixelToAngle(i);
 		double d = angular_distance(ang[0], ang[1], ahi, delta) +
@@ -441,7 +441,7 @@ G3SkyMap::QueryAlphaEllipse(double alpha, double delta, double a, double b) cons
 }
 
 static double
-skymap_getitem(const G3SkyMap &skymap, int i)
+skymap_getitem(const G3SkyMap &skymap, ssize_t i)
 {
 
 	if (i < 0)
@@ -455,7 +455,7 @@ skymap_getitem(const G3SkyMap &skymap, int i)
 }
 
 static void
-skymap_setitem(G3SkyMap &skymap, int i, double val)
+skymap_setitem(G3SkyMap &skymap, ssize_t i, double val)
 {
 
 	if (i < 0)
diff --git a/maps/src/HealpixSkyMap.cxx b/maps/src/HealpixSkyMap.cxx
index 3d006770..515ec867 100644
--- a/maps/src/HealpixSkyMap.cxx
+++ b/maps/src/HealpixSkyMap.cxx
@@ -8,7 +8,6 @@
 
 #include <G3Units.h>
 #include <maps/HealpixSkyMap.h>
-#include <maps/chealpix.h>
 #include <maps/G3SkyMapMask.h>
 
 #include "mapdata.h"
@@ -79,10 +78,10 @@ HealpixSkyMap::HealpixSkyMap(boost::python::object v, bool weighted,
 			PyBuffer_Release(&view);
 			log_fatal("Only 1-D maps supported");
 		}
+		ssize_t npix = view.shape[0];
 		PyBuffer_Release(&view);
 
-		nside = npix2nside(view.shape[0]);
-		info_.initialize(nside, nested, shift_ra);
+		info_.initialize(npix, nested, shift_ra, true);
 
 		FillFromArray(v);
 
@@ -284,7 +283,7 @@ HealpixSkyMap::ConvertToDense()
 	dense_ = new std::vector<double>(size(), 0);
 
 	if (ring_sparse_) {
-		long idx;
+		size_t idx;
 		for (auto i = ring_sparse_->begin(); i != ring_sparse_->end();
 		    i++) {
 			idx = info_.RingToPixel(i.x, i.y);
@@ -341,7 +340,7 @@ HealpixSkyMap::ConvertToIndexedSparse()
 	indexed_sparse_ = new std::unordered_map<uint64_t, double>;
 
 	if (ring_sparse_) {
-		long idx;
+		size_t idx;
 		for (auto i = ring_sparse_->begin(); i != ring_sparse_->end();
 		    i++) {
 			if ((*i) == 0)
@@ -728,7 +727,6 @@ G3SkyMap &HealpixSkyMap::operator /=(const G3SkyMap &rhs) {
 			} else
 				zero = true;
 		} else if (ring_sparse_) {
-			long i = 0;
 			if (b.dense_ || b.ring_sparse_ || b.indexed_sparse_) {
 				for (size_t i = 0; i < size(); i++) {
 					double valb = b.at(i);
@@ -1009,7 +1007,7 @@ HealpixSkyMap::PixelToAngle(size_t pixel) const
 	return info_.PixelToAngle(pixel);
 }
 
-void HealpixSkyMap::GetRebinAngles(long pixel, size_t scale,
+void HealpixSkyMap::GetRebinAngles(size_t pixel, size_t scale,
     std::vector<double> & alphas, std::vector<double> & deltas) const
 {
 	info_.GetRebinAngles(pixel, scale, alphas, deltas);
@@ -1017,12 +1015,12 @@ void HealpixSkyMap::GetRebinAngles(long pixel, size_t scale,
 
 void
 HealpixSkyMap::GetInterpPixelsWeights(double alpha, double delta,
-    std::vector<long> & pixels, std::vector<double> & weights) const
+    std::vector<size_t> & pixels, std::vector<double> & weights) const
 {
 	info_.GetInterpPixelsWeights(alpha, delta, pixels, weights);
 }
 
-std::vector<long>
+std::vector<size_t>
 HealpixSkyMap::QueryDisc(double alpha, double delta, double radius) const
 {
 	return info_.QueryDisc(alpha, delta, radius);
@@ -1048,18 +1046,12 @@ G3SkyMapPtr HealpixSkyMap::Rebin(size_t scale, bool norm) const
 	else
 		return out;
 
-	const size_t scale2 = scale * scale;
-	const double sqscal = norm ? scale2 : 1.0;
+	const double sqscal = norm ? scale * scale : 1.0;
 
 	for (auto i : *this) {
 		if (i.second == 0)
 			continue;
-		long ip = i.first;
-		if (!nested())
-			ring2nest(nside(), ip, &ip);
-		ip /= scale2;
-		if (!nested())
-			nest2ring(out->nside(), ip, &ip);
+		size_t ip = info_.RebinPixel(i.first, scale);
 		(*out)[ip] += i.second / sqscal;
 	}
 
@@ -1149,7 +1141,7 @@ HealpixSkyMap_nonzeropixels(const HealpixSkyMap &m)
 }
 
 static double
-skymap_getitem(const G3SkyMap &skymap, int i)
+skymap_getitem(const G3SkyMap &skymap, ssize_t i)
 {
 
 	if (i < 0)
@@ -1163,7 +1155,7 @@ skymap_getitem(const G3SkyMap &skymap, int i)
 }
 
 static void
-skymap_setitem(G3SkyMap &skymap, int i, double val)
+skymap_setitem(G3SkyMap &skymap, ssize_t i, double val)
 {
 
 	if (i < 0)
@@ -1262,7 +1254,7 @@ static PyBufferProcs healpixskymap_bufferprocs;
 
 
 static void
-HealpixSkyMap_setitem_1d(G3SkyMap &skymap, size_t i, double val)
+HealpixSkyMap_setitem_1d(G3SkyMap &skymap, ssize_t i, double val)
 {
 
 	if (i < 0)
diff --git a/maps/src/HealpixSkyMapInfo.cxx b/maps/src/HealpixSkyMapInfo.cxx
index 3059b553..471bbdfb 100644
--- a/maps/src/HealpixSkyMapInfo.cxx
+++ b/maps/src/HealpixSkyMapInfo.cxx
@@ -6,7 +6,7 @@
 #include <G3Units.h>
 
 #include <maps/HealpixSkyMapInfo.h>
-#include <maps/chealpix.h>
+#include "chealpix.h"
 
 static const double twothird = 2.0 / 3.0;
 static const double twopi = 2 * M_PI;
@@ -27,9 +27,13 @@ HealpixSkyMapInfo::HealpixSkyMapInfo(const HealpixSkyMapInfo & other)
 }
 
 void
-HealpixSkyMapInfo::initialize(size_t nside, bool nested, bool shifted)
+HealpixSkyMapInfo::initialize(size_t nside_or_npix, bool nested, bool shifted,
+   bool is_npix)
 {
-	SetNSide(nside);
+	if (is_npix)
+		SetNPix(nside_or_npix);
+	else
+		SetNSide(nside_or_npix);
 	SetNested(nested);
 	SetShifted(shifted);
 }
@@ -86,6 +90,16 @@ HealpixSkyMapInfo::res() const
 	return sqrt(4.0 * M_PI / npix_) * G3Units::rad;
 }
 
+void
+HealpixSkyMapInfo::SetNPix(size_t npix)
+{
+	ssize_t nside = npix2nside64(npix);
+	if (nside < 0)
+		log_fatal("Invalid npix %zu", npix);
+
+	SetNSide(nside);
+}
+
 void
 HealpixSkyMapInfo::SetNSide(size_t nside)
 {
@@ -152,12 +166,12 @@ HealpixSkyMapInfo::PixelToRing(size_t pixel) const
 {
 	static const std::pair<size_t, size_t> bad = {(size_t) -1, (size_t)-1};
 
-	if (pixel < 0 || pixel >= npix_)
+	if (pixel >= npix_)
 		return bad;
 
 	if (nested_) {
-		long pix = pixel;
-		nest2ring(nside_, pix, &pix);
+		ssize_t pix = pixel;
+		nest2ring64(nside_, pix, &pix);
 		pixel = pix;
 	}
 
@@ -168,13 +182,13 @@ HealpixSkyMapInfo::PixelToRing(size_t pixel) const
 		iring = (size_t)((pixel - ncap_) / nring_ + nside_);
 	else /* South Polar cap */
 		iring = nring_ - (size_t)(0.5 * (1 + sqrt(2 * (npix_ - pixel) - 0.5)));
-	if (iring < 0 || iring >= nring_)
+	if (iring >= nring_)
 		return bad;
 
 	const HealpixRingInfo & ring = rings_[iring];
 
 	size_t ringpix = pixel - ring.pix0;
-	if (ringpix < 0 || ringpix >= ring.npix)
+	if (ringpix >= ring.npix)
 		return bad;
 
 	if (shifted_)
@@ -186,21 +200,24 @@ HealpixSkyMapInfo::PixelToRing(size_t pixel) const
 size_t
 HealpixSkyMapInfo::RingToPixel(size_t iring, size_t ringpix) const
 {
-	if (iring < 0 || iring >= nring_)
+	if (iring >= nring_)
 		return -1;
 	const HealpixRingInfo & ring = rings_[iring];
 
 	if (shifted_)
 		ringpix = (ringpix + ring.npix / 2) % ring.npix;
-	if (ringpix < 0 || ringpix >= ring.npix)
+	if (ringpix >= ring.npix)
 		return -1;
 
-	long pixel = ring.pix0 + ringpix;
-	if (pixel < 0 || pixel >= npix_)
+	size_t pixel = ring.pix0 + ringpix;
+	if (pixel >= npix_)
 		return -1;
 
-	if (nested_)
-		ring2nest(nside_, pixel, &pixel);
+	if (nested_) {
+		ssize_t pix = pixel;
+		ring2nest64(nside_, pix, &pix);
+		pixel = pix;
+	}
 
 	return pixel;
 }
@@ -219,11 +236,11 @@ HealpixSkyMapInfo::AngleToPixel(double alpha, double delta) const
 	while (alpha < 0)
 		alpha += twopi;
 
-	long outpix;
+	ssize_t outpix;
 	if (nested_)
-		ang2pix_nest(nside_, theta, alpha, &outpix);
+		ang2pix_nest64(nside_, theta, alpha, &outpix);
 	else
-		ang2pix_ring(nside_, theta, alpha, &outpix);
+		ang2pix_ring64(nside_, theta, alpha, &outpix);
 
 	if (outpix < 0 || outpix >= npix_)
 		return -1;
@@ -239,9 +256,9 @@ HealpixSkyMapInfo::PixelToAngle(size_t pixel) const
 
 	double alpha, delta;
 	if (nested_)
-		pix2ang_nest(nside_, pixel, &delta, &alpha);
+		pix2ang_nest64(nside_, pixel, &delta, &alpha);
 	else
-		pix2ang_ring(nside_, pixel, &delta, &alpha);
+		pix2ang_ring64(nside_, pixel, &delta, &alpha);
 
 	while(alpha > M_PI)
 		alpha -= twopi;
@@ -255,26 +272,43 @@ HealpixSkyMapInfo::PixelToAngle(size_t pixel) const
 	return {alpha, delta};
 }
 
+size_t
+HealpixSkyMapInfo::RebinPixel(size_t pixel, size_t scale) const
+{
+	ssize_t pix = pixel;
+
+	if (!nested_)
+		ring2nest64(nside_, pix, &pix);
+	pix /= scale * scale;
+	if (!nested_)
+		nest2ring64(nside_ / scale, pix, &pix);
+
+	return pix;
+}
+
 void
-HealpixSkyMapInfo::GetRebinAngles(long pixel, size_t scale,
+HealpixSkyMapInfo::GetRebinAngles(size_t pixel, size_t scale,
     std::vector<double> & alphas, std::vector<double> & deltas) const
 {
 	if (nside_ % scale != 0)
 		log_fatal("Nside must be a multiple of rebinning scale");
 
-	if (!nested_)
-		ring2nest(nside_, pixel, &pixel);
+	if (!nested_) {
+		ssize_t pix = pixel;
+		ring2nest64(nside_, pix, &pix);
+		pixel = pix;
+	}
 
 	alphas = std::vector<double>(scale * scale);
 	deltas = std::vector<double>(scale * scale);
 
 	size_t nside_rebin = nside_ * scale;
-	long pixmin = pixel * scale * scale;
+	size_t pixmin = pixel * scale * scale;
 
 	for (size_t i = 0; i < (scale * scale); i++) {
-		long p = pixmin + i;
+		size_t p = pixmin + i;
 		double theta, phi;
-		pix2ang_nest(nside_rebin, p, &theta, &phi);
+		pix2ang_nest64(nside_rebin, p, &theta, &phi);
 
 		while (phi > M_PI)
 			phi -= twopi;
@@ -296,12 +330,12 @@ HealpixSkyMapInfo::RingAbove(double z) const
 
 void
 HealpixSkyMapInfo::GetInterpPixelsWeights(double alpha, double delta,
-    std::vector<long> & pixels, std::vector<double> & weights) const
+    std::vector<size_t> & pixels, std::vector<double> & weights) const
 {
 	alpha /= G3Units::rad;
 	delta /= G3Units::rad;
 
-	pixels = std::vector<long>(4, -1);
+	pixels = std::vector<size_t>(4, (size_t) -1);
 	weights = std::vector<double>(4, 0);
 
 	double theta = M_PI_2 - delta;
@@ -319,11 +353,11 @@ HealpixSkyMapInfo::GetInterpPixelsWeights(double alpha, double delta,
 		const HealpixRingInfo & ring = rings_[ir1];
 		z1 = ring.z;
 		double tmp = phi / ring.dphi - ring.shift;
-		long i1 = (tmp < 0) ? tmp - 1 : tmp;
+		ssize_t i1 = (tmp < 0) ? tmp - 1 : tmp;
 		double w1 = (phi - (i1 + ring.shift) * ring.dphi) / ring.dphi;
 		if (i1 < 0)
 			i1 += ring.npix;
-		long i2 = i1 + 1;
+		ssize_t i2 = i1 + 1;
 		if (i2 >= ring.npix)
 			i2 -= ring.npix;
 		pixels[0] = ring.pix0 + i1;
@@ -336,11 +370,11 @@ HealpixSkyMapInfo::GetInterpPixelsWeights(double alpha, double delta,
 		const HealpixRingInfo & ring = rings_[ir2];
 		z2 = ring.z;
 		double tmp = phi / ring.dphi - ring.shift;
-		long i1 = (tmp < 0) ? tmp - 1 : tmp;
+		ssize_t i1 = (tmp < 0) ? tmp - 1 : tmp;
 		double w1 = (phi - (i1 + ring.shift) * ring.dphi) / ring.dphi;
 		if (i1 < 0)
 			i1 += ring.npix;
-		long i2 = i1 + 1;
+		ssize_t i2 = i1 + 1;
 		if (i2 >= ring.npix)
 			i2 -= ring.npix;
 		pixels[2] = ring.pix0 + i1;
@@ -380,21 +414,22 @@ HealpixSkyMapInfo::GetInterpPixelsWeights(double alpha, double delta,
 	}
 
 	if (nested_) {
-		ring2nest(nside_, pixels[0], &pixels[0]);
-		ring2nest(nside_, pixels[1], &pixels[1]);
-		ring2nest(nside_, pixels[2], &pixels[2]);
-		ring2nest(nside_, pixels[3], &pixels[3]);
+		for (size_t i = 0; i < pixels.size(); i++) {
+			ssize_t pix = pixels[i];
+			ring2nest64(nside_, pix, &pix);
+			pixels[i] = pix;
+		}
 	}
 }
 
-std::vector<long>
+std::vector<size_t>
 HealpixSkyMapInfo::QueryDisc(double alpha, double delta, double radius) const
 {
-	auto pixels = std::vector<long>();
+	auto pixels = std::vector<size_t>();
 
 	radius /= G3Units::rad;
 	if (radius >= M_PI) {
-		for (long i = 0; i < npix_; i++)
+		for (ssize_t i = 0; i < npix_; i++)
 			pixels.push_back(i);
 		return pixels;
 	}
@@ -416,7 +451,7 @@ HealpixSkyMapInfo::QueryDisc(double alpha, double delta, double radius) const
 	if ((rlat1 <= 0) && (irmin > 1)) {
 		// north pole in the disk
 		const HealpixRingInfo & ring = rings_[irmin - 1];
-		for (long i = 0; i < ring.pix0 + ring.npix; i++)
+		for (ssize_t i = 0; i < ring.pix0 + ring.npix; i++)
 			pixels.push_back(i);
 	}
 
@@ -435,10 +470,10 @@ HealpixSkyMapInfo::QueryDisc(double alpha, double delta, double radius) const
 		double dphi = atan2(sqrt(ysq), x);
 
 		// highest pixel number in the ring
-		long ipix2 = ring.pix0 + ring.npix - 1;
+		ssize_t ipix2 = ring.pix0 + ring.npix - 1;
 
-		long ip_lo = (long)floor((phi - dphi) / ring.dphi - ring.shift) + 1;
-		long ip_hi = (long)floor((phi + dphi) / ring.dphi - ring.shift);
+		ssize_t ip_lo = (ssize_t)floor((phi - dphi) / ring.dphi - ring.shift) + 1;
+		ssize_t ip_hi = (ssize_t)floor((phi + dphi) / ring.dphi - ring.shift);
 		if (ip_lo > ip_hi)
 			continue;
 
@@ -447,12 +482,12 @@ HealpixSkyMapInfo::QueryDisc(double alpha, double delta, double radius) const
 			ip_hi -= ring.npix;
 		}
 		if (ip_lo < 0) {
-			for (long i = ring.pix0; i < ring.pix0 + ip_hi + 1; i++)
+			for (ssize_t i = ring.pix0; i < ring.pix0 + ip_hi + 1; i++)
 				pixels.push_back(i);
-			for (long i = ring.pix0 + ip_lo + ring.npix; i < ipix2 + 1; i++)
+			for (ssize_t i = ring.pix0 + ip_lo + ring.npix; i < ipix2 + 1; i++)
 				pixels.push_back(i);
 		} else {
-			for (long i = ring.pix0 + ip_lo; i < ring.pix0 + ip_hi + 1; i++)
+			for (ssize_t i = ring.pix0 + ip_lo; i < ring.pix0 + ip_hi + 1; i++)
 				pixels.push_back(i);
 		}
 	}
@@ -460,13 +495,16 @@ HealpixSkyMapInfo::QueryDisc(double alpha, double delta, double radius) const
 	if ((rlat2 >= M_PI) && (irmax + 1 < nring_)) {
 		// south pole in the disk
 		const HealpixRingInfo & ring = rings_[irmax + 1];
-		for (long i = ring.pix0; i < npix_; i++)
+		for (ssize_t i = ring.pix0; i < npix_; i++)
 			pixels.push_back(i);
 	}
 
 	if (nested_) {
-		for (size_t i = 0; i < pixels.size(); i++)
-			ring2nest(nside_, pixels[i], &pixels[i]);
+		for (size_t i = 0; i < pixels.size(); i++) {
+			ssize_t pix = pixels[i];
+			ring2nest64(nside_, pix, &pix);
+			pixels[i] = pix;
+		}
 		std::sort(pixels.begin(), pixels.end());
 	}
 
diff --git a/maps/src/MapMockObserver.cxx b/maps/src/MapMockObserver.cxx
index 38b0774b..ba81012d 100644
--- a/maps/src/MapMockObserver.cxx
+++ b/maps/src/MapMockObserver.cxx
@@ -215,7 +215,7 @@ MapMockObserver::Process(G3FramePtr frame, std::deque<G3FramePtr> &out)
 			    pcoupling, U_->GetPolConv());
 			for (size_t i = 0; i < det.size(); i++) {
 				if (interp_) {
-					std::vector<long> pixels;
+					std::vector<size_t> pixels;
 					std::vector<double> weights;
 					T_->GetInterpPixelsWeights(alpha[i], delta[i], pixels, weights);
 					det[i] =
diff --git a/maps/src/SingleDetectorBoresightBinner.cxx b/maps/src/SingleDetectorBoresightBinner.cxx
index dbf81b61..6c9b5d4f 100644
--- a/maps/src/SingleDetectorBoresightBinner.cxx
+++ b/maps/src/SingleDetectorBoresightBinner.cxx
@@ -184,7 +184,7 @@ SingleDetectorBoresightBinner::Process(G3FramePtr frame,
 	
 	// Conjugate pointing rotation with boresight vector
 	auto pointing_vec = (*pointing)*quat(0,1,0,0)/(*pointing);
-	std::vector<size_t> pixels = template_->QuatsToPixels(pointing_vec);
+	auto pixels = template_->QuatsToPixels(pointing_vec);
 
 	for (size_t i = 0; i < pixels.size(); i++)
 		(*map_weights_->TT)[pixels[i]] += 1;
diff --git a/maps/src/chealpix.c b/maps/src/chealpix.c
index 05bd9b82..12baac79 100644
--- a/maps/src/chealpix.c
+++ b/maps/src/chealpix.c
@@ -30,7 +30,7 @@
 #include <stdlib.h>
 #include <string.h>
 #include <math.h>
-#include <maps/chealpix.h>
+#include "chealpix.h"
 
 static const double twothird=2.0/3.0;
 static const double pi=3.141592653589793238462643383279502884197;
diff --git a/maps/include/maps/chealpix.h b/maps/src/chealpix.h
similarity index 97%
rename from maps/include/maps/chealpix.h
rename to maps/src/chealpix.h
index 3aefb07c..4b835e70 100644
--- a/maps/include/maps/chealpix.h
+++ b/maps/src/chealpix.h
@@ -106,9 +106,8 @@ void pix2vec_ring(long nside, long ipix, double *vec);
 /* operations on Nside values up to 2^29 */
 
 /*! 64bit integer type
-    \note We are not using \c int64_t, since this type is not part of the C++
-    standard, and we want the header to be usable from C++. */
-typedef long long hpint64;
+    \note We use ssize_t here since this is properly defined on all platforms. */
+typedef ssize_t hpint64;
 
 /*! Sets \a *ipix to the pixel number in NEST scheme at resolution \a nside,
     which contains the position \a theta, \a phi. */


Test patch:
diff --git a/maps/tests/healpix_maps.py b/maps/tests/healpix_maps.py
index 3857c4f7..0e141820 100755
--- a/maps/tests/healpix_maps.py
+++ b/maps/tests/healpix_maps.py
@@ -221,3 +221,9 @@
 lon, lat = maps.azel.convert_radec_to_gal(alpha, delta)
 pix = xeq.angles_to_pixels(lon, lat)
 assert(np.allclose(np.asarray(pix), np.asarray(xeq)))
+
+# high-res
+hx = maps.HealpixSkyMap(16384)
+pix1 = hx.angle_to_pixel(0, -45 * core.G3Units.deg)
+pix2 = hp.ang2pix(16384, 3 * np.pi / 4, 0)
+assert(pix1 == pix2)