Evaluating determinants

[raymond]

I recently had to use determinants to solve a problem but didn't have a function readily available to resolve determinants. I have thus prepared one which I will eventually include in my next revision of the FPULIB. Meanwhile, if you have a need for such a function, you could use the following base procedure.

The usual precautions for preserving FPU data and making sure enough FPU registers are available to perform the function will be added for the FPULIB version. This would be YOUR responsibility with this base procedure.

The algo is based on the simplification of the given determinant by subtracting multiples of lines from other lines. A 3x3 example would be:

 a  b  c      a  b  c
 d  e  f  ->  0  j  k
 g  h  i      0  l  m

by subtracting the proper multiple of the first line from the second and third lines.

 a  b  c      a  b  c
 0  j  k  ->  0  j  k
 0  l  m      0  0  n

by subtracting the proper multiple of the new second line from the new third line.

The final result would be: a x j x n

The only restriction is that it needs the Windows API to allocate a block of memory. The final version will require the user to specify the address of a sufficient memory block so that it could also be used with some other OS.

Raymond

Code Select Expand

; #########################################################################
;
;                         FPU Determinants
;
;##########################################################################
;
; Procedure to evaluate a NxN determinant
;
; lpSrc = pointer to array of the determinant's data
; uSize = size N of the determinant
; uType = type of the data
;           0 = DWORD
;           1 = QWORD
;           2 = REAL4
;           3 = REAL8
;           4 = REAL10
;
; The result is returned in ST(0) on the FPU.
;
; The original determinant's data is not modified.
;
; -----------------------------------------------------------------------

    .386
    .model flat, stdcall  ; 32 bit memory model
    option casemap :none  ; case sensitive

      include \masm32\include\windows.inc

      include \masm32\include\kernel32.inc
      includelib \masm32\lib\kernel32.lib

; The above is simply to indicate the minimum requirements.
; Do not repeat if you paste the proc in your app.

; #########################################################################

FpuDet proc public uses ebx esi edi lpSrc:DWORD, uSize:DWORD, uType:DWORD

LOCAL hMem  :DWORD

      mov   eax,uSize
      mul   eax
      push  eax
      imul  eax,10
      invoke LocalAlloc,LPTR,eax
      mov   hMem,eax
      mov   edi,eax
      mov   esi,lpSrc
      pop   ecx
      .if   uType == 0
       @@:
          fild  dword ptr[esi]
          fstp  tbyte ptr[edi]
          add   esi,4
          add   edi,10
          dec   ecx
          jnz   @B
      .elseif uType == 1
       @@:
          fild  qword ptr[esi]
          fstp  tbyte ptr[edi]
          add   esi,8
          add   edi,10
          dec   ecx
          jnz   @B
      .elseif uType == 2
       @@:
          fld  dword ptr[esi]
          fstp  tbyte ptr[edi]
          add   esi,4
          add   edi,10
          dec   ecx
          jnz   @B
      .elseif uType == 3
       @@:
          fld  qword ptr[esi]
          fstp  tbyte ptr[edi]
          add   esi,8
          add   edi,10
          dec   ecx
          jnz   @B
      .elseif uType == 4
            imul  ecx,5
            rep   movsw
      .endif

      fld1
      mov   ebx,uSize
det0:
      mov   eax,uSize
      sub   eax,ebx
      mul   uSize
      imul  eax,10
      mov   edx,uSize
      sub   edx,ebx
      imul  edx,10
      add   eax,edx
      mov   esi,hMem
      add   esi,eax
      mov   edi,uSize
      imul  edi,10
      add   edi,esi
     
      fld   tbyte ptr[esi]
      fmul  st(1),st
      mov   ecx,ebx
      dec   ecx
      jnz   @F
      fstp  st
      invoke LocalFree,hMem
      ret

   @@:
      mov   edx,ecx
det1:
      fld   tbyte ptr[edi]
      fdiv  st,st(1)
      push  edi
      push  esi
      push  edx
det2:
      add   esi,10
      add   edi,10
      fld   tbyte ptr[esi]
      fmul  st,st(1)
      fld   tbyte ptr[edi]
      fsubr
      fstp  tbyte ptr[edi]
      dec   edx
      jnz   det2
      pop   edx
      pop   esi
      pop   edi
      mov   eax,uSize
      imul  eax,10
      add   edi,eax
      fstp  st
      dec   ecx
      jnz   det1
      fstp  st
      dec   ebx
      jmp   det0

FpuDet endp


Edit: Corrected the "dword" to "qword" for the uType == 1 and == 3

[roticv]

That's way too complicated to calculate det of a 3x3 matrix. Say Matrix A (3x3) is


a b c
d e f
g h i


Then det A = a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g)

[Eóin]

roticv you are correct, though I believe Raymond's method is the one used when the matrices get larger.

[raymond]

roticv

The 3x3 determinant was given as an example of the algo.

I do agree that a simpler code could be written to handle that particular size. However, would YOU also use a similar approach for a 10x10 or a 20x20 determinant as Eóin suspected? :dazzled:

(BTW my problem involved determinants up to the 10x10 size.)

Raymond

[joe]

Great work! Effective code in a few lines. It's good start to create inverse matrix, solve linear equations ...

[raymond]

Thanks joe. Actually, my need was for solving sets of linear equations. I was also pleasantly surprised when my little program which included a lot of other computations spew out the answer in less than 1 millisec.

Raymond

[roticv]

Hmmm I was thinking if your neeed is to solve linear equations, why not just get the augmented matrix into reduced row echelon form and you will have the solutions, just a bit more steps since you have already been gotten row echelon form when finding the det.

[joe]

I tested this procedure with 3x3 matrices only & didn't observe any problems. But didn't needed to use pivot (max. element)? Without it may lose precision.