Раздел «Язык Ruby».InverseMatrix:

Inverting Matrix and the problem "Get rid of Irrationality"

I'd like to present you method for solving linear equations exemplified by solution of the following problem:

One has array of primes p and irrational number u:

u = 1 / (sqrt(p[0]) + sqrt(p[1]) + .. + sqrt(p[n]))

Write the program, that outputs pairs [a[i], q[i]] (a[i] — rational, q[i] — square-free positive integer > 1) such that

u == Sum_i (  a[i]*sqrt(q[i]) )

require 'set'

class Hash
  # add values of the second hash h to corresponding values of self
  def add(h)
    h.each_key do |k|
      self[k] ||= 0
      self[k] += h[k]
      self.delete(k) if self[k] == 0
    end
    self
  end
  alias :+ :add

  # multiply all values by v
  def times(v)
    self.each_key {|i|
      self[i] *= v 
    }
    self
  end
  alias :* :times
end

# Solves matrix equation m * x = b
#
# * m is Hash of Hashes : row_key -> column_key -> value;
# * b is Hash row_key -> value;   
# 
# default value is 0;
#
# It finds such Hash x col_key->value that:
#  calculations
#     answer = {}
#     m.each_key { |k| x.each_key { |c| answer[k] += m[k][c] * x[c] } }
#  brings to
#     answer == b
# 
#  OR 
#   
#   sum_c ( m[k][c] * x[c] ) == b[k]  for any k
# 
def matrix_solve(m, b)
  m = Marshal.load(Marshal.dump(m))
  col_keys = Set.new
  m.each {|k, row| col_keys += row.keys }

  m.each_key do |k|
    m[k][:res] = b[k] 
  end
  
  row_keys = m.keys 
  rest_row_keys = row_keys.dup
  row_main_key = {}
  
  row_keys.each do |k|
    row = m[k]
    c = row.keys.find {|c| col_keys.include?(c) && row[c] != 0}
    if c.nil?
      m.delete(c)
      next
    end
    col_keys.delete(c)
    row_main_key[k] = c
    rest_row_keys.delete(k)
    rest_row_keys.each do |k2|
      z = - m[k2][c] / row[c]
      # puts "k=#{k2} c=#{c} z=#{z}"
      m[k2].+( row.dup.*( z ) ) if z != 0
    end
  end
  
  done_row_keys = Set.new
  row_keys = m.keys
  row_keys.each do |k|
    row = m[k]
    c = row_main_key[k]
    next unless c
    done_row_keys.each do |k2|
      z = - m[k2][c] / row[c]
      # puts "k=#{k2} c=#{c} z=#{z}"
      m[k2].+( row.dup.*( z ) ) if z
    end
    done_row_keys << k
  end
  
  # puts m.inspect

  x = {}
  row_keys.each {|k| 
    row = m[k]
    next unless row_main_key[k] # this row should have only zeros
                                # no solution if it is not the case
    next if row[:res] == 0
    x[ row_main_key[k] ] = row[:res] / row[row_main_key[k]]
  }
  x
end

p = [2, 3, 5, 7, 11]

m = p.size

matrix = Hash.new

(2**m).times do |a|
  matrix[a] = Hash.new(0)
end

(2**m).times do |a|
  m.times {|k|
    if (a & 2**k) != 0
      matrix[a][a ^ 2**k] += Rational(1, 1)
    else
      matrix[a][a ^ 2**k] += Rational(p[k], 1)
    end
  }
end

v = Hash.new(0)
v[0] = 1

x = matrix_solve(matrix, v)

#~ puts x.inspect

#~ answer = Hash.new(0)
#~ matrix.each_key {|k| 
  #~ x.each_key {|c| answer[k] += m[k][c] * x[c]  }  
#~ }
#~ answer.delete_if {|k,v| v == 0}
#~ puts answer.inspect


x.each {|a,v|
  next if v == 0
  x  = 1
  a.to_s(2).split('').reverse.each_with_index {|c,i|
    x *= p[i] if c == '1'
  }
  puts "#{v}  #{x}"
}

Output

For p = [2, 3, 5, 7, 11] we have the following answer:

-1061/43684  11
1105/43684  165
1823/43684  42
73/43684  385
3767/21842  2
-3863/43684  30
4883/43684  3
-611/43684  70
7399/43684  7
-533/43684  66
-583/43684  105
-1577/43684  154
147/21842  2310
-5231/43684  5
-1183/43684  231
1233/43684  110

For p = [2, 3, 5, 7, 11, 13, 17] we have the following answer:

-945544276436469543/228263782235476891276  1547
5937473619309020545/228263782235476891276  195
-896653270908170123/228263782235476891276  14586
-1574307149774503251/114131891117738445638  286
6897175279415436589/684791346706430673828  561
-541971010737977319/228263782235476891276  13090
14759133188785854189/114131891117738445638  11
-662454149240556559/228263782235476891276  4290
-3163069258783021867/228263782235476891276  455
449962608194253617/228263782235476891276  34034
-2604642531419964613/228263782235476891276  1309
7921707362713183465/684791346706430673828  165
665920702988310091/114131891117738445638  42
-806811056258973599/342395673353215336914  23205
-2341313676396337153/228263782235476891276  385
88813429590514381/228263782235476891276  10010
-51175043790150726/57065945558869222819  19635
31601067011694235/228263782235476891276  510510
-588353868840492181/228263782235476891276  429
130150413554256224/171197836676607668457  51051
1157053533198150823/228263782235476891276  1001
2092645848147654423/228263782235476891276  1105
10096628561798003/342395673353215336914  102
-7404241398752938521/57065945558869222819  2
308209841071143197/171197836676607668457  15015
1409267941341249957/228263782235476891276  935
-117945279993497829/228263782235476891276  24310
1699769205220083361/114131891117738445638  238
-25836076827031874593/342395673353215336914  30
1035661130382208929/228263782235476891276  3570
5950238305401736333/114131891117738445638  70
-1099662493030098627/228263782235476891276  2431
-25453081784846770819/114131891117738445638  3
3726356605381096381/114131891117738445638  78
-2163848886010771285/684791346706430673828  9282
-159781020136952847/114131891117738445638  36465
8225562984004498157/342395673353215336914  66
-3344926344876935899/114131891117738445638  182
4961647667308648387/228263782235476891276  357
16466840581643829853/114131891117738445638  7
-4108851071566082775/228263782235476891276  715
-500125974117969973/228263782235476891276  7854
5678092069761145383/228263782235476891276  105
-3623315758862814747/114131891117738445638  170
-706006947279665473/684791346706430673828  2730
-1710849584132074109/114131891117738445638  154
207226889708284983/114131891117738445638  85085
-3358377496875766513/228263782235476891276  273
-6454357975943625757/114131891117738445638  17
1695351194574441861/114131891117738445638  442
-1475748604255059689/228263782235476891276  2310
-2145349183467471725/114131891117738445638  5
-9941487793521642281/684791346706430673828  255
-4888251239928649387/228263782235476891276  231
4009681137113992301/684791346706430673828  6006
2139744463550751069/114131891117738445638  130
493761179381222801/114131891117738445638  374
499336508596894659/228263782235476891276  6630
4841583417501407971/114131891117738445638  13
1276070762737375901/228263782235476891276  5610
-153168318982594399/228263782235476891276  595
3174357322245559621/114131891117738445638  110
-388116727504945521/228263782235476891276  15470
1565337485411589237/228263782235476891276  663

C++ version

If you have C++ solution for this problem, you are welcome. Send it to me ( artem.voroztsov // gmail.com ) and I'll place it here with your name.

-- ArtemVoroztsov - 18 Oct 2007