(Sorry if this is a double post, I sent this message the first time before completing registration on the mailing list.) Sometime ago I wrote a small set of functions for construction and use transformation matrices, including matrix inversion. Here is the code: https://github.com/lhartmann/openscad_m4lib/blob/master/m4.scad Matrix inversion code works most of the time, except I just detected a bug I can not solve. A simple dual rotation matrix, 90º around Y then 90º around Z, is yielding nan+inf on all elements of the inverse. Works nicely for other angles, though... M = m4rz(90) * m4ry(90); // Rotates around Y first, then Z Minv = m4inv(M); // Fails I managed to narrow it down to m4inv_solve(...) which does diagonal zero detection. This is the testing code I used to isolate the error: M = m4rz(90)*m4ry(90); // Test matrix echo( // Failing m4inv_solve(0, // Solve row 0: Detect diagonal zero, swap rows, calls solve2(...). m4inv_mx(M) // Expand to Row-Reduced Echelon Format ) ); echo( // Working by bypassing m4inv_solve(...) m4inv_solve2(0, // Solve row 0 (normalize, scale and subtract from all other rows) m4inv_rowswap(0,2, // Swap rows 0 and 2 by hand m4inv_mx(M) // Expand to Row-Reduced Echelon Format ) ) ); This is the code for the misbehaving function, it should check move zeroes away from the diagonal of row i: // Solve one step, ensure nonzero at diagonal element. function m4inv_zero(x) = abs(x) < 1e-5; function m4inv_solve(i,mx) = !m4inv_zero(mx[i][i ]) ? m4inv_solve2(i, mx ) : i+1<4 && !m4inv_zero(mx[i][i+1]) ? m4inv_solve2(i, m4inv_rowswap(i,i+1,mx)) : i+2<4 && !m4inv_zero(mx[i][i+2]) ? m4inv_solve2(i, m4inv_rowswap(i,i+2,mx)) : i+3<4 && !m4inv_zero(mx[i][i+3]) ? m4inv_solve2(i, m4inv_rowswap(i,i+3,mx)) : m4identity(); // Singular matrix, reset to identity I looked at the code for m4inv_solve(...) above for a while, but don't get why it is failing... Any ideas? -- Lucas Vinicius Hartmann Dizem que se você rodar o CD do Windows ao contrário ele mostra uma mensagem demoníaca... Mas isso nem é o pior, se você rodar ele normal ele instala o Windows!

Are you sure this code ever run properly on a matrix with 0 in diagonal? look at this trace: > mx = m4rx(90)*m4ry(90); > i=0; > echo(m4inv_solve(0,m4inv_mx(mx))); // first call evaluated in m4inv(mx) > echo(i+2<4 && !m4inv_zero(mx[i][i+2]) ? m4inv_solve2(i, > m4inv_rowswap(i,i+2,mx)):"not taken"); // echo(m4inv_rowswap(i,i+2,mx)); > // d > echo((i==0) ? mx[i]/mx[i][i]:"not taken"); > > Answer: > ECHO: [[0, 0, 1, 0], [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1]] // mx > ECHO: [[nan, inf, nan, nan, nan, nan, inf, nan], [nan, -inf, nan, nan, > nan, nan, -inf, nan], [nan, nan, nan, nan, nan, nan, nan, nan], [nan, nan, > nan, nan, nan, nan, nan, nan]] > ECHO: [[nan, inf, nan, nan], [nan, -inf, nan, nan], [nan, nan, nan, nan], > [nan, nan, nan, nan]] > ECHO: [[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]] // row swap > ECHO: [nan, nan, inf, nan] So, you test for a column !m4inv_zero(mx[i][i+2]) then do a row swap m4inv_rowswap(i,i+2,mx)) and then divide by 0 mx[0]/mx[0][0] You don't know, what you are "buying" by the row swap. You should test for the row with first element non-zero and swap this one -- View this message in context: http://forum.openscad.org/Matrix-inversion-tp17647p17648.html Sent from the OpenSCAD mailing list archive at Nabble.com.

Ups I pressed the wrong button ... I'd say this works better: > function m4inv_solve(i,mx) = > !m4inv_zero(mx[i][i ]) ? m4inv_solve2(i, mx > ) : > i+1<4 && !m4inv_zero(mx[i+1][i]) ? m4inv_solve2(i, > m4inv_rowswap(i,i+1,mx)) : > i+2<4 && !m4inv_zero(mx[i+2][i]) ? m4inv_solve2(i, > m4inv_rowswap(i,i+2,mx)) : > i+3<4 && !m4inv_zero(mx[i+3][i]) ? m4inv_solve2(i, > m4inv_rowswap(i,i+3,mx)) : > m4identity(); // Singular matrix, reset to identity -- View this message in context: http://forum.openscad.org/Matrix-inversion-tp17647p17649.html Sent from the OpenSCAD mailing list archive at Nabble.com.

  mx

You are absolutely right. Thanks! Somehow I ended up swapping rows and columns there... Fix commited to github. -- Lucas Vinicius Hartmann Dizem que se você rodar o CD do Windows ao contrário ele mostra uma mensagem demoníaca... Mas isso nem é o pior, se você rodar ele normal ele instala o Windows! 2016-06-11 21:14 GMT-03:00 Parkinbot <rudolf@parkinbot.com>: > Ups I pressed the wrong button ... > > I'd say this works better: > > > > function m4inv_solve(i,mx) = > > !m4inv_zero(mx[i][i ]) ? m4inv_solve2(i, > mx > > ) : > > i+1<4 && !m4inv_zero(mx[i+1][i]) ? m4inv_solve2(i, > > m4inv_rowswap(i,i+1,mx)) : > > i+2<4 && !m4inv_zero(mx[i+2][i]) ? m4inv_solve2(i, > > m4inv_rowswap(i,i+2,mx)) : > > i+3<4 && !m4inv_zero(mx[i+3][i]) ? m4inv_solve2(i, > > m4inv_rowswap(i,i+3,mx)) : > > m4identity(); // Singular matrix, reset to identity > > > > > > -- > View this message in context: > http://forum.openscad.org/Matrix-inversion-tp17647p17649.html > Sent from the OpenSCAD mailing list archive at Nabble.com. > > _______________________________________________ > OpenSCAD mailing list > Discuss@lists.openscad.org > http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org >