View unanswered posts | View active topics It is currently Wed Jul 30, 2014 11:08 am






Reply to topic  [ 9 posts ] 
recursive calculation of matrix functions (eg. determinants) 
Author Message
Guru
User avatar

Joined: Sat Mar 01, 2008 12:52 pm
Posts: 1030
Post 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
Profile
Moderator
Moderator
User avatar

Joined: Wed Mar 05, 2008 8:14 am
Posts: 3165
Location: Rotterdam, The Netherlands
Post 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
Profile WWW
Guru
User avatar

Joined: Sat Mar 01, 2008 12:52 pm
Posts: 1030
Post 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
Profile
Moderator
Moderator
User avatar

Joined: Wed Mar 05, 2008 8:14 am
Posts: 3165
Location: Rotterdam, The Netherlands
Post 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
Profile WWW
Moderator
Moderator
User avatar

Joined: Wed Mar 05, 2008 8:14 am
Posts: 3165
Location: Rotterdam, The Netherlands
Post 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
Profile WWW
Guru
User avatar

Joined: Sat Mar 01, 2008 12:52 pm
Posts: 1030
Post 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:

Code:
  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;

_________________
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
Profile
Guru
User avatar

Joined: Sat Mar 01, 2008 12:52 pm
Posts: 1030
Post 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
Profile
Moderator
Moderator
User avatar

Joined: Wed Mar 05, 2008 8:14 am
Posts: 3165
Location: Rotterdam, The Netherlands
Post 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
Profile WWW
Guru
User avatar

Joined: Sat Mar 01, 2008 12:52 pm
Posts: 1030
Post 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
Profile
Display posts from previous:  Sort by  
Reply to topic   [ 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

Search for:
Jump to:  



Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group.
Designed by ST Software for PTF.