Skip to content

Instantly share code, notes, and snippets.

@MrSmith33
Created January 31, 2014 13:59
Show Gist options
  • Save MrSmith33/8732520 to your computer and use it in GitHub Desktop.
Save MrSmith33/8732520 to your computer and use it in GitHub Desktop.
determinant.d
import rational;
import std.conv : to;
import std.algorithm;
import std.stdio;
import std.traits;
import matrix;
import vector;
ElementType pow(ElementType, I)(ElementType number, I power) if (isIntegral!I)
{
ElementType result = number;
foreach(_; 0..power-1) number *= number;
return number;
}
struct CoeffMatrix(ElementType, CoeffType)
{
Matrix!ElementType matrix;
CoeffType coefficient = to!CoeffType(1);
string toString()
{
static if (isRational!CoeffType)
return (coefficient.denominator==1 ? to!string(coefficient.numerator) : to!string(coefficient.toString)) ~"\n" ~ matrix.toString();
}
}
alias TempMatrix(ElementType) = CoeffMatrix!(Matrix!ElementType, ElementType);
ElementType solve2x2(ElementType)(Matrix!ElementType matrix)
{
assert(matrix.size == 2);
return matrix[0, 0] * matrix[1, 1] - matrix[0, 1] * matrix[1, 0];
}
CoeffMatrix!(Matrix!ElementType, ElementType) makeTemp(T : CoeffMatrix!(ElementType, ElementType), ElementType)(T matrix)
{
auto result = CoeffMatrix!(Matrix!ElementType, ElementType)(Matrix!(Matrix!ElementType)(matrix.matrix.size-1), matrix.coefficient);
foreach(x; 0..result.matrix.size)
{
foreach(y; 0..result.matrix.size)
{
result.matrix[x, y] = Matrix!ElementType(2, [matrix.matrix[0, 0], matrix.matrix[x+1, 0],
matrix.matrix[0, y+1], matrix.matrix[x+1, y+1]]);
}
}
result.coefficient = matrix.coefficient * 1/(pow(matrix.matrix[0, 0], matrix.matrix.size-2));
return result;
}
CoeffMatrix!(ElementType, ElementType) solveTemp(T : CoeffMatrix!(Matrix!ElementType, ElementType), ElementType)(T cmatrix)
{
writeln("solveTemp start");
auto result = CoeffMatrix!(ElementType, ElementType)( Matrix!ElementType(cmatrix.matrix.size), cmatrix.coefficient );
foreach(ind; 0..(result.matrix.size^^2))
{
result.matrix.elements[ind] = solve2x2!ElementType(cmatrix.matrix.elements[ind]);
}
writeln("before solveTemp return");
return result;
}
void solveAndPrint(T : CoeffMatrix!(ElementType, ElementType), ElementType)(T _matrix)
{
auto matrix = _matrix;
TempMatrix!ElementType temp;
auto num = 3;
foreach(_; 0..num)
{
writeln("foreach start");
temp = makeTemp(matrix);
//writeln(temp);
writeln("after makeTemp");
matrix = solveTemp(temp);
//writeln(matrix);
writeln("foreach end");
stdout.flush();
}
writeln("done foreach");
stdout.flush();
}
void main()
{
alias CoeffMatrixType = CoeffMatrix!(Rational!int, Rational!int);
alias TempMatrixType = CoeffMatrix!(CoeffMatrix!(Rational!int, Rational!int), Rational!int);
CoeffMatrixType input = CoeffMatrixType(Matrix!(Rational!int)(3, [2, -1, 1, 1, 2, 3, 1, -3, -2]));
// 5, 5, -5, -5. Det = 0
solveAndPrint(input);
writeln("solved");
}
import std.conv : to;
import std.exception : enforce;
import rational;
struct Matrix(ElementType)
{
ElementType[] elements;
uint size;
//@disable this();
this(uint size)
{
this.size = size;
elements.length = size*size;
}
this(Element)(uint size, Element[] initialValues)
{
enforce(initialValues.length == size*size);
this.size = size;
static if(is(Element == ElementType))
elements = initialValues;
else
foreach(ind, val; initialValues)
elements ~= to!ElementType(val);
}
ElementType opIndex(size_t x, size_t y)
{
return elements[x + y*size];
}
ElementType opIndexAssign(ElementType value, size_t x, size_t y)
{
return elements[x + y*size] = value;
}
ElementType opIndexOpAssign(string op)(ElementType value, size_t x, size_t y)
{
return mixin("elements[x + y*size] "~op~"= value");
}
void swapRows(uint a, uint b)
{
auto temp = elements[a*size..a*size+size];
elements[a*size..a*size+size] = elements[b*size..b*size+1];
elements[b*size..b*size+1] = temp;
}
void cycleRows()
{
elements = elements[size..$]~elements[0..size];
}
string toString()
{
string result;
static if(isRational!(ElementType))
{
foreach(ind, elem; elements)
{
result ~= (elem.denominator==1 ? to!string(elem.numerator) : elem.toString()) ~ ' ';
if ((ind + 1) % size == 0)
result ~= "\n";
}
}
else
{
foreach(ind, elem; elements)
{
result ~= to!string(elem) ~ ' ';
if ((ind + 1) % size == 0)
result ~= "\n";
}
}
return result;
}
}
import std.conv : to;
struct Vector(ElementType)
{
ElementType[] elements;
uint size;
@disable this();
this(uint size)
{
this.size = size;
elements.length = size;
}
this(Element)(uint size, Element[] initialValues)
{
this.size = size;
static if(is(Element == ElementType))
elements = initialValues;
else
foreach(ind, val; initialValues)
elements ~= ElementType(val);
}
ElementType opIndex(size_t index)
{
return elements[index];
}
ElementType opIndexAssign(ElementType value, size_t index)
{
return elements[index] = value;
}
ElementType opIndexOpAssign(string op)(ElementType value, size_t index)
{
return mixin("elements[index] "~op~"= value");
}
void swapElements(uint a, uint b)
{
auto temp = elements[a];
elements[a] = elements[b];
elements[b] = temp;
}
void cycleElements()
{
elements = elements[1..$] ~ elements[0];
}
string toString()
{
string result;
foreach(ind, elem; elements)
{
result ~= (elem.denominator==1 ? to!string(elem.numerator) : elem.toString()) ~ ' ';
}
return result;
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment