Re: [Maxima-discuss] complete_square

[Mike V. <mic...@gm...>]

Eventually, after this becomes good enough, I'm going to add it to my
github collection.


/* This file is released under the GNU General Public License 3 */
/* http://www.gnu.org/licenses/gpl-3.0.en.html */

/* This is the basic outline for completing the square in N dimensions */
/* I would appreciate it if someone can do some better error checking */
/* known issue when c_mat is non-invertible:
 complete_square(expand(x^2+y^2), [x,y]) */

"Use: complete_square(Quadratic expression, [vars])";
"Option: complete_square_alt_form: 0, 1, or 2";
" 0) This returns 3 elements [v,C,a] where the completed square is
transpose(v) . C . v + a";
" 1) This mode carries out the matrix multiplication in 0";
" 2) This mode attempts to carry out the matrix multiplication in a
non-standard way to produce a nicer form";

complete_square_alt_form:0$
complete_square(expr, vars):=block([a_const:expr, b_vec, c_mat, c_mat_inv,
new_const, offset, x],
  /* Local vars */
  local(new_const, offset, b_vec, c_mat, c_mat_inv, x),

  /* Needs proper error checking */
  if(not listp(vars)) then error("complete_square: the variable vars should
be be a list of variables"),
  for i: 1 thru length(vars) do
    if(not atom(vars[i]) or subvarp(vars[i])) then
      error("complete_square: the variable vars should be a list of atoms
or indexed variables"),

  /* Allocate space */
  b_vec: genmatrix(lambda([i,j], 0), length(vars), 1),
  c_mat: genmatrix(lambda([i,j], 0), length(vars), length(vars)),

  /* Build the vector of coefficients of variables with degree 1, purge
cross-terms, and isolate the constant */
  for i: 1 thru length(vars) do block(
    b_vec[i,1]: coeff(expr, vars[i]),
    b_vec[i,1]: subst(map(lambda([avar], avar=0), vars),b_vec[i,1]),
    a_const: ratcoef(a_const, vars[i], 0)
  ),

  /* Build the matrix of coefficients for degree 2 */
  for i: 1 thru length(vars) do
    for j: 1 thru length(vars) do
      c_mat[i,j]: ratcoef(expr, vars[i]*vars[j])/(2-kron_delta(i,j)),

  /* Build the offset and new constant*/
  c_mat_inv:invert(c_mat),
  offset: ratsimp( -c_mat_inv . b_vec / 2),
  new_const: ratsimp(a_const - transpose(b_vec) . c_mat_inv . b_vec/4),
  x: ''(vars - offset),

  /* Output in the most pretty sense */
  if(length(vars)=1) then ((vars - offset)^2 + ratsimp(new_const))[1,1]
  elseif(complete_square_alt_form=2) then matrix_sum_elements((x .
transpose(x))*c_mat) + ratsimp(new_const)
  elseif(complete_square_alt_form=1) then transpose(x) . c_mat . x +
ratsimp(new_const)
  else [x, c_mat, ratsimp(new_const)]
)$

matrix_sum_elements(amat):= block([res_sum:0],
  for i: 1 thru length(amat) do res_sum: res_sum + lsum(x,x,amat[i]),
  res_sum
)$