diff --git a/libsharp/sharp_core.c b/libsharp/sharp_core.c
index 2dcee90..690bf2a 100644
--- a/libsharp/sharp_core.c
+++ b/libsharp/sharp_core.c
@@ -280,6 +280,58 @@ NOINLINE static void calc_alm2map (sharp_job * restrict job,
   alm2map_kernel(d, rf, alm, l, lmax, nv2);
   }
 
+NOINLINE static void calc_alm2map_alt (sharp_job * restrict job,
+  const sharp_Ylmgen_C * restrict gen, s0data_v * restrict d, int nth)
+  {
+  int l=gen->m, lmax=gen->lmax;
+  int nv2 = (nth+VLEN-1)/VLEN;
+  job->opcnt += (lmax+1-l) * 6*nth;
+
+const double inv_sqrt4pi = 0.2820947917738781434740397257803862929220;
+  Tv mfac = vload(gen->mfac[gen->m]);
+  for (int i=0; i<nv2; ++i)
+    {
+    d->lam1[i]=vzero;
+    mypow(d->sth[i],l,gen->powlimit,&d->lam2[i],&d->scale[i]);
+    d->lam2[i] *= mfac;
+    Tvnormalize(&d->lam2[i],&d->scale[i],sharp_ftol);
+    }
+
+  const dcmplx * restrict alm=job->almtmp;
+  dcmplx * restrict alm2=RALLOC(dcmplx, gen->lmax+5);
+  {
+  for (int il=0, l=gen->m; l<=gen->lmax; ++il,l+=2)
+    {
+    dcmplx al = alm[l];
+    dcmplx al1 = (l+1>gen->lmax) ? 0. : alm[l+1];
+    dcmplx al2 = (l+2>gen->lmax) ? 0. : alm[l+2];
+    alm2[l  ] = gen->alpha[il]*(gen->eps[l+1]*al + gen->eps[l+2]*al2);
+    alm2[l+1] = gen->alpha[il]*al1;
+    }
+  }
+  for (int i=0; i<nv2; ++i)
+    getCorfac(d->scale[i], &d->corfac[i], gen->cf);
+
+  for (int il=0, l=gen->m; l<=lmax; ++il, l+=2)
+    {
+    Tv ar1=vload(creal(alm2[l  ])), ai1=vload(cimag(alm2[l  ]));
+    Tv ar2=vload(creal(alm2[l+1])), ai2=vload(cimag(alm2[l+1]));
+    for (int i=0; i<nv2; ++i)
+      {
+      d->p1r[i] += d->lam2[i]*d->corfac[i]*ar1;
+      d->p1i[i] += d->lam2[i]*d->corfac[i]*ai1;
+      d->p2r[i] += d->cth[i]*d->lam2[i]*d->corfac[i]*ar2;
+      d->p2i[i] += d->cth[i]*d->lam2[i]*d->corfac[i]*ai2;
+      Tv tmp = (gen->a[il]*d->cth[i]*d->cth[i] + gen->b[il])*d->lam2[i] + d->lam1[i];
+      d->lam1[i] = d->lam2[i];
+      d->lam2[i] = tmp;
+      if (rescale(&d->lam1[i], &d->lam2[i], &d->scale[i], vload(sharp_ftol)))
+        getCorfac(d->scale[i], &d->corfac[i], gen->cf);
+      }
+    }
+  DEALLOC(alm2);
+  }
+
 NOINLINE static void map2alm_kernel(s0data_v * restrict d,
   const sharp_ylmgen_dbl2 * restrict rf, dcmplx * restrict alm, int l,
   int lmax, int nv2)
@@ -812,7 +864,7 @@ NOINLINE static void inner_loop_a2m(sharp_job *job, const int *ispair,
               d.s.sth[i]=d.s.sth[nth-1];
               d.s.p1r[i]=d.s.p1i[i]=d.s.p2r[i]=d.s.p2i[i]=0.;
               }
-            calc_alm2map (job, gen, &d.v, nth);
+            calc_alm2map_alt (job, gen, &d.v, nth);
             for (int i=0; i<nth; ++i)
               {
               int tgt=itgt[i];
diff --git a/libsharp/sharp_ylmgen_c.c b/libsharp/sharp_ylmgen_c.c
index 9cb9dbb..e44d88b 100644
--- a/libsharp/sharp_ylmgen_c.c
+++ b/libsharp/sharp_ylmgen_c.c
@@ -74,13 +74,17 @@ 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+5);
-    gen->iroot = RALLOC(double,2*gen->lmax+5);
-    for (int m=0; m<2*gen->lmax+5; ++m)
+    gen->root = RALLOC(double,2*gen->lmax+6);
+    gen->iroot = RALLOC(double,2*gen->lmax+6);
+    for (int m=0; m<2*gen->lmax+6; ++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->a=RALLOC(double, gen->lmax/2+10);
+gen->b=RALLOC(double, gen->lmax/2+10);
     }
   else
     {
@@ -137,6 +141,10 @@ void sharp_Ylmgen_destroy (sharp_Ylmgen_C *gen)
     DEALLOC(gen->mfac);
     DEALLOC(gen->root);
     DEALLOC(gen->iroot);
+DEALLOC(gen->eps);
+DEALLOC(gen->alpha);
+DEALLOC(gen->a);
+DEALLOC(gen->b);
     }
   else
     {
@@ -165,6 +173,18 @@ void sharp_Ylmgen_prepare (sharp_Ylmgen_C *gen, int m)
       gen->rf[l].f[0] = tmp*gen->root[2*l+1];
       gen->rf[l].f[1] = tmp*gen->root[l+m]*gen->root[l-m]*gen->iroot[2*l-1];
       }
+for (int l=m; l<gen->lmax+10; ++l)
+  gen->eps[l] = sqrt((l*l-m*m)/(4.*l*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)
+  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)
+  {
+  gen->a[il] = ((il&1) ? -1 : 1)*gen->alpha[il]*gen->alpha[il];
+  double t1 = gen->eps[l+2], t2 = gen->eps[l+1];
+  gen->b[il] = -gen->a[il]*(t1*t1+t2*t2);
+  }
     }
   else
     {
diff --git a/libsharp/sharp_ylmgen_c.h b/libsharp/sharp_ylmgen_c.h
index 63b23cd..2e9e5ff 100644
--- a/libsharp/sharp_ylmgen_c.h
+++ b/libsharp/sharp_ylmgen_c.h
@@ -58,6 +58,8 @@ typedef struct
   double *mfac;
   sharp_ylmgen_dbl2 *rf;
 
+double *eps, *alpha, *a, *b;
+
 /* used if s!=0 */
   int sinPow, cosPow, preMinus_p, preMinus_m;
   double *prefac;