Commits

Barry Schwartz committed 594944f

Polynomial multiplication and related functions in Pure for imatrix and smatrix.

Comments (0)

Files changed (1)

pure/ebeno/curve.pure

 
 namespace ebeno;
 
-public mo_of_sb sb_of_mo sr_of_sb sb_of_sr;
+public
+  be_of_sb sb_of_be
+  mo_of_sb sb_of_mo
+  sr_of_sb sb_of_sr
+  split_sr unsplit_sr
+  dropterms_sr;
 
 //-------------------------------------------------------------------------
 
     ebeno_convolve_pa da db (pack a) (pack b) c;
   end
   if rowvectorp a && rowvectorp b;
+convolve a::imatrix b::imatrix |
+convolve a::smatrix b::smatrix = {term j | j = 0..dc}
+  with
+    term j::int = foldl (+) 0 (zipwith (*) a_slice b_slice)
+      when
+        a_slice = a_rev!!(da-j..da);
+        b_slice = b!!(j-da..db);
+      end;
+  end
+  when
+    da::int = deg a;
+    db::int = deg b;
+    dc::int = da + db;
+    a_rev = reverse a;
+  end;
 
 // Compose polynomials.
 private compose_np_ compose_pa_;
     ebeno_mul_be_pa da db (pack a) (pack b) c;
   end
   if rowvectorp a && rowvectorp b;
+mul_be a::imatrix b::imatrix |
+mul_be a::smatrix b::smatrix =
+  be_of_sb $ convolve (sb_of_be a) (sb_of_be b)
+  if rowvectorp a && rowvectorp b;
 
 // Multiplication of polynomials in Sánchez-Reyes (symmetric power)
 // basis.
     ebeno_mul_sr_pa da db (pack a) (pack b) c;
   end
   if rowvectorp a && rowvectorp b;
+mul_sr a::imatrix b::imatrix |
+mul_sr a::smatrix b::smatrix =
+  dropterms_sr (unsplit_sr (c1,c2)) (deg a + deg b)
+  when
+    a1,a2 = split_sr a;
+    b1,b2 = split_sr b;
+    d = {0, convolve (a2 - a1) (b2 - b1)};
+    c1 = {convolve a1 b1, 0} - d;
+    c2 = {convolve a2 b2, 0} - d;
+  end
+  if rowvectorp a && rowvectorp b;
 
 // Derivatives.
 private deriv_np_ deriv_pa_;
 sr_of_mo b::imatrix |
 sr_of_mo b::smatrix = sr_of_sb $ sb_of_mo b; // FIXME: Do this directly.
 
+// Drop terms of degree higher than n.
+dropterms_mo b::dmatrix n::int =
+  if n < 0 then {0.0} else c
+  when
+    count::int = (dim b)!0;
+    db::int = deg b;
+    dc::int = max 0 (min n db);
+    c::dmatrix = {0.0 | _ = 0..dc};
+    ebeno_dropterms_mo_np_many db count (pack b) n c;
+  end;
+dropterms_mo b::cmatrix n::int =
+  if n < 0 then {0.0+:0.0} else c
+  when
+    count::int = (dim b)!0;
+    db::int = deg b;
+    dc::int = max 0 (min n db);
+    c::cmatrix = {0.0+:0.0 | _ = 0..dc};
+    ebeno_dropterms_mo_pa_many db count (pack b) n c;
+  end;
+dropterms_mo b::imatrix n::int |
+dropterms_mo b::smatrix n::int =
+  if n < 0 then
+    colvector $ repeatn (dim b!0) 0
+  else
+    b!!(0..(dim b!0)-1, 0..n);
+dropterms_sr b::dmatrix n::int =
+  if n < 0 then {0.0} else c
+  when
+    count::int = (dim b)!0;
+    db::int = deg b;
+    dc::int = max 0 (min n db);
+    c::dmatrix = {0.0 | _ = 0..dc};
+    ebeno_dropterms_sr_np_many db count (pack b) n c;
+  end;
+dropterms_sr b::cmatrix n::int =
+  if n < 0 then {0.0+:0.0} else c
+  when
+    count::int = (dim b)!0;
+    db::int = deg b;
+    dc::int = max 0 (min n db);
+    c::cmatrix = {0.0+:0.0 | _ = 0..dc};
+    ebeno_dropterms_sr_pa_many db count (pack b) n c;
+  end;
+dropterms_sr b::imatrix n::int |
+dropterms_sr b::smatrix n::int =
+  if n < 0 then
+    colvector $ repeatn (dim b!0) 0
+  else if db <= n then
+    b
+  else
+    (redim (dim b!0,n+1) {do_row v | v = rows b}
+     with
+       do_row v =
+         case n_mod_2 of
+           0 = {v!!(0..q-1), (v!q + v!(db-q))%2, v!!(db-q+1..db)};
+           1 = {v!!(0..q-1), v!!(db-q+1..db)};
+         end;
+     end
+     when
+       n_mod_2::int = n mod 2;
+       q::int = n div 2 + n_mod_2;
+     end)
+  when
+    db::int = deg b;
+  end;
+
 // Evaluation of a Bézier curve at time t.
 eval_be b::dmatrix t::double = ebeno_eval_be_np (deg b) (pack b) t
   if rowvectorp b;
     ebeno_split_sr_pa db (pack b) c;
   end
   if rowvectorp b;
+split_sr b::imatrix |
+split_sr b::smatrix =
+  case db_mod_2 of
+    0 = b!!(0..q), b!!(db:db-1..q);
+    1 = b!!(0..q-1), b!!(db:db-1..q);
+  end
+  when
+    db::int = deg b;
+    db_mod_2::int = db mod 2;
+    q::int = db div 2 + db_mod_2;
+  end
+  if rowvectorp b;
 unsplit_sr c::dmatrix = b
   when
     c_size::int = (dim c)!1;
     ebeno_unsplit_sr_pa db (pack c) b;
   end
   if (dim c)!0 == 2;
+unsplit_sr (a::imatrix, b::imatrix) |
+unsplit_sr (a::smatrix, b::smatrix) = {a, reverse b}
+  if rowvectorp a && rowvectorp b && #a == #b;
 
 //-------------------------------------------------------------------------
 //
 type bern_dvector (Bern a::dmatrix) = rowvectorp a;
 type bern_cvector (Bern a::cmatrix) = rowvectorp a;
 type bern_ivector (Bern a::imatrix) = rowvectorp a;
-type bern_ivector (Bern a::smatrix) = rowvectorp a;
+type bern_svector (Bern a::smatrix) = rowvectorp a;
 type bern_matrix m::bern_dmatrix;
 type bern_matrix m::bern_cmatrix;
 type bern_matrix m::bern_imatrix;
 a::sanch_dvector ::- b::sanch_dvector = Sanch $ minus_sr (vector a) (vector b);
 
 a::bern_cvector ::* b::bern_cvector |
-a::bern_dvector ::* b::bern_dvector = Bern $ mul_be (vector a) (vector b);
+a::bern_dvector ::* b::bern_dvector |
+a::bern_ivector ::* b::bern_ivector |
+a::bern_svector ::* b::bern_svector = Bern $ mul_be (vector a) (vector b);
 a::sbern_cvector ::* b::sbern_cvector |
-a::sbern_dvector ::* b::sbern_dvector = SBern $ convolve (vector a) (vector b);
+a::sbern_dvector ::* b::sbern_dvector |
+a::sbern_ivector ::* b::sbern_ivector |
+a::sbern_svector ::* b::sbern_svector = SBern $ convolve (vector a) (vector b);
 a::mono_cvector ::* b::mono_cvector |
-a::mono_dvector ::* b::mono_dvector = Mono $ convolve (vector a) (vector b);
+a::mono_dvector ::* b::mono_dvector |
+a::mono_ivector ::* b::mono_ivector |
+a::mono_svector ::* b::mono_svector = Mono $ convolve (vector a) (vector b);
 a::sanch_cvector ::* b::sanch_cvector |
-a::sanch_dvector ::* b::sanch_dvector = Sanch $ mul_sr (vector a) (vector b);
+a::sanch_dvector ::* b::sanch_dvector |
+a::sanch_ivector ::* b::sanch_ivector |
+a::sanch_svector ::* b::sanch_svector = Sanch $ mul_sr (vector a) (vector b);
 
 find_roots v::bern_dvector = findroots_be $ vector v;
 find_roots v::sbern_dvector = findroots_sb $ vector v;
 intersect a::bern_cvector b::bern_cvector =
   intersect_be (vector a) (vector b);
 
-drop_terms b::mono_dmatrix n::int =
-  if n < 0 then Mono {0.0} else Mono c
+drop_terms b::mono_matrix n::int = Mono $ dropterms_mo (matrix b) n;
+drop_terms b::sanch_matrix n::int = Sanch $ dropterms_sr (matrix b) n;
+
+split b::sanch_vector = Mono v, Mono w
   when
-    m::dmatrix = matrix b;
-    count::int = (dim m)!0;
-    db::int = deg b;
-    dc::int = max 0 (min n db);
-    c::dmatrix = {0.0 | _ = 0..dc};
-    ebeno_dropterms_mo_np_many db count (pack m) n c;
+    v,w = split_sr $ vector b;
   end;
-drop_terms b::mono_cmatrix n::int =
-  if n < 0 then Mono {0.0+:0.0} else Mono c
-  when
-    m::cmatrix = matrix b;
-    count::int = (dim m)!0;
-    db::int = deg b;
-    dc::int = max 0 (min n db);
-    c::cmatrix = {0.0+:0.0 | _ = 0..dc};
-    ebeno_dropterms_mo_pa_many db count (pack m) n c;
-  end;
-drop_terms b::sanch_dmatrix n::int =
-  if n < 0 then Sanch {0.0} else Sanch c
-  when
-    m::dmatrix = matrix b;
-    count::int = (dim m)!0;
-    db::int = deg b;
-    dc::int = max 0 (min n db);
-    c::dmatrix = {0.0 | _ = 0..dc};
-    ebeno_dropterms_sr_np_many db count (pack m) n c;
-  end;
-drop_terms b::sanch_cmatrix n::int =
-  if n < 0 then Sanch {0.0+:0.0} else Sanch c
-  when
-    m::cmatrix = matrix b;
-    count::int = (dim m)!0;
-    db::int = deg b;
-    dc::int = max 0 (min n db);
-    c::cmatrix = {0.0+:0.0 | _ = 0..dc};
-    ebeno_dropterms_sr_pa_many db count (pack m) n c;
-  end;
-
-split b::sanch_vector = split_sr $ vector b;
-unsplit c::dmatrix |
-unsplit c::cmatrix = Sanch $ unsplit_sr c if (dim c)!0 == 2;
+unsplit (a::mono_vector, b::mono_vector) =
+  Sanch {vector a, reverse (vector b)}
+  if deg a == deg b;
 
 //-------------------------------------------------------------------------