Accelerated multiply_poly

src/ao_one_e_ints/pot_ao_ints.irp.f

@@ -455,10 +455,12 @@ recursive subroutine I_x1_pol_mult_one_e(a,c,R1x,R1xp,R2x,d,nd,n_pt_in)
     do ix=0,nx
       X(ix) *= dble(c)
     enddo
-    call multiply_poly(X,nx,R2x,2,d,nd)
+!    call multiply_poly(X,nx,R2x,2,d,nd)
+    call multiply_poly_c2(X,nx,R2x,d,nd)
     ny=0
     call I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,Y,ny,n_pt_in)
-    call multiply_poly(Y,ny,R1x,2,d,nd)
+!    call multiply_poly(Y,ny,R1x,2,d,nd)
+    call multiply_poly_c2(Y,ny,R1x,d,nd)
   else
     do ix=0,n_pt_in
       X(ix) = 0.d0
@@ -469,7 +471,8 @@ recursive subroutine I_x1_pol_mult_one_e(a,c,R1x,R1xp,R2x,d,nd,n_pt_in)
     do ix=0,nx
       X(ix) *= dble(a-1)
     enddo
-    call multiply_poly(X,nx,R2x,2,d,nd)
+!    call multiply_poly(X,nx,R2x,2,d,nd)
+    call multiply_poly_c2(X,nx,R2x,d,nd)
 
     nx = nd
     do ix=0,n_pt_in
@@ -479,10 +482,12 @@ recursive subroutine I_x1_pol_mult_one_e(a,c,R1x,R1xp,R2x,d,nd,n_pt_in)
     do ix=0,nx
       X(ix) *= dble(c)
     enddo
-    call multiply_poly(X,nx,R2x,2,d,nd)
+!    call multiply_poly(X,nx,R2x,2,d,nd)
+    call multiply_poly_c2(X,nx,R2x,d,nd)
     ny=0
     call I_x1_pol_mult_one_e(a-1,c,R1x,R1xp,R2x,Y,ny,n_pt_in)
-    call multiply_poly(Y,ny,R1x,2,d,nd)
+!    call multiply_poly(Y,ny,R1x,2,d,nd)
+    call multiply_poly_c2(Y,ny,R1x,d,nd)
   endif
 end
 
@@ -519,7 +524,8 @@ recursive subroutine I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,dim)
     do ix=0,nx
       X(ix) *= dble(c-1)
     enddo
-    call multiply_poly(X,nx,R2x,2,d,nd)
+!    call multiply_poly(X,nx,R2x,2,d,nd)
+    call multiply_poly_c2(X,nx,R2x,d,nd)
     ny = 0
     do ix=0,dim
       Y(ix) = 0.d0
@@ -527,7 +533,8 @@ recursive subroutine I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,dim)
 
     call I_x1_pol_mult_one_e(0,c-1,R1x,R1xp,R2x,Y,ny,dim)
     if(ny.ge.0)then
-      call multiply_poly(Y,ny,R1xp,2,d,nd)
+!      call multiply_poly(Y,ny,R1xp,2,d,nd)
+      call multiply_poly_c2(Y,ny,R1xp,d,nd)
     endif
   endif
 end

src/ao_two_e_ints/two_e_integrals.irp.f

@@ -975,18 +975,7 @@ recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt
 
 !  !DIR$ FORCEINLINE
 !  call multiply_poly(X,nx,B_10,2,d,nd)
-  if (nx >= 0) then
+  call multiply_poly_c2(X,nx,B_10,d,nd)
-    integer :: ib
-    do ib=0,nx
-      d(ib  ) = d(ib  ) + B_10(0) * X(ib)
-      d(ib+1) = d(ib+1) + B_10(1) * X(ib)
-      d(ib+2) = d(ib+2) + B_10(2) * X(ib)
-    enddo
-
-    do nd = nx+2,0,-1
-      if (d(nd) /= 0.d0) exit
-    enddo
-  endif
 
   nx = nd
   !DIR$ LOOP COUNT(8)
@@ -1009,17 +998,7 @@ recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt
     endif
 !    !DIR$ FORCEINLINE
 !    call multiply_poly(X,nx,B_00,2,d,nd)
-    if (nx >= 0) then
+    call multiply_poly_c2(X,nx,B_00,d,nd)
-       do ib=0,nx
-           d(ib  ) = d(ib  ) + B_00(0) * X(ib)
-           d(ib+1) = d(ib+1) + B_00(1) * X(ib)
-           d(ib+2) = d(ib+2) + B_00(2) * X(ib)
-       enddo
-
-       do nd = nx+2,0,-1
-         if (d(nd) /= 0.d0) exit
-       enddo
-    endif
   endif
 
   ny=0
@@ -1038,17 +1017,7 @@ recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt
 
 !  !DIR$ FORCEINLINE
 !  call multiply_poly(Y,ny,C_00,2,d,nd)
-   if (ny >= 0) then
+  call multiply_poly_c2(Y,ny,C_00,d,nd)
-     do ib=0,ny
-         d(ib  ) = d(ib  ) + C_00(0) * Y(ib)
-         d(ib+1) = d(ib+1) + C_00(1) * Y(ib)
-         d(ib+2) = d(ib+2) + C_00(2) * Y(ib)
-     enddo
-
-   do nd = ny+2,0,-1
-     if (d(nd) /= 0.d0) exit
-   enddo
-  endif
 end
 
 recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
@@ -1088,18 +1057,7 @@ recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
 
 !  !DIR$ FORCEINLINE
 !  call multiply_poly(X,nx,B_00,2,d,nd)
-  if (nx >= 0) then
+  call multiply_poly_c2(X,nx,B_00,d,nd)
-    integer                        :: ib
-    do ib=0,nx
-      d(ib  ) = d(ib  ) + B_00(0) * X(ib)
-      d(ib+1) = d(ib+1) + B_00(1) * X(ib)
-      d(ib+2) = d(ib+2) + B_00(2) * X(ib)
-    enddo
-
-    do nd = nx+2,0,-1
-      if (d(nd) /= 0.d0) exit
-    enddo
-  endif
 
   ny=0
 
@@ -1111,17 +1069,7 @@ recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
 
 !  !DIR$ FORCEINLINE
 !  call multiply_poly(Y,ny,C_00,2,d,nd)
-  if (ny >= 0) then
+  call multiply_poly_c2(Y,ny,C_00,d,nd)
-    do ib=0,ny
-      d(ib  ) = d(ib  ) + C_00(0) * Y(ib)
-      d(ib+1) = d(ib+1) + C_00(1) * Y(ib)
-      d(ib+2) = d(ib+2) + C_00(2) * Y(ib)
-    enddo
-
-    do nd = ny+2,0,-1
-      if (d(nd) /= 0.d0) exit
-    enddo
-  endif
 
 end
 
@@ -1150,18 +1098,7 @@ recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
 
 !  !DIR$ FORCEINLINE
 !  call multiply_poly(X,nx,B_10,2,d,nd)
-  if (nx >= 0) then
+  call multiply_poly_c2(X,nx,B_10,d,nd)
-    integer :: ib
-    do ib=0,nx
-      d(ib  ) = d(ib  ) + B_10(0) * X(ib)
-      d(ib+1) = d(ib+1) + B_10(1) * X(ib)
-      d(ib+2) = d(ib+2) + B_10(2) * X(ib)
-    enddo
-
-    do nd = nx+2,0,-1
-      if (d(nd) /= 0.d0) exit
-    enddo
-  endif
 
   nx = nd
   !DIR$ LOOP COUNT(8)
@@ -1181,17 +1118,7 @@ recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
 
 !  !DIR$ FORCEINLINE
 !  call multiply_poly(X,nx,B_00,2,d,nd)
-  if (nx >= 0) then
+  call multiply_poly_c2(X,nx,B_00,d,nd)
-    do ib=0,nx
-      d(ib  ) = d(ib  ) + B_00(0) * X(ib)
-      d(ib+1) = d(ib+1) + B_00(1) * X(ib)
-      d(ib+2) = d(ib+2) + B_00(2) * X(ib)
-    enddo
-
-    do nd = nx+2,0,-1
-      if (d(nd) /= 0.d0) exit
-    enddo
-  endif
 
   ny=0
   !DIR$ LOOP COUNT(8)
@@ -1203,17 +1130,7 @@ recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
 
 !  !DIR$ FORCEINLINE
 !  call multiply_poly(Y,ny,C_00,2,d,nd)
-  if (ny >= 0) then
+  call multiply_poly_c2(Y,ny,C_00,d,nd)
-    do ib=0,ny
-        d(ib  ) = d(ib  ) + C_00(0) * Y(ib)
-        d(ib+1) = d(ib+1) + C_00(1) * Y(ib)
-        d(ib+2) = d(ib+2) + C_00(2) * Y(ib)
-    enddo
-
-    do nd = ny+2,0,-1
-      if (d(nd) /= 0.d0) exit
-    enddo
-  endif
 end
 
 recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
@@ -1262,18 +1179,7 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
 
 !      !DIR$ FORCEINLINE
 !      call multiply_poly(Y,ny,D_00,2,d,nd)
-      if (ny >= 0) then
+      call multiply_poly_c2(Y,ny,D_00,d,nd)
-        integer :: ib
-        do ib=0,ny
-            d(ib  ) = d(ib  ) + D_00(0) * Y(ib)
-            d(ib+1) = d(ib+1) + D_00(1) * Y(ib)
-            d(ib+2) = d(ib+2) + D_00(2) * Y(ib)
-        enddo
-
-        do nd = ny+2,0,-1
-          if (d(nd) /= 0.d0) exit
-        enddo
-      endif
 
       return
 
@@ -1293,17 +1199,7 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
 
 !      !DIR$ FORCEINLINE
 !      call multiply_poly(X,nx,B_01,2,d,nd)
-      if (nx >= 0) then
+      call multiply_poly_c2(X,nx,B_01,d,nd)
-        do ib=0,nx
-          d(ib  ) = d(ib  ) + B_01(0) * X(ib)
-          d(ib+1) = d(ib+1) + B_01(1) * X(ib)
-          d(ib+2) = d(ib+2) + B_01(2) * X(ib)
-        enddo
-
-        do nd = nx+2,0,-1
-          if (d(nd) /= 0.d0) exit
-        enddo
-      endif
 
       ny = 0
       !DIR$ LOOP COUNT(6)
@@ -1314,17 +1210,7 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
 
 !      !DIR$ FORCEINLINE
 !      call multiply_poly(Y,ny,D_00,2,d,nd)
-      if (ny >= 0) then
+      call multiply_poly_c2(Y,ny,D_00,d,nd)
-        do ib=0,ny
-            d(ib  ) = d(ib  ) + D_00(0) * Y(ib)
-            d(ib+1) = d(ib+1) + D_00(1) * Y(ib)
-            d(ib+2) = d(ib+2) + D_00(2) * Y(ib)
-        enddo
-
-        do nd = ny+2,0,-1
-          if (d(nd) /= 0.d0) exit
-        enddo
-      endif
 
   end select
 end

src/bi_ort_ints/three_body_ints_bi_ort.irp.f

@@ -79,28 +79,30 @@ subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral)
   integer,          intent(in)  :: n, l, k, m, j, i
   double precision, intent(out) :: integral
   integer                       :: ipoint
-  double precision              :: weight
+  double precision              :: weight, tmp
 
   PROVIDE mo_l_coef mo_r_coef
   PROVIDE int2_grad1_u12_bimo_t
 
   integral = 0.d0
   do ipoint = 1, n_points_final_grid
-    weight = final_weight_at_r_vector(ipoint)                                                                          
 
-    integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & 
+    tmp =     mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) &
               * ( int2_grad1_u12_bimo_t(ipoint,1,n,m) * int2_grad1_u12_bimo_t(ipoint,1,l,j)    &
                 + int2_grad1_u12_bimo_t(ipoint,2,n,m) * int2_grad1_u12_bimo_t(ipoint,2,l,j)    &
                 + int2_grad1_u12_bimo_t(ipoint,3,n,m) * int2_grad1_u12_bimo_t(ipoint,3,l,j) )
-    integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & 
+
+    tmp = tmp + mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) &
               * ( int2_grad1_u12_bimo_t(ipoint,1,n,m) * int2_grad1_u12_bimo_t(ipoint,1,k,i)    &
                 + int2_grad1_u12_bimo_t(ipoint,2,n,m) * int2_grad1_u12_bimo_t(ipoint,2,k,i)    &
                 + int2_grad1_u12_bimo_t(ipoint,3,n,m) * int2_grad1_u12_bimo_t(ipoint,3,k,i) )
-    integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) &
+
+    tmp = tmp + mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) &
               * ( int2_grad1_u12_bimo_t(ipoint,1,l,j) * int2_grad1_u12_bimo_t(ipoint,1,k,i)    &
                 + int2_grad1_u12_bimo_t(ipoint,2,l,j) * int2_grad1_u12_bimo_t(ipoint,2,k,i)    &
                 + int2_grad1_u12_bimo_t(ipoint,3,l,j) * int2_grad1_u12_bimo_t(ipoint,3,k,i) )
 
+    integral = integral + tmp * final_weight_at_r_vector(ipoint)
   enddo
 
 end subroutine give_integrals_3_body_bi_ort

src/utils/integration.irp.f

@@ -468,114 +468,6 @@ end subroutine
 
 
 
-subroutine multiply_poly_0c(b,c,nc,d,nd)
-  implicit none
-  BEGIN_DOC
-  ! Multiply two polynomials
-  ! D(t) += B(t)*C(t)
-  END_DOC
-
-  integer, intent(in)            :: nc
-  integer, intent(out)           :: nd
-  double precision, intent(in)   :: b(0:0), c(0:nc)
-  double precision, intent(inout) :: d(0:0+nc)
-
-  integer                        :: ic
-
-  do ic = 0,nc
-    d(ic) = d(ic) + c(ic) * b(0)
-  enddo
-
-  do nd = nc,0,-1
-    if (d(nd) /= 0.d0) exit
-  enddo
-
-end
-
-subroutine multiply_poly_1c(b,c,nc,d,nd)
-  implicit none
-  BEGIN_DOC
-  ! Multiply two polynomials
-  ! D(t) += B(t)*C(t)
-  END_DOC
-
-  integer, intent(in)            :: nc
-  integer, intent(out)           :: nd
-  double precision, intent(in)   :: b(0:1), c(0:nc)
-  double precision, intent(inout) :: d(0:1+nc)
-
-  integer                        :: ic, id
-  if(nc < 0) return
-
-  do ic = 0,nc
-    d(  ic) = d(  ic) + c(ic) * b(0)
-    d(1+ic) = d(1+ic) + c(ic) * b(1)
-  enddo
-
-  do nd = nc+1,0,-1
-    if (d(nd) /= 0.d0) exit
-  enddo
-
-end
-
-
-subroutine multiply_poly_2c(b,c,nc,d,nd)
-  implicit none
-  BEGIN_DOC
-  ! Multiply two polynomials
-  ! D(t) += B(t)*C(t)
-  END_DOC
-
-  integer, intent(in)            :: nc
-  integer, intent(out)           :: nd
-  double precision, intent(in)   :: b(0:2), c(0:nc)
-  double precision, intent(inout) :: d(0:2+nc)
-
-  integer                        :: ic, id, k
-  if (nc <0) return
-
-  do ic = 0,nc
-    d(  ic) = d(  ic) + c(ic) * b(0)
-    d(1+ic) = d(1+ic) + c(ic) * b(1)
-    d(2+ic) = d(2+ic) + c(ic) * b(2)
-  enddo
-
-  do nd = nc+2,0,-1
-    if (d(nd) /= 0.d0) exit
-  enddo
-
-end
-
-subroutine multiply_poly_3c(b,c,nc,d,nd)
-  implicit none
-  BEGIN_DOC
-  ! Multiply two polynomials
-  ! D(t) += B(t)*C(t)
-  END_DOC
-
-  integer, intent(in)            :: nc
-  integer, intent(out)           :: nd
-  double precision, intent(in)   :: b(0:3), c(0:nc)
-  double precision, intent(inout) :: d(0:3+nc)
-
-  integer                        :: ic, id
-  if (nc <0) return
-
-  do ic = 1,nc
-    d(  ic) = d(1+ic) + c(ic) * b(0)
-    d(1+ic) = d(1+ic) + c(ic) * b(1)
-    d(2+ic) = d(1+ic) + c(ic) * b(2)
-    d(3+ic) = d(1+ic) + c(ic) * b(3)
-  enddo
-
-  do nd = nc+3,0,-1
-    if (d(nd) /= 0.d0) exit
-  enddo
-
-end
-
-
-
 subroutine multiply_poly(b,nb,c,nc,d,nd)
   implicit none
   BEGIN_DOC
@@ -604,6 +496,254 @@ subroutine multiply_poly(b,nb,c,nc,d,nd)
 
 end
 
+
+subroutine multiply_poly_b0(b,c,nc,d,nd)
+  implicit none
+  BEGIN_DOC
+  ! Multiply two polynomials
+  ! D(t) += B(t)*C(t)
+  END_DOC
+
+  integer, intent(in)            :: nc
+  integer, intent(out)           :: nd
+  double precision, intent(in)   :: b(0:0), c(0:nc)
+  double precision, intent(inout) :: d(0:nc)
+
+  integer                        :: ndtmp
+  integer                        :: ic, id, k
+  if(nc < 0) return !False if nc>=0
+
+  do ic = 0,nc
+    d(ic) = d(ic) + c(ic) * b(0)
+  enddo
+
+  do nd = nc,0,-1
+    if (d(nd) /= 0.d0) exit
+  enddo
+
+end
+
+subroutine multiply_poly_b1(b,c,nc,d,nd)
+  implicit none
+  BEGIN_DOC
+  ! Multiply two polynomials
+  ! D(t) += B(t)*C(t)
+  END_DOC
+
+  integer, intent(in)            :: nc
+  integer, intent(out)           :: nd
+  double precision, intent(in)   :: b(0:1), c(0:nc)
+  double precision, intent(inout) :: d(0:1+nc)
+
+  integer                        :: ndtmp
+  integer                        :: ib, ic, id, k
+  if(nc < 0) return !False if nc>=0
+
+
+  select case (nc)
+    case (0)
+      d(0) = d(0) + c(0) * b(0)
+      d(1) = d(1) + c(0) * b(1)
+
+    case (1)
+      d(0) = d(0) + c(0) * b(0)
+      d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
+      d(2) = d(2) + c(1) * b(1)
+
+    case default
+      d(0) = d(0) + c(0) * b(0)
+      do ic = 1,nc
+        d(ic) = d(ic) + c(ic) * b(0) + c(ic-1) * b(1)
+      enddo
+      d(nc+1) = d(nc+1) + c(nc) * b(1)
+
+  end select
+
+  do nd = 1+nc,0,-1
+    if (d(nd) /= 0.d0) exit
+  enddo
+
+end
+
+
+subroutine multiply_poly_b2(b,c,nc,d,nd)
+  implicit none
+  BEGIN_DOC
+  ! Multiply two polynomials
+  ! D(t) += B(t)*C(t)
+  END_DOC
+
+  integer, intent(in)            :: nc
+  integer, intent(out)           :: nd
+  double precision, intent(in)   :: b(0:2), c(0:nc)
+  double precision, intent(inout) :: d(0:2+nc)
+
+  integer                        :: ndtmp
+  integer                        :: ib, ic, id, k
+  if(nc < 0) return !False if nc>=0
+
+  select case (nc)
+    case (0)
+      d(0) = d(0) + c(0) * b(0)
+      d(1) = d(1) + c(0) * b(1)
+      d(2) = d(2) + c(0) * b(2)
+
+    case (1)
+      d(0) = d(0) + c(0) * b(0)
+      d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
+      d(2) = d(2) + c(0) * b(2) + c(1) * b(1)
+      d(3) = d(3) + c(1) * b(2)
+
+    case (2)
+      d(0) = d(0) + c(0) * b(0)
+      d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
+      d(2) = d(2) + c(0) * b(2) + c(1) * b(1) + c(2) * b(0)
+      d(3) = d(3) + c(2) * b(1) + c(1) * b(2)
+      d(4) = d(4) + c(2) * b(2)
+
+    case default
+
+      d(0) = d(0) + c(0) * b(0)
+      d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
+      do ic = 2,nc
+        d(ic) = d(ic) + c(ic) * b(0) + c(ic-1) * b(1) + c(ic-2) * b(2)
+      enddo
+      d(nc+1) = d(nc+1) + c(nc) * b(1) + c(nc-1) * b(2)
+      d(nc+2) = d(nc+2) + c(nc) * b(2)
+
+  end select
+
+  do nd = 2+nc,0,-1
+    if (d(nd) /= 0.d0) exit
+  enddo
+
+end
+
+
+subroutine multiply_poly_c0(b,nb,c,d,nd)
+  implicit none
+  BEGIN_DOC
+  ! Multiply two polynomials
+  ! D(t) += B(t)*C(t)
+  END_DOC
+
+  integer, intent(in)            :: nb
+  integer, intent(out)           :: nd
+  double precision, intent(in)   :: b(0:nb), c(0:0)
+  double precision, intent(inout) :: d(0:nb)
+
+  integer                        :: ndtmp
+  integer                        :: ib, ic, id, k
+  if(nb < 0) return !False if nb>=0
+
+  do ib=0,nb
+      d(ib) = d(ib) + c(0) * b(ib)
+  enddo
+
+  do nd = nb,0,-1
+    if (d(nd) /= 0.d0) exit
+  enddo
+
+end
+
+
+subroutine multiply_poly_c1(b,nb,c,d,nd)
+  implicit none
+  BEGIN_DOC
+  ! Multiply two polynomials
+  ! D(t) += B(t)*C(t)
+  END_DOC
+
+  integer, intent(in)            :: nb
+  integer, intent(out)           :: nd
+  double precision, intent(in)   :: b(0:nb), c(0:1)
+  double precision, intent(inout) :: d(0:nb+1)
+
+  integer                        :: ndtmp
+  integer                        :: ib, ic, id, k
+  if(nb < 0) return !False if nb>=0
+
+  select case (nb)
+    case (0)
+      d(0) = d(0) + c(0) * b(0)
+      d(1) = d(1) + c(1) * b(0)
+
+    case (1)
+      d(0) = d(0) + c(0) * b(0)
+      d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
+      d(2) = d(2) + c(1) * b(1)
+
+    case default
+      d(0) = d(0) + c(0) * b(0)
+      do ib=1,nb
+        d(ib) = d(ib) + c(0) * b(ib) + c(1) * b(ib-1)
+      enddo
+      d(nb+1) = d(nb+1) + c(1) * b(nb)
+
+  end select
+
+  do nd = nb+1,0,-1
+    if (d(nd) /= 0.d0) exit
+  enddo
+
+end
+
+
+subroutine multiply_poly_c2(b,nb,c,d,nd)
+  implicit none
+  BEGIN_DOC
+  ! Multiply two polynomials
+  ! D(t) += B(t)*C(t)
+  END_DOC
+
+  integer, intent(in)            :: nb
+  integer, intent(out)           :: nd
+  double precision, intent(in)   :: b(0:nb), c(0:2)
+  double precision, intent(inout) :: d(0:nb+2)
+
+  integer                        :: ndtmp
+  integer                        :: ib, ic, id, k
+  if(nb < 0) return !False if nb>=0
+
+  select case (nb)
+    case (0)
+      d(0) = d(0) + c(0) * b(0)
+      d(1) = d(1) + c(1) * b(0)
+      d(2) = d(2) + c(2) * b(0)
+
+    case (1)
+      d(0) = d(0) + c(0) * b(0)
+      d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
+      d(2) = d(2) + c(1) * b(1) + c(2) * b(0)
+      d(3) = d(3) + c(2) * b(1)
+
+    case (2)
+      d(0) = d(0) + c(0) * b(0)
+      d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
+      d(2) = d(2) + c(0) * b(2) + c(1) * b(1) + c(2) * b(0)
+      d(3) = d(3) + c(1) * b(2) + c(2) * b(1)
+      d(4) = d(4) + c(2) * b(2)
+
+    case default
+      d(0) = d(0) + c(0) * b(0)
+      d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
+      do ib=2,nb
+        d(ib) = d(ib) + c(0) * b(ib) + c(1) * b(ib-1) + c(2) * b(ib-2)
+      enddo
+      d(nb+1) = d(nb+1) + c(1) * b(nb) + c(2) * b(nb-1)
+      d(nb+2) = d(nb+2) + c(2) * b(nb)
+
+  end select
+
+  do nd = nb+2,0,-1
+    if (d(nd) /= 0.d0) exit
+  enddo
+
+end
+
+
+
+
 subroutine multiply_poly_v(b,nb,c,nc,d,nd,n_points)
   implicit none
   BEGIN_DOC