{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a43dfc3e",
   "metadata": {},
   "source": [
    "# Matrix generators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ae9a526f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sgen(i,c,t):\n",
    "    r\"\"\"\n",
    "    returns the generator s_i(c)\n",
    "    \"\"\"\n",
    "    M=MatrixSpace(QQ,t,t)\n",
    "    I=identity_matrix(QQ, t)\n",
    "    E=M.basis()\n",
    "    return( I-E[i-1,i-1]-E[i,i]+E[i-1,i]+E[i,i-1]+c*E[i-1,i-1] )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b14bb6ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sinv(i,c,t):\n",
    "    r\"\"\"\n",
    "    returns  s_i(c)^{-1}\n",
    "    \"\"\"\n",
    "    M=MatrixSpace(QQ,t,t)\n",
    "    E=M.basis()\n",
    "    I=identity_matrix(QQ, t)\n",
    "    return( I-E[i-1,i-1]-E[i,i]+E[i-1,i]+E[i,i-1]-c*E[i,i] )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3e9b32a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def hgen(i,c,t):\n",
    "    r\"\"\"\n",
    "    returns the generator h_i(c)\n",
    "    \"\"\"\n",
    "    M=MatrixSpace(QQ,t,t)\n",
    "    E=M.basis()\n",
    "    I=identity_matrix(QQ, t)\n",
    "    return( I+(c-1)*E[i-1,i-1] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e5daf997",
   "metadata": {},
   "outputs": [],
   "source": [
    "def xgen(i,j,c,t):\n",
    "    r\"\"\"\n",
    "    returns the generator x_{ij}i(c)\n",
    "    \"\"\"\n",
    "    M=MatrixSpace(QQ,t,t)\n",
    "    E=M.basis()\n",
    "    I=identity_matrix(QQ, t)\n",
    "    if i < j:\n",
    "        return( I + c*E[i-1,j-1] )\n",
    "    else:\n",
    "        print(\"must have i not equal to j in x(i,j,c) so returning 1\")\n",
    "        return(I)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e0241ff0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1 0 0 0]\n",
       "[0 1 0 0]\n",
       "[0 0 2 1]\n",
       "[0 0 1 0]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sgen(3,2,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a5c9f4e0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[ 1  0  0  0]\n",
       "[ 0  0  1  0]\n",
       "[ 0  1 -5  0]\n",
       "[ 0  0  0  1]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sinv(2,5,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "108d50db",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1 0 0 0]\n",
       "[0 5 0 0]\n",
       "[0 0 1 0]\n",
       "[0 0 0 1]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hgen(2,5,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b739186a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1 0 0 0]\n",
       "[0 1 7 0]\n",
       "[0 0 1 0]\n",
       "[0 0 0 1]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xgen(2,3,7,4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6374f63d-fa50-4560-b862-33cde2db6ae7",
   "metadata": {},
   "source": [
    "# Some matrices for examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "eed00307-8dd7-4066-8e8a-5a06be3412a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "M46=MatrixSpace(QQ,4,6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "7aa94e25-3fa4-4be4-b8e5-11543bdcf835",
   "metadata": {},
   "outputs": [],
   "source": [
    "M55=MatrixSpace(QQ,5,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "bdd1034b",
   "metadata": {},
   "outputs": [],
   "source": [
    "A=M46([-3,-36,-39,-42,-46,-84, -9,-108,-233/2,-114,-247/2,-167, 3,36,39,42,45,48, 0,0,1/4,6,25/4,13/2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "399bf122",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[    -3    -36    -39    -42    -46    -84]\n",
       "[    -9   -108 -233/2   -114 -247/2   -167]\n",
       "[     3     36     39     42     45     48]\n",
       "[     0      0    1/4      6   25/4   13/2]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "46f6c313-7e44-45fd-a685-4d8c81db4dc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "B0=M46([3,36,39,42,45,48, 0,0,1/4,6,25/4,13/2, 0,0,0,0,-1,-36, 0,0,0,0,0,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "19366c72-35b0-4eb6-bfd2-cec286fbadd6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[   3   36   39   42   45   48]\n",
       "[   0    0  1/4    6 25/4 13/2]\n",
       "[   0    0    0    0   -1  -36]\n",
       "[   0    0    0    0    0    0]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "d2561f0d-8094-44e9-88a1-36847d551360",
   "metadata": {},
   "outputs": [],
   "source": [
    "C0=M46([1,12,13,14,15,16, 0,0,1,24,25,26, 0,0,0,0,1,36, 0,0,0,0,0,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "bfb3145f-2e4b-4c39-a397-4f0e105f7bb2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[ 1 12 13 14 15 16]\n",
       "[ 0  0  1 24 25 26]\n",
       "[ 0  0  0  0  1 36]\n",
       "[ 0  0  0  0  0  0]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "c2ef8b1f-fafe-4c6d-8503-bf66d8b14257",
   "metadata": {},
   "outputs": [],
   "source": [
    "V=M55([4,7,9,10,1, 3,6,8,1,0, 2,5,1,0,0, 1,1,0,0,0, 1,0,0,0,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "59e4e67e-79db-4f38-98b1-df72ad262795",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[ 4  7  9 10  1]\n",
       "[ 3  6  8  1  0]\n",
       "[ 2  5  1  0  0]\n",
       "[ 1  1  0  0  0]\n",
       "[ 1  0  0  0  0]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "d5150a08-806d-4e43-b2a8-8259625aa334",
   "metadata": {},
   "outputs": [],
   "source": [
    "w0=M55([0,0,0,0,1, 0,0,0,1,0, 0,0,1,0,0, 0,1,0,0,0, 1,0,0,0,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "9ba12aee-5737-4d3b-b71f-d8f3d29763c1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0 0 0 0 1]\n",
       "[0 0 0 1 0]\n",
       "[0 0 1 0 0]\n",
       "[0 1 0 0 0]\n",
       "[1 0 0 0 0]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "8244699e-3c3a-4592-b79a-24e0b3145875",
   "metadata": {},
   "outputs": [],
   "source": [
    "U=M55([1,10,9,7,4, 0,1,8,6,3, 0,0,1,5,2, 0,0,0,1,1, 0,0,0,0,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "14933581-ccc9-4fef-8f42-f69dad8c2c05",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[ 1 10  9  7  4]\n",
       "[ 0  1  8  6  3]\n",
       "[ 0  0  1  5  2]\n",
       "[ 0  0  0  1  1]\n",
       "[ 0  0  0  0  1]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9724cfde-b9a7-486e-b787-9e31e3a7d0e6",
   "metadata": {},
   "source": [
    "# First nonzero sequence, rank and status checks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "895c5ce3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fnz(A):\n",
    "    r\"\"\"\n",
    "    Constructs the first non-zero sequence of A, i.e.\n",
    "    the sequence recording the position of the first nonzero\n",
    "    entry in each row of A\n",
    "    \"\"\"\n",
    "    f = [-1]\n",
    "    t = A.nrows()\n",
    "    s = A.ncols()\n",
    "    for i in range(t):\n",
    "        nonzero_list = [j for j in range(s) if A[i,j] != 0]\n",
    "        if len([j for j in range(s) if A[i,j] != 0])>0:\n",
    "            fi = min([j for j in range(s) if A[i,j] != 0])\n",
    "        else:\n",
    "            fi = None\n",
    "        f.append(fi)\n",
    "    return(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "08f91f17",
   "metadata": {},
   "outputs": [],
   "source": [
    "def finccheck(f):\n",
    "    r\"\"\"\n",
    "    Checks if the sequence f is strictly increasing\n",
    "    \"\"\"\n",
    "    for i in range(1,len(f)):\n",
    "        if f[i-1] == None and f[i] != None:\n",
    "            return(False)\n",
    "        if f[i-1] != None and f[i] != None and f[i] <= f[i-1]:\n",
    "            return(False)\n",
    "    else:\n",
    "        return(True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "933b1e24",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ifionecheck(A):\n",
    "    r\"\"\"\n",
    "    Checks that A is in echelon form and that the first nonzero entry in each onzero row is 1\n",
    "    \"\"\"\n",
    "    f = fnz(A)\n",
    "    if finccheck(f) == True:\n",
    "        r = erank(A)\n",
    "        for i in range(1, r+1):\n",
    "            if A[i-1, f[i]] != 1:\n",
    "                return(False)\n",
    "            else:\n",
    "                return(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "ae9b772f-33b0-49d2-8354-c0899307c8a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def xredcheck(A):\n",
    "    r\"\"\"\n",
    "    This checks that A is in echelon form and that the first nonzero entry in each nonzero row is 1 and that the entries\n",
    "    above the first nonzero in each row are all 0\n",
    "    \"\"\"\n",
    "    f = fnz(A)\n",
    "    if finccheck(f) == False:\n",
    "        print(\"the matrix is not in echelon form\")\n",
    "        return(False)\n",
    "    elif ifionecheck(A) == False:\n",
    "        print(\"the matrix is in echelon form but the pivots are not all 1\")\n",
    "        return(False)\n",
    "    else:\n",
    "        r = erank(A)\n",
    "        #print(f\"r is {r} and f is {f}\")\n",
    "        for i in range(1,r):\n",
    "            for j in range(i):\n",
    "                #print(f\"i is {i} and j is {j} and f[i+1] is {f[i+1]} and the [j,f[i+1]] entry is {A[j,f[i+1]]}\")\n",
    "                if A[j,f[i+1]] != 0:\n",
    "                    return(False)\n",
    "    return(True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "0c0bc7b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def erank(A):\n",
    "    r\"\"\"\n",
    "    This returns the rank of a matrix in echelon form\n",
    "    \"\"\"\n",
    "    f = fnz(A)\n",
    "    if finccheck(f) == True:\n",
    "        r = max([i for i in range(len(f)) if f[i] != None])\n",
    "        return(r)\n",
    "    else:\n",
    "        print(\"A is not in echelon form, first apply step one\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "c997b51f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-1, 0, 0, 0, 2]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fnz(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "4748a424",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "finccheck(fnz(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "7221f24e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "finccheck([-1,3,6,9,None,None])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "559e81c6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "finccheck([-1,3,6,9,9,12])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "cdd96b3c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ifionecheck(B0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "4efafa2d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ifionecheck(C0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "099f37c7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A is not in echelon form, first apply step one\n"
     ]
    }
   ],
   "source": [
    "erank(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "dde622fa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "erank(B0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a17fadc-ec02-4728-aecc-9ae0ba04449b",
   "metadata": {},
   "source": [
    "# Factoring steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "246d8d15",
   "metadata": {},
   "outputs": [],
   "source": [
    "def stepone(A, slist=[]):\n",
    "    r\"\"\"\n",
    "    Iterates step one to produce an echelon matrix and an slist\n",
    "    The slist should default to [] if not specified in the call\n",
    "    \"\"\"\n",
    "    t=A.nrows()\n",
    "    s=A.ncols()\n",
    "    f = fnz(A)\n",
    "    #print(f\"the increasing sequence is{f}\")\n",
    "    if finccheck(f):\n",
    "        #print(\"the output matrix is now in echelon form, execute steptwo\")\n",
    "        return(A,slist)\n",
    "    else:\n",
    "        m = min([f[k] for k in range(1, len(f)) if f[k]!= None and (f[k-1]==None or f[k-1]>=f[k])])\n",
    "        #print(f\"the min column with a pivot is {m}\")\n",
    "        i= max([j for j in range(1, len(f)) if f[j]==m])\n",
    "        #print(f\"the max row with a pivot in column {m} is {i}\")\n",
    "        c = A[i-2,f[i]] / A[i-1,f[i]]\n",
    "        B = sinv(i-1,c,t) * A\n",
    "        slist.append([i-1,c])\n",
    "        #print(f\"B is {B} and slist is {slist}\")\n",
    "        return(stepone(B, slist))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "14e9faa1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def steptwo(B):\n",
    "    r\"\"\"\n",
    "    This takes the matrix B and makes the pivots all 1.\n",
    "    The output is the new matrix and the hlist\n",
    "    \"\"\"\n",
    "    t=B.nrows()\n",
    "    s=B.ncols()\n",
    "    f = fnz(B)\n",
    "    if finccheck(f): \n",
    "        hlist=[] \n",
    "        for i in range(1,len(f)):\n",
    "            if f[i] == None:\n",
    "                B=B\n",
    "                hlist.append(1)\n",
    "            else:\n",
    "                hlist.append(B[i-1,f[i]])\n",
    "                B = hgen(i,1/B[i-1,f[i]],t) * B\n",
    "        return(B,hlist)\n",
    "    else:\n",
    "        print(\"A is not in echelon form, execute step one\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "f92f4927",
   "metadata": {},
   "outputs": [],
   "source": [
    "def stepthree(C):\n",
    "    t=C.nrows()\n",
    "    s=C.ncols()\n",
    "    f = fnz(C)\n",
    "    if finccheck(f):\n",
    "        if ifionecheck(C):\n",
    "            #print(f)\n",
    "            r = max([k for k in range(1,len(f)) if f[k]!=None])\n",
    "            xlist = []\n",
    "            for i in range(r,1,-1):\n",
    "                #print(f\"r and i and f[i] are {r} and {i} and {f[i]}\")\n",
    "                for j in range(i-1,0,-1):\n",
    "                    #print(f\"j and f[i] {j} and {f[i]}\")\n",
    "                    #print(f\"calling xgen({j},{i},{C[j-1,f[i]]},{t})\")\n",
    "                    if C[j-1,f[i]] != 0:\n",
    "                        xlist.append([j,i,C[j-1,f[i]]])\n",
    "                    C = xgen(j,i,-C[j-1,f[i]],t) * C\n",
    "            return(C,xlist)\n",
    "        else:\n",
    "            print(\"A is in echelon form but pivots are not all 1, execute steptwo\")\n",
    "    else:\n",
    "        print(\"A is not in echelon form, execute stepone\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "ebbb8f65-6181-415d-b323-2d8228847a60",
   "metadata": {},
   "outputs": [],
   "source": [
    "def stepfour(R):\n",
    "    t=R.nrows()\n",
    "    s=R.ncols()\n",
    "    f = fnz(R)\n",
    "    if finccheck(f):\n",
    "        if ifionecheck(R):\n",
    "            if xredcheck(R):\n",
    "                #print(f)\n",
    "                r = erank(R)\n",
    "                rxlist = []\n",
    "                for i in range(1,r+1):\n",
    "                    for j in range(f[i]+2,s+1):\n",
    "                        #print(f\"rmult by xgen({f[i]+1},{j},{-R[i-1,j-1]},{s})\")\n",
    "                        if R[i-1,j-1] != 0:\n",
    "                            #This next line was getting the list in the wrong order\n",
    "                            #rxlist.append([f[i]+1,j,R[i-1,j-1]])\n",
    "                            #Lets prepend instead\n",
    "                            rxlist.insert(0,[f[i]+1,j,R[i-1,j-1]])\n",
    "                            R = R * xgen(f[i]+1,j,-R[i-1,j-1],s)\n",
    "                            #print(f\"new R is {R}\")\n",
    "                return(R,rxlist)\n",
    "            else:\n",
    "                print(\"A is in echelon form with pivots all 1 but is not reduced, execute stepthree\")\n",
    "        else:\n",
    "            print(\"A is in echelon form but pivots are not all 1, execute steptwo\")\n",
    "    else:\n",
    "        print(\"A is not in echelon form, execute stepone\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "76626c09",
   "metadata": {},
   "outputs": [],
   "source": [
    "def stepfive(L):\n",
    "    t = L.nrows()\n",
    "    s = L.ncols()\n",
    "    f = fnz(L)\n",
    "    if finccheck(f):\n",
    "        if ifionecheck(L):\n",
    "            if xredcheck(L):\n",
    "                r = erank(L)\n",
    "                #print(f\"r is {r}\")\n",
    "                rslist = []\n",
    "                for i in range(1,r+1):\n",
    "                    #print(f\"i-1 is {i-1} and f[i] is {f[i]} and f[i]>i-1 is {f[i]>i-1}\")\n",
    "                    if f[i]>i-1:\n",
    "                        for j in range(f[i]+1,i,-1):\n",
    "                            #print(f\"j is {j}\")\n",
    "                            L = L * sgen(j-1,0,s)\n",
    "                            #print(f\"the new L is {L}\")\n",
    "                            #This next line was getting the list in the wrong order\n",
    "                            #rslist.append([j-1,0])\n",
    "                            #Lets prepend instead\n",
    "                            rslist.insert(0,[j-1,0])\n",
    "#                    else:\n",
    "#                        print(\"Go on please\")\n",
    "                return(L,rslist)\n",
    "            else:\n",
    "                print(\"A is in echelon form with pivots all 1 but is not reduced, execute stepthree\")\n",
    "        else:\n",
    "            print(\"A is in echelon form but pivots are not all 1, execute steptwo\")\n",
    "    else:\n",
    "        print(\"A is not in echelon form, execute stepone\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "f4013a35-2a18-4539-9768-cbb952c76ff2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def computePQfact(A):\n",
    "    r\"\"\"\n",
    "    Combine steponeiter, steptwo and stepthree to produce\n",
    "    the factorization A = P1_rQ\n",
    "    \"\"\"\n",
    "    t = A.nrows()\n",
    "    s = A.ncols()\n",
    "    B=stepone(A,[])\n",
    "    C=steptwo(B[(0)])\n",
    "    D=stepthree(C[0])\n",
    "    R=stepfour(D[0])\n",
    "    F=stepfive(R[0])\n",
    "    slist = B[1]\n",
    "    slen = len(slist)\n",
    "    S =identity_matrix(QQ, t)\n",
    "    for i in range(slen):\n",
    "        sfactor = slist[i]\n",
    "        S = S*sgen(sfactor[0],sfactor[1],t)\n",
    "    hlist = C[1]\n",
    "    hlen = len(hlist)\n",
    "    H =identity_matrix(QQ, t)\n",
    "    for i in range(hlen):\n",
    "        hfactor = hlist[i]\n",
    "        H = H*hgen(i+1,hfactor,t)\n",
    "    xlist = D[1]\n",
    "    xlen = len(xlist)\n",
    "    X =identity_matrix(QQ, t)\n",
    "    for i in range(xlen):\n",
    "        xfactor = xlist[i]\n",
    "        X = X*xgen(xfactor[0],xfactor[1],xfactor[2],t)\n",
    "    P = S*H*X\n",
    "    rxlist = R[1]\n",
    "    #print(f\"rxlist is {rxlist}\")\n",
    "    rxlen = len(rxlist)\n",
    "    RX =identity_matrix(QQ, s)\n",
    "    for i in range(rxlen):\n",
    "        rxfactor = rxlist[i]\n",
    "        #print(f\"i is {i} and rxfactor is {rxfactor}\")\n",
    "        RX = xgen(rxfactor[0],rxfactor[1],rxfactor[2],s)*RX\n",
    "    rslist = F[1]\n",
    "    rslen = len(rslist)\n",
    "    RS =identity_matrix(QQ, s)\n",
    "    for i in range(rslen):\n",
    "        rsfactor = rslist[i]\n",
    "        RS = RS*sgen(rsfactor[0],rsfactor[1],s)\n",
    "    Q=RS*RX\n",
    "    #checkone = S*B[0]\n",
    "    #print(\"checkone: S*B[0] is\")\n",
    "    #print(f\"{checkone}\")\n",
    "    #checktwo = H*C[0]\n",
    "    #print(\"checktwo: H*C[0] is\")\n",
    "    #print(f\"{checktwo}\")\n",
    "    #checkthree = X*D[0]\n",
    "    #print(\"checkthree: X*D[0] is\")\n",
    "    #print(f\"{checkthree}\")\n",
    "    #checkfour = R[0]*RX\n",
    "    #print(f\"checkfour: R[0]*RX is {checkfour}\")\n",
    "    #P=S*H*X\n",
    "    #checkfive = P*D[0]\n",
    "    #print(\"checkfive: P*D[0] is\")\n",
    "    #print(f\"{checkfive}\")\n",
    "    #checksix = F[0]*RS\n",
    "    #print(\"checksix: F[0]*RS is\")\n",
    "    #print(f\"{checksix}\")\n",
    "    #checkseven = F[0]*Q\n",
    "    #print(\"checkseven: F[0]*Q is\")\n",
    "    #print(f\"{checkseven}\")\n",
    "    print(f\"Computing the factorization of A=\")\n",
    "    print(f\"{A}\")\n",
    "    print(f\"The s(i,c,{t}) factors are {slist}\")\n",
    "    print(f\"The h(i,c,{t}) factors are {hlist}\")\n",
    "    print(f\"The x(i,j,c,{t}) factors are {xlist}\")\n",
    "    print(f\"The right s(i,0,{s}) factors are {rslist}\")\n",
    "    print(f\"The right x(i,j,c.{s}) factors are {rxlist}\")\n",
    "    print(f\"The left orbit representative in reduced row echelon form is R=\")\n",
    "    print(f\"{D[0]}\")\n",
    "    print(f\"The rank of A is r={erank(D[0])}\")\n",
    "    checkeight = P*F[0]*Q\n",
    "    print(\"check: P*F[0]*Q is\")\n",
    "    print(f\"{checkeight}\")\n",
    "    print(\"A,P,1_r and Q are\")\n",
    "    return(A,P,F[0],Q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "11b201aa-2064-49ca-a65f-2df7654a0c79",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing the factorization of A=\n",
      "[    -3    -36    -39    -42    -46    -84]\n",
      "[    -9   -108 -233/2   -114 -247/2   -167]\n",
      "[     3     36     39     42     45     48]\n",
      "[     0      0    1/4      6   25/4   13/2]\n",
      "The s(i,c,4) factors are [[2, -3], [1, -1], [3, 2], [2, 0], [3, 1]]\n",
      "The h(i,c,4) factors are [3, 1/4, -1, 1]\n",
      "The x(i,j,c,4) factors are [[2, 3, 25], [1, 3, 15], [1, 2, 13]]\n",
      "The right s(i,0,6) factors are [[3, 0], [4, 0], [2, 0]]\n",
      "The right x(i,j,c.6) factors are [[5, 6, 36], [3, 6, -874], [3, 4, 24], [1, 6, 10838], [1, 4, -298], [1, 2, 12]]\n",
      "The left orbit representative in reduced row echelon form is R=\n",
      "[    1    12     0  -298     0 10838]\n",
      "[    0     0     1    24     0  -874]\n",
      "[    0     0     0     0     1    36]\n",
      "[    0     0     0     0     0     0]\n",
      "The rank of A is r=3\n",
      "check: P*F[0]*Q is\n",
      "[    -3    -36    -39    -42    -46    -84]\n",
      "[    -9   -108 -233/2   -114 -247/2   -167]\n",
      "[     3     36     39     42     45     48]\n",
      "[     0      0    1/4      6   25/4   13/2]\n",
      "A,P,1_r and Q are\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(\n",
       "[    -3    -36    -39    -42    -46    -84]\n",
       "[    -9   -108 -233/2   -114 -247/2   -167]\n",
       "[     3     36     39     42     45     48]\n",
       "[     0      0    1/4      6   25/4   13/2],\n",
       "\n",
       "[    -3    -39    -46      1]  [1 0 0 0 0 0]\n",
       "[    -9 -233/2 -247/2      0]  [0 1 0 0 0 0]\n",
       "[     3     39     45      0]  [0 0 1 0 0 0]\n",
       "[     0    1/4   25/4      0], [0 0 0 0 0 0],\n",
       "\n",
       "[    1    12     0  -298     0 10838]\n",
       "[    0     0     1    24     0  -874]\n",
       "[    0     0     0     0     1    36]\n",
       "[    0     1     0     0     0     0]\n",
       "[    0     0     0     1     0     0]\n",
       "[    0     0     0     0     0     1]\n",
       ")"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "computePQfact(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1dcf1f5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 10.9",
   "language": "sage",
   "name": "sagemath-10.9"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.14.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}