loadflint[] := 
 With[{file1 = 
    FileNameJoin[{$UserBaseDirectory, "mylib", "flintsolver.dll"}], 
   file2 = FileNameJoin[{$UserBaseDirectory, "mylib", 
      "libwinpthread-1.dll"}]}, 
  If[! FileExistsQ[file1], 
   With[{dir = FileNameJoin[{$UserBaseDirectory, "mylib"}]}, 
    If[! DirectoryQ[dir], CreateDirectory[dir]];]; 
   URLDownload["https://files.doubilee.com/lib/win/flintsolver.dll", 
    file1];];
  If[! FileExistsQ[file2], 
   With[{dir = FileNameJoin[{$UserBaseDirectory, "mylib"}]}, 
    If[! DirectoryQ[dir], CreateDirectory[dir]];]; 
   URLDownload[
    "https://files.doubilee.com/lib/win/libwinpthread-1.dll", file2];];
  LibraryLoad[file2];
  preflintrref = 
   LibraryFunctionLoad[file1, 
    "wll_flint_rref_safe", {{Integer, 1}, {Integer, 2}, UTF8String, 
     UTF8String}, UTF8String];
  ]
unloadflint[] := 
 With[{file2 = 
    FileNameJoin[{$UserBaseDirectory, "mylib", 
      "libwinpthread-1.dll"}]}, LibraryUnload[file2]; 
  LibraryFunctionUnload /@ {preflintrref};]
loadflint[];

Flintrref[mat_]:=With[{dims=Dimensions[mat],posmat=mat["ExplicitPositions"],strs={StringTake[ToString[Numerator@#,InputForm],{2,-2}],StringTake[ToString[Denominator@#,InputForm],{2,-2}]}&@mat["ExplicitValues"]},With[{nums=First[strs],dens=Last[strs]},With[{result=preflintrref[dims,posmat,nums,dens]},SparseArray[ToExpression[result],dims]]]]
FlintSolver[eqs_,vars_]:=With[{mat=Last@CoefficientArrays[eqs,vars]},With[{dims=Dimensions[mat],posmat=mat["ExplicitPositions"],strs={StringTake[ToString[Numerator@#,InputForm],{2,-2}],StringTake[ToString[Denominator@#,InputForm],{2,-2}]}&@mat["ExplicitValues"]},With[{nums=First[strs],dens=Last[strs]},With[{result=ToExpression@preflintrref[dims,posmat,nums,dens]},Function[row,vars[[row[[1,1,2]]]]->(-(Last/@#)) . (vars[[#[[1,2]]]]&/@#)&@(Rest[row])]/@GatherBy[result,#[[1,1]]&]]]]]

(*linear algebra*)
mycoefficientarrays[exp_,vars_,max_:Infinity]:=With[{len=Length[exp]/max},If[max===Infinity||len<=1,Normal[Prepend[#[[2]] vars,#[[1]]]]&@CoefficientArrays[exp,vars],Block[{temp=0},Do[temp=temp+mycoefficientarrays[exp[[1+max i;;Min[(i+1) max,Length[exp]]]],vars,Infinity],{i,0,len}];temp]]]
mycollect[exp_,head_,max_:Infinity]:=DeleteCases[With[{tensors=Cases[1+exp,head,Infinity]//Union},mycoefficientarrays[exp,tensors,max]],0]
simplecollect[exp_,ss_]:=With[{hh=Cases[1+exp,ss,Infinity]//Union},Normal@Last@CoefficientArrays[exp,hh]]
(*mysolveceqs[ceqs_]:=With[{cs=Union@Cases[ceqs,_c,Infinity]},With[{fp=First@FirstPosition[#,1]},cs[[fp]]->-ReplacePart[#,fp->0].Append[cs,1]]&/@DeleteCases[Rationalize[#,10^-4]&@Chop[#,10^-4]&@RowReduce[EchoFunction["size of mat: ",Dimensions]@MapThread[Append,Normal@N@Reverse@CoefficientArrays[ceqs,cs]]],{0..}]]*)
mysolveceqs[ceqs_]:=With[{cs=Union@Cases[ceqs,_c,Infinity]},With[{fp=First@FirstPosition[#,1]},cs[[fp]]->-ReplacePart[#,fp->0].Append[cs,1]]&/@DeleteCases[Normal@Flintrref[EchoFunction["size of mat: ",Dimensions]@SparseArray@MapThread[Append,Normal@Reverse@CoefficientArrays[ceqs,cs]]],{0..}]]

sizeofceqs[ceqs_]:={Length[ceqs],Length[Cases[ceqs,_c,Infinity]//Union]}
mycoefficientarraysp[exp_List,vars_,max_:Infinity]:=With[{len=Length[exp]/max},If[max===Infinity||len<=1,Last@CoefficientArrays[exp,vars],Block[{temp={}},Do[temp=Join[temp,mycoefficientarraysp[exp[[1+max i;;Min[(i+1) max,Length[exp]]]],vars,Infinity]],{i,0,len}];temp]]]
(*mysolveceqsp[ceqs_,max_:Infinity]:=With[{cs=Union@Cases[ceqs,_c,Infinity]},With[{fp=First@FirstPosition[#,1]},cs[[fp]]->-ReplacePart[#,fp->0].cs]&/@With[{jj=DeleteCases[Rationalize[#,10^-4]&@Chop[#,10^-4]&@RowReduce[N[#]&@(Union@Normal@mycoefficientarraysp[ceqs,cs,max])],{0..}]},(*bigmat=jj;*)jj]]*)
mysolveceqsp[ceqs_,max_:Infinity]:=With[{cs=Union@Cases[ceqs,_c,Infinity]},With[{fp=First@FirstPosition[#,1]},cs[[fp]]->-ReplacePart[#,fp->0].cs]&/@With[{jj=DeleteCases[Normal@Flintrref[mycoefficientarraysp[ceqs,cs,max]],{0..}]},(*bigmat=jj;*)jj]]
mysolveceqspq[preceqs_,max_:Infinity]:=With[{ceqs=(preceqs/. _c:>0) c[$$last$$]+(preceqs-(preceqs/. _c:>0))},mysolveceqsp[ceqs,max]/. c[$$last$$]->1]
blockceqs[ceqs_,max_:Infinity]:=With[{cs=Union@Cases[ceqs,_c,Infinity]},With[{amat=mycoefficientarraysp[Union@ceqs,cs,max]},ceqs[[#]]&/@(((Select[#,Head[#]===$eq$&]/. $eq$->Identity)&/@ConnectedComponents[Graph[UndirectedEdge[$eq$[#[[1]]],$var$[#[[2]]]]&/@First/@Most[If[#[[2]]=!=0,#[[1]]->1,#]&/@ArrayRules[amat]]]]))]]
myblockslowsolve[ceqs_,max_:Infinity]:=With[{bceqs=EchoFunction["size of subsystems:",sizeofceqs/@#&]@blockceqs[ceqs,max]},Flatten[mysolveceqsp[#,max]&/@bceqs]]

findzeroc[tceqs_,set_]:=With[{hh=Cases[Select[tceqs,Count[1+#,_c,Infinity]===1&],_c,Infinity]//Union},If[hh==={},{tceqs,set},findzeroc[Union[tceqs//.Dispatch[#->0&/@hh]],Join[set,hh]]]]
solveshortc[tceqs_,set_]:=With[{hh=findzeroc[tceqs,{}]},With[{kk=Select[First[hh],Length[Variables[#]]<3&]},If[Length[kk]===0,{tceqs,set},With[{jj=Join[#->0&/@Last[hh],First@Solve[kk==0]]},solveshortc[DeleteCases[Union[First[hh]/. Dispatch[jj]],0],Join[set/. Dispatch[jj],jj]]]]]]
myslowsolve[ceqs_,max_:Infinity]:=If[Length[ceqs]==0,{},With[{hh=findzeroc[ceqs,{}]},Join[#->0&/@hh[[2]],myblockslowsolve[hh[[1]],max]]]]

myblockslowsolvequiet[ceqs_,max_:Infinity]:=With[{bceqs=blockceqs[ceqs,max]},Flatten[mysolveceqsp[#,max]&/@bceqs]]
myslowsolvequiet[ceqs_,max_:Infinity]:=If[Length[ceqs]==0,{},With[{hh=findzeroc[ceqs,{}]},Join[#->0&/@hh[[2]],myblockslowsolvequiet[hh[[1]],max]]]]

rrefNullSpace[mat_]:=With[{smat=SparseArray[mat]},With[{pivots=#[[1,1,2]]&/@GatherBy[Select[ArrayRules[smat],Last[#]===1&],#[[1,1]]&],dim=Dimensions[smat][[2]]},SparseArray[Function[hh,SparseArray[Select[Join[{hh->1},pivots[[#[[1,1]]]]->-Last[#]&/@Most[ArrayRules@(Transpose@smat)[[hh]]]],First[#]<=hh&],dim]]/@Complement[Range[dim],pivots]]]]
