PureBoard

Das deutsche PureBasic-Forum
https://www.purebasic.fr/german/

Ein paar Matrix-Funktionen

https://www.purebasic.fr/german/viewtopic.php?t=15311

Seite 1 von 2

Ein paar Matrix-Funktionen

Verfasst: 02.01.2008 17:58
von NicTheQuick
Falls jemand ein paar Matrix-Funktionen gebrauchen kann, die noch nicht ganz
fertig sind und wo z.B. die 'refMatrix()'-Procedure im numerischen Aspekt weit
hinten liegt, kann sich hier bedienen:

///Edit 1:
Neuste Version ab jetzt immer hier im ersten Post.

///Edit 2:
Fehler in 'fromStringMatrix_cp()' bei negativen Realteil beseitigt.
Es folgen Möglichkeiten komplexe Zahlen darzustellen. 'x' sei der reele Teil
und 'y' der imaginäre.
- x+y
- -x+y
- x-y
- -x-y
Vor dem 'x' darf höchstens ein '-' stehen, aber kein '+'. Vor dem 'y' muss
ein '-' oder '+' stehen. Beispiele, wie man es falsch macht:
- +4-3 ('4' wird als Imaginärteil erkannt, '3' wird ignoriert)
- 4i-3 ('4' wird als Realteil erkannt, '-3' als Imaginärteil, nicht andersrum)
- -4euihas+3asdfkhj ('-4' wird als Realteil, '3' als Imaginärteil erkannt)

///Edit 3:
Wieder ein Update.
Der Code ist nun wesentlich kürzer geworden, da ich ein Include gebastelt
habe, dass mit einheitlichen Aufrufen gleichzeitig mit rationalen, reellen
und komplexen Zahlen umgehen kann und alles zur Verfügung stellt, was
ich bei meinen Matrix-Funktionen brauche. Die Include "Value.pbi" ist dazu
nötig.
Es sind momentan leider nicht so viele Beispiel wie vorher dabei, aber die
Funktionsweise sollte klar sein.
Es sind auch noch ein paar Bugs drin, wenn man versucht Matrizen eines
Typs mit Skalaren oder Matrizen anderen Typs benutzt. Also bevor man
das macht, besser die Matrizen mit convertMatrix() auf den selben Typ
bringen.

Code: Alles auswählen

XIncludeFile "Value.pbi"

EnableExplicit
;works fine:
; - transMatrix
; - add, subtract, multiply, product
; - copy
; - cutted matrices
; - determinant
; - invert
; - trace
; - convert
;
;ToDo:
;
;in progress:
;

Structure Matrix
  type.l      ;
  x.l         ;x-size
  y.l         ;y-size
  *m
  *pm.Matrix   ;linked Matrix
  xpos.l       ;left corner
  ypos.l       ;upper corner
EndStructure
Structure Matrix_b
  b.b[0]
EndStructure

Procedure newMatrix(x.l, y.l, type.l = #Value_Real) ;creates a new matrix
  Protected *m.Matrix, xp.l, yp.l, *v.Value_S
  
  *m = AllocateMemory(SizeOf(Matrix))
  If Not *m : ProcedureReturn #False : EndIf
  
  *m\x = x
  *m\y = y
  *m\type = type
  *m\m = AllocateMemory(x * y * SizeOf(Value_S))
  If Not *m\m : FreeMemory(*m) : ProcedureReturn #False : EndIf
  
  *v = *m\m
  For yp = 1 To y
    For xp = 1 To x
      Value_New(type, *v)
      *v + SizeOf(Value_S)
    Next
  Next
  
  ProcedureReturn *m
EndProcedure

Procedure freeMatrix(*m.Matrix) ;destroys the matrix and frees its memory
  If *m\m : FreeMemory(*m\m) : EndIf
  FreeMemory(*m)
  ProcedureReturn #True
EndProcedure

Procedure transMatrix(*m.Matrix) ;translates the matrix without temporary memory allocation
  Protected x.l, y.l, *d1.Value, *d2.Value, n1.l, n2.l, n3.l, *b.Matrix_b, max.l
 
  If *m\x = 1 Or *m\y = 1 ;if only one row or column: O(1)
    If *m\pm And *m\x <> *m\y : ProcedureReturn #False : EndIf
    Swap *m\x, *m\y
  ElseIf *m\x = *m\y ;if it is a square matrix, mirror it: O(n(n-1)/2)
    If *m\pm
      For y = 2 To *m\y
        *d1 = *m\pm\m + (*m\xpos + y - 1 + (*m\ypos + y - 2) * *m\pm\x) * SizeOf(Value_S)
        *d2 = *m\pm\m + (*m\xpos + y - 2 + (*m\ypos + y - 1) * *m\pm\x) * SizeOf(Value_S)
        For x = y To *m\x
          *d1\Flip(*d2)
          *d1 + SizeOf(Value_S)
          *d2 + *m\pm\x * SizeOf(Value_S)
        Next
      Next
    Else
      *d1 = *m\m + SizeOf(Value_S)         ;(col 2, row 1)
      *d2 = *m\m + SizeOf(Value_S) * *m\x  ;(col 1, row 2)
      For y = 2 To *m\y
        For x = y To *m\x
          *d1\Flip(*d2)
          *d1 + SizeOf(Value_S)         ;(col+1)
          *d2 + SizeOf(Value_S) * *m\x  ;(row+1)
        Next
        *d1 + SizeOf(Value_S) * y ;(col+y)
        *d2 + SizeOf(Value_S) * (1 - *m\y * (*m\y - y))
      Next
    EndIf
  Else ;if is something else: O(max(x,y)*(min(x,y)-1))
    If *m\pm : ProcedureReturn #False : EndIf
    *d1 = *m\m
    max = *m\x * *m\y - 2
    *b = AllocateMemory(SizeOf(Byte) * max)
    If Not *b : ProcedureReturn #False : EndIf
    For n1 = 1 To max
      If *b\b[n1 - 1] = 0
        n2 = n1
        Repeat
          *b\b[n2 - 1] = 1
          n3 = n2
          n2 = n2 / *m\y + (n2 % *m\y) * *m\x
          If n2 = n1: Break : EndIf
          *d1 = *m\m + n3 * SizeOf(Value_S)
          *d2 = *m\m + n2 * SizeOf(Value_S)
          *d1\Flip(*d2)
          ;*m\ma\v[n3]\Flip(*m\ma\v[n2])
        ForEver
      EndIf
    Next
    Swap *m\x, *m\y
    FreeMemory(*b)
  EndIf
 
  ProcedureReturn *m
EndProcedure 

Procedure mulRowMatrix(*m.Matrix, y.l, *mul.Value) ;multiplies row y with mul
  Protected *d.Value, x.l
  
  If *m\pm
    *d = *m\pm\m + (*m\xpos + (*m\ypos + y - 1) * *m\pm\x) * SizeOf(Value_S)
  Else
    *d = *m\m + (y - 1) * *m\x * SizeOf(Value_S)
  EndIf
  For x = 1 To *m\x
    *d\Mul(*mul)
    *d + SizeOf(Value_S)
  Next
  
  ProcedureReturn *m
EndProcedure

Procedure mulColMatrix(*m.Matrix, x.l, *mul.Value) ;multiplies column x with mul
  Protected *d.Value, y.l, mx.l
  
  If *m\pm
    *d = *m\pm\m + (*m\xpos + x - 1 + *m\ypos * *m\pm\x) * SizeOf(Value_S)
    mx = *m\pm\x
  Else
    *d = *m\m + (x - 1) * SizeOf(Value_S)
    mx = *m\x
  EndIf
  For y = 1 To *m\y
    *d\Mul(*mul)
    *d + SizeOf(Value_S) * mx
  Next
  
  ProcedureReturn *m
EndProcedure

Procedure addmulRowMatrix(*m.Matrix, y1.l, y2.l, *mul.Value) ;adds the product of row y1 with mul to row y2
  Protected *d1.Value, *d2.Value, x.l, *d3.Value = Value_New(*m\type)
  
  If *m\pm
    *d1 = *m\pm\m + (*m\xpos + (*m\ypos + y1 - 1) * *m\pm\x) * SizeOf(Value_S)
    *d2 = *m\pm\m + (*m\xpos + (*m\ypos + y2 - 1) * *m\pm\x) * SizeOf(Value_S)
  Else
    *d1 = *m\m + (y1 - 1) * SizeOf(Value_S) * *m\x
    *d2 = *m\m + (y2 - 1) * SizeOf(Value_S) * *m\x
  EndIf
  For x = 1 To *m\x
    *d3\Set(*d1)
    *d3\Mul(*mul)
    *d2\Add(*d3)
    *d1 + SizeOf(Value_S)
    *d2 + SizeOf(Value_S)
  Next
  
  *d3\Free()
  
  ProcedureReturn *m
EndProcedure

Procedure addmulColMatrix(*m.Matrix, x1.l, x2.l, *mul.Value) ;adds the product of col x1 with mul to col x2
  Protected *d1.Value, *d2.Value, y.l, *d3.Value = Value_New(*m\type), mx.l
  
  If *m\pm
    *d1 = *m\pm\m + ((*m\xpos + x1 - 1) + *m\ypos * *m\pm\x) * SizeOf(Value_S)
    *d2 = *m\pm\m + ((*m\xpos + x2 - 1) + *m\ypos * *m\pm\x) * SizeOf(Value_S)
    mx = *m\pm\x
  Else
    *d1 = *m\m + (x1 - 1) * SizeOf(Value_S)
    *d2 = *m\m + (x2 - 1) * SizeOf(Value_S)
    mx = *m\x
  EndIf
  For y = 1 To *m\y
    *d3\Set(*d1)
    *d3\Mul(*mul)
    *d2\Add(*d3)
    *d1 + SizeOf(Value_S) * mx
    *d2 + SizeOf(Value_S) * mx
  Next
  
  *d3\Free()
  
  ProcedureReturn *m
EndProcedure

CompilerIf #PB_OS_Linux = #PB_Compiler_OS
Procedure.s toStringMatrix(*m.Matrix, nbDecimals.l = 4) ;creates a string from a matrix (eg. "1"+tab+"2"+tab+"3"+return+"4"+tab+"5"+tab+"6")
  Protected *d.Value, x.l, y.l, s.s
  
  If *m\pm
    *d = *m\pm\m + (*m\xpos + *m\ypos * *m\pm\x) * SizeOf(Value_S)
  Else 
    *d = *m\m
  EndIf
  
  For y = 1 To *m\y
    For x = 1 To *m\x
      s + *d\Str(nbDecimals)
      If x < *m\x : s + Chr(9) : EndIf
      *d + SizeOf(Value_S)
    Next
    If y < *m\y : s + Chr(10) : EndIf
    If *m\pm : *d + (*m\pm\x - *m\x) * SizeOf(Value_S) : EndIf
  Next
  
  ProcedureReturn s
EndProcedure
CompilerElse
Procedure.s toStringMatrix(*m.Matrix, nbDecimals.l = 4) ;creates a string from a matrix (eg. "1,2,3|4,5,6")
  Protected *d.Value_S, x.l, y.l, s.s
  
  If *m\pm
    *d = *m\pm\m + (*m\xpos + *m\ypos * *m\pm\x) * SizeOf(Value_S)
  Else 
    *d = *m\m
  EndIf
  
  For y = 1 To *m\y
    For x = 1 To *m\x
      s + *d\Str(nbDecimals)
      If x < *m\x : s + "," : EndIf
      *d + SizeOf(Value_S)
    Next
    If y < *m\y : s + "|" : EndIf
    If *m\pm : *d + (*m\pm\x - *m\x) * SizeOf(Value_S) : EndIf
  Next
  
  ProcedureReturn s
EndProcedure
CompilerEndIf

Procedure fromStringMatrix(String.s, type.l) ;creates a matrix from a string (eg. "1,2,3|4,5,6|7,8,9")
  Protected *c.Character, sx.l, sy.l, n.l
  Protected *m.Matrix, *d.Value, s.s
  
  sy = 1
  *c = @String
  n = 1
  While *c\c
    If *c\c = ','
      n + 1
    ElseIf *c\c = '|'
      If n > sx : sx = n : EndIf
      n = 1
      sy + 1
    EndIf
    *c + SizeOf(Character)
  Wend
  If n > sx : sx = n : EndIf
  *m = newMatrix(sx, sy, type)
  If Not *m : ProcedureReturn #False : EndIf
  *c = @String
  *d = *m\m
  n = 0
  While *c\c
    If *c\c = ','
      *d\Val(s)
      *d + SizeOf(Value_S)
      n + 1
      s = ""
    ElseIf *c\c = '|'
      *d\Val(s)
      *d + SizeOf(Value_S) * (sx - n)
      n = 0
      s = ""
    Else
      s + Chr(*c\c)
    EndIf
    *c + SizeOf(Character)
  Wend
  *d\Val(s)
  
  ProcedureReturn *m
EndProcedure

Procedure swapRowsMatrix(*m.Matrix, y1.l, y2.l) ;swaps two rows of a matrix
  Protected *d1.Value, *d2.Value, x.l
  
  If y1 = y2 : ProcedureReturn *m : EndIf
  If *m\pm
    *d1 = *m\pm\m + (*m\xpos + (*m\ypos + y1 - 1) * *m\pm\x) * SizeOf(Value_S)
    *d2 = *m\pm\m + (*m\xpos + (*m\ypos + y2 - 1) * *m\pm\x) * SizeOf(Value_S)
  Else
    *d1 = *m\m + (y1 - 1) * SizeOf(Value_S) * *m\x
    *d2 = *m\m + (y2 - 1) * SizeOf(Value_S) * *m\x
  EndIf
  For x = 1 To *m\x
    *d1\Flip(*d2)
    *d1 + SizeOf(Value_S)
    *d2 + SizeOf(Value_S)
  Next
  
  ProcedureReturn *m
EndProcedure

Procedure swapColsMatrix(*m.Matrix, x1.l, x2.l) ;swaps two columns of a matrix
  Protected *d1.Value, *d2.Value, y.l, mx.l
  
  If x1 = x2 : ProcedureReturn *m : EndIf
  If *m\pm
    *d1 = *m\pm\m + (*m\xpos + x1 - 1 + *m\ypos * *m\pm\x) * SizeOf(Value_S)
    *d2 = *m\pm\m + (*m\xpos + x2 - 1 + *m\ypos * *m\pm\x) * SizeOf(Value_S)
    mx = *m\pm\x
  Else
    *d1 = *m\m + (x1 - 1) * SizeOf(Value_S) * *m\x
    *d2 = *m\m + (x2 - 1) * SizeOf(Value_S) * *m\x
    mx = *m\x
  EndIf
  For y = 1 To *m\y
    *d1\Flip(*d2)
    *d1 + SizeOf(Value_S) * mx
    *d2 + SizeOf(Value_S) * mx
  Next
  
  ProcedureReturn *m
EndProcedure

Procedure refMatrix(*m.Matrix, *swaps.Long = 0) ;returns the reduced echelon form of a matrix
  Protected mx.l = 1, my.l = 1, *d.Value
  Protected *dtmp.Value = Value_New(*m\type), *dtmp2.Value = Value_New(*m\type)
  Protected x.l, y.l, y2.l, mxs.l, *rm, xp.l, yp.l
  
  If *m\pm
    mxs = *m\pm\x
    *rm = *m\pm\m
    xp = *m\xpos
    yp = *m\ypos
  Else
    mxs = *m\x
    *rm = *m\m
  EndIf
  Repeat
    *d = *rm + ((yp + my - 1) * mxs + xp + mx - 1) * SizeOf(Value_S)
    For y = my To *m\y
      If Not *d\Null()
        Break
      EndIf
      *d + mxs * SizeOf(Value_S)
    Next
    If y <= *m\y
      *dtmp\Set(*d)
      *d = *rm + ((yp + my - 1) * mxs + xp + mx - 1) * SizeOf(Value_S)
      For y2 = my To *m\y
        *dtmp2\Set(*d)
        *dtmp2\Div(*dtmp)
        *dtmp2\Neg()
        If y2 <> y
          addmulRowMatrix(*m, y, y2, *dtmp2)
        EndIf
        *d + SizeOf(Value_S) * mxs
      Next
      If *swaps And y <> my : *swaps\l + 1 : EndIf
      swapRowsMatrix(*m, y, my)
      mx + 1
      my + 1
    Else
      mx + 1
    EndIf
    If mx > *m\x : Break : EndIf
    ;Debug toStringMatrix(*m)
  ForEver
  
  *dtmp\Free()
  *dtmp2\Free()
  
  ProcedureReturn *m
EndProcedure

Procedure diaMatrix(*m.Matrix) ;cuts the main diagonal of the matrix
  Protected *d.Value, x.l, y.l, xm.l
  
  If *m\pm
    *d = *m\pm\m
    xm = *m\pm\x - *m\x
  Else
    *d = *m\m
    xm = 0
  EndIf
  
  For y = 1 To *m\y
    For x = 1 To *m\x
      If x <> y
        *d\SetNull()
      EndIf
      *d + SizeOf(Value_S)
    Next
    *d + xm * SizeOf(Value_S)
  Next
EndProcedure

Procedure addMatrix(*m1.Matrix, *m2.Matrix) ;adds *m2 to *m1
  Protected x.l, y.l, *d1.Value, *d2.Value, xd1.l = 0, xd2.l = 0
  
  If *m1\x <> *m2\x Or *m1\y <> *m2\y : ProcedureReturn #False : EndIf
  If *m1\pm
    *d1 = *m1\pm\m + (*m1\xpos + *m1\ypos * *m1\pm\x) * SizeOf(Value_S)
    xd1 = *m1\pm\x - *m1\x
  Else
    *d1 = *m1\m
  EndIf
  If *m2\pm
    *d2 = *m2\pm\m + (*m2\xpos + *m2\ypos * *m2\pm\x) * SizeOf(Value_S)
    xd2 = *m2\pm\x - *m2\x
  Else
    *d2 = *m2\m
  EndIf
  For y = 1 To *m1\y
    For x = 1 To *m1\x
      *d1\Add(*d2)
      *d1 + SizeOf(Value_S)
      *d2 + SizeOf(Value_S)
    Next
    *d1 + xd1 * SizeOf(Value_S)
    *d2 + xd2 * SizeOf(Value_S)
  Next
  
  ProcedureReturn *m1
EndProcedure

Procedure subMatrix(*m1.Matrix, *m2.Matrix) ;subtracts *m2 from *m1
  Protected x.l, y.l, *d1.Value, *d2.Value, xd1.l = 0, xd2.l = 0
  
  If *m1\x <> *m2\x Or *m1\y <> *m2\y : ProcedureReturn #False : EndIf
  If *m1\pm
    *d1 = *m1\pm\m + (*m1\xpos + *m1\ypos * *m1\pm\x) * SizeOf(Value_S)
    xd1 = *m1\pm\x - *m1\x
  Else
    *d1 = *m1\m
  EndIf
  If *m2\pm
    *d2 = *m2\pm\m + (*m2\xpos + *m2\ypos * *m2\pm\x) * SizeOf(Value_S)
    xd2 = *m2\pm\x - *m2\x
  Else
    *d2 = *m2\m
  EndIf
  For y = 1 To *m1\y
    For x = 1 To *m1\x
      *d1\Sub(*d2)
      *d1 + SizeOf(Value_S)
      *d2 + SizeOf(Value_S)
    Next
    *d1 + xd1 * SizeOf(Value_S)
    *d2 + xd2 * SizeOf(Value_S)
  Next
  
  ProcedureReturn *m1
EndProcedure

Procedure mulMatrix(*m.Matrix, *mul.Value) ;multiplies *m with a Value
  Protected x.l, y.l, *d.Value
  
  If *m\pm
    *d = *m\pm\m + (*m\xpos + *m\ypos * *m\pm\x) * SizeOf(Value_S)
    For y = 1 To *m\y
      For x = 1 To *m\x
        *d\Mul(*mul)
        *d + SizeOf(Value_S)
      Next
      *d + (*m\pm\x - *m\x) * SizeOf(Value_S)
    Next
  Else
    x = *m\x * *m\y
    *d = *m\m
    While x
      *d\Mul(*mul)
      *d + SizeOf(Value_S)
      x - 1
    Wend
  EndIf
  ProcedureReturn *m
EndProcedure

Procedure productMatrix(*m1.Matrix, *m2.Matrix) ;returns the matrix which is the product of m1 and m2
  Protected *mr.Matrix, x.l, y.l, *dr.Value
  Protected *d1.Value, *d2.Value, m.l, m2x.l, *dtmp.Value = Value_New(*m1\type)
  
  If *m1\x <> *m2\y : ProcedureReturn #False : EndIf
  If *m1\type <> *m2\type : ProcedureReturn #False : EndIf
  
  *mr = newMatrix(*m2\x, *m1\y, *m1\type)
  If Not *mr : ProcedureReturn #False : EndIf
  *dr = *mr\m
  
  If *m2\pm : m2x = *m2\pm\x : Else : m2x = *m2\x : EndIf
  For y = 1 To *m1\y
    For x = 1 To *m2\x
      *dr\SetNull()
      If *m1\pm
        *d1 = *m1\pm\m + (*m1\xpos + (*m1\ypos + y - 1) * *m1\pm\x) * SizeOf(Value_S)
      Else
        *d1 = *m1\m + (y - 1) * *m1\x * SizeOf(Value_S)
      EndIf
      If *m2\pm
        *d2 = *m2\pm\m + (*m2\xpos + x - 1 + *m2\ypos * *m2\pm\x) * SizeOf(Value_S)
      Else
        *d2 = *m2\m + (x - 1) * SizeOf(Value_S)
      EndIf
      For m = 1 To *m1\x
        *dtmp\Set(*d1)
        *dtmp\Mul(*d2)
        *dr\Add(*dtmp)
        *d1 + SizeOf(Value_S)
        *d2 + SizeOf(Value_S) * m2x
      Next
      *dr + SizeOf(Value_S)
    Next
  Next
  
  *dtmp\Free()
  ProcedureReturn *mr
EndProcedure

Procedure copyMatrix(*m.Matrix) ;returns a copy of *m
  Protected *cm.Matrix, *d1.Value, *d2.Value, y.l, x.l, mx.l
  
  *cm = newMatrix(*m\x, *m\y, *m\type)
  If Not *cm : ProcedureReturn #False : EndIf
  
  *d2 = *cm\m
  If *m\pm
    *d1 = *m\pm\m + (*m\xpos + *m\ypos * *m\pm\x) * SizeOf(Value_S)
    mx = *m\pm\x - *m\x
  Else
    *d1 = *m\m
    mx = 0
  EndIf
  For y = 1 To *m\y
    For x = 1 To *m\x
      *d2\Set(*d1)
      *d1 + SizeOf(Value_S)
      *d2 + SizeOf(Value_S)
    Next
    *d1 + mx * SizeOf(Value_S)
  Next
  
  ProcedureReturn *cm
EndProcedure

Procedure convertMatrix(*m.Matrix, type.l, nbDecimals.l = 10)
  Protected *d.Value, x.l, y.l
  
  If *m\pm : ProcedureReturn #False : EndIf
  
  *m\type = type
  *d = *m\m
  For y = 1 To *m\y
    For x = 1 To *m\x
      *d\Convert(type, nbDecimals)
      *d + SizeOf(Value_S)
    Next
  Next
  
  ProcedureReturn #True
EndProcedure

;cutMatrix(*m, 1, 1, 2, 2)
; /1 2 3 4\
;| 5 6 7 8 | ___\ / 6 7 \
;| 9 0 1 2 |    / \ 0 1 /
; \3 4 5 6/
Procedure cutMatrix(*m.Matrix, xpos.l, ypos.l, xsize.l, ysize.l) ;ready to use for all types of matrices
  Protected *cm.Matrix
  
  If xsize = 0 Or ysize = 0 : ProcedureReturn #False : EndIf
  If *m\x < xpos + xsize
    If xpos => *m\x : ProcedureReturn #False : EndIf
    xsize = *m\x - xpos
  EndIf
  If *m\y < ypos + ysize
    If ypos >= *m\y : ProcedureReturn #False : EndIf
    ysize = *m\y - ypos
  EndIf
  
  *cm = AllocateMemory(SizeOf(Matrix))
  If Not *cm : ProcedureReturn #False : EndIf
  
  *cm\x = xsize
  *cm\y = ysize
  *cm\type = *m\type
  If *m\pm
    *cm\pm = *m\pm
    *cm\xpos = *m\xpos + xpos
    *cm\ypos = *m\ypos + ypos
  Else
    *cm\pm = *m
    *cm\xpos = xpos
    *cm\ypos = ypos
  EndIf 
  
  ProcedureReturn *cm
EndProcedure
Procedure isCutMatrix(*m.Matrix) ;returns 0 if *m is not a cutted matrix, else pointer to the parent matrix
  ProcedureReturn *m\pm
EndProcedure

Procedure detMatrix(*m.Matrix, *det.Value = 0) ;calculates the determinant of a matrix
  Protected *mc.Matrix, *d.Value, n.l, swaps.l
  
  If *m\x <> *m\y : ProcedureReturn #False : EndIf
  
  *mc = copyMatrix(*m)
  If Not *mc : ProcedureReturn #False : EndIf
  
  If Not *det
    *det = Value_New(*mc\type)
    If Not *det : ProcedureReturn #False : EndIf
  EndIf
  
  refMatrix(*mc, @swaps)
  
  *d = *mc\m
  *det\Set(*d)
  For n = 2 To *mc\x
    *d + (1 + *mc\x) * SizeOf(Value_S)
    *det\Mul(*d)
  Next
  
  freeMatrix(*mc)
  
  If swaps & 1 : *det\Neg() : EndIf
  
  ProcedureReturn *det
EndProcedure

Procedure traceMatrix(*m.Matrix, *trace.Value = 0) ;calculates the trace of a matrix
  Protected *d.Value, n.l
  
  If *m\x <> *m\y : ProcedureReturn #False : EndIf
  
  If Not *trace
    *trace = Value_New(*m\type)
    If Not *trace : ProcedureReturn #False : EndIf
  EndIf
  
  *trace\SetNull()
  If *m\pm
    *d = *m\pm\m
    For n = 1 To *m\x
      *trace\Add(*d)
      *d + (1 + *m\pm\x) * SizeOf(Value_S)
    Next
  Else
    *d = *m\m
    For n = 1 To *m\x
      *trace\Add(*d)
      *d + (1 + *m\x) * SizeOf(Value_S)
    Next
  EndIf
  
  ProcedureReturn *trace
EndProcedure

Procedure invertMatrix(*m.Matrix)
  Protected *mi.Matrix, *d.Value, x.l, y.l, mx.l = 1, my.l = 1, *di.Value
  Protected *dtmp.Value = Value_New(*m\type), *dtmp2.Value = Value_New(*m\type)
  Protected y2.l, mxs.l, *rm, xp.l, yp.l
  
  If *m\x <> *m\y : ProcedureReturn #False : EndIf
  
  *mi = newMatrix(*m\x, *m\y, *m\type)
  If Not *mi : ProcedureReturn #False : EndIf
  
  *d = *mi\m
  For x = 1 To *mi\x
    *d\SetOne()
    *d + (*mi\x + 1) * SizeOf(Value_S)
  Next
  
  If *m\pm
    mxs = *m\pm\x
    *rm = *m\pm\m
    xp = *m\xpos
    yp = *m\ypos
  Else
    mxs = *m\x
    *rm = *m\m
  EndIf
;   Debug toStringMatrix(*m)
;   Debug toStringMatrix(*mi)
;   Debug ""
  Repeat ;Zeilenstufenform
    *d = *rm + ((yp + my - 1) * mxs + xp + mx - 1) * SizeOf(Value_S)
    For y = my To *m\y
      If Not *d\Null()
        Break
      EndIf
      *d + mxs * SizeOf(Value_S)
    Next
    If y <= *m\y
      *dtmp\Set(*d)
      *d  = *rm   + ((yp + my - 1) * mxs   + xp + mx - 1) * SizeOf(Value_S)
      *di = *mi\m + ((my - 1)      * *mi\x + mx - 1)      * SizeOf(Value_S)
      For y2 = my To *m\y
        *dtmp2\Set(*d)
        *dtmp2\Div(*dtmp)
        *dtmp2\Neg()
        If y2 <> y
          addmulRowMatrix(*m, y, y2, *dtmp2)
          addmulRowMatrix(*mi, y, y2, *dtmp2)
        EndIf
        *d + SizeOf(Value_S) * mxs
      Next
      swapRowsMatrix(*m, y, my)
      swapRowsMatrix(*mi, y, my)
      mx + 1
      my + 1
    Else
      mx + 1
    EndIf
    If mx > *m\x : Break : EndIf
;     Debug toStringMatrix(*m)
;     Debug toStringMatrix(*mi)
;     Debug ""
  ForEver
  
;   Debug "2. Schritt"
  For mx = 1 To *m\x
    *d = *rm + (mx - 1) * (1 + mxs) * SizeOf(Value_S)
    If *d\Null()
      freeMatrix(*mi)
      *dtmp\Free()
      *dtmp2\Free()
      ProcedureReturn #False
    EndIf
    *dtmp\Set(*d)
    *d = *rm + (mx - 1) * SizeOf(Value_S)
    For y = 1 To mx - 1
      *dtmp2\Set(*d)
      *dtmp2\Div(*dtmp)
      *dtmp2\Neg()
      addmulRowMatrix(*m, mx, y, *dtmp2)
      addmulRowMatrix(*mi, mx, y, *dtmp2)
      *d + mxs * SizeOf(Value_S)
    Next
    *dtmp\Invert()
    mulRowMatrix(*m, mx, *dtmp)
    mulRowMatrix(*mi, mx, *dtmp)
;     Debug toStringMatrix(*m)
;     Debug toStringMatrix(*mi)
;     Debug ""
  Next
    
  *dtmp\Free()
  *dtmp2\Free()
  
  If *m\pm
  Else
    Swap *m\m, *mi\m
  EndIf
  
  freeMatrix(*mi)
  
  ProcedureReturn *m
EndProcedure
DisableExplicit
Kleines Beispiel:

Code: Alles auswählen

Define *m1.Matrix, *m2.Matrix, *m3.Matrix, *m4.Matrix, *m5.Matrix, *m6.Matrix, *det.Value, *trace.Value, *a.Value
Define n.l = 200, *d.Double, a.l, time.l, s.s
*a = Value_New(#Value_Real)
*a\Val("2.138")

*m1 = fromStringMatrix("0,3,2,-1,2|1,-1,0,0,0|2,6,-8,4,1|-4,0,3,3,-1|0,2,-2,-4,-2", #Value_Rational)

*m2 = copyMatrix(*m1)
invertMatrix(*m2)

Debug toStringMatrixEx(*m1)
Debug toStringMatrixEx(*m2)

*m3 = productMatrix(*m1, *m2)
Debug toStringMatrixEx(*m3)

mulMatrix(*m1, *a)

Debug toStringMatrixEx(*m1)

End
Getestet unter Linux

Verfasst: 02.01.2008 18:37
von STARGÅTE
Nicht schlecht, aber viele wichtige Dinge fehlen noch, aber du hast ja geschrieben, das sie noch nicht ganz fertig sind:
- Determinate einer quadratischen Matrix
- Transponierte Matrix zu einer Matrix
- inverse Matrix zu einer Matrix
Und was noch wichtig wäre:
- die Benutzung von komplexen Zahlen für einträge in der Matrix

...

Dann könnte man daraus ein nettes kleines Anwendungsprogramm für den Bereich "Matrizen" schreiben, was bestimmt gut ankommt.

EDIT:
ich erhalte von debugger die Nachricht:

Code: Alles auswählen

---------------------------
PureBasic - Assembler error
---------------------------
PureBasic.asm [2680]:

 dd SM_l

error: undefined symbol 'SM_l'.


---------------------------
OK   
---------------------------

Verfasst: 02.01.2008 19:07
von NicTheQuick
- Determinate einer quadratischen Matrix
Die letzte Procedure, die da steht, ist dafür gedacht. :wink:
- Transponierte Matrix zu einer Matrix
Die Transponierte einer transponierten Matrix A ist die Matrix A.
- inverse Matrix zu einer Matrix
Die Inverse eine inversen Matrix A ist die Matrix A.
Und was noch wichtig wäre:
- die Benutzung von komplexen Zahlen für einträge in der Matrix
Das wäre wohl das geringste Problem, ist aber eine gute Idee, die ich
noch einbauen könnte.
EDIT:
ich erhalte von debugger die Nachricht:

Code: Alles auswählen

---------------------------
PureBasic - Assembler error
---------------------------
PureBasic.asm [2680]:

 dd SM_l

error: undefined symbol 'SM_l'.


---------------------------
OK   
---------------------------
Tja, hmm... Da weiß ich jetzt nicht, woran es liegt.
Hab's grad nochmal unter PB 4.10 (Windows) getestet und da funktioniert es.

Unter Linux werden übrigens die Matrizen im Debugger richtig dargestellt,
also nicht alles nebeneinander, sondern Zeile für Zeile.

Verfasst: 02.01.2008 19:22
von STARGÅTE
gut dann liegt es daran das ich 4.20 habe -.-
naja ^^ ist halt doch nicht immer gut auf dem neusten Stand zu sein :D

Verfasst: 15.01.2008 15:16
von NicTheQuick
Hab mal meine Matrix-Funktionen etwas erweitert.
Die Inverse und die Determinante sind zwar immer noch nicht drin, aber
dafür kann jetzt mit komplexen Zahlen gerechnet werden.
Außerdem ist es möglich aus einer Matrix eine kleinere Matrix
'auszuschneiden', die allerdings nur eine Referenz zur großen Matrix ist.
Das heißt, wenn man etwas an der kleinen Matrix ändert, ändert sich der
Teil in der großen automatisch mit.

Ein paar Beispiele sind enthalten.

Hinweise:
- 'transMatrix()' funktioniert bei referenzierten Matrizen nur, wenn diese quadratisch sind
- Vorsicht mit referenzierten Matrizen und der Transponierung ihrer Referenz: Im Allgemeinen wird die ref. Matrix keinen Sinn mehr ergeben, aber auch keinen "Invalid Memory Access" liefern.
- Vorsicht bei referenzierten Matrizen und der Anwendung von 'freeMatrix()' auf ihre Referenz: Die referenzierten Matrizen werden dann ins leere führen.

Neuer Code: Siehe ersten Post.

Verfasst: 15.01.2008 15:45
von STARGÅTE
man sieht das viel aus ^^

Wenn du willst kann ich dir mal meine Matrix Determinatenberechnung schicken.

Verfasst: 15.01.2008 16:08
von NicTheQuick
Ja, das kannst du mir gerne mal schicken.
Aber wenn du Laplace oder das Gaußverfahren nutzt, brauch ich es nicht.
Diese beiden Verfahren könnte ich selbst schnell implementieren, habe ich
aber bisher noch nicht getan, weil sie für große Matrizen zu langsam sind.

Verfasst: 15.01.2008 16:31
von STARGÅTE
jo, ich habe es halt Rekurvis gemacht, also nach Laplace.

ich meine Rekursiv musst du es doch so oder so machen, egal ob du es direkt an der Matrix machst, oder ob du dir n Formel rekursiv zusammen stellst und dann die mega formel ausrechnest.

Obwohl mir gerade eine Idee kommt das es vllt noch schneller geht.....


Hier erst mal meine Det-Variante
Code in PureBasic hat geschrieben:; Berechnet die Determinate einer quadratischen Matrix
Procedure.f MatrixDeterminante(*Matrix)
 xx = PeekL(*Matrix+0) : yy = PeekL(*Matrix+4)
 Zeilen = xx+1 : Spalten = yy+1
 If Zeilen = Spalten
  If Zeilen = 1
   ProcedureReturn PeekF(*Matrix+#MatrixKopf)
  Else 
   Determinante.f = 0
   For y = 0 To yy
    TempEintrag.f = PeekF(*Matrix+y*Zeilen*#MatrixEintrag+0*#MatrixEintrag+#MatrixKopf)
    *TempMatrix = CreateMatrix(Zeilen-1, Spalten-1)
    TempZeilen = Zeilen-1 
    ty = 0
    For t = 0 To yy
     If t <> y
      Laenge = 4*xx
      CopyMemory(*Matrix+t*Zeilen*#MatrixEintrag+1*#MatrixEintrag+#MatrixKopf, *TempMatrix+ty*TempZeilen*#MatrixEintrag+#MatrixKopf, Laenge)
      ty + 1
     EndIf
    Next t
    TempDeterminante.f = MatrixDeterminante(*TempMatrix)
    DeleteMatrix(*TempMatrix)
    Determinante = Determinante + Pow(-1, 1+(y+1))*TempEintrag*TempDeterminante
   Next y
   ProcedureReturn Determinante
  EndIf
 Else
  ProcedureReturn #False 
 EndIf
EndProcedure

Verfasst: 15.01.2008 19:05
von NicTheQuick
Hab mal die Determinante nach Gauß eingebaut.

Der Code befindet sich im ersten Post.

Verfasst: 15.01.2008 20:32
von STARGÅTE
ist vllt n dumme Frage, aber wieso hast du jetzt immer alle Proceduren "Doppelt" für Reelle und Komplexe Zahlen.

Du kannst das doch in einer vereinigen, denn die Komplexe Zahlen benhalten ja die Reellen, nur halt ohne den Imaginärteil.
Dann könntest du vllt Macros draus machen, sodass der Code je nach Zahlenart angepasst wird, das erspart dir n menge Code.

[OffTopic]
Ich hasse es z.B. wenn ich mehrere proceduren haben die fast das gleiche machen, denn dann muss man bei änderungen immer beide anpassen, und igrnedwann sind sie dann doch nicht mehr gleich.
Dann mache ich lieber n If-Abfrage mehr rein und schreibe die in eine gemeinsamme Procedure.
Aber ist ja dein Code *sich also raushalt*
[/OffTopic]
Alle Zeiten sind UTC+02:00
Seite 1 von 2

Powered by phpBB® Forum Software © phpBB Limited

Deutsche Übersetzung durch phpBB.de