title: Linear system solution with fractions in numpy
body: <p>I have matrix <code>A</code> and a right-hand-side vector <code>y</code> expressed in terms of <code>fractions.Fraction</code> objects:</p>
  
  <pre><code>import random, fractions, numpy as np
  
  A = np.zeros((3, 3), dtype=fractions.Fraction)
  y = np.zeros((3, 1), dtype=fractions.Fraction)
  for i in range(3):
      for j in range(3):
          A[i, j] = fractions.Fraction(np.random.randint(0, 4), np.random.randint(1, 6))
      y[i] = fractions.Fraction(np.random.randint(0, 4), np.random.randint(1, 6))
  </code></pre>
  
  <p>I would like to solve the system <code>A*x = y</code> using the provided functions in <code>numpy</code> and get a result expressed in fraction objects, but unfortunately the basic <code>x = np.linalg.solve(A, y)</code> returns the result in standard floating point values:</p>
  
  <pre><code>&gt;&gt;&gt; np.linalg.solve(A, y)
  
  array([[-1.5245283 ],
         [ 2.36603774],
         [ 0.56352201]])
  </code></pre>
  
  <p>Is there a way of getting the exact result with fraction objects?</p>
  
  <hr>
  
  <p>EDIT</p>
  
  <p>What I would like to do is just not feasible with the built-in functionalities of numpy (as of version 1.10 - see Mad Physicist's answer). What one could do is implementing his/her own linear solver based on Gauss elimination, which relies on sum, subtraction, multiplication and division, all of which are well-defined and executed exactly with fraction objects (as long as the numerators and denominators fit in the data type, which I think is arbitrarily long).</p>
  
  <p>If you are really interested in having this, just implement a solver yourself, it will be easy and fast to do (follow one of the many tutorials online). I'm not that much interested, so I will stick to the floating point result.</p>