diff --git a/libsharp/sharp_core.c b/libsharp/sharp_core.c
index 00ff6b3..b9489a9 100644
--- a/libsharp/sharp_core.c
+++ b/libsharp/sharp_core.c
@@ -219,15 +219,11 @@ NOINLINE static void alm2map_kernel(s0data_v * restrict d,
     Tv a2=vload(ab[il+1].f[0]), b2=vload(ab[il+1].f[1]);
     for (int i=0; i<nv2; ++i)
       {
-      d->p1r[i] += d->lam2[i]*ar1;
-      d->p1i[i] += d->lam2[i]*ai1;
-      d->p2r[i] += d->lam2[i]*ar2;
-      d->p2i[i] += d->lam2[i]*ai2;
       d->lam1[i] = (a1*d->csq[i] + b1)*d->lam2[i] + d->lam1[i];
-      d->p1r[i] += d->lam1[i]*ar3;
-      d->p1i[i] += d->lam1[i]*ai3;
-      d->p2r[i] += d->lam1[i]*ar4;
-      d->p2i[i] += d->lam1[i]*ai4;
+      d->p1r[i] += d->lam2[i]*ar1 + d->lam1[i]*ar3;
+      d->p1i[i] += d->lam2[i]*ai1 + d->lam1[i]*ai3;
+      d->p2r[i] += d->lam2[i]*ar2 + d->lam1[i]*ar4;
+      d->p2i[i] += d->lam2[i]*ai2 + d->lam1[i]*ai4;
       d->lam2[i] = (a2*d->csq[i] + b2)*d->lam1[i] + d->lam2[i];
       }
     }
diff --git a/libsharp/sharp_testsuite.c b/libsharp/sharp_testsuite.c
index f712ae4..fdc842a 100644
--- a/libsharp/sharp_testsuite.c
+++ b/libsharp/sharp_testsuite.c
@@ -241,13 +241,14 @@ static int good_fft_size(int n)
   }
 
 static void get_infos (const char *gname, int lmax, int *mmax, int *gpar1,
-  int *gpar2, sharp_geom_info **ginfo, sharp_alm_info **ainfo)
+  int *gpar2, sharp_geom_info **ginfo, sharp_alm_info **ainfo, int verbose)
   {
   UTIL_ASSERT(lmax>=0,"lmax must not be negative");
   if (*mmax<0) *mmax=lmax;
   UTIL_ASSERT(*mmax<=lmax,"mmax larger than lmax");
 
-  if (mytask==0) printf ("lmax: %d, mmax: %d\n",lmax,*mmax);
+  verbose &= (mytask==0);
+  if (verbose) printf ("lmax: %d, mmax: %d\n",lmax,*mmax);
 
   sharp_make_triangular_alm_info(lmax,*mmax,1,ainfo);
 #ifdef USE_MPI
@@ -259,14 +260,14 @@ static void get_infos (const char *gname, int lmax, int *mmax, int *gpar1,
     if (*gpar1<1) *gpar1=lmax/2;
     if (*gpar1==0) ++(*gpar1);
     sharp_make_healpix_geom_info (*gpar1, 1, ginfo);
-    if (mytask==0) printf ("HEALPix grid, nside=%d\n",*gpar1);
+    if (verbose) printf ("HEALPix grid, nside=%d\n",*gpar1);
     }
   else if (strcmp(gname,"gauss")==0)
     {
     if (*gpar1<1) *gpar1=lmax+1;
     if (*gpar2<1) *gpar2=2*(*mmax)+1;
     sharp_make_gauss_geom_info (*gpar1, *gpar2, 0., 1, *gpar2, ginfo);
-    if (mytask==0)
+    if (verbose)
       printf ("Gauss-Legendre grid, nlat=%d, nlon=%d\n",*gpar1,*gpar2);
     }
   else if (strcmp(gname,"fejer1")==0)
@@ -274,21 +275,21 @@ static void get_infos (const char *gname, int lmax, int *mmax, int *gpar1,
     if (*gpar1<1) *gpar1=2*lmax+1;
     if (*gpar2<1) *gpar2=2*(*mmax)+1;
     sharp_make_fejer1_geom_info (*gpar1, *gpar2, 0., 1, *gpar2, ginfo);
-    if (mytask==0) printf ("Fejer1 grid, nlat=%d, nlon=%d\n",*gpar1,*gpar2);
+    if (verbose) printf ("Fejer1 grid, nlat=%d, nlon=%d\n",*gpar1,*gpar2);
     }
   else if (strcmp(gname,"fejer2")==0)
     {
     if (*gpar1<1) *gpar1=2*lmax+1;
     if (*gpar2<1) *gpar2=2*(*mmax)+1;
     sharp_make_fejer2_geom_info (*gpar1, *gpar2, 0., 1, *gpar2, ginfo);
-    if (mytask==0) printf ("Fejer2 grid, nlat=%d, nlon=%d\n",*gpar1,*gpar2);
+    if (verbose) printf ("Fejer2 grid, nlat=%d, nlon=%d\n",*gpar1,*gpar2);
     }
   else if (strcmp(gname,"cc")==0)
     {
     if (*gpar1<1) *gpar1=2*lmax+1;
     if (*gpar2<1) *gpar2=2*(*mmax)+1;
     sharp_make_cc_geom_info (*gpar1, *gpar2, 0., 1, *gpar2, ginfo);
-    if (mytask==0)
+    if (verbose)
       printf("Clenshaw-Curtis grid, nlat=%d, nlon=%d\n",*gpar1,*gpar2);
     }
   else if (strcmp(gname,"smallgauss")==0)
@@ -318,7 +319,7 @@ static void get_infos (const char *gname, int lmax, int *mmax, int *gpar1,
         ofs+=pring;
         }
       }
-    if (mytask==0)
+    if (verbose)
       {
       ptrdiff_t npix=get_npix(*ginfo);
       printf("Small Gauss grid, nlat=%d, npix=%ld, savings=%.2f%%\n",
@@ -485,6 +486,16 @@ static void check_accuracy (sharp_geom_info *ginfo, sharp_alm_info *ainfo,
   DEALLOC(err_abs);
   }
 
+static void run(int lmax, int mmax, int nlat, int nlon, int spin)
+  {
+  sharp_geom_info *ginfo;
+  sharp_alm_info *ainfo;
+  get_infos ("gauss", lmax, &mmax, &nlat, &nlon, &ginfo, &ainfo, 0);
+  check_accuracy(ginfo,ainfo,spin);
+  sharp_destroy_alm_info(ainfo);
+  sharp_destroy_geom_info(ginfo);
+  }
+
 static void sharp_acctest(void)
   {
   if (mytask==0) sharp_module_startup("sharp_acctest",1,1,"",1);
@@ -495,17 +506,15 @@ static void sharp_acctest(void)
 
   if (mytask==0) printf("Testing map analysis accuracy.\n");
 
-  sharp_geom_info *ginfo;
-  sharp_alm_info *ainfo;
-  int lmax=127, mmax=127, nlat=128, nlon=256;
-  get_infos ("gauss", lmax, &mmax, &nlat, &nlon, &ginfo, &ainfo);
-  check_accuracy(ginfo,ainfo,0);
-  check_accuracy(ginfo,ainfo,1);
-  check_accuracy(ginfo,ainfo,2);
-  check_accuracy(ginfo,ainfo,3);
-  check_accuracy(ginfo,ainfo,30);
-  sharp_destroy_alm_info(ainfo);
-  sharp_destroy_geom_info(ginfo);
+  run(127, 127, 128, 256, 0);
+  run(127, 127, 128, 256, 1);
+  run(127, 127, 128, 256, 2);
+  run(127, 127, 128, 256, 3);
+  run(127, 127, 128, 256, 30);
+  run(5, 0, 6, 1, 0);
+  run(5, 0, 7, 2, 0);
+  run(8, 8, 9, 17, 0);
+  run(8, 8, 9, 17, 2);
   if (mytask==0) printf("Passed.\n\n");
   }
 
@@ -524,7 +533,7 @@ static void sharp_test (int argc, const char **argv)
 
   sharp_geom_info *ginfo;
   sharp_alm_info *ainfo;
-  get_infos (argv[2], lmax, &mmax, &gpar1, &gpar2, &ginfo, &ainfo);
+  get_infos (argv[2], lmax, &mmax, &gpar1, &gpar2, &ginfo, &ainfo, 1);
 
   int ncomp = (spin==0) ? 1 : 2;
   double t_a2m=1e30, t_m2a=1e30;
diff --git a/libsharp/sharp_ylmgen_c.c b/libsharp/sharp_ylmgen_c.c
index f0672b0..f9de49c 100644
--- a/libsharp/sharp_ylmgen_c.c
+++ b/libsharp/sharp_ylmgen_c.c
@@ -73,16 +73,16 @@ void sharp_Ylmgen_init (sharp_Ylmgen_C *gen, int l_max, int m_max, int spin)
     gen->mfac[0] = inv_sqrt4pi;
     for (int m=1; m<=gen->mmax; ++m)
       gen->mfac[m] = gen->mfac[m-1]*sqrt((2*m+1.)/(2*m));
-    gen->root = RALLOC(double,2*gen->lmax+7);
-    gen->iroot = RALLOC(double,2*gen->lmax+7);
-    for (int m=0; m<2*gen->lmax+7; ++m)
+    gen->root = RALLOC(double,2*gen->lmax+8);
+    gen->iroot = RALLOC(double,2*gen->lmax+8);
+    for (int m=0; m<2*gen->lmax+8; ++m)
       {
       gen->root[m] = sqrt(m);
       gen->iroot[m] = (m==0) ? 0. : 1./gen->root[m];
       }
-    gen->eps=RALLOC(double, gen->lmax+10);
-    gen->alpha=RALLOC(double, gen->lmax/2+10);
-    gen->ab=RALLOC(sharp_ylmgen_dbl2, gen->lmax/2+10);
+    gen->eps=RALLOC(double, gen->lmax+4);
+    gen->alpha=RALLOC(double, gen->lmax/2+2);
+    gen->ab=RALLOC(sharp_ylmgen_dbl2, gen->lmax/2+2);
     }
   else
     {
@@ -161,15 +161,16 @@ void sharp_Ylmgen_prepare (sharp_Ylmgen_C *gen, int m)
 
   if (gen->s==0)
     {
-    for (int l=m; l<gen->lmax+10; ++l)
+    gen->eps[m] = 0.;
+    for (int l=m+1; l<gen->lmax+4; ++l)
       gen->eps[l] = gen->root[l+m]*gen->root[l-m]
                    *gen->iroot[2*l+1]*gen->iroot[2*l-1];
     gen->alpha[0] = 1./gen->eps[m+1];
     gen->alpha[1] = gen->eps[m+1]/(gen->eps[m+2]*gen->eps[m+3]);
-    for (int il=1, l=m+2; l<gen->lmax+5; ++il, l+=2)
+    for (int il=1, l=m+2; l<gen->lmax+1; ++il, l+=2)
       gen->alpha[il+1]= ((il&1) ? -1 : 1)
                        /(gen->eps[l+2]*gen->eps[l+3]*gen->alpha[il]);
-    for (int il=0, l=m; l<gen->lmax+5; ++il, l+=2)
+    for (int il=0, l=m; l<gen->lmax+2; ++il, l+=2)
       {
       gen->ab[il].f[0] = ((il&1) ? -1 : 1)*gen->alpha[il]*gen->alpha[il];
       double t1 = gen->eps[l+2], t2 = gen->eps[l+1];