diff --git a/libsharp/sharp_core.c b/libsharp/sharp_core.c
index 86a40c2..afeded5 100644
--- a/libsharp/sharp_core.c
+++ b/libsharp/sharp_core.c
@@ -283,54 +283,113 @@ NOINLINE static void calc_alm2map (sharp_job * restrict job,
 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 lmax=gen->lmax;
   int nv2 = (nth+VLEN-1)/VLEN;
-  job->opcnt += (lmax+1-l) * 6*nth;
-
-const double inv_sqrt4pi = 0.2820947917738781434740397257803862929220;
+  job->opcnt += (lmax+1-gen->m) * 6*nth;
+  int full_ieee;
+Tv csq[nv0];
   Tv mfac = vload(gen->mfac[gen->m]);
   if (gen->m&1) mfac=-mfac;
+  int below_limit = 1;
   for (int i=0; i<nv2; ++i)
     {
+csq[i] = d->cth[i]*d->cth[i];
     d->lam1[i]=vzero;
-    mypow(d->sth[i],l,gen->powlimit,&d->lam2[i],&d->scale[i]);
+    mypow(d->sth[i],gen->m,gen->powlimit,&d->lam2[i],&d->scale[i]);
     d->lam2[i] *= mfac;
     Tvnormalize(&d->lam2[i],&d->scale[i],sharp_ftol);
+    below_limit &= vallTrue(vlt(d->scale[i],vload(sharp_limscale)));
+    }
+  int l=gen->m, il=0;
+  while (below_limit)
+    {
+    if (l+2>gen->lmax) {l=gen->lmax+1;return;}
+    below_limit=1;
+    for (int i=0; i<nv2; ++i)
+      {
+      Tv tmp = (gen->a[il]*csq[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)))
+        below_limit &= vallTrue(vlt(d->scale[i],vload(sharp_limscale)));
+      }
+    l+=2; ++il;
     }
 
-  const dcmplx * restrict alm=job->almtmp;
-  dcmplx * restrict alm2=RALLOC(dcmplx, gen->lmax+5);
+   const dcmplx * restrict alm=job->almtmp;
+   full_ieee=1;
+   for (int i=0; i<nv2; ++i)
+     {
+     getCorfac(d->scale[i], &d->corfac[i], gen->cf);
+     full_ieee &= vallTrue(vge(d->scale[i],vload(sharp_minscale)));
+     }
   {
-  for (int il=0, l=gen->m; l<=gen->lmax; ++il,l+=2)
+ // int il=0, l=gen->m;
+  for (; (l<=lmax) && (!full_ieee); ++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]));
+    Tv ar1=vload(creal(alm[l  ])), ai1=vload(cimag(alm[l  ]));
+    Tv ar2=vload(creal(alm[l+1])), ai2=vload(cimag(alm[l+1]));
+    full_ieee=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->p2r[i] += d->lam2[i]*d->corfac[i]*ar2;
+      d->p2i[i] += d->lam2[i]*d->corfac[i]*ai2;
+      Tv tmp = (gen->a[il]*csq[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);
+        full_ieee &= vallTrue(vge(d->scale[i],vload(sharp_minscale)));
+        }
+   //   full_ieee=0;
       }
     }
-  DEALLOC(alm2);
+  for (int i=0; i<nv2; ++i)
+    {
+    d->lam1[i] *= d->corfac[i];
+    d->lam2[i] *= d->corfac[i];
+    }
+  for (; l<=lmax-2; il+=2, l+=4)
+    {
+    Tv ar1=vload(creal(alm[l  ])), ai1=vload(cimag(alm[l  ]));
+    Tv ar2=vload(creal(alm[l+1])), ai2=vload(cimag(alm[l+1]));
+    Tv ar3=vload(creal(alm[l+2])), ai3=vload(cimag(alm[l+2]));
+    Tv ar4=vload(creal(alm[l+3])), ai4=vload(cimag(alm[l+3]));
+    Tv a1=vload(gen->a[il]), a2=vload(gen->a[il+1]),
+       b1=vload(gen->b[il]), b2=vload(gen->b[il+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*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->lam2[i] = (a2*csq[i] + b2)*d->lam1[i] + d->lam2[i];
+     }
+    }
+  for (; l<=lmax; ++il, l+=2)
+    {
+    Tv ar1=vload(creal(alm[l  ])), ai1=vload(cimag(alm[l  ]));
+    Tv ar2=vload(creal(alm[l+1])), ai2=vload(cimag(alm[l+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;
+      Tv tmp = (gen->a[il]*csq[i] + gen->b[il])*d->lam2[i] + d->lam1[i];
+      d->lam1[i] = d->lam2[i];
+      d->lam2[i] = tmp;
+      }
+    }
+  }
   }
 
 NOINLINE static void map2alm_kernel(s0data_v * restrict d,
@@ -833,6 +892,17 @@ NOINLINE static void inner_loop_a2m(sharp_job *job, const int *ispair,
       {
       if (job->spin==0)
         {
+//adjust the a_lm for the new algorithm
+dcmplx * restrict alm=job->almtmp;
+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];
+  alm[l  ] = gen->alpha[il]*(gen->eps[l+1]*al + gen->eps[l+2]*al2);
+  alm[l+1] = gen->alpha[il]*al1;
+  }
+
         const int nval=nv0*VLEN;
         int ith=0;
         int itgt[nval];
@@ -868,6 +938,9 @@ NOINLINE static void inner_loop_a2m(sharp_job *job, const int *ispair,
             calc_alm2map_alt (job, gen, &d.v, nth);
             for (int i=0; i<nth; ++i)
               {
+//adjust for new algorithm
+d.s.p2r[i]*=d.s.cth[i];
+d.s.p2i[i]*=d.s.cth[i];
               int tgt=itgt[i];
               int phas_idx = tgt*job->s_th + mi*job->s_m;
               complex double r1 = d.s.p1r[i] + d.s.p1i[i]*_Complex_I,