[Propose] Add new method for Bunch-Kaufman factorization

[yoshoku]

I would like to add new method for LDLt or Bunch-Kaufman factorization like scipy.linalg.ldl method. Bunch-Kaufman factorization decomposes a symmetric/Hermitian matrix into the product of an upper / lower triangular matrix and block diagonal matrix such as A = U D U^T.

> a=Numo::DFloat.new(3,3).rand
> a=0.5*(a+a.transpose)
=> Numo::DFloat#shape=[3,3]
[[0.0617545, 0.287054, 0.667382],
 [0.287054, 0.116041, 0.540924],
 [0.667382, 0.540924, 0.165089]]
> u,d,p=Numo::Linalg.ldl(a)
> u.dot(d.dot(u.transpose))
=> Numo::DFloat#shape=[3,3]
[[0.0617545, 0.287054, 0.667382],
 [0.287054, 0.116041, 0.540924],
 [0.667382, 0.540924, 0.165089]]
> u[p,true]
=> Numo::DFloat(view)#shape=[3,3]
[[1, 0.720613, 0.36344],
 [0, 1, 0],
 [0, 0, 1]]