<p>I need to programmatically solve a system of linear equations in C, Objective C, or (if needed) C++.</p>

<p>Here's an example of the equations:</p>

<pre><code>-44.3940 = a * 50.0 + b * 37.0 + tx
-45.3049 = a * 43.0 + b * 39.0 + tx
-44.9594 = a * 52.0 + b * 41.0 + tx
</code></pre>

<p>From this, I'd like to get the best approximation for <code>a</code>, <code>b</code>, and <code>tx</code>.</p>

[{'answer_id': 774, 'author': 'nlucaroni', 'author_id': 157, 'author_profile': 'https://Stackoverflow.com/users/157', 'pm_score': 2, 'selected': False, 'text': '<p>Are you looking for a software package that\'ll do the work or actually doing the matrix operations and such and do each step?</p>\n\n<p>The the first, a coworker of mine just used <a href="http://ocaml-glpk.sourceforge.net/" rel="nofollow noreferrer" title="Ocaml GLPK">Ocaml GLPK</a>. It is just a wrapper for the <a href="http://www.gnu.org/software/glpk/" rel="nofollow noreferrer" title="GNU Linear Programming Kit">GLPK</a>, but it removes a lot of the steps of setting things up. It looks like you\'re going to have to stick with the GLPK, in C, though. For the latter, thanks to delicious for saving an old article I used to learn LP awhile back, <a href="http://www.math.ucla.edu/~tom/LP.pdf" rel="nofollow noreferrer" title="Linear Programming">PDF</a>. If you need specific help setting up further, let us know and I\'m sure, me or someone will wander back in and help, but, I think it\'s fairly straight forward from here. Good Luck! </p>\n'}, {'answer_id': 781, 'author': 'Brian Jorgensen', 'author_id': 229, 'author_profile': 'https://Stackoverflow.com/users/229', 'pm_score': 5, 'selected': True, 'text': '<p><a href="http://en.wikipedia.org/wiki/Cramers_rule" rel="noreferrer">Cramer\'s Rule</a>\nand\n<a href="http://en.wikipedia.org/wiki/Gaussian_elimination" rel="noreferrer">Gaussian Elimination</a>\nare two good, general-purpose algorithms (also see <a href="http://en.wikipedia.org/wiki/Simultaneous_linear_equations" rel="noreferrer">Simultaneous Linear Equations</a>). If you\'re looking for code, check out <a href="http://en.wikipedia.org/wiki/GiNaC" rel="noreferrer">GiNaC</a>, <a href="http://maxima.sourceforge.net/" rel="noreferrer">Maxima</a>, and <a href="http://issc.uj.ac.za/symbolic/symbolic.html" rel="noreferrer">SymbolicC++</a> (depending on your licensing requirements, of course).</p>\n\n<p>EDIT: I know you\'re working in C land, but I also have to put in a good word for <a href="http://code.google.com/p/sympy/" rel="noreferrer">SymPy</a> (a  computer algebra system in Python). You can learn a lot from its algorithms (if you can read a bit of python). Also, it\'s under the new BSD license, while most of the free math packages are GPL.</p>\n'}, {'answer_id': 6206, 'author': 'Community', 'author_id': -1, 'author_profile': 'https://Stackoverflow.com/users/-1', 'pm_score': 3, 'selected': False, 'text': '<p>For a 3x3 system of linear equations I guess it would be okay to roll out your own algorithms. </p>\n\n<p>However, you might have to worry about accuracy, division by zero or really small numbers and what to do about infinitely many solutions. My suggestion is to go with a standard numerical linear algebra package such as <a href="http://en.wikipedia.org/wiki/Lapack" rel="nofollow noreferrer">LAPACK</a>. </p>\n'}, {'answer_id': 26578, 'author': 'Frank Krueger', 'author_id': 338, 'author_profile': 'https://Stackoverflow.com/users/338', 'pm_score': 2, 'selected': False, 'text': '<p>Personally, I\'m partial to the algorithms of <a href="http://www.nr.com/" rel="nofollow noreferrer">Numerical Recipes</a>. (I\'m fond of the C++ edition.)</p>\n\n<p>This book will teach you why the algorithms work, plus show you some pretty-well debugged implementations of those algorithms.</p>\n\n<p>Of course, you could just blindly use <a href="http://www.netlib.org/clapack/" rel="nofollow noreferrer">CLAPACK</a> (I\'ve used it with great success), but I would first hand-type a Gaussian Elimination algorithm to at least have a faint idea of the kind of work that has gone into making these algorithms stable.</p>\n\n<p>Later, if you\'re doing more interesting linear algebra, looking around the source code of <a href="http://www.gnu.org/software/octave/" rel="nofollow noreferrer">Octave</a> will answer a lot of questions.</p>\n'}, {'answer_id': 26637, 'author': 'Baltimark', 'author_id': 1179, 'author_profile': 'https://Stackoverflow.com/users/1179', 'pm_score': 2, 'selected': False, 'text': '<p><a href="http://math.nist.gov/tnt/download.html" rel="nofollow noreferrer">Template Numerical Toolkit</a> from NIST has tools for doing that. </p>\n\n<p>One of the more reliable ways is to use a <a href="http://en.wikipedia.org/wiki/QR_decomposition" rel="nofollow noreferrer">QR Decomposition</a>. </p>\n\n<p>Here\'s an example of a wrapper so that I can call "GetInverse(A, InvA)" in my code and it will put the inverse into InvA.</p>\n\n<pre><code>void GetInverse(const Array2D&lt;double&gt;&amp; A, Array2D&lt;double&gt;&amp; invA)\n   {\n   QR&lt;double&gt; qr(A);  \n   invA = qr.solve(I); \n   }\n</code></pre>\n\n<p>Array2D is defined in the library. </p>\n'}, {'answer_id': 83758, 'author': 'David Nehme', 'author_id': 14167, 'author_profile': 'https://Stackoverflow.com/users/14167', 'pm_score': 2, 'selected': False, 'text': '<p>From the wording of your question, it seems like you have more equations than unknowns and you want to minimize the inconsistencies.  This is typically done with linear regression, which minimizes the sum of the squares of the inconsistencies.  Depending on the size of the data, you can do this in a spreadsheet or in a statistical package.  R is a high-quality, free package that does linear regression, among a lot of other things.  There is a lot to linear regression (and a lot of gotcha\'s), but as it\'s straightforward to do for simple cases.  Here\'s an R example using your data.  Note that the "tx" is the intercept to your model.</p>\n\n<pre><code>&gt; y &lt;- c(-44.394, -45.3049, -44.9594)\n&gt; a &lt;- c(50.0, 43.0, 52.0)\n&gt; b &lt;- c(37.0, 39.0, 41.0)\n&gt; regression = lm(y ~ a + b)\n&gt; regression\n\nCall:\nlm(formula = y ~ a + b)\n\nCoefficients:\n(Intercept)            a            b  \n  -41.63759      0.07852     -0.18061  \n</code></pre>\n'}, {'answer_id': 99339, 'author': 'Thelema', 'author_id': 12874, 'author_profile': 'https://Stackoverflow.com/users/12874', 'pm_score': 2, 'selected': False, 'text': '<p>In terms of run-time efficiency, others have answered better than I.  If you always will have the same number of equations as variables, I like <a href="http://en.wikipedia.org/wiki/Cramer&#39;s_rule" rel="nofollow noreferrer">Cramer\'s rule</a> as it\'s easy to implement.  Just write a function to calculate determinant of a matrix (or use one that\'s already written, I\'m sure you can find one out there), and divide the determinants of two matrices.</p>\n'}, {'answer_id': 590060, 'author': 'paxdiablo', 'author_id': 14860, 'author_profile': 'https://Stackoverflow.com/users/14860', 'pm_score': 4, 'selected': False, 'text': '<p>You can solve this with a program exactly the same way you solve it by hand (with multiplication and subtraction, then feeding results back into the equations). This is pretty standard secondary-school-level mathematics.</p>\n\n<pre><code>-44.3940 = 50a + 37b + c (A)\n-45.3049 = 43a + 39b + c (B)\n-44.9594 = 52a + 41b + c (C)\n\n(A-B): 0.9109 =  7a -  2b (D)\n(B-C): 0.3455 = -9a -  2b (E)\n\n(D-E): 1.2564 = 16a (F)\n\n(F/16):  a = 0.078525 (G)\n\nFeed G into D:\n       0.9109 = 7a - 2b\n    =&gt; 0.9109 = 0.549675 - 2b (substitute a)\n    =&gt; 0.361225 = -2b (subtract 0.549675 from both sides)\n    =&gt; -0.1806125 = b (divide both sides by -2) (H)\n\nFeed H/G into A:\n       -44.3940 = 50a + 37b + c\n    =&gt; -44.3940 = 3.92625 - 6.6826625 + c (substitute a/b)\n    =&gt; -41.6375875 = c (subtract 3.92625 - 6.6826625 from both sides)\n</code></pre>\n\n<p>So you end up with:</p>\n\n<pre><code>a =   0.0785250\nb =  -0.1806125\nc = -41.6375875\n</code></pre>\n\n<p>If you plug these values back into A, B and C, you\'ll find they\'re correct.</p>\n\n<p>The trick is to use a simple 4x3 matrix which reduces in turn to a 3x2 matrix, then a 2x1 which is "a = n", n being an actual number. Once you have that, you feed it into the next matrix up to get another value, then those two values into the next matrix up until you\'ve solved all variables.</p>\n\n<p>Provided you have N distinct equations, you can always solve for N variables. I say distinct because these two are not:</p>\n\n<pre><code> 7a + 2b =  50\n14a + 4b = 100\n</code></pre>\n\n<p>They are the <em>same</em> equation multiplied by two so you cannot get a solution from them - multiplying the first by two then subtracting leaves you with the true but useless statement:</p>\n\n<pre><code>0 = 0 + 0\n</code></pre>\n\n<hr>\n\n<p>By way of example, here\'s some C code that works out the simultaneous equations that you\'re placed in your question. First some necessary types, variables, a support function for printing out an equation, and the start of <code>main</code>:</p>\n\n<pre><code>#include &lt;stdio.h&gt;\n\ntypedef struct { double r, a, b, c; } tEquation;\ntEquation equ1[] = {\n    { -44.3940,  50, 37, 1 },      // -44.3940 = 50a + 37b + c (A)\n    { -45.3049,  43, 39, 1 },      // -45.3049 = 43a + 39b + c (B)\n    { -44.9594,  52, 41, 1 },      // -44.9594 = 52a + 41b + c (C)\n};\ntEquation equ2[2], equ3[1];\n\nstatic void dumpEqu (char *desc, tEquation *e, char *post) {\n    printf ("%10s: %12.8lf = %12.8lfa + %12.8lfb + %12.8lfc (%s)\\n",\n        desc, e-&gt;r, e-&gt;a, e-&gt;b, e-&gt;c, post);\n}\n\nint main (void) {\n    double a, b, c;\n</code></pre>\n\n<p>Next, the reduction of the three equations with three unknowns to two equations with two unknowns:</p>\n\n<pre><code>    // First step, populate equ2 based on removing c from equ.\n\n    dumpEqu ("&gt;", &amp;(equ1[0]), "A");\n    dumpEqu ("&gt;", &amp;(equ1[1]), "B");\n    dumpEqu ("&gt;", &amp;(equ1[2]), "C");\n    puts ("");\n\n    // A - B\n    equ2[0].r = equ1[0].r * equ1[1].c - equ1[1].r * equ1[0].c;\n    equ2[0].a = equ1[0].a * equ1[1].c - equ1[1].a * equ1[0].c;\n    equ2[0].b = equ1[0].b * equ1[1].c - equ1[1].b * equ1[0].c;\n    equ2[0].c = 0;\n\n    // B - C\n    equ2[1].r = equ1[1].r * equ1[2].c - equ1[2].r * equ1[1].c;\n    equ2[1].a = equ1[1].a * equ1[2].c - equ1[2].a * equ1[1].c;\n    equ2[1].b = equ1[1].b * equ1[2].c - equ1[2].b * equ1[1].c;\n    equ2[1].c = 0;\n\n    dumpEqu ("A-B", &amp;(equ2[0]), "D");\n    dumpEqu ("B-C", &amp;(equ2[1]), "E");\n    puts ("");\n</code></pre>\n\n<p>Next, the reduction of the two equations with two unknowns to one equation with one unknown:</p>\n\n<pre><code>    // Next step, populate equ3 based on removing b from equ2.\n\n    // D - E\n    equ3[0].r = equ2[0].r * equ2[1].b - equ2[1].r * equ2[0].b;\n    equ3[0].a = equ2[0].a * equ2[1].b - equ2[1].a * equ2[0].b;\n    equ3[0].b = 0;\n    equ3[0].c = 0;\n\n    dumpEqu ("D-E", &amp;(equ3[0]), "F");\n    puts ("");\n</code></pre>\n\n<p>Now that we have a formula of the type <code>number1 = unknown * number2</code>, we can simply work out the unknown value with <code>unknown &lt;- number1 / number2</code>. Then, once you\'ve figured that value out, substitute it into one of the equations with two unknowns and work out the second value. Then substitute both those (now-known) unknowns into one of the original equations and you now have the values for all three unknowns:</p>\n\n<pre><code>    // Finally, substitute values back into equations.\n\n    a = equ3[0].r / equ3[0].a;\n    printf ("From (F    ), a = %12.8lf (G)\\n", a);\n\n    b = (equ2[0].r - equ2[0].a * a) / equ2[0].b;\n    printf ("From (D,G  ), b = %12.8lf (H)\\n", b);\n\n    c = (equ1[0].r - equ1[0].a * a - equ1[0].b * b) / equ1[0].c;\n    printf ("From (A,G,H), c = %12.8lf (I)\\n", c);\n\n    return 0;\n}\n</code></pre>\n\n<p>The output of that code matches the earlier calculations in this answer:</p>\n\n<pre><code>         &gt;: -44.39400000 =  50.00000000a +  37.00000000b +   1.00000000c (A)\n         &gt;: -45.30490000 =  43.00000000a +  39.00000000b +   1.00000000c (B)\n         &gt;: -44.95940000 =  52.00000000a +  41.00000000b +   1.00000000c (C)\n\n       A-B:   0.91090000 =   7.00000000a +  -2.00000000b +   0.00000000c (D)\n       B-C:  -0.34550000 =  -9.00000000a +  -2.00000000b +   0.00000000c (E)\n\n       D-E:  -2.51280000 = -32.00000000a +   0.00000000b +   0.00000000c (F)\n\nFrom (F    ), a =   0.07852500 (G)\nFrom (D,G  ), b =  -0.18061250 (H)\nFrom (A,G,H), c = -41.63758750 (I)\n</code></pre>\n'}, {'answer_id': 770947, 'author': 'Bobby Ortiz', 'author_id': 25843, 'author_profile': 'https://Stackoverflow.com/users/25843', 'pm_score': 3, 'selected': False, 'text': '<p>Take a look at the <a href="http://code.msdn.microsoft.com/solverfoundation" rel="noreferrer">Microsoft Solver Foundation</a>.</p>\n\n<p>With it you could write code like this:</p>\n\n<pre><code>  SolverContext context = SolverContext.GetContext();\n  Model model = context.CreateModel();\n\n  Decision a = new Decision(Domain.Real, "a");\n  Decision b = new Decision(Domain.Real, "b");\n  Decision c = new Decision(Domain.Real, "c");\n  model.AddDecisions(a,b,c);\n  model.AddConstraint("eqA", -44.3940 == 50*a + 37*b + c);\n  model.AddConstraint("eqB", -45.3049 == 43*a + 39*b + c);\n  model.AddConstraint("eqC", -44.9594 == 52*a + 41*b + c);\n  Solution solution = context.Solve();\n  string results = solution.GetReport().ToString();\n  Console.WriteLine(results); \n</code></pre>\n\n<p><strong>Here is the output:</strong><br/>\n===Solver Foundation Service Report===<br/>\nDatetime: 04/20/2009 23:29:55<br/>\nModel Name: Default<br/>\nCapabilities requested: LP<br/>\nSolve Time (ms): 1027<br/>\nTotal Time (ms): 1414<br/>\nSolve Completion Status: Optimal<br/>\nSolver Selected: Microsoft.SolverFoundation.Solvers.SimplexSolver<br/>\nDirectives:<br/>\nMicrosoft.SolverFoundation.Services.Directive<br/>\nAlgorithm: Primal<br/>\nArithmetic: Hybrid<br/>\nPricing (exact): Default<br/>\nPricing (double): SteepestEdge<br/>\nBasis: Slack<br/>\nPivot Count: 3<br/>\n===Solution Details===<br/>\nGoals:<br/>\n<br/>\nDecisions:<br/>\na: 0.0785250000000004<br/>\nb: -0.180612500000001<br/>\nc: -41.6375875<br/></p>\n'}, {'answer_id': 27708081, 'author': 'Bulent S.', 'author_id': 4405846, 'author_profile': 'https://Stackoverflow.com/users/4405846', 'pm_score': 1, 'selected': False, 'text': '<pre><code>function x = LinSolve(A,y)\n%\n% Recursive Solution of Linear System Ax=y\n% matlab equivalent: x = A\\y \n% x = n x 1\n% A = n x n\n% y = n x 1\n% Uses stack space extensively. Not efficient.\n% C allows recursion, so convert it into C. \n% ----------------------------------------------\nn=length(y);\nx=zeros(n,1);\nif(n&gt;1)\n    x(1:n-1,1) = LinSolve( A(1:n-1,1:n-1) - (A(1:n-1,n)*A(n,1:n-1))./A(n,n) , ...\n                           y(1:n-1,1) - A(1:n-1,n).*(y(n,1)/A(n,n))); \n    x(n,1) = (y(n,1) - A(n,1:n-1)*x(1:n-1,1))./A(n,n); \nelse\n    x = y(1,1) / A(1,1);\nend\n</code></pre>\n'}, {'answer_id': 69017457, 'author': 'Hussain Chitalwala', 'author_id': 11280050, 'author_profile': 'https://Stackoverflow.com/users/11280050', 'pm_score': -1, 'selected': False, 'text': '<p>For general cases, you could use python along with numpy for Gaussian elimination. And then plug in values and get the remaining values.</p>\n'}]
['https://Stackoverflow.com/questions/769', 'https://Stackoverflow.com', 'https://Stackoverflow.com/users/79/']