// Copyright (c) the JPEG XL Project Authors. All rights reserved.

//

// Use of this source code is governed by a BSD-style

// license that can be found in the LICENSE file.

#include "lib/jxl/enc_linalg.h"

#include <cmath>

#include "lib/jxl/base/status.h"

namespace jxl {

void ConvertToDiagonal(const Matrix2x2& A, Vector2& diag, Matrix2x2& U) {

  // Check A is symmetric.

  JXL_DASSERT(std::abs(A[0][1] - A[1][0]) < 1e-15);

  double b = -(A[0][0] + A[1][1]);

  double c = A[0][0] * A[1][1] - A[0][1] * A[0][1];

  double d = b * b - 4.0 * c;

  if (std::abs(A[0][1]) < 1e-10 || d < 0) {

    // Already diagonal.

    diag[0] = A[0][0];

    diag[1] = A[1][1];

    U[0][0] = U[1][1] = 1.0;

    U[0][1] = U[1][0] = 0.0;

    return;

  }

  double sqd = std::sqrt(d);

  double l1 = (-b - sqd) * 0.5;

  double l2 = (-b + sqd) * 0.5;

  Vector2 v1 = {A[0][0] - l1, A[1][0]};

  double v1n = 1.0 / std::hypot(v1[0], v1[1]);

  v1[0] = v1[0] * v1n;

  v1[1] = v1[1] * v1n;

  diag[0] = l1;

  diag[1] = l2;

  U[0][0] = v1[1];

  U[0][1] = -v1[0];

  U[1][0] = v1[0];

  U[1][1] = v1[1];

}

}  // namespace jxl