Sparse Matrix (lab) 深刻なビジネス です

{$r-,q-}
// Das Feen Reich!
// Fus Roh Dah!~

unit SparseMatrix;

interface

uses
  Math, SysUtils;

type
  sint32 = longint;
  uint32 = longword;
  sint = sint32;
  uint = uint32;       pUint = ^uint;
  float = extended;
  PtrUint = uint32; // unsigned integer with same size as pointer

const
  FloatEps = 1e-17;

type
  tSMHashValue = type uint32;

  pSMHashRec = ^tSMHashRec;
  tSMHashRec = object
    h: tSMHashValue;
    x, y: uint;
    value: float;
    next: pSMHashRec;
    procedure UpdateHash;
  end;

  tSparseMatrixElemProc = procedure(var x, y: uint; var value: float; param: pointer);

  CSparseMatrix = class
  private
    _nh: uint;
    _h: array of pSMHashRec;
    _sizeX, _sizeY: uint;
    _relocating: boolean;
    procedure _AddHashRec(r: pSMHashRec);
    function _FindHashRec(x, y: uint): pSMHashRec;
    procedure _RemHashRec(x, y: uint);
    procedure _RelocateHash(force: boolean = false);
    function _GetCell(x, y: uint): float;
    procedure _SetCell(x, y: uint; const value: float);
    procedure _SetSizeX(newSizeX: uint);
    procedure _SetSizeY(newSizeY: uint);
  public
    constructor Create(newSizeX, newSizeY: uint); overload;
    constructor Create(m: CSparseMatrix); overload;
    destructor Destroy; override;
    procedure ForEach(proc: tSparseMatrixElemProc; param: pointer);
    procedure RemoveCol(x: uint);
    procedure RemoveRow(y: uint);
    procedure SwapRows(y1, y2: uint);
    procedure AddRow(target, row: uint; factor: float);
    function GetMinor(x, y: uint): float;
    function GetAlgebraicComplement(x, y: uint): float;
    function Det: float;
    procedure DoForwardGauss;

    property SizeX: uint read _sizeX write _SetSizeX;
    property SizeY: uint read _sizeY write _SetSizeY;
    property Cells[x, y: uint]: float read _GetCell write _SetCell; default;
  end;

  function SparseMul(a, b: CSparseMatrix): CSparseMatrix;


implementation


  {function HashUint32(x: uint32): tSMHashValue;
  begin
    x := (x + $7ed55d16) + (x shl 12);
    x := (x xor $c761c23c) xor (x shr 19);
    x := (x + $165667b1) + (x shl 5);
    x := (x + $d3a2646c) xor (x shl 9);
    x := (x + $fd7046c5) + (x shl 3);
    x := (x xor $b55a4f09) xor (x shr 16);
    result := x;
  end;}

  function HashUint32(x: uint32): tSMHashValue;
  begin
    x := x + not (x shl 15);
    x := x xor (x shr 10);
    x := x + (x shl 3);
    x := x xor (x shr 6);
    x := x + not (x shl 11);
    x := x xor (x shr 16);
    result := x;
  end;

  function HashUint32XY(x, y: uint32): tSMHashValue;
  begin
    result := HashUint32(x) xor HashUint32(y);
  end;

  function IsZero(x: float): boolean;
  begin
    result := abs(x) < FloatEps;
  end;

  function NotZero(x: float): boolean;
  begin
    result := abs(x) >= FloatEps;
  end;

  function IsNaN(x: float): boolean;
  begin
    result := x <> x;
  end;

  function IsInf(x: float): boolean;
  begin
    result := x > x;
  end;

  function IsNegInf(x: float): boolean;
  begin
    result := x < x;
  end;


  procedure tSMHashRec.UpdateHash;
  begin
    h := HashUint32XY(x, y);
  end;

  procedure CSparseMatrix._AddHashRec(r: pSMHashRec);
  var
    i: uint;
  begin
    if not Assigned(r) then exit;
    _RemHashRec(r^.x, r^.y);
    i := r^.h mod uint(length(_h));
    r^.next := _h[i];
    _h[i] := r;
    inc(_nh);
    _RelocateHash;
  end;

  function CSparseMatrix._FindHashRec(x, y: uint): pSMHashRec;
  var
    i: uint;
  begin
    i := HashUint32XY(x, y) mod uint(length(_h));
    result := _h[i];
    while (Assigned(result)) and ((result^.x <> x) or (result^.y <> y)) do
      result := result^.next;
  end;

  procedure CSparseMatrix._RemHashRec(x, y: uint);
  var
    r, cur, prev: pSMHashRec;
    i: uint;
  begin
    r := _FindHashRec(x, y);
    if not Assigned(r) then exit;
    i := r^.h mod uint(length(_h));
    prev := nil;
    cur := _h[i];
    while Assigned(cur) do
    begin
      if cur = r then
      begin
        if Assigned(prev) then
          prev^.next := cur^.next
        else
          _h[i] := cur^.next;
        break;
      end;
      prev := cur;
      cur := cur^.next;
    end;
    dispose(r);
    dec(_nh);
    _RelocateHash;
  end;

  procedure CSparseMatrix._RelocateHash(force: boolean = false);
  var
    store, c, t: pSMHashRec;
    i: uint;
  begin
    if (_nh >= uint(length(_h) div 2)) and (_nh <= uint(2*length(_h))) and
      (not force) then
      exit;
    if _relocating then exit;
    _relocating := true;
    store := nil;
    for i := 0 to High(_h) do
    begin
      c := _h[i];
      while Assigned(c) do
      begin
        t := c^.next;
        c^.next := store;
        store := c;
        c := t;
      end;
    end;
    if _nh > 0 then SetLength(_h, _nh) else SetLength(_h, 1);
    for i := 0 to High(_h) do _h[i] := nil;
    _nh := 0;
    while Assigned(store) do
    begin
      t := store;
      store := store^.next;
      _AddHashRec(t);
    end;
    _relocating := false;
  end;

  function CSparseMatrix._GetCell(x, y: uint): float;
  var
    c: pSMHashRec;
  begin
    if (x >= _sizeX) or (y >= _sizeY) then
    begin
      result := 0.0;
      exit;
    end;
    c := _FindHashRec(x, y);
    if Assigned(c) then
      result := c^.value
    else
      result := 0.0;
  end;

  procedure CSparseMatrix._SetCell(x, y: uint; const value: float);
  var
    c: pSMHashRec;
  begin
    if (x >= _sizeX) or (y >= _sizeY) then exit;
    if IsZero(value) then
      _RemHashRec(x, y)
    else
    begin
      new(c);
      c^.x := x;
      c^.y := y;
      c^.value := value;
      c^.UpdateHash;
      _AddHashRec(c);
    end;
  end;

  procedure __elemfits(var x, y: uint; var value: float; param: pointer);
  var
    m: CSparseMatrix absolute param;
  begin
    if (x >= m.SizeX) or (y >= m.SizeY) then value := 0.0;
  end;

  procedure CSparseMatrix._SetSizeX(newSizeX: uint);
  begin
    if _sizeX <> newSizeX then
    begin
      _sizeX := newSizeX;
      ForEach(__elemfits, self);
    end;
  end;

  procedure CSparseMatrix._SetSizeY(newSizeY: uint);
  begin
    if _sizeY <> newSizeY then
    begin
      _sizeY := newSizeY;
      ForEach(__elemfits, self);
    end;
  end;

  constructor CSparseMatrix.Create(newSizeX, newSizeY: uint);
  begin
    inherited Create;
    _relocating := false;
    _nh := 0;
    SetLength(_h, 1);
    _h[0] := nil;
    _sizeX := 0;
    _sizeY := 0;
    SizeX := newSizeX;
    SizeY := newSizeY;
  end;

  procedure __AddElem(var x, y: uint; var value: float; param: pointer);
  begin
    CSparseMatrix(param)[x, y] := value;
  end;

  constructor CSparseMatrix.Create(m: CSparseMatrix);
  begin
    Create(m.SizeX, m.SizeY);
    m.ForEach(__AddElem, self);
  end;

  destructor CSparseMatrix.Destroy;
  var
    c, t: pSMHashRec;
    i: sint;
  begin
    for i := 0 to High(_h) do
    begin
      c := _h[i];
      while Assigned(c) do
      begin
        t := c;
        c := c^.next;
        dispose(t);
      end;
    end;
    inherited Destroy;
  end;

  procedure CSparseMatrix.ForEach(proc: tSparseMatrixElemProc; param: pointer);
  var
    i: sint;
    prev, c, t: pSMHashRec;
  begin
    for i := 0 to High(_h) do
    begin
      prev := nil;
      c := _h[i];
      while Assigned(c) do
      begin
        t := c;
        c := c^.next;
        proc(t^.x, t^.y, t^.value, param);
        if NotZero(t^.value) then
        begin
          t^.UpdateHash;
          prev := t;
        end else
        begin
          if Assigned(prev) then
            prev^.next := c
          else
            _h[i] := c;
          dispose(t);
          dec(_nh);
        end;
      end;
    end;
    _RelocateHash(true);
  end;

  procedure __RemCol(var x, y: uint; var value: float; param: pointer);
  begin
    if x = PtrUint(param) then value := 0.0 else
    if x > PtrUint(param) then dec(x);
  end;

  procedure __RemRow(var x, y: uint; var value: float; param: pointer);
  begin
    if y = PtrUint(param) then value := 0.0 else
    if y > PtrUint(param) then dec(y);
  end;

  procedure CSparseMatrix.RemoveCol(x: uint);
  begin
    if x >= _sizeX then exit;
    ForEach(__RemCol, pointer(PtrUint(x)));
    SizeX := SizeX - 1;
  end;

  procedure CSparseMatrix.RemoveRow(y: uint);
  begin
    if y >= _sizeX then exit;
    ForEach(__RemRow, pointer(PtrUint(y)));
    SizeY := SizeY - 1;
  end;

  procedure CSparseMatrix.SwapRows(y1, y2: uint);
  var
    t: float;
    i: sint;
  begin
    if (y1 >= _sizeY) or (y2 >= _sizeY) or (y1 = y2) then exit;
    for i := 0 to _sizeX-1 do
    begin
      t := self[i, y1];
      self[i, y1] := self[i, y2];
      self[i, y2] := t;
    end;
  end;

  procedure CSparseMatrix.AddRow(target, row: uint; factor: float);
  var
    i: sint;
  begin
    if (target >= _sizeY) or (row >= _sizeY) then exit;
    for i := 0 to _sizeX-1 do
      self[i, target] := self[i, target] + self[i, row] * factor;
  end;

  function CSparseMatrix.GetMinor(x, y: uint): float;
  var
    tm: CSparseMatrix;
  begin
    if (x >= _sizeX) or (y >= _sizeY) or (_sizeX <> _sizeY) or (_sizeX <= 1) then
    begin
      result := NaN;
      exit;
    end;
    // terrible! :D
    tm := CSparseMatrix.Create(self);
    tm.RemoveCol(x);
    tm.RemoveRow(y);
    result := tm.Det;
    tm.Free;
  end;

  function CSparseMatrix.GetAlgebraicComplement(x, y: uint): float;
  begin
    result := GetMinor(x, y);
    if (not IsNaN(result)) and odd(x + y) then result := -result;
  end;

  function CSparseMatrix.Det: float;
  var
    i: sint;
    v: float;
  begin
    if (_sizeX <> _sizeY) or (_sizeX = 0) then
    begin
      result := NaN;
      exit;
    end;
    if _sizeX = 1 then
    begin
      result := self[0, 0];
      exit;
    end;
    result := 0.0;
    for i := 0 to _sizeX-1 do
    begin
      v := self[i, 0];
      if NotZero(v) then
        result := result + v * GetAlgebraicComplement(i, 0);
    end;
  end;

  procedure CSparseMatrix.DoForwardGauss;
  var
    i, j, sy, firstnonzero: sint;
    x: float;
  begin
    sy := 0;
    for i := 0 to _sizeX-1 do
    begin
      firstnonzero := -1;
      for j := sy to _sizeY-1 do
        if NotZero(self[i, j]) then
        begin
          firstnonzero := j;
          break;
        end;
      if firstnonzero = -1 then continue;
      SwapRows(sy, firstnonzero);
      x := self[i, sy];
      for j := sy+1 to _sizeY-1 do
        AddRow(j, sy, -self[i, j] / x);
      inc(sy);
    end;
  end;

  function SparseMul(a, b: CSparseMatrix): CSparseMatrix;
  var
    x, y, k: sint;
    v: float;
  begin
    if a.SizeX <> b.SizeY then
    begin
      result := nil;
      exit;
    end;
    result := CSparseMatrix.Create(b.SizeX, a.SizeY);
    for x := 0 to result.SizeX-1 do
      for y := 0 to result.SizeY-1 do
      begin
        v := 0.0;
        for k := 0 to a.SizeX-1 do
          v := v + a[k, y] * b[x, k];
        result[x, y] := v;
      end;
  end;

end.