|
Page 1 of 1
|
[ 9 posts ] |
|
recursive calculation of matrix functions (eg. determinants)
Author |
Message |
Ford Prefect
Guru
Joined: Sat Mar 01, 2008 12:52 pm Posts: 1030
|
 recursive calculation of matrix functions (eg. determinants)
hey, I'm curious how to port the following recursive C code into RobotC (exemplarily: the determinant of square matrices):  |  |  |  | Code: // matrix determinant
double MatrixDet(int N, double A[N][N]) // or float instead { int i,j,i_count,j_count, count=0; double Asub[N-1][N-1], det=0;
if(N==1) return A[0][0]; if(N==2) return (A[0][0]*A[1][1] - A[0][1]*A[1][0]);
for(count=0; count<N; count++) { i_count=0; for(i=1; i<N; i++) { j_count=0; for(j=0; j<N; j++) { if(j == count) continue; Asub[i_count][j_count] = A[i][j]; j_count++; } i_count++; } det += pow(-1, count) * A[0][count] * MatrixDet(N-1,Asub); } return det; }
int main() { double A[6][6]; // or float instead A[0][0]=1; A[0][1]=2; A[0][2]=3; A[0][3]=4; A[0][4]=5; A[0][5]=6;
A[1][0]=10; A[1][1]=11; A[1][2]=12; A[1][3]=13; A[1][4]=14; A[1][5]=15;
A[2][0]=19; A[2][1]=20; A[2][2]=21; A[2][3]=22; A[2][4]=23; A[2][5]=24;
A[3][0]=32; A[3][1]=33; A[3][2]=34; A[3][3]=35; A[3][4]=36; A[3][5]=37;
A[4][0]=39; A[4][1]=40; A[4][2]=41; A[4][3]=42; A[4][4]=43; A[4][5]=44;
A[5][0]=51; A[5][1]=52; A[5][2]=53; A[5][3]=54; A[5][4]=55; A[5][5]=56;
int n=6; double det=MatrixDet(n, A);
printf("%9.3f\n", det); return 1; }
|  |  |  |  |
_________________ regards, HaWe aka Ford #define S sqrt(t+2*i*i)<2 #define F(a,b) for(a=0;a<b;++a) float x,y,r,i,s,j,t,n;task main(){F(y,64){F(x,99){r=i=t=0;s=x/33-2;j=y/32-1;F(n,50&S){t=r*r-i*i;i=2*r*i+j;r=t+s;}if(S){PutPixel(x,y);}}}while(1)}
|
Tue Dec 31, 2013 4:47 am |
|
 |
mightor
Site Admin
Joined: Wed Mar 05, 2008 8:14 am Posts: 3654 Location: Rotterdam, The Netherlands
|
 Re: recursive calculation of matrix functions (eg. determina
You need to find a way to get rid of the runtime sized array and then it would be fine.
= Xander
_________________| Professional Conduit of Reasonableness| (Title bestowed upon on the 8th day of November, 2013) | My Blog: I'd Rather Be Building Robots| ROBOTC 3rd Party Driver Suite: [ Project Page]
|
Tue Dec 31, 2013 5:43 am |
|
 |
Ford Prefect
Guru
Joined: Sat Mar 01, 2008 12:52 pm Posts: 1030
|
 Re: recursive calculation of matrix functions (eg. determina
yes, that's the point. IIRC: you already tried it, too, didn't you? But you can't "get rid of it", it's simply the way the methematical algorithm works.
ps, the recursive computation of submatrices is not only an issue for determinants but also for matrix adjugates or inverse matrices !
_________________ regards, HaWe aka Ford #define S sqrt(t+2*i*i)<2 #define F(a,b) for(a=0;a<b;++a) float x,y,r,i,s,j,t,n;task main(){F(y,64){F(x,99){r=i=t=0;s=x/33-2;j=y/32-1;F(n,50&S){t=r*r-i*i;i=2*r*i+j;r=t+s;}if(S){PutPixel(x,y);}}}while(1)}
|
Tue Dec 31, 2013 5:49 am |
|
 |
mightor
Site Admin
Joined: Wed Mar 05, 2008 8:14 am Posts: 3654 Location: Rotterdam, The Netherlands
|
 Re: recursive calculation of matrix functions (eg. determina
Well, I didn't try very hard  Mostly because I don't know enough about the algorithm or the actual operation that this is. I think there may be a way to do this, I'll see what I can come up with. = Xander
_________________| Professional Conduit of Reasonableness| (Title bestowed upon on the 8th day of November, 2013) | My Blog: I'd Rather Be Building Robots| ROBOTC 3rd Party Driver Suite: [ Project Page]
|
Tue Dec 31, 2013 9:01 am |
|
 |
mightor
Site Admin
Joined: Wed Mar 05, 2008 8:14 am Posts: 3654 Location: Rotterdam, The Netherlands
|
 Re: recursive calculation of matrix functions (eg. determina
Does this do what it's supposed to?  |  |  |  | Code: #define MAX_N 6 float MatrixDet(int N, float *A) // or float instead { writeDebugStreamLine("MatrixDet: N=%d, A: %p", N, A); int i,j,i_count,j_count, count=0; float Asub[MAX_N * MAX_N]; float det = 0;
if(N==1) return A[0]; if(N==2) return (A[0*N+0]*A[1*N+1] - A[0*N+1]*A[1*N+0]);
for(count=0; count<N; count++) { i_count=0; for(i=1; i<N; i++) { j_count=0; for(j=0; j<N; j++) { if(j == count) continue; Asub[i*(N-1)+j] = A[i*N+j]; j_count++; } i_count++; } det += pow(-1, count) * A[0*N+count] * MatrixDet(N-1, &Asub[0]); } return det; }
|  |  |  |  |
= Xander
_________________| Professional Conduit of Reasonableness| (Title bestowed upon on the 8th day of November, 2013) | My Blog: I'd Rather Be Building Robots| ROBOTC 3rd Party Driver Suite: [ Project Page]
|
Tue Dec 31, 2013 9:27 am |
|
 |
Ford Prefect
Guru
Joined: Sat Mar 01, 2008 12:52 pm Posts: 1030
|
 Re: recursive calculation of matrix functions (eg. determina
Do I see it correctly that you have to linearize the matrix? what result do you get with the 6x6 matrix above? --- and then as a test this one:
_________________ regards, HaWe aka Ford #define S sqrt(t+2*i*i)<2 #define F(a,b) for(a=0;a<b;++a) float x,y,r,i,s,j,t,n;task main(){F(y,64){F(x,99){r=i=t=0;s=x/33-2;j=y/32-1;F(n,50&S){t=r*r-i*i;i=2*r*i+j;r=t+s;}if(S){PutPixel(x,y);}}}while(1)}
|
Tue Dec 31, 2013 10:44 am |
|
 |
Ford Prefect
Guru
Joined: Sat Mar 01, 2008 12:52 pm Posts: 1030
|
 Re: recursive calculation of matrix functions (eg. determina
for test reasons: this is the program for the EV3:  |  |  |  | Code: // hw MatrixLib
#include <stdio.h> #include <math.h> #include <fcntl.h> #include <string.h> #include <sys/mman.h> #include <sys/ioctl.h>
#include "ev3_button.h" #include "ev3_lcd.h" #include "ev3_constants.h" #include "ev3_command.h" #include "ev3_timer.h" #include "ev3_lcd.h" #include "ev3_sound.h" #include "ev3_output.h"
double MatrixDet(int N, double A[N][N]) // or float instead { int i,j,i_count,j_count, count=0; double Asub[N-1][N-1], det=0;
if(N==1) return A[0][0]; if(N==2) return (A[0][0]*A[1][1] - A[0][1]*A[1][0]);
for(count=0; count<N; count++) { i_count=0; for(i=1; i<N; i++) { j_count=0; for(j=0; j<N; j++) { if(j == count) continue; Asub[i_count][j_count] = A[i][j]; j_count++; } i_count++; } det += pow(-1, count) * A[0][count] * MatrixDet(N-1,Asub); } return det; }
int main() {
ButtonLedInit(); LcdInit(); LcdClean();
double A[6][6]; // or float instead
A[0][0]=-1; A[0][1]= 2; A[0][2]=-3; A[0][3]= 0; A[0][4]= 0; A[0][5]= 0;
A[1][0]= 2; A[1][1]= 1; A[1][2]= 0; A[1][3]= 0; A[1][4]= 0; A[1][5]= 0;
A[2][0]= 4; A[2][1]=-2; A[2][2]= 0; A[2][3]= 0; A[2][4]=23; A[2][5]=24;
A[3][0]= 0; A[3][1]= 0; A[3][2]= 0; A[3][3]= 2; A[3][4]= 1; A[3][5]=-1;
A[4][0]= 0; A[4][1]= 0; A[4][2]=10; A[4][3]= 2; A[4][4]= 1; A[4][5]= 0;
A[5][0]= 0; A[5][1]= 0; A[5][2]= 0; A[5][3]= 5; A[5][4]= 2; A[5][5]=-3;
int n=6; double det=MatrixDet(n, A); char buf[100];
sprintf(buf, "%9.3f", det); LcdText(1, 0,20, buf);
while(ButtonWaitForAnyPress(100) != BUTTON_ID_LEFT);
LcdExit(); ButtonLedExit(); return 1; }
|  |  |  |  |
this matrix has a det of 74, the det of the matrix of the TOP is 0. (very exact results!)
_________________ regards, HaWe aka Ford #define S sqrt(t+2*i*i)<2 #define F(a,b) for(a=0;a<b;++a) float x,y,r,i,s,j,t,n;task main(){F(y,64){F(x,99){r=i=t=0;s=x/33-2;j=y/32-1;F(n,50&S){t=r*r-i*i;i=2*r*i+j;r=t+s;}if(S){PutPixel(x,y);}}}while(1)}
|
Tue Dec 31, 2013 11:14 am |
|
 |
mightor
Site Admin
Joined: Wed Mar 05, 2008 8:14 am Posts: 3654 Location: Rotterdam, The Netherlands
|
 Re: recursive calculation of matrix functions (eg. determina
This program seems to produce the required results:  |  |  |  | Code: // matrix determinant
#define MAX_N 6
// taken from float Determinant(int n, float *a) { int i,j,j1,j2; float det = 0; float m[MAX_N * MAX_N];
if (n < 1) { /* Error */ } else if (n == 1) { /* Shouldn't get used */ det = a[0*n+0]; } else if (n == 2) { det = a[0*n+0] * a[1*n+1] - a[1*n+0] * a[0*n+1]; } else { det = 0; for (j1=0;j1<n;j1++) { for (i=1;i<n;i++) { j2 = 0; for (j=0;j<n;j++) { if (j == j1) continue; m[(i-1)*(n-1)+j2] = a[i*n+j]; j2++; } } det += pow(-1.0,1.0+j1+1.0) * a[0*n+j1] * Determinant(n-1, m); } } return(det); }
task main() { float A[6][6];
A[0][0]=1; A[0][1]=2; A[0][2]=3; A[0][3]=4; A[0][4]=5; A[0][5]=6;
A[1][0]=10; A[1][1]=11; A[1][2]=12; A[1][3]=13; A[1][4]=14; A[1][5]=15;
A[2][0]=19; A[2][1]=20; A[2][2]=21; A[2][3]=22; A[2][4]=23; A[2][5]=24;
A[3][0]=32; A[3][1]=33; A[3][2]=34; A[3][3]=35; A[3][4]=36; A[3][5]=37;
A[4][0]=39; A[4][1]=40; A[4][2]=41; A[4][3]=42; A[4][4]=43; A[4][5]=44;
A[5][0]=51; A[5][1]=52; A[5][2]=53; A[5][3]=54; A[5][4]=55; A[5][5]=56;
float det1 = Determinant(6, A);
A[0][0]=-1; A[0][1]= 2; A[0][2]=-3; A[0][3]= 0; A[0][4]= 0; A[0][5]= 0;
A[1][0]= 2; A[1][1]= 1; A[1][2]= 0; A[1][3]= 0; A[1][4]= 0; A[1][5]= 0;
A[2][0]= 4; A[2][1]=-2; A[2][2]= 0; A[2][3]= 0; A[2][4]=23; A[2][5]=24;
A[3][0]= 0; A[3][1]= 0; A[3][2]= 0; A[3][3]= 2; A[3][4]= 1; A[3][5]=-1;
A[4][0]= 0; A[4][1]= 0; A[4][2]=10; A[4][3]= 2; A[4][4]= 1; A[4][5]= 0;
A[5][0]= 0; A[5][1]= 0; A[5][2]= 0; A[5][3]= 5; A[5][4]= 2; A[5][5]=-3;
float det2 = Determinant(6, A);
writeDebugStreamLine("%9.3f\n", det1); writeDebugStreamLine("%9.3f\n", det2); } |  |  |  |  |
_________________| Professional Conduit of Reasonableness| (Title bestowed upon on the 8th day of November, 2013) | My Blog: I'd Rather Be Building Robots| ROBOTC 3rd Party Driver Suite: [ Project Page]
|
Tue Dec 31, 2013 1:19 pm |
|
 |
Ford Prefect
Guru
Joined: Sat Mar 01, 2008 12:52 pm Posts: 1030
|
 Re: recursive calculation of matrix functions (eg. determina
yes, great, thanks, and that's actually almost exactly my code except this original ** thing (and your linearized adaption)!  happy new year everybody
_________________ regards, HaWe aka Ford #define S sqrt(t+2*i*i)<2 #define F(a,b) for(a=0;a<b;++a) float x,y,r,i,s,j,t,n;task main(){F(y,64){F(x,99){r=i=t=0;s=x/33-2;j=y/32-1;F(n,50&S){t=r*r-i*i;i=2*r*i+j;r=t+s;}if(S){PutPixel(x,y);}}}while(1)}
|
Tue Dec 31, 2013 2:28 pm |
|
|
|
Page 1 of 1
|
[ 9 posts ] |
|
Who is online |
Users browsing this forum: No registered users and 2 guests |
|
You cannot post new topics in this forum You cannot reply to topics in this forum You cannot edit your posts in this forum You cannot delete your posts in this forum You cannot post attachments in this forum
|
|