1_solve.rb

def solve(a)
  i = 0
  last_error = 0

  begin
    i += 1
    errors1 = [
      2182449 - (a[0][0] + a[0][1] + a[0][2]),
      893841 - (a[1][0] + a[1][1] + a[1][2]),
      1307065 - (a[2][0] + a[2][1] + a[2][2])
    ]

    errors2 = [
      3466855 - (a[0][0] + a[1][0] + a[2][0]),
      873224 - (a[0][1] + a[1][1] + a[2][1]),
      52389 - (a[0][2] + a[1][2] + a[2][2])
    ]

    errors = Array.new

    max = 0
    coords = nil

    [0, 1, 2].repeated_permutation(2).each do |i, j|
      errors[i] ||= Array.new
      errors[i][j] = errors1[i] + errors2[j]
      if max.abs <= errors[i][j].abs
        max = errors[i][j]
        coords = [i, j]
      end
    end

    a[coords[0]][coords[1]] += max * 0.05
    a[coords[0]][coords[1]] = 0 if a[coords[0]][coords[1]] < 0 # clip

    error = errors.flatten.map(&:abs).inject(&:+)
    delta = (last_error - error).abs
    last_error = error

  end while delta > 0.0001

  [error, i, a]
end

100000.times do
  a = [
    [Random.rand(3000000), Random.rand(600000), Random.rand(30000)],
    [Random.rand(3000000), Random.rand(600000), Random.rand(30000)],
    [Random.rand(3000000), Random.rand(600000), Random.rand(30000)],
  ]

  error, iterations, solution = solve(a)

  next if error > 0.1
  print solution.flatten.map(&:to_i).join(',')
  print ','
  print error
  print ','
  print iterations
  puts
end

2_plot.R

solutions <- read.csv("~/solutions.csv", header=FALSE)

options(scipen=5)

par(mfrow=c(3,3))
hist(solutions$V1, xlim=c(0,2500000), prob=T, xlab=NULL, ylab=NULL, main=NULL)
hist(solutions$V2, xlim=c(0,1000000), prob=T, xlab=NULL, ylab=NULL, main=NULL)
hist(solutions$V3, xlim=c(0,50000), prob=T, xlab=NULL, ylab=NULL, main=NULL)

hist(solutions$V4, xlim=c(0,2500000), prob=T, xlab=NULL, ylab=NULL, main=NULL)
hist(solutions$V5, xlim=c(0,1000000), prob=T, xlab=NULL, ylab=NULL, main=NULL)
hist(solutions$V6, xlim=c(0,50000), prob=T, xlab=NULL, ylab=NULL, main=NULL)

hist(solutions$V7, xlim=c(0,2500000), prob=T, xlab=NULL, ylab=NULL, main=NULL)
hist(solutions$V8, xlim=c(0,1000000), prob=T, xlab=NULL, ylab=NULL, main=NULL)
hist(solutions$V9, xlim=c(0,50000), prob=T, xlab=NULL, ylab=NULL, main=NULL)