function balance (arr, N, prec, type, findVectors) {
    const big = type === 'BigNumber'
    const cplx = type === 'Complex'

    const realzero = big ? bignumber(0) : 0
    const one = big ? bignumber(1) : cplx ? complex(1) : 1
    const realone = big ? bignumber(1) : 1

    // base of the floating-point arithmetic
    const radix = big ? bignumber(10) : 2
    const radixSq = multiplyScalar(radix, radix)

    // the diagonal transformation matrix R
    let Rdiag
    if (findVectors) {
      Rdiag = Array(N).fill(one)
    }

    // this isn't the only time we loop thru the matrix...
    let last = false

    while (!last) {
      // ...haha I'm joking! unless...
      last = true

      for (let i = 0; i < N; i++) {
        // compute the taxicab norm of i-th column and row
        // TODO optimize for complex numbers
        let colNorm = realzero
        let rowNorm = realzero

        for (let j = 0; j < N; j++) {
          if (i === j) continue
          colNorm = addScalar(colNorm, abs(arr[j][i]))
          rowNorm = addScalar(rowNorm, abs(arr[i][j]))
        }

        if (!equal(colNorm, 0) && !equal(rowNorm, 0)) {
          // find integer power closest to balancing the matrix
          // (we want to scale only by integer powers of radix,
          // so that we don't lose any precision due to round-off)

          let f = realone
          let c = colNorm

          const rowDivRadix = divideScalar(rowNorm, radix)
          const rowMulRadix = multiplyScalar(rowNorm, radix)

          while (smaller(c, rowDivRadix)) {
            c = multiplyScalar(c, radixSq)
            f = multiplyScalar(f, radix)
          }
          while (larger(c, rowMulRadix)) {
            c = divideScalar(c, radixSq)
            f = divideScalar(f, radix)
          }

          // check whether balancing is needed
          // condition = (c + rowNorm) / f < 0.95 * (colNorm + rowNorm)
          const condition = smaller(divideScalar(addScalar(c, rowNorm), f), multiplyScalar(addScalar(colNorm, rowNorm), 0.95))

          // apply balancing similarity transformation
          if (condition) {
            // we should loop once again to check whether
            // another rebalancing is needed
            last = false

            const g = divideScalar(1, f)

            for (let j = 0; j < N; j++) {
              if (i === j) {
                continue
              }
              arr[i][j] = multiplyScalar(arr[i][j], g)
              arr[j][i] = multiplyScalar(arr[j][i], f)
            }

            // keep track of transformations
            if (findVectors) {
              Rdiag[i] = multiplyScalar(Rdiag[i], g)
            }
          }
        }
      }
    }

    // return the diagonal row transformation matrix
    return findVectors ? diag(Rdiag) : null
  }
@param {number[][]} arr
@param {number} N
@param {number} prec
@param {'number'|'BigNumber'|'Complex'} type
@returns {number[][]}