Re: [Maxima-discuss] complete_square

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

Thanks for pointing that out. At one point I factored a 1/2 out of the
c_mat, but undid that change. So I had to fix division elsewhere.
complete_square_alt_form:0$
complete_square(2*x^2-y^2-4*x+4*y+2,[x,y]);
complete_square_alt_form:1$
complete_square(2*x^2-y^2-4*x+4*y+2,[x,y]);
complete_square_alt_form:2$
complete_square(2*x^2-y^2-4*x+4*y+2,[x,y]);

[matrix([x-1],[y-2]),matrix([2,0],[0,-1]),4]
(2-y)*(y-2)+2*(x-1)^2+4
-(y-2)^2+2*(x-1)^2+4

So complete_square_alt_form:2$, produces output most closely like what you
wanted. I'm actually looking at alternative factoring ways to reduce the
number of times a variable appears. I was thinking of trying factoring
using the largest symmetric truncated SVD of the coefficients. Also there
exist techniques for completing the cube, etc. The main purpose of this
work is improve interval estimates by reducing dependencies.


/* 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],
  /* 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]: ev(b_vec[i,1], map(lambda([avar], avar=0), vars)),
    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
)$