{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SurfP2P1:\n",
    "    import sympy \n",
    "    ##from sage.symbolic.ring            \n",
    "    import numpy\n",
    "    import numpy as np\n",
    "    from scipy import linalg\n",
    "    from scipy.linalg import eig\n",
    "    from sage.plot.line import Line\n",
    "    from sage.plot.point import Point\n",
    "    from sage.rings.polynomial.polydict import PolyDict\n",
    "    from time import gmtime,strftime,time\n",
    "    from datetime import datetime, date, time,timedelta\n",
    "    \n",
    "    def __init__(self,uvDegree,list,string,field):\n",
    "        assert uvDegree > 0;\n",
    "        self._uvDegree = uvDegree;\n",
    "        self._list=list;\n",
    "        self._string=string;\n",
    "        self._field=field;        \n",
    "\n",
    "        R.<u,v>= self._field['u','v'];\n",
    "        R.gens();          \n",
    "        \n",
    " \n",
    "\n",
    "\n",
    "        if bool(self._list==[]):\n",
    "            \n",
    "            base=((1+u+v)^(self._uvDegree));\n",
    "            Poly0=0;\n",
    "            Poly1=0;\n",
    "            Poly2=0;\n",
    "            Poly3=0;\n",
    "                \n",
    "            for i in range(len(base.monomials())):\n",
    "\n",
    "                if(self._field.cardinality()<infinity):\n",
    "                            \n",
    "                    Poly0=Poly0+(-1)^(field.random_element())*(ZZ.random_element(1,field.cardinality()-1))*base.monomials()[i];\n",
    "                    Poly1=Poly1+(-1)^(field.random_element())*(ZZ.random_element(1,field.cardinality()-1))*base.monomials()[i];\n",
    "                    Poly2=Poly2+(-1)^(field.random_element())*(ZZ.random_element(1,field.cardinality()-1))*base.monomials()[i];\n",
    "                    Poly3=Poly3+(-1)^(field.random_element())*(ZZ.random_element(1,field.cardinality()-1))*base.monomials()[i];\n",
    "                    \n",
    "                else:\n",
    "                    \n",
    "                    ##if we take field = QQ then -1^QQ.random_element() not  good idea to write\n",
    "                                       \n",
    "                    Poly0=Poly0+(-1)^(ZZ.random_element())*(self.randomElement())*base.monomials()[i];\n",
    "                    Poly1=Poly1+(-1)^(ZZ.random_element())*(self.randomElement())*base.monomials()[i];\n",
    "                    Poly2=Poly2+(-1)^(ZZ.random_element())*(self.randomElement())*base.monomials()[i];\n",
    "                    Poly3=Poly3+(-1)^(ZZ.random_element())*(self.randomElement())*base.monomials()[i];    \n",
    "\n",
    "            p=vector(((self._field).random_element(),(self._field).random_element(),(self._field).random_element())); #print \"random point is:\",p;\n",
    "            \n",
    "            \n",
    "            ##In order to keep Polyi's in ZZ, otherwise then are seen in rational field.\n",
    "            if(self._field == ZZ):\n",
    "                Poly0=Poly0.numerator();\n",
    "                Poly1=Poly1.numerator();\n",
    "                Poly2=Poly2.numerator();\n",
    "                Poly3=Poly3.numerator();          \n",
    "    \n",
    "                \n",
    "            if (self._string==\"rational\"):\n",
    "                self._list.append(Poly0);\n",
    "                self._list.append(Poly1);\n",
    "                self._list.append(Poly2);\n",
    "                self._list.append(Poly3);\n",
    "                self._list.append(p);\n",
    "                    \n",
    "            elif (self._string==\"nonrational\"):    \n",
    "                Poly0=R(1);\n",
    "                self._list.append(Poly0);\n",
    "                self._list.append(Poly1);\n",
    "                self._list.append(Poly2);\n",
    "                self._list.append(Poly3);\n",
    "                self._list.append(p); \n",
    "            else:\n",
    "                print \"You can choose either 'rational' or 'nonrational' ! \";  \n",
    "                return \"You can choose either 'rational' or 'nonrational' ! \";    \n",
    "                \n",
    "               \n",
    "    def randomElement(self):\n",
    "        \"\"\"\n",
    "        INPUT: no input.\n",
    "        \n",
    "        OUTPUT : non-zero random element of the field where the class is defined.\n",
    "        \"\"\"\n",
    "        a = 0;\n",
    "        while a==self._field(0):\n",
    "            a =(self._field).random_element();  \n",
    "            #print \"parent of random element is :\",parent(a);\n",
    "        return a;       \n",
    "        \n",
    "        \n",
    "    def uvDegree(self,listOfPolynomials=[]):\n",
    "        \"\"\"\n",
    "        INPUT(optional): listOfPolynomials - list. listOfPolynomials=[] by default.\n",
    "        \n",
    "        OUTPUT : integer - Degree of non-homogeneous parameterization.\n",
    "        \"\"\"\n",
    "        bound=var('bound');\n",
    "        uvMax=[];  \n",
    "        \n",
    "        if listOfPolynomials==[]:\n",
    "            listOfPolynomials=self._list; \n",
    "            bound=4;\n",
    "            if bool(listOfPolynomials[0] in RR)==True:\n",
    "                j=1;       \n",
    "            else:\n",
    "                j=0;\n",
    "        else:  \n",
    "            bound=len(listOfPolynomials);\n",
    "            j=0;\n",
    "        #print \"listOfPolynomials is :\",listOfPolynomials;     \n",
    "        for i in range(j,bound):\n",
    "            #print \"listOfPolynomials[i] is :\",listOfPolynomials[i];\n",
    "            list1=listOfPolynomials[i].exponents();\n",
    "            #print list1;\n",
    "            uvMax.append(max(list1[i][0]+list1[i][1] for i in range (len(list1))));\n",
    "            #print list1[0];\n",
    "            #print list1[0][0];\n",
    "        return max(uvMax[i] for i in range(len(uvMax)));         \n",
    "        \n",
    "    def uvwDegree(self,poly=[]): \n",
    "\n",
    "        \"\"\"\n",
    "        INPUT(optional): poly - list. poly=[] by default.\n",
    "        \n",
    "        OUTPUT : integer- Degree of homogeneous parameterization.\n",
    "        \"\"\"\n",
    "        R.<u,v>=self._field[];R.gens();\n",
    "        R1.<u,v,w,la,mu>=self._field[];R1.gens();\n",
    "        R1.gen(0);\n",
    "        R.gen(0);\n",
    "        u=R1.gen(0);  \n",
    "        uvwMax=[];\n",
    "        #print \"poly0 is :\",poly[0];\n",
    "        #print \"type\",type(poly[0]);\n",
    "        #print \"parent\",parent(poly[0]);\n",
    "        #print \"bakalim neymis\",((R1(poly[0])).exponents());\n",
    "        #print \"bakalim bu neymis\",((R1(poly[0])).exponents())[0];\n",
    "        exponent=((poly[0]).exponents())[0];   \n",
    "        #print \"exponent is :\",exponent;          \n",
    "        uvwMax = exponent[0]+exponent[1]+exponent[2];\n",
    "        return uvwMax;         \n",
    "\n",
    "    \n",
    "    #is only used for the homogenization of the parametrization of the given surface\n",
    "    def homogenizationInR(self,l=[]):\n",
    "        \"\"\"\n",
    "        INPUT(optional): l - list. l=[] by default. \n",
    "        \n",
    "        OUTPUT : list - Gives homogenization of the given list of polynomials in P^2, all \n",
    "                  the polynomials are homogenized in their own degree.\n",
    "        \"\"\"\n",
    "        k=self._field;\n",
    "        R.<u,v>=self._field[];R.gens();\n",
    "        R1.<u,v,w>=self._field[];R1.gens();\n",
    "        R1.gen(0);\n",
    "        R.gen(0);\n",
    "        u=R1.gen(0); \n",
    "        length=0;\n",
    "  \n",
    "        if l==[]:\n",
    "            l=[self._list[0],self._list[1],self._list[2],self._list[3]];\n",
    "            length=3;\n",
    "            uvDegree=self._uvDegree;\n",
    "\n",
    "        else:\n",
    "            \n",
    "            var('deg');\n",
    "            if l[0] in R:\n",
    "                deg=self.uvDegree(l);\n",
    "            if l[0] in R1 :\n",
    "                deg=self.uvwDegree(l);\n",
    "                \n",
    "        homogeneousL=[0 for i in range(len(l))];     \n",
    "\n",
    "        for i in range(len(l)):\n",
    "            monome=l[i].monomials();\n",
    "            coef=l[i].coefficients();\n",
    "            #print \"monomes are :\",monome;\n",
    "            for j in range(len(monome)):\n",
    "                degressOfMonome=monome[j].exponents();\n",
    "                \n",
    "                #print \"degreesOfMonome\",degressOfMonome;\n",
    "\n",
    "                sumDeg=degressOfMonome[0][0]+degressOfMonome[0][1];\n",
    "                \n",
    "                if(sumDeg<uvDegree):\n",
    "                    monome[j]=monome[j]*w^(uvDegree-sumDeg);\n",
    "                homogeneousL[i]=homogeneousL[i]+monome[j]*coef[j];        \n",
    "        return homogeneousL;  \n",
    "\n",
    "\n",
    "    def HomogeneousParametrization(self):\n",
    "        \"\"\"\n",
    "        INPUT: no input. \n",
    "        \n",
    "        OUTPUT : list - Homogenization of the polynomials defining the triangular surface.\n",
    "        \"\"\"\n",
    "        k=self._field;\n",
    "        R.<u,v>= self._field['u','v'];\n",
    "        R1.<u,v,w>=self._field[];R1.gens();\n",
    "        R1.gen(0);\n",
    "        R.gen(0);\n",
    "        u=R1.gen(0);\n",
    "        HomFi=[];\n",
    "        \n",
    "        if self._field==RR:\n",
    "            return self.homogenizationInR();\n",
    "        else:\n",
    "                    \n",
    "            for i in range(4):\n",
    "                a=self._list[i];\n",
    "                #print a.exponents()\n",
    "                expo=a.exponents();\n",
    "                lenExpo=len(a.exponents());\n",
    "                newA=0;\n",
    "                coefA=a.coefficients();\n",
    "                #print \"coef\",coefA;\n",
    "                for i in range(lenExpo):\n",
    "                    #print \"expo[i][0]\",expo[i][0];\n",
    "                    difference=self._uvDegree-(expo[i][0]+expo[i][1]);\n",
    "                    #print difference;\n",
    "                    if difference==0:\n",
    "                        newA=newA+coefA[i]*u^(expo[i][0])*v^(expo[i][1]);\n",
    "                    if difference!=0:\n",
    "                        #print \"coefA[i]\",coefA[i];\n",
    "                        newA=newA+coefA[i]*u^(expo[i][0])*v^(expo[i][1])*w^(difference);\n",
    "                        #print \"newA\",newA; \n",
    "                HomFi.append(newA);        \n",
    "                #this does not work with finite fields in SageMath 8.6        \n",
    "                #print expand((self._list[i]*(w^self._uvDegree)).subs(u=u/w,v=v/w));\n",
    "                #HomFi.append(expand((self._list[i]*(w^self._uvDegree)).subs(u=u/w,v=v/w)));\n",
    "            return HomFi \n",
    "\n",
    "        \n",
    "    #NormalVectorsViaJacobian function works only for the objects of this class.... because we take s,u Degree of the object inside of the function \n",
    "    def NormalVectorsViaJacobian(self,Fi=[]):\n",
    "        \"\"\"\n",
    "        INPUT(optional): Fi-list, Fi=[] by default, Fi: list of homogeneous polynomials \n",
    "                         defining the triangular surface. \n",
    "        \n",
    "        OUTPUT : list - Homogeneous normal vectors with respect to the given homogeneous parameterization.\n",
    "        \"\"\"               \n",
    "        k=self._field;\n",
    "        R.<u,v>= self._field['u','v'];\n",
    "        R.gens(); \n",
    "        R1.<u,v,w>=self._field[];R1.gens();\n",
    "        R1.gen(0);\n",
    "        R.gen(0);\n",
    "        #u=R1.gen(0); \n",
    "        field=RR;\n",
    "        if k==RR:\n",
    "            field=CDF;\n",
    "        else:\n",
    "            field=k;\n",
    "        R2.<T0,T1,T2,T3,u,v,w>=field[];R2.gens();   \n",
    "        u=R2.gen(4); \n",
    "                       \n",
    "        string1=self._string;\n",
    "        d=var('d');minor0=var('minor0');minor1=var('minor1');minor2=var('minor2');minor3=var('minor3');\n",
    "\n",
    "        K=R1.fraction_field();\n",
    "        \n",
    "            \n",
    "        M=matrix(R2,4,4,[[(Fi[0]).gradient()[0],(Fi[1]).gradient()[0],\n",
    "                              (Fi[2]).gradient()[0],(Fi[3]).gradient()[0]],\n",
    "                             [(Fi[0]).gradient()[1],(Fi[1]).gradient()[1],\n",
    "                              (Fi[2]).gradient()[1],(Fi[3]).gradient()[1]],\n",
    "                             [(Fi[0]).gradient()[2],(Fi[1]).gradient()[2],\n",
    "                              (Fi[2]).gradient()[2],(Fi[3]).gradient()[2]],[T0,T1,T2,T3]]);    \n",
    "       \n",
    "\n",
    "\n",
    "        \n",
    "        #print \"M is :\",M;\n",
    "        \n",
    "        d=det(M);\n",
    "        #print \"d is :\",d;\n",
    "                \n",
    "        minor0=d.subs(T0=1,T1=0,T2=0,T3=0);\n",
    "        minor1=d.subs(T1=1,T0=0,T2=0,T3=0);\n",
    "        minor2=d.subs(T2=1,T0=0,T1=0,T3=0);\n",
    "        minor3=d.subs(T3=1,T0=0,T1=0,T2=0);\n",
    "        \n",
    "        #print \"0 :\",minor0;     \n",
    "        #print \"1 :\",minor1;              \n",
    "        #print \"2 :\",minor2;              \n",
    "        #print \"3 :\",minor3; \n",
    "        \n",
    "        normalVector=[0 for i in range(3)];\n",
    "        \n",
    "        if string1==\"nonrational\":\n",
    " \n",
    "            dw=self._uvDegree;\n",
    "            #print \"dw is :\",dw;\n",
    "            u=R1.gen(0);  \n",
    "            \n",
    "            minors=[R1(minor1),R1(minor2),R1(minor3)];                \n",
    "            for i in range(3):\n",
    "                monomes=minors[i].monomials();\n",
    "                #print \"monomes\",monomes;\n",
    "                coef=minors[i].coefficients();\n",
    "                #print \"coefficients are:\",coef;                    \n",
    "                for j in range(len(monomes)):  \n",
    "                    degreesOfMonomial=monomes[j].exponents();\n",
    "                    #print \"degreesOfMonomial\",degreesOfMonomial;\n",
    "                    if degreesOfMonomial[0][2]>=dw-1:\n",
    "                        normalVector[i]=normalVector[i]+coef[j]*u^(degreesOfMonomial[0][0])*v^(degreesOfMonomial[0][1])*w^(degreesOfMonomial[0][2]-dw+1);\n",
    "            return [R1(normalVector[0]),R1(normalVector[1]),R1(normalVector[2])] \n",
    "\n",
    "        else:\n",
    "            \n",
    "            dw=self._uvDegree;   \n",
    "                        \n",
    "            return [R1(minor1),R1(minor2),R1(minor3)];       \n",
    "\n",
    " \n",
    "            \n",
    "            \n",
    "    def NormalVectorsViaCrossProduct(self,writeDetails=False):  \n",
    "        #here we add string1 as input in case if we want to use the function for a polynomial different from of the instance of the class\n",
    "        \"\"\"\n",
    "        INPUT(optional): writeDetails-bool, writeDetails=False by default. \n",
    "        \n",
    "        OUTPUT : vector - non homogeneous normal vectors computed via the cross product of the partial\n",
    "                  derivatives of the parameterization of the given triangular surface.\n",
    "                  If writeDetails==True, it writes the two partial derivatives and their cross product.\n",
    "        \"\"\"         \n",
    "        \n",
    "        string1=self._string;\n",
    "        K=R.fraction_field();\n",
    "        \n",
    "        if self._field==RR:\n",
    "            \n",
    "            vec=[];\n",
    "            print \"type(self._list[0])\",type(self._list[0]);\n",
    "            if self._string==\"nonrational\":\n",
    "                dsf0=0;\n",
    "                duf0=0;\n",
    "            else :\n",
    "                dsf0=self._list[0].gradient()[0];\n",
    "                duf0=self._list[0].gradient()[1];                \n",
    "                \n",
    "\n",
    "            dsf1=self._list[1].gradient()[0];\n",
    "            dsf2=self._list[2].gradient()[0];\n",
    "            dsf3=self._list[3].gradient()[0];\n",
    "            \n",
    "\n",
    "            duf1=self._list[1].gradient()[1];\n",
    "            duf2=self._list[2].gradient()[1];\n",
    "            duf3=self._list[3].gradient()[1];\n",
    "            \n",
    "            f0=self._list[0];\n",
    "            f1=self._list[1];\n",
    "            f2=self._list[2];\n",
    "            f3=self._list[3];\n",
    "            \n",
    "            vec.append(expand(-dsf3*duf2*f0 + dsf2*duf3*f0 + dsf3*duf0*f2 - dsf0*duf3*f2 -dsf2*duf0*f3 + dsf0*duf2*f3));\n",
    "            vec.append(expand(dsf3*duf1*f0 - dsf1*duf3*f0 - dsf3*duf0*f1 + dsf0*duf3*f1 + dsf1*duf0*f3- dsf0*duf1*f3));\n",
    "            vec.append(expand(-dsf2*duf1*f0 + dsf1*duf2*f0 + dsf2*duf0*f1 - dsf0*duf2*f1 -dsf1*duf0*f2 + dsf0*duf1*f2));\n",
    "            \n",
    "            return vec;\n",
    "        \n",
    "        if (string1==\"rational\"):\n",
    "            \n",
    "            Phi=vector([self._list[1]/self._list[0], self._list[2]/self._list[0], self._list[3]/self._list[0]]);            \n",
    "          \n",
    "        else:\n",
    "            \n",
    "            Phi=vector([self._list[1], self._list[2], self._list[3]]);                    \n",
    "        \n",
    "\n",
    "        delta_u = vector([simplify(Phi[0].derivative(u)), simplify(Phi[1].derivative(u)), simplify(Phi[2].derivative(u))]);\n",
    "        delta_v = vector([simplify(Phi[0].derivative(v)), simplify(Phi[1].derivative(v)), simplify(Phi[2].derivative(v))]);\n",
    "        Delta=simplify((delta_u).cross_product(delta_v)); #by saying simplify, we simplify the cross product by the factor common f0 \n",
    "        Num=[simplify(Delta[0].numerator()),simplify(Delta[1].numerator()),simplify(Delta[2].numerator())];                                      \n",
    "                       \n",
    "        if (writeDetails==True):\n",
    "            print \"Partial derivatives of the parametrization wrt u  are :\",delta_u;\n",
    "            print \"Partial derivatives of the parametrization wrt v  are :\",delta_v;            \n",
    "            #print \"Cross product of partial derivatives wrt u and v is:\",Delta;\n",
    "            return R(Num[0]),R(Num[1]),R(Num[2])          \n",
    "        elif (writeDetails==False):                            \n",
    "            return R(Num[0]),R(Num[1]),R(Num[2])                \n",
    "        else:\n",
    "            print \"You can choose either 'True' or 'False' \"; \n",
    "            return \"You can choose either 'True' or 'False' \"; \n",
    "            \n",
    "            \n",
    "    def HomogeneousCongruenceOfNormalLines(self,writeDetails=False,method=\"cross_product\"):\n",
    "        \"\"\"\n",
    "        INPUT(optional): writeDetails-bool, writeDetails=False by default. \n",
    "                         method-string : either jacobian matrix or cross_product computation for the \n",
    "                         congruence of normal lines.\n",
    "        \n",
    "        OUTPUT : numpy.array - homogeneous parameterization of normal lines to the given triangular surface\n",
    "                  If writeDetails==True, it writes at which degree it computes the homogenization.\n",
    "        \"\"\"            \n",
    "        if method!=\"cross_product\" and method!=\"jacobian\":\n",
    "            print \"You can choose either 'cross_product' or 'jacobian' \";  \n",
    "            return \"You can choose either 'True' or 'False' \";    \n",
    "\n",
    "            \n",
    "        c=[];\n",
    "        b=var('b,uvDeg');\n",
    "        nV=[];\n",
    "        if method==\"cross_product\" :\n",
    "            R.<u,v>=self._field[];R.gens();\n",
    "            u=R.gen(0);\n",
    "            R1.<u,v,w,la,mu>=self._field[];R1.gens();\n",
    "            R1.gen(0);\n",
    "            u=R1.gen(0); \n",
    "            K1=R1.fraction_field();K1.gens();  \n",
    "            \n",
    "            nV=self.NormalVectorsViaCrossProduct(False);\n",
    "                        \n",
    "                                   \n",
    "            c.append(simplify(mu*(R1(self._list[0]))));\n",
    "            c.append(simplify(mu*R1(self._list[1]))+(la*R1(nV[0])));\n",
    "            c.append(simplify(mu*R1(self._list[2]))+(la*R1(nV[1])));\n",
    "            c.append(simplify(mu*R1(self._list[3]))+(la*R1(nV[2])));\n",
    "\n",
    "            uvDeg=(self.uvDegree(nV));\n",
    "\n",
    "            #print \"FOR CROSS_PRODUCT - normal vector is of degree :\",uvDeg;\n",
    "            \n",
    "            if self._field==RR:\n",
    "                b=[0,0,0,0];    \n",
    "                for i in range(len(c)):\n",
    "                    monome=c[i].monomials();\n",
    "                    coef=c[i].coefficients();\n",
    "                    for j in range(len(monome)):\n",
    "                        exponents=monome[j].exponents();\n",
    "                        sumDeg=exponents[0][0]+exponents[0][1];\n",
    "                        uvDeg=max(uvDeg,sumDeg);\n",
    "                        if sumDeg<uvDeg:\n",
    "                            monome[j]=monome[j]*w^(uvDeg-sumDeg);\n",
    "                        b[i]=b[i]+monome[j]*R1(coef[j]);    \n",
    "                                    \n",
    "            else:\n",
    "                uvDeg=max(uvDeg,self._uvDegree);                 \n",
    "                b= numpy.array([((w^uvDeg*c[0]).subs(u=u/w,v=v/w)).numerator(),((w^uvDeg*c[1]).subs(u=u/w,v=v/w)).numerator(),((w^uvDeg*c[2]).subs(u=u/w,v=v/w)).numerator(),((w^uvDeg*c[3]).subs(u=u/w,v=v/w)).numerator()]);\n",
    "            #print parent(b[1]);\n",
    "            \n",
    "            c=[];  \n",
    "            nV=[];\n",
    "            \n",
    "        else: #string2==\"jacobian\":\n",
    "             #R.<u,v>=self._field[];R.gens();\n",
    "             #s=R.gen(0);             \n",
    "            R1.<u,v,w,la,mu>=self._field[];R1.gens();   \n",
    "            R1.gen(0);\n",
    "            u=R1.gen(0); \n",
    "            HomFi=self.HomogeneousParametrization();\n",
    "            nV=self.NormalVectorsViaJacobian(HomFi); \n",
    "            \n",
    "            c.append(simplify(mu*(self._list[0])));\n",
    "            c.append(simplify(mu*(self._list[1])));  #+la*nV[0])));\n",
    "            c.append(simplify(mu*(self._list[2])));  #+la*nV[1])));\n",
    "            c.append(simplify(mu*(self._list[3])));  #+la*nV[2])));\n",
    "             \n",
    "            uvDeg=(self.uvwDegree(nV));      \n",
    "            #print \"FOR JACOBIAN - normal vector is of degree :\",uvDeg;\n",
    "              \n",
    "            T=True;\n",
    "            if T==True:\n",
    "            #if self._field==RR:\n",
    "                b=[0,0,0,0];\n",
    "                \n",
    "                #print \"REELLER\";\n",
    "                #print \"uvDeg is :\",uvDeg;\n",
    "                \n",
    "                for i in range(len(c)):\n",
    "                    monome=c[i].monomials();\n",
    "                    coef=c[i].coefficients();\n",
    "                    for j in range(len(monome)):\n",
    "                        exponents=monome[j].exponents();\n",
    "                        sumDeg=exponents[0][0]+exponents[0][1];\n",
    "                        m=max(uvDeg,self._uvDegree);\n",
    "                        if sumDeg<=m:\n",
    "                            monome[j]=monome[j]*w^(m-sumDeg);\n",
    "                        b[i]=b[i]+monome[j]*R1(coef[j]);  \n",
    "                \n",
    "                    \n",
    "                if uvDeg>=self._uvDegree:            \n",
    "                    b[1]=b[1]+la*nV[0];\n",
    "                    b[2]=b[2]+la*nV[1];\n",
    "                    b[3]=b[3]+la*nV[2];\n",
    "                else:\n",
    "                    b[1]=b[1]+la*nV[0]*w;\n",
    "                    b[2]=b[2]+la*nV[1]*w;\n",
    "                    b[3]=b[3]+la*nV[2]*w;\n",
    "            else:\n",
    "            \n",
    "                #print \"uvDeg is :\",uvDeg;\n",
    "                uvDeg=max(uvDeg,self._uvDegree); \n",
    "                #print \"uvDeg is :\",uvDeg;\n",
    "                 \n",
    "                if self._uvDegree==1:\n",
    "                                      \n",
    "                    b = numpy.array([((w^uvDeg*(c[0])).subs(u=u/w,v=v/w)).numerator(),((w^uvDeg*(c[1])).subs(u=u/w,v=v/w)).numerator(),((w^uvDeg*(c[2])).subs(u=u/w,v=v/w)).numerator(),((w^uvDeg*(c[3])).subs(u=u/w,v=v/w)).numerator()]);   \n",
    "                    b[1]=b[1]+la*nV[0]*w;\n",
    "                    b[2]=b[2]+la*nV[1]*w;\n",
    "                    b[3]=b[3]+la*nV[2]*w;                 \n",
    "                else:\n",
    "                     \n",
    "                    b = numpy.array([((w^uvDeg*c[0]).subs(u=u/w,v=v/w)).numerator(),((w^uvDeg*c[1]).subs(u=u/w,v=v/w)).numerator(),((w^uvDeg*c[2]).subs(u=u/w,v=v/w)).numerator(),((w^uvDeg*c[3]).subs(u=u/w,v=v/w)).numerator()]);                            \n",
    "                    b[1]=b[1]+la*nV[0];\n",
    "                    b[2]=b[2]+la*nV[1];\n",
    "                    b[3]=b[3]+la*nV[2];\n",
    "                                       \n",
    "                c=[];\n",
    "                nV=[];\n",
    "                     \n",
    "        if (writeDetails==True):\n",
    "            print \"We homogenize in degree :\",uvDeg;\n",
    "            return b;\n",
    "        elif (writeDetails==False):\n",
    "            return b;  \n",
    "        else: \n",
    "            print \"You can choose either 'True' or 'False' \";   \n",
    "            return \"You can choose either 'True' or 'False' \";   \n",
    "             \n",
    "            \n",
    "    def regularity(self,d): #it takes a and the homogeneous parametrization as an input\n",
    "        \"\"\"\n",
    "        INPUT: 1)d - integer, d:degree of the non-homogeneous parameterization of the choosen surface.\n",
    "        \n",
    "        OUTPUT : integer - the threshold degree for the regularity\n",
    "        \"\"\"    \n",
    "        if self._string==\"rational\":\n",
    "            return 9*d-11;\n",
    "        if self._string==\"nonrational\":\n",
    "            return 6*d-8;\n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "\n",
    "    def SNu(self,Nu,D,writeDetails=False,method=\"cross_product\",Psi=[]):\n",
    "        \"\"\"\n",
    "        INPUT: 1)Nu-integer:regularity \n",
    "               2)D-integer:degree of the homogeneous parameterization of the congruence of normal lines.\n",
    "               3)(optional) writeDetails-bool\n",
    "               4)(optional) method-string : either \"cross_product\" or \"jacobian\".\n",
    "               5)(optional) given homogeneous parameterization in degree D in variables u,v,w, in degree 1 in la,mu\n",
    "                 over the field of the class.\n",
    "        \n",
    "        OUTPUT : matrix - intermediary matrix whose kernel is the matrix representation of the homogeneous parameterization\n",
    "                  of the normal lines to the surface.\n",
    "                  If writeDetails==True, then it writes rank and the size of the intermediary matrix.\n",
    "        \"\"\" \n",
    "        count=0;\n",
    "        R.<u,v>= self._field['u','v']; \n",
    "        R1.<u,v,w,la,mu>=self._field[];R1.gens();\n",
    "        R1.gen(0);\n",
    "        R.gen(0);\n",
    "        u=R1.gen(0);\n",
    "        \n",
    "        if Psi==[]:\n",
    "              \n",
    "            if(method==\"jacobian\"):\n",
    "                Psi=(self.HomogeneousCongruenceOfNormalLines(False,\"jacobian\")); \n",
    "            elif(method==\"cross_product\"):              \n",
    "                Psi=(self.HomogeneousCongruenceOfNormalLines(False));\n",
    "            else:\n",
    "                return \"You can choose among jacobian or cross_product\";\n",
    "        \n",
    "        temp1 = var('temp1');        \n",
    "        \n",
    "        #print \"where :\",parent(Psi[1]);\n",
    "        \n",
    "        S_Nu = matrix(self._field,((Nu+D+1)*(Nu+D+2)),(2*(Nu+1)*(Nu+2)),0);\n",
    "        a=(D+Nu+1)*(D+Nu+2)/2;\n",
    "        for l in range(4): # counts Psi_i s.\n",
    "            #print \"Psi no :\",l;\n",
    "            for i in range(Nu+1):# counts the degree of u in the monominal bases that we multiply Psi_i s with in terms of u,v,w\n",
    "                for k in range(i+1): # counts degree of v in the monomial bases that we multiplly Psi_i s with in terms of u,v,w\n",
    "                    temp1=((Psi[l])*u^((Nu)-i)*v^(k)*w^(i-k));\n",
    "                 \n",
    "                    #print temp1;\n",
    "                    #t1=PolyDict(temp1.dict());\n",
    "                    list1=temp1.exponents(); # gives the degrees respectively in u,v,w,la,mu \n",
    "                    #print \"list1 is :\",list1;\n",
    "                    coef1=temp1.coefficients(); #exponents and coefficients give the outputs in the same order\n",
    "                    #print \"coef1 is :\",coef1;\n",
    "                    #print \"length of coef1 is:\",len(coef1);\n",
    "                    #print \"coef1 is :\",coef1;\n",
    "                    #print \"length list1 is :\",len(list1);\n",
    "                 \n",
    "                    for h in range(len(coef1)):\n",
    "                        #print \"h is\",h;\n",
    "                        if (list1[h][3]==0): #if the degree of lambda=0, we compute the mu block.\n",
    "                            #print \"list in elemani ve coef 1 mu block:\",list1[h],coef1[h];\n",
    "                            S_Nu[a+(D+Nu-list1[h][0])*(D+Nu-list1[h][0]+1)/2+list1[h][1],count] = coef1[h];\n",
    "                            #S_Nu[a+(D+Nu-(list1[h][0]+1)+1)*(D+Nu-(list1[h][0]+1)+2)/2+k,l*(Nu +1)*(Nu+2)/2+i*(i+1)/2+k] = coef1[h];\n",
    "                            #print \"row is :\",a+(D+Nu-list1[h][0])*(D+Nu-list1[h][0]+1)/2+list1[h][1];\n",
    "                            #print \"column is :\",2+i*(i+1)/2+k,count;\n",
    "                            #print \"coef is:\",coef1[h];\n",
    "                            #print (S_Nu.transpose()[l*(Nu +1)*(Nu+2)/2+i*(i+1)/2+list1[h][1]]);\n",
    "                         \n",
    "                        if (list1[h][4]==0): #if the degree of mu=0, we compute the lambda block.\n",
    "                            #print \"list in elemani ve coef 1 lambda block:\",list1[h],coef1[h];\n",
    "                            S_Nu[(D+Nu-(list1[h][0]))*(D+Nu-(list1[h][0])+1)/2+list1[h][1],count] = coef1[h];\n",
    "                            #S_Nu[(D+Nu-(list1[h][0]+1)+1)*(D+Nu-(list1[h][0]+1)+2)/2+k,l*(Nu +1)*(Nu+2)/2+i*(i+1)/2+k] = coef1[h]; \n",
    "                            #print \"row is :\",(D+Nu-(list1[h][0]))*(D+Nu-(list1[h][0])+1)/2+list1[h][1]; #(D+Nu-list1[h][0])*(D+Nu-list1[h][0]+1)/2+list1[h][1];\n",
    "                            #print \"column is :\",2+i*(i+1)/2+k,count; \n",
    "                            #print \"coef is:\",coef1[h];  \n",
    "                            #print S_Nu.transpose()[l*(Nu +1)*(Nu+2)/2+i*(i+1)/2+list1[h][1]];\n",
    "                    list1=0;\n",
    "                    coef1=0; \n",
    "                    count=count+1; \n",
    "        \n",
    "        if (writeDetails==True):\n",
    "            print \"rank of S_Nu is :\",((S_Nu)).rank();\n",
    "            print \"size of S_nu is :\",S_Nu.nrows(),S_Nu.ncols();\n",
    "            return S_Nu\n",
    "        elif (writeDetails==False):\n",
    "            return S_Nu\n",
    "        else:\n",
    "            return \"You can choose either 'True' or 'False' \";                \n",
    "                \n",
    "    def matrixRepresentation(self,Nu,D,Psi,writeDetails=False,p=[],M=[],method=\"cross_product\"):\n",
    "        \"\"\"\n",
    "        INPUT: 1)Nu-integer: regularity,\n",
    "               2)D-integer: degree of the homogeneous parameterization of congruence of normal lines,\n",
    "               3)Psi-list : list of homogeneous parameterization of congruence of normal lines,\n",
    "               3)(optional) writeDetails-bool, writeDetails=False by default,\n",
    "               4)(optional) p-[float,float,float], p=[] by default,\n",
    "               5)(optional) M-list : intermediary matrix, can be computed seperatedly via the command .Snu \n",
    "               in this module,\n",
    "               M-matrix : intermediary matrix,\n",
    "               6)method-string, either method=\"cross_product\" or method=\"jacobian\", method=\"cross_product\" \n",
    "               by default.\n",
    "               \n",
    "        OUTPUT : numpy.matrix,integer,(float,float,float),timedate() - matrix representation, denoted by Mrep, of the congruence of normal lines to the\n",
    "                 given triangular surface evaluated at point p, the corank of Mrep, point p in R^3,\n",
    "                 time of computation.\n",
    "        \"\"\" \n",
    "        if p==[]:\n",
    "            p=self._list[4];\n",
    "        if M==[]:\n",
    "            M=self.SNu(Nu,D,writeDetails,method,Psi) ;\n",
    "        (mrepG,time1)=self.foncMrepG(M,writeDetails);\n",
    "        # [m0,m1,m2,m3],MrepG3-MrepG1\n",
    "        mrep=self.foncMrepEval(mrepG,writeDetails);\n",
    "        #M_rep,M_rep.shape[0]-RankMrepatP,p;\n",
    "        return mrep,time1;\n",
    "                \n",
    "    def verificationSNu(self,m,D,Nu,method=\"cross_product\"):\n",
    "        \"\"\"\n",
    "        INPUT: 1)m-matrix : intermediary matrix\n",
    "               2)D-integer: degree of the parameterization\n",
    "               3)Nu-integer: regularity\n",
    "               4)(optional) method-string: either \"cross_product\" or \"jacobian\", \"cross_product\" by default.\n",
    "               \n",
    "        OUTPUT : bool - construction of the intermediary matrix is whether True.\n",
    "        \"\"\" \n",
    "        \n",
    "        R1.<u,v,w,la,mu>= ZZ['u','v','w','la','mu'];R1.gens();\n",
    "        u=R1.gen(0);\n",
    "        \n",
    "        ColList=[];LineList=[];\n",
    "        \n",
    "        if(method==\"jacobian\"):\n",
    "            Psi=(self.HomogeneousCongruenceOfNormalLines(False,\"jacobian\"));   \n",
    "        elif(method==\"cross_product\"): \n",
    "            R.<u,v>=self._field[];R.gens();\n",
    "            u=R.gen(0);\n",
    "            Psi=(self.HomogeneousCongruenceOfNormalLines(False));\n",
    "        else:\n",
    "            print \"You can choose either jacobian or cross_product\";\n",
    "            return \"You can choose either jacobian or cross_product\";\n",
    "        \n",
    "        #print \"parent of Psi is :\",parent(psi);\n",
    "      \n",
    "        for l in range(4): # counts Psi_i s.\n",
    "            for i in range((Nu)+1):# counts the degree of u in the monominal bases that we multiply Psi_i s with in terms of u,v,w\n",
    "                for k in range(i+1): # counts degree of v in the monomial bases that we multiplly Psi_i s with in terms of u,v,w\n",
    "                    temp1=((Psi[l])*u^((Nu)-i)*v^(k)*w^(i-k));\n",
    "                    ColList.append(temp1);\n",
    "        #print \"Collist is:\",ColList;\n",
    "        for l in range(2): #for lambda and for mu\n",
    "            for i in range((Nu+D)+1):# counts the degree of u in the monominal bases that we multiply Psi_i s with in terms of u,v,w\n",
    "                for k in range(i+1): # counts degree of v in the monomial bases that we multiplly Psi_i s with in terms of u,v,w\n",
    "                    if l==0:\n",
    "                        temp2=u^((Nu+D)-i)*v^(k)*w^(i-k)*la;\n",
    "                    if l==1:\n",
    "                        temp2=u^((Nu+D)-i)*v^(k)*w^(i-k)*mu;                \n",
    "                    LineList.append(temp2);    \n",
    "                    \n",
    "                    \n",
    "                    \n",
    "        #print \"LineList is:\",LineList;\n",
    "        print \"lenght of LineList is :\",len(LineList);\n",
    "        print \"lenght of ColList is :\",len(ColList);        \n",
    "        #print \"parent of LHS:\", parent(vector(ColList));\n",
    "        #print \"parent of RHS:\", parent(Snu.T*vector(LineList));\n",
    "        \n",
    "        return vector(ColList)==m.T*vector(LineList);     \n",
    "        \n",
    "        \n",
    "    def foncMrepGExact(self,SN,writeDetails=False):\n",
    "        \"\"\"\n",
    "        INPUT: 1)SN-matrix : intermediary matrix\n",
    "               2)(optional) writeDetails-bool : writeDetails=False by default.\n",
    "               \n",
    "        OUTPUT : list of 4 matrices - corresponding to a point in P^3.\n",
    "                If writeDetails==True, then it writes, if it is exact computation, the rank and size of \n",
    "                the intermediary matrix,the size of its kernel. Also, it writes the size of the output matrices,\n",
    "                the computation time.\n",
    "                This method only works for numeric computations, i.e, for the computations over RR.\n",
    "        \"\"\"         \n",
    "        MrepG1=datetime.now();        \n",
    "        \n",
    "        kernel=SN.right_kernel(basis='pivot').basis_matrix().transpose();\n",
    "        #print \"kernel is size:\",kernel.nrows(),kernel.ncols();\n",
    "        MrepG2=datetime.now();       \n",
    "        \n",
    "        nCs=kernel.nrows()/4;\n",
    "        \n",
    "        m0=kernel[[i for i in range(0,nCs)],:]; #print \"m0 is:\",m0;\n",
    "        m1=kernel[[i for i in range(nCs,2*nCs)],:]; #print \"m1 is:\",m1;\n",
    "        m2=kernel[[i for i in range(2*nCs,3*nCs)],:]; #print \"m2 is:\",m2;\n",
    "        m3=kernel[[i for i in range(3*nCs,4*nCs)],:]; #print \"m3 is:\",m3;  \n",
    "        \n",
    "\n",
    "        MrepG3=datetime.now();        \n",
    "        if (writeDetails==True):\n",
    "            print \"Exact computation !\"\n",
    "            print \"SNu is of rank:\", SN.rank();\n",
    "            print \"Kernel of SNu is of size :\",kernel.nrows(),kernel.ncols();\n",
    "            print \"Size of M-rep is:\",m1.nrows(),m1.ncols();  \n",
    "            print \"Time to compute the kernel of Snu is :\",MrepG2-MrepG1;\n",
    "            print \"Time to compute the Mrep, without computing the kernel of Snu is :\",MrepG3-MrepG2;\n",
    "            print \"TotalComputationTime is:\",MrepG3-MrepG1;            \n",
    "            return [m0,m1,m2,m3],MrepG3-MrepG1\n",
    "        elif (writeDetails==False): \n",
    "            return [m0,m1,m2,m3],MrepG3-MrepG1\n",
    "        else:\n",
    "            print \"You can choose either 'True' or 'False' !\";\n",
    "            return \"You can choose either 'True' or 'False' !\";   \n",
    "            \n",
    "            \n",
    "            \n",
    "    def foncMrepG(self,SN,writeDetails=False):\n",
    "        \"\"\"\n",
    "        INPUT: 1)SN-matrix : intermediary matrix\n",
    "               2)(optional) writeDetails-bool : writeDetails=False by default.\n",
    "               \n",
    "        OUTPUT : list of 4 matrices - corresponding to a point in P^3.\n",
    "                If writeDetails==True, then it writes, if it is exact computation, the rank and size of \n",
    "                the intermediary matrix,the size of its kernel. Also, it writes the size of the output matrices,\n",
    "                the computation time.\n",
    "                This method only works for EXACT computations, i.e, for the computations over QQ or finite fields.\n",
    "        \"\"\"           \n",
    "        #T=True;\n",
    "        #if T=True:\n",
    "        if(self._field.cardinality()<infinity) or self._field==QQ:\n",
    "            #print \"Exact computation\";\n",
    "            return self.foncMrepGExact(SN,string);\n",
    "        #if T=False:\n",
    "        else:    \n",
    "            \n",
    "\n",
    "            #print strftime(\"INSIDE 1, %a, %d %b %Y %H:%M:%S \", gmtime());\n",
    "        \n",
    "            RankSN=self.rankCorankViaNorm2(SN,False)[0];\n",
    "            MrepG1=datetime.now();        \n",
    "            #print strftime(\"INSIDE 2, %a, %d %b %Y %H:%M:%S \", gmtime());\n",
    "        \n",
    "            kernel=self.numericalKernel(SN,RankSN,False);\n",
    "            MrepG2=datetime.now();        \n",
    "            #print strftime(\"INSIDE 3, %a, %d %b %Y %H:%M:%S \", gmtime());\n",
    "        \n",
    "            nCs=kernel.shape[0]/4;\n",
    "        \n",
    "            m0=kernel[[i for i in range(0,nCs)],:]; #print \"m0 is:\",m0;\n",
    "            m1=kernel[[i for i in range(nCs,2*nCs)],:]; #print \"m1 is:\",m1;\n",
    "            m2=kernel[[i for i in range(2*nCs,3*nCs)],:]; #print \"m2 is:\",m2;\n",
    "            m3=kernel[[i for i in range(3*nCs,4*nCs)],:]; #print \"m3 is:\",m3;  \n",
    "        \n",
    "            MrepG3=datetime.now();        \n",
    "        if (writeDetails==True):\n",
    "            print \"SNu is of rank:\", RankSN;\n",
    "            print \"Kernel of SNu is of size :\",kernel.shape[0],kernel.shape[1];\n",
    "            print \"Size of M-rep is:\",nCs,kernel.shape[1];  \n",
    "            print \"Time to compute the kernel of Snu is :\",MrepG2-MrepG1;\n",
    "            print \"Time to compute the Mrep, without computing the kernel of Snu is :\",MrepG3-MrepG2;\n",
    "            print \"TotalComputationTime is:\",MrepG3-MrepG1;            \n",
    "            return [m0,m1,m2,m3],MrepG3-MrepG1\n",
    "        elif (writeDetails==False): \n",
    "            return [m0,m1,m2,m3],MrepG3-MrepG1\n",
    "        else:\n",
    "            print \"You can choose either 'True' or 'False' !\";\n",
    "            return \"You can choose either 'True' or 'False' !\";             \n",
    "            \n",
    "                     \n",
    "    def rankCorank(self,M,writeDetails=False): #output : rank, corank\n",
    "        \"\"\"\n",
    "        INPUT: 1)M-matrix,\n",
    "               2)(optional) writeDetails-bool : writeDetails=False by default.\n",
    "               \n",
    "        OUTPUT : (integer,integer) = rank, corank of the input matrix.\n",
    "                If writeDetails==True, then if writes the singular values of the given matrix, the size of the\n",
    "                given matrix.\n",
    "                This method uses SVD for numeric computations, computations over RR,\n",
    "                .rank, and .right_kernel() commands of SageMath for exact computations, i.e, the computations over QQ.\n",
    "                It uses norm infinity* 10^(-8) as tolerence.\n",
    "        \"\"\"   \n",
    "        if self._field==RR:\n",
    "            M=numpy.matrix(M);\n",
    "            Sn= linalg.svd(M,compute_uv= False); #here we use a package of python, to compute SVD of given matrix M           \n",
    "            rank=0;    \n",
    "            if (len(Sn)==1):\n",
    "                count=1;\n",
    "                rank=1;\n",
    "            if (len(Sn)>1):\n",
    "                delta=self.normInfinity(M);\n",
    "                #print \"delta is :\",delta;\n",
    "                count=0;\n",
    "                for i in range(len(Sn)):\n",
    "                    if Sn[len(Sn)-1-i]<delta:\n",
    "                        count=count+1;                 \n",
    "                    else :\n",
    "                        break; \n",
    "                rank=len(Sn)-count;            \n",
    "            if (writeDetails==True):\n",
    "                #print \"length of singular values:\",len(Sn);\n",
    "                print \"All singular values are:\",Sn;  \n",
    "                print \"Size of the matrix is :\",M.shape[0],M.shape[1];  \n",
    "                return rank, M.shape[0]-(rank)\n",
    "            elif (writeDetails==False):\n",
    "                return rank, M.shape[0]-(rank)\n",
    "            else: \n",
    "                print \"You should either choose True or False!\" \n",
    "                return \"You should either choose True or False!\"        \n",
    "        else:\n",
    "            #print \"Exact computation!\";\n",
    "            rank=M.rank();\n",
    "            corank=M.nrows()-rank;\n",
    "            if writeDetails==True:\n",
    "                un,sn,vn=linalg.svd(M);\n",
    "                print \"All singular values are:\",sn;\n",
    "                print \"size of the matrix is :\",M.nrows(),M.ncols();\n",
    "                return rank,corank;\n",
    "            elif writeDetails==False:\n",
    "                return rank,corank;\n",
    "            else:\n",
    "                print \"You should either choose True or False!\" \n",
    "                return \"You should either choose True or False!\"  \n",
    "            \n",
    "    # theoretical bound with respect to the article of Sturmfels: euclidean distance degree        \n",
    "    def theoreticalBound(self):\n",
    "        \"\"\"\n",
    "        INPUT: no input.\n",
    "        OUTPUT : integer - number of the expected orthogonal projections, with respect to the article, \n",
    "                Euclidean distance degree.\n",
    "        \"\"\"  \n",
    "        typeOfParametrization=self._string;\n",
    "        bound=0;\n",
    "        uvDegree=self._uvDegree;\n",
    "        if typeOfParametrization==\"rational\":\n",
    "            bound=7*uvDegree^2-9*uvDegree+3;\n",
    "        else:\n",
    "            bound=(2*uvDegree-1)^2;\n",
    "\n",
    "    def rankCorankViaNorm2(self,MM,writeDetails=False): #output : rank, corank\n",
    "        \"\"\"\n",
    "        INPUT: 1)M-matrix,\n",
    "               2)(optional) writeDetails-bool : writeDetails=False by default.\n",
    "               \n",
    "        OUTPUT : (integer,integer) = rank, corank of the input matrix.\n",
    "                If writeDetails==True, then if writes the singular values of the given matrix, the size of the\n",
    "                given matrix.\n",
    "                This method uses SVD for numeric computations, computations over RR,\n",
    "                .rank, and .right_kernel() commands of SageMath for exact computations, i.e, the computations over QQ.\n",
    "                It uses norm2 * 10^(-8) as tolerence.\n",
    "        \"\"\"      \n",
    "        if self._field==RR:\n",
    "            M=numpy.matrix(MM);\n",
    "            Sn= linalg.svd(M,compute_uv= False); #here we use a package of python, to compute SVD of given matrix M           \n",
    "            rank=0;    \n",
    "            if (len(Sn)==1):\n",
    "                count=1;\n",
    "                rank=1;\n",
    "            if (len(Sn)>1):\n",
    "                delta=self.norm2(M);\n",
    "                #print \"delta is :\",delta;\n",
    "                count=0;\n",
    "                for i in range(len(Sn)):\n",
    "                    if Sn[len(Sn)-1-i]<delta:\n",
    "                        count=count+1;                 \n",
    "                    else :\n",
    "                        break; \n",
    "                rank=len(Sn)-count;            \n",
    "            if (writeDetails==True):\n",
    "                #print \"length of singular values:\",len(Sn);\n",
    "                print \"All singular values are:\",Sn;  \n",
    "                print \"Size of the matrix is :\",M.shape[0],M.shape[1];  \n",
    "                return rank, M.shape[0]-(rank)\n",
    "            elif (writeDetails==False):\n",
    "                return rank, M.shape[0]-(rank)\n",
    "            else: \n",
    "                print \"You should either choose True or False!\"\n",
    "                return \"You should either choose True or False!\"        \n",
    "\n",
    "\n",
    "        else:\n",
    "            #print \"Exact computation!\";\n",
    "            rank=M.rank();\n",
    "            corank=M.nrows()-rank;\n",
    "            if writeDetails==True:\n",
    "                un,sn,vn=linalg.svd(M);\n",
    "                print \"All singular values are:\",sn;\n",
    "                print \"size of the matrix is :\",M.nrows(),M.ncols();\n",
    "                return rank,corank;\n",
    "            elif writeDetails==False:\n",
    "                return rank,corank;\n",
    "            else:\n",
    "                print \"You should either choose True or False!\"\n",
    "                return \"You should either choose True or False!\"  \n",
    "\n",
    "\n",
    "\n",
    "    def numericalKernel(self,M,r,writeDetails=False):\n",
    "        \"\"\"\n",
    "        INPUT: 1)M-matrix.\n",
    "               2)rank of the given matrix -integer.\n",
    "               3)(optional) writeDetails-string\n",
    "               \n",
    "        OUTPUT : Kernel of the input matrix- numpy.matrix, if it is numerical computation, i.e. over RR. \n",
    "                 Then it is computed via SVD.\n",
    "                 \n",
    "                 Kernel of the input matrix- matrix (SageMath), if it is exact computation. Then,it is \n",
    "                 computed via, .right_kernel(), command of SageMath.\n",
    "        \"\"\"         \n",
    "        if self._field==RR:\n",
    "            M=numpy.matrix(M);               \n",
    "            Un,Sn,Vn = linalg.svd(M, full_matrices=True,compute_uv= True)\n",
    "            V=numpy.matrix(Vn).transpose();\n",
    "        \n",
    "            #print \"min(Row,Col) = \",min(M.shape[0],M.shape[1]);\n",
    "            \n",
    "\n",
    "            kernel=V[:, [i for i in range(r,M.shape[1])]]; \n",
    "                #kernel=U[[i for i in range(r,M.shape[0])],:];\n",
    "            if (writeDetails==True):\n",
    "                print \"verification of the right null space :\";\n",
    "                print M*Vn.T;\n",
    "                return kernel\n",
    "            elif (writeDetails==False):\n",
    "                return kernel;\n",
    "            else:\n",
    "                return \"You should either choose True or False\";   \n",
    "        else:\n",
    "            if writeDetails!=True and writeDetails!=False:\n",
    "                return \"You should either choose True or False\";\n",
    "                \n",
    "            kernel=M.right_kernel(basis='pivot').basis_matrix().transpose();\n",
    "            print \"kernel is size:\",kernel.nrows(),kernel.ncols();\n",
    "            return kernel;   \n",
    "            \n",
    "    def normInfinity(self,M):\n",
    "        \"\"\"\n",
    "        INPUT: 1)M-matrix.\n",
    "               \n",
    "        OUTPUT : norm infinity of the given matrix - float\n",
    "        \"\"\"  \n",
    "        \n",
    "        M=numpy.matrix(M);\n",
    "        sumL=[];\n",
    "        #print M.shape[0],M.shape[1];\n",
    "        for i in range (M.shape[0]):            \n",
    "            sumL.append(sum(numpy.absolute(M[i,j]) for j in range(M.shape[1])));   \n",
    "        return max(sumL[i] for i in range(M.shape[0]))*10^(-5).n()\n",
    "         \n",
    "    def norm2(self,M):\n",
    "        \"\"\"\n",
    "        INPUT: 1)M-matrix.\n",
    "               \n",
    "        OUTPUT : norm2 of the given matrix - float\n",
    "        \"\"\"          \n",
    "        \n",
    "        M=numpy.matrix(M);\n",
    "        Sn= linalg.svd(M,compute_uv= False);\n",
    "        return (Sn[0]*10^(-8));\n",
    "\n",
    "    def foncMrepExactEval(self,m,writeDetails=False,p=[]):\n",
    "        \"\"\"\n",
    "        INPUT: 1)m-matrix or numpy.matrix.\n",
    "               2)writeDetails- bool, writeDetails=False, by default.\n",
    "               3)p=(float,float,float), point in R^3, p=[] by default.\n",
    "               \n",
    "        OUTPUT : list[numpy.matrix,integer,p), matrix representation of the congurence of the normal lines to the given surface,\n",
    "                 at point p, corank of the matrix representation at point p and the coordinates of the point p.\n",
    "                 This method is called only for EXACT computations by the other methods in this module.\n",
    "                 If writeDetails==True, then it writes it is exact computation, size, rank and corank of the \n",
    "                 input matrix.\n",
    "        \"\"\"     \n",
    "        \n",
    "        if p==[]:\n",
    "            p=self._list[4];\n",
    "        M_rep=((m[0])+ ((p[0]))*(m[1])+((p[1]))*(m[2])+((p[2]))*(m[3]));\n",
    "        RankMrepatP=M_rep.rank();\n",
    "        coRankMrepatP=M_rep.nrows()-RankMrepatP;\n",
    "        \n",
    "        if (writeDetails==True): \n",
    "            print \"Exact computation !\"\n",
    "            print \"Size of M-rep is:\",M_rep.nrows(),M_rep.ncols();\n",
    "            print \"Rank/CoRank of Mrep is at a randomly choosen point is :\",RankMrepatP,coRankMrepatP; \n",
    "            return M_rep,coRankMrepatP,p;\n",
    "        elif (writeDetails==False): \n",
    "            return M_rep,coRankMrepatP,p;\n",
    "        else:\n",
    "            print \"You can choose either 'True' or 'False' \";  \n",
    "            return \"You can choose either 'True' or 'False' \";  \n",
    "             \n",
    "        \n",
    "    def foncMrepEval(self,m,writeDetails=False,p=[]):\n",
    "        \"\"\"\n",
    "        INPUT: 1)m-matrix or numpy.matrix.\n",
    "               2)writeDetails- bool, writeDetails=False, by default.\n",
    "               3)p=(float,float,float), point in R^3, p=[] by default.\n",
    "               \n",
    "        OUTPUT : list[numpy.matrix,integer,p), matrix representation of the congurence of the normal lines to the given surface,\n",
    "                 at point p, corank of the matrix representation at point p and the coordinates of the point p.\n",
    "                 This method is called only for NUMERIC computations by the other methods in this module.\n",
    "                 If writeDetails==True, then it writes it is exact computation, size, rank and corank of the \n",
    "                 input matrix.\n",
    "        \"\"\" \n",
    "        \n",
    "        if p==[]:\n",
    "            p=self._list[4];\n",
    "        \n",
    "        #T=True;\n",
    "        #if T==True:\n",
    "        if(self._field.cardinality()<infinity) or self._field==QQ:\n",
    "            #print \"calcul exact\";\n",
    "            return self.foncMrepExactEval(m,string,p);\n",
    "        #if T==False:\n",
    "        else: \n",
    "        \n",
    "            M_rep=((m[0])+ (QQ(p[0]))*(m[1])+(QQ(p[1]))*(m[2])+(QQ(p[2]))*(m[3]));\n",
    "            (RankMrepatP,coRankMrepatP)=self.rankCorankViaNorm2(M_rep,False);\n",
    "        \n",
    "        if (writeDetails==True): \n",
    "            print \"Size of M-rep is:\",M_rep.shape[0],M_rep.shape[1];\n",
    "            print \"Rank/CoRank of Mrep is at a randomly choosen point is :\",RankMrepatP,coRankMrepatP; \n",
    "            return M_rep,M_rep.shape[0]-RankMrepatP,p;\n",
    "        elif (writeDetails==False): \n",
    "            return M_rep,M_rep.shape[0]-RankMrepatP,p;\n",
    "        else:\n",
    "            print \"You can choose either 'True' or 'False' \";\n",
    "            return \"You can choose either 'True' or 'False' \";\n",
    "            \n",
    "            \n",
    "\n",
    "    def partialDerivatives(self,writeDetails=False):   \n",
    "        \"\"\"\n",
    "        INPUT : 1)writeDetails=bool, writeDetails=False by default.\n",
    "                                    \n",
    "        OUTPUT : 1)delta_u - vector : partial derivative with respect to u,\n",
    "                 2)delta_v - vector : partial derivative with respect to v.\n",
    "                 If writeDetails==True, then the function prints the partial derivatives\n",
    "                 of the parameterization of given triangular surface and also the cross product\n",
    "        \"\"\"\n",
    "        R.<u,v>= self._field['u','v'];\n",
    "        R.gens();\n",
    "        string1=self._string;\n",
    "        K=R.fraction_field();\n",
    "        \n",
    "        if (string1==\"rational\"):\n",
    "            \n",
    "            Phi=vector([self._list[1]/self._list[0], self._list[2]/self._list[0], self._list[3]/self._list[0]]);            \n",
    "            \n",
    "        else:\n",
    "            \n",
    "            Phi=vector([self._list[1], self._list[2], self._list[3]]);            \n",
    "            \n",
    "            \n",
    "        delta_u = vector([simplify(Phi[0].derivative(u)), simplify(Phi[1].derivative(u)), simplify(Phi[2].derivative(u))]);\n",
    "        delta_v = vector([simplify(Phi[0].derivative(v)), simplify(Phi[1].derivative(v)), simplify(Phi[2].derivative(v))]);\n",
    "                                 \n",
    "                \n",
    "        \n",
    "        if (writeDetails==True):\n",
    "            print \"Partial derivatives of the parametrization wrt u  are :\",delta_u;\n",
    "            print \"Partial derivatives of the parametrization wrt v  are :\",delta_v;            \n",
    "            print \"Cross product of partial derivatives wrt u and v is:\",Delta;\n",
    "            return delta_u,delta_v\n",
    "            \n",
    "        elif (writeDetails==False):              \n",
    "            return delta_u,delta_v\n",
    "                  \n",
    "        else:\n",
    "            return \"You can choose either 'True' or 'False' \";            \n",
    "            \n",
    "             \n",
    "                                              \n",
    "    def homogeneousPartialDerivatives(self,Fi=[]):\n",
    "        \"\"\"\n",
    "        INPUT : 1)Fi-list, list Polynomials in variables u,v,w over the chosen field, Fi=[] by default.\n",
    "                                    \n",
    "        OUTPUT : 1)delta_u - vector : partial derivative with respect to u,\n",
    "                 2)delta_v - vector : partial derivative with respect to v,\n",
    "                 2)delta_w - vector : partial derivative with respect to w.\n",
    "                 It only computes the partial derivatives of : (Fi[1]/Fi[0],Fi[2]/Fi[0],Fi[3]/Fi[0]).\n",
    "                 It does not compute the partial derivatives for the rest of the elements of the fiven list Fi.\n",
    "        \"\"\"\n",
    "        \n",
    "        \n",
    "        R1.<u,v,w>=self._field[];\n",
    "        R1.gens();\n",
    "        R1.gen(0);\n",
    "        s=R1.gen(0); \n",
    "        #K=R1.fraction_field();\n",
    "        if Fi==[]:\n",
    "            Fi=self.HomogeneousParametrization();\n",
    "        RFi=[Fi[1]/Fi[0],Fi[2]/Fi[0],Fi[3]/Fi[0]];\n",
    " \n",
    "\n",
    "        delta_u = vector([(RFi[0].derivative(u)).numerator(), \n",
    "                           (RFi[1].derivative(u)).numerator(), \n",
    "                            (RFi[2].derivative(u)).numerator()]);\n",
    "                           \n",
    "        delta_v = vector([(RFi[0].derivative(v)).numerator(),\n",
    "                          (RFi[1].derivative(v)).numerator(),\n",
    "                          (RFi[2].derivative(v)).numerator()]); \n",
    "                            \n",
    "        delta_w = vector([(RFi[0].derivative(w)).numerator(), \n",
    "                          (RFi[1].derivative(w)).numerator(), \n",
    "                          (RFi[2].derivative(w)).numerator()]);    \n",
    "                            \n",
    "        return delta_u,delta_v,delta_w\n",
    "\n",
    "\n",
    "    def verificationHomogeneousNormalVectorViaJacobian(self,Fi): \n",
    "        ##We verify only with all partial derivatives wrt s,t,u,v but only one normal vector (without considering dv)\n",
    "        ##(d_ F1/F0,d_ F2/F0,d_ F3/F0)*(N1,N1,N2) == 0 ?        \n",
    "        \"\"\"\n",
    "        INPUT : 1)Fi-list, list Polynomials in variables u,v,w.\n",
    "                                    \n",
    "        OUTPUT : bool, if the computed homogeneous normal vectors are perpendicular to the surface, simply does\n",
    "                 dot product computation.\n",
    "        \"\"\"        \n",
    "        HomogeneousNormalVector = self.NormalVectorsViaJacobian(Fi);\n",
    "        HomogeneousDer = self.homogeneousPartialDerivatives(Fi);\n",
    "        \n",
    "\n",
    "\n",
    "        v1=vector((HomogeneousDer[0][0],HomogeneousDer[0][1],HomogeneousDer[0][2]))*vector((HomogeneousNormalVector[0],HomogeneousNormalVector[1],HomogeneousNormalVector[2]));\n",
    "        v2=vector((HomogeneousDer[1][0],HomogeneousDer[1][1],HomogeneousDer[1][2]))*vector((HomogeneousNormalVector[0],HomogeneousNormalVector[1],HomogeneousNormalVector[2]));        \n",
    "        v3=vector((HomogeneousDer[2][0],HomogeneousDer[2][1],HomogeneousDer[2][2]))*vector((HomogeneousNormalVector[0],HomogeneousNormalVector[1],HomogeneousNormalVector[2]));                   \n",
    "        \n",
    "    \n",
    "        if (v1==0 and v2==0 and v3==0):  #(v1<10^(-3)  and  v2<10^(-3) and v3<10^(-3)):\n",
    "            return true;\n",
    "        vList=[v1,v2,v3];\n",
    "        totalLength=0;\n",
    "        counter=0;\n",
    "        for i in range(3):\n",
    "            coef=vList[i].coefficients();\n",
    "            length=len(coef);\n",
    "            totalLength=totalLength+length;\n",
    "            for j in range(length):\n",
    "                if coef<10^(-3) or coef<e:\n",
    "                    counter=counter+1;\n",
    "            if counter==totalLength:\n",
    "                return true;\n",
    "        else:\n",
    "            print v1,v2,v3\n",
    "            return false;  \n",
    "\n",
    "            \n",
    "            \n",
    "            \n",
    "        \n",
    "    def leftKernel(self,M,r,writeDetails=False):\n",
    "        \"\"\"\n",
    "        INPUT : 1)M-matrix or numpy.matrix\n",
    "                2)r - integer:rank of the input matrix,\n",
    "                3)writeDetails-bool, writeDetails=False by default.\n",
    "                                    \n",
    "        OUTPUT : numpy.matrix, numerical left kernel of the input matrix by SVD computation.\n",
    "        \"\"\"         \n",
    "        M=numpy.matrix(M);               \n",
    "        Un,Sn,Vn = linalg.svd(M, full_matrices=True,compute_uv= True) \n",
    "        U=numpy.matrix(Un).transpose();\n",
    "        \n",
    "        kernel=U[[i for i in range(r,M.shape[0])],:];\n",
    "        if (writeDetails==True):\n",
    "            print \"verification of the left null space :\";\n",
    "            print Un*M;\n",
    "            return kernel\n",
    "        elif (writeDetails==False):\n",
    "            return kernel;\n",
    "        else:\n",
    "            print \"You should either choose True or False\";\n",
    "            return \"You should either choose True or False\";        \n",
    "\n",
    "               \n",
    "        if ((writeDetails==True) and len(Sn)-count>1):\n",
    "            nR=V.shape[1]; \n",
    "            nC=V.shape[1]-count;            \n",
    "            print \"all singular values are:\",Sn;\n",
    "            #print \"type of the singular values are :\",type(Sn);\n",
    "            #print \"type of one element among all of sing values is :\",type(Sn[0]);\n",
    "            print \"lenght of Singular values list:\",len(Sn);\n",
    "            #print \"Size of the given matrix:\",len(Sn),V.shape[1];\n",
    "            print \"number of the last columns of V defining kernel of S_Nu are:\",V.shape[1]-count;    \n",
    "            #print \"type of V is:\",type(Vn);\n",
    "            #print \"Approximation of the kernel is :\",V;\n",
    "            return kernel,nR,nC,count \n",
    "        elif ((writeDetails==False) and len(Sn)-count>1):\n",
    "            nR=V.shape[1]; \n",
    "            nC=V.shape[1]-count;        \n",
    "            return kernel,nR,nC,count \n",
    "        elif (len(Sn)-count==1):\n",
    "            print \"The given matrix is a column matrix, not enough to compute the kernel\";\n",
    "            print \"ATTENTION\"; \n",
    "            return V[:,len(Sn)-1],V.shape[1],1,1 \n",
    "        elif (len(Sn)-count==0):\n",
    "            print \"It is full rank!!!\";\n",
    "            return 0,V.shape[1],0,min(V.shape[1],len(Sn))                 \n",
    "        else:\n",
    "            print \"You should either choose True or False\";\n",
    "            return \"You should either choose True or False\";    \n",
    "            \n",
    "            \n",
    "    def foncIsFullRank(self,M):\n",
    "        \"\"\"\n",
    "        INPUT : 1)M-matrix or numpy.matrix\n",
    "                                    \n",
    "        OUTPUT : bool, if the input matrix is full rank.\n",
    "        \"\"\"   \n",
    "        \n",
    "        M=numpy.matrix(M);\n",
    "        Sn= linalg.svd(M,compute_uv=False );\n",
    "        delta=Sn[0]*10^(-8);\n",
    "        #delta=self.normInfinity(M);\n",
    "        #print \"delta is :\",delta;\n",
    "        if Sn[len(Sn)-1]<delta :\n",
    "            return false;\n",
    "        else:\n",
    "            return true;                \n",
    "            \n",
    "\n",
    "\n",
    "    def IsMrepBigEnoughToInvert(self,M,nu,a=0):\n",
    "        \"\"\"\n",
    "        INPUT : 1)M-matrix or numpy.matrix\n",
    "                2)nu-integer, the regularity\n",
    "                3)a-integer, a=0 by default. a is the number how many times this method has been used before.\n",
    "                                    \n",
    "        OUTPUT : bool, if the input matrix is big enough to compute the parameter values corresponding to the\n",
    "                 orthogonal projection on the the given triangular surface.\n",
    "        \"\"\"\n",
    "        \n",
    "        if nu==1:\n",
    "            return False;\n",
    "        \n",
    "        corankM=M[1];\n",
    "        MM=numpy.matrix(M[0]);\n",
    "        if a==0:\n",
    "            if nu==1:\n",
    "                return False;\n",
    "            elif nu>1 and corankM+(nu+1)<=MM.shape[0]:\n",
    "                return True;\n",
    "            else:\n",
    "                return False;\n",
    "        elif a==1:\n",
    "            return False;        \n",
    "        else:\n",
    "            return \"ERROR: a should be either 0 or 1, use 1 for when you get ERROR from inversion command!\";\n",
    "            \n",
    "    def PreInversion(self,M,Deg,a=0,writeDetails=False,method=\"cross_product\"):\n",
    "        \"\"\"\n",
    "        INPUT : 1)M-matrix or numpy.matrix\n",
    "                2)Deg-integer, degree of the homogeneous parameterization of the normal lines to the surface.\n",
    "                3)a-integer, a=0 by default. a is the number how many times this method has been used before.\n",
    "                4)writeDetails-bool, writeDetails=False, by default.\n",
    "                5)method-string : either \"cross_product\" or \"jacobian\", method=\"cross_product\" by default.\n",
    "        OUTPUT : numpy.matrix,integer,integer : new matrix representation, new regularity, new rank.\n",
    "                 If writeDetails==True, then it writes the size and corank of the new matrix, and \n",
    "                 the computation time.\n",
    "        \"\"\"        \n",
    "        if string2!=\"cross_product\" and string2!=\"jacobian\":\n",
    "            return \"You should choose either jacobian or cross_product!\"\n",
    "            \n",
    "        m=M[0];    \n",
    "        nu=self.regularity(Deg);\n",
    "        #print \"nu is:\",nu;\n",
    "        \n",
    "        if self.IsMrepBigEnoughToInvert(M,nu,a)==True:\n",
    "            #print \"Mrep is big enough to invert\";\n",
    "            MM=numpy.matrix((M[0]));\n",
    "            return M,nu,MM.shape[0]-M[1];\n",
    "        else:\n",
    "            print \"Second Mrep is being computed!\";\n",
    "            m=M[0];\n",
    "            coRankOfFirstMrep=M[1];\n",
    "        \n",
    "            p=M[2];\n",
    "        \n",
    "            #print \"Size of first Mrep is :\",m.shape[0],m.shape[1];\n",
    "            #print \"corank of fist Mrep is:\",coRankOfm;\n",
    "            \n",
    "            #print \"PREINVERSION1\";\n",
    "            preInv1= datetime.now();        \n",
    "\n",
    "            nu= nu+1;\n",
    "            \n",
    "            print \"new nu is :\",nu;              \n",
    "            NewSnu=self.SNu(nu,Deg,False,method,[]) ;\n",
    "            \n",
    "            #print \"PREINVERSION2\";            \n",
    "            #Notice here we only make computations for RR so no problem with writing NewSnu.n();\n",
    "            \n",
    "            NewMrepG=self.foncMrepG(NewSnu,False)[0];\n",
    "            \n",
    "            NewMrep=self.foncMrepEval(NewMrepG);\n",
    "            coRankOfNewMrep=NewMrep[1];\n",
    "            m=NewMrep[0];  \n",
    "            \n",
    "            preInv1bis = datetime.now()     \n",
    "            Rankm=m.shape[0]-coRankOfFirstMrep;   \n",
    "\n",
    "            #print \"PREINVERSION3\";  \n",
    "            print \"Coranks of first and second Mrep are:\",coRankOfFirstMrep,coRankOfNewMrep;          \n",
    "            if coRankOfFirstMrep!=coRankOfNewMrep:\n",
    "                print \"Coranks of first and second Mrep are:\",coRankOfFirstMrep,coRankOfNewMrep;\n",
    "                print \"ERROR : corank of old Mrep and corank of new Mrep are different!\";\n",
    "                return \"ERROR : corank of old Mrep and corank of new Mrep are different!\";\n",
    "                \n",
    "            if writeDetails==False:    \n",
    "\n",
    "                return  NewMrep,nu,Rankm;  \n",
    "                \n",
    "            elif writeDetails==True:\n",
    "                print \"Size of New Mrep is :\",m.shape[0],m.shape[1];\n",
    "                print \"Corank of New Mrep is :\",coRankOfNewMrep;\n",
    "                print \"TimeOfComputationOfNewMrepFrom0:\",preInv1bis-preInv1;\n",
    "                return  NewMrep,nu,Rankm; \n",
    "           \n",
    "            else:\n",
    "                print \"You can either choose True or False\"; \n",
    "                return \"You can either choose True or False\";                 \n",
    "                                              \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "            \n",
    "            \n",
    "            \n",
    "    def inversion(self,M,nu,rankM,writeDetails=False): #inversion step1:  #inversion step1: Divide into 2 submatrix of size=max rank; then compute the eigen values, i.e. output is s coordinates at point p. It takes directly Mrep of m where m[0] is its Mrep and m[1] is its corank at point p and m[2] is p.\n",
    "        \"\"\"\n",
    "        INPUT : 1)M-matrix or numpy.matrix,\n",
    "                2)nu-integer, the regularity,\n",
    "                3)rankM-integer,rank of the input matrix,\n",
    "                4)writeDetails-bool, writeDetails=False, by default.\n",
    "        OUTPUT : list,list,integer,list,vector,numpy.matrix,numpy.matrix,numpy.array,numpy.array,timedate :\n",
    "                 real u parameter values, real v parameter values, number of real parameter values,\n",
    "                 complex v parameter values, point that we project onto the surface orthogonaly,\n",
    "                 cokernel of the input matrix, submatrix of the input matrix, submatrix of the input matrix, \n",
    "                 generalized eigen vectors with respect to the two submatrix of the two submatrix of the input matrix,\n",
    "                 generalized eigen values with respect to the two submatrix of the two submatrix of the input matrix\n",
    "                 computation time in timedate() for the inversion.\n",
    "                 If writeDetails==True, it writes the computation times at each step.\n",
    "\n",
    "        \"\"\"        \n",
    "    \n",
    "        invNumber=1;     \n",
    "        timeOfComputation=[];\n",
    "        Coords=[];\n",
    "        EV=[];\n",
    "        Sols=[]; #Has real-valued eigen values\n",
    "        lSols=[];\n",
    "        Nlist=[];\n",
    "        uCoords=[]; \n",
    "        uVecs=[];\n",
    "         \n",
    "        m=M[0];\n",
    "        coRankOfMrep=M[1];\n",
    "        #print \"coRankOfMrep is :\",coRankOfMrep;\n",
    "        p=M[2];\n",
    "    \n",
    "        inv1=datetime.now();\n",
    "        \n",
    "        \n",
    "          \n",
    "        k1=self.leftKernel(m,Integer(rankM),False).transpose();\n",
    "        #print \"k1 is of size:\",k1.shape[0],k1.shape[1];        \n",
    "        inv2 = datetime.now();   \n",
    "        \n",
    "        #bound=floor(sqrt(2*coRankOfMrep));\n",
    "        #bound=(nu+1)*(nu+2)/2;                     \n",
    " \n",
    "        A=k1[[i for i in range(coRankOfMrep)],:];\n",
    "        counter=0;\n",
    "        index=var('index');\n",
    "        indiceList=[];\n",
    "\n",
    "        for i in range(1,nu+1):\n",
    "            for j in range(i):\n",
    "                #print i;\n",
    "                #print \"index no :\",counter+i;\n",
    "                index=counter+i;\n",
    "                indiceList.append(index);\n",
    "                counter=counter+1;\n",
    "                if counter==coRankOfMrep:    \n",
    "                    break; \n",
    "            if len(indiceList)==coRankOfMrep:\n",
    "                 break;       \n",
    "        #print \"indiceList is :\", indiceList;\n",
    "        B=k1[[indiceList[i] for i in range(len(indiceList))],:];\n",
    "        \n",
    "        #print \"A is of size :\", A.shape[0],A.shape[1];\n",
    "        #print \"B is of size :\", B.shape[0],B.shape[1];\n",
    "\n",
    "\n",
    "        \n",
    "        f=False;\n",
    "        #if self.foncIsFullRank(A)==False or self.foncIsFullRank(B)==False:\n",
    "        if f==True:  \n",
    "              \n",
    "            #print \"A and B are choosen LIGNE BY LIGNE!\";\n",
    "      \n",
    "            firstRowVariable = 0;\n",
    "            squareFound = false;    \n",
    "        \n",
    "            counter=1;\n",
    "            IndexList=[];\n",
    "            IndexList.append(firstRowVariable);\n",
    "            ligneVar=k1.shape[0]-Nu_2-1;  # we do not take into account the last block of k1\n",
    "            #print \"ligne Var is :\",ligneVar;\n",
    "            AIsDone=False;\n",
    "       \n",
    "\n",
    "            k2=k1[[i for i in range(ligneVar+1)],:];   \n",
    "            #print \"k2 is of size:\", k2.shape[0],k2.shape[1];    \n",
    "            testMatrix=np.vstack((k2,k2));\n",
    "            #print \"Size of testMatrix is :\",testMatrix.shape[0],testMatrix.shape[1];\n",
    "               \n",
    "            while(squareFound == false):\n",
    "                A=testMatrix[firstRowVariable,:];\n",
    "                for j in range(firstRowVariable+1,firstRowVariable+ligneVar):\n",
    "                    if self.foncIsFullRank(np.vstack((A,testMatrix[j,:])))==true:\n",
    "                        if j<ligneVar:                       \n",
    "                            #print \"j<=ligneVar\",j,len(IndexList),IndexList;\n",
    "                            IndexList.append(j);\n",
    "                        else:\n",
    "                            IndexList.append(j-ligneVar); \n",
    "                            #print j-ligneVar-1,len(IndexList);\n",
    "                        A=np.vstack((A,testMatrix[j,:]));\n",
    "                        #print A;\n",
    "                        if (A.shape[0] == A.shape[1]):\n",
    "                            squareFound = true;\n",
    "                            break;             \n",
    "                        counter=counter+1;\n",
    "                #print \"arada A is of size\",j,A.shape[0],A.shape[1];  \n",
    "                if(A.shape[0]!=A.shape[1]):  \n",
    "                    #print \"A is of size\",A.shape[0],A.shape[1];                          \n",
    "                    firstRowVariable = firstRowVariable + 1;\n",
    "                    #print \"firstRowVariable\";\n",
    "                    IndexList=[];\n",
    "                    IndexList.append(firstRowVariable);\n",
    "                    counter=1;\n",
    "                    #print \"firstRowVariable\",firstRowVariable;\n",
    "           \n",
    "            #print \"Last firstRowVariable is :\",firstRowVariable;\n",
    "             \n",
    "             \n",
    "            if (squareFound == false):\n",
    "                print \"ERROR: Square matrix could not be found in input matrix!\";\n",
    "                return \"ERROR: Square matrix could not be found in input matrix!\";\n",
    " \n",
    "                        \n",
    "            #print \"Rank of new A is :\",self.rankCorankViaNorm2(A)[0];\n",
    "            #print \"The size of A is : \",A.shape[0],A.shape[1];  \n",
    "         \n",
    "            B=k1[firstRowVariable+(Nu_2+1),:];\n",
    "\n",
    "            #print \"new Nu_2 is :\",Nu_2;\n",
    "            #print \"lenght of IndexList is :\",len(IndexList),IndexList;\n",
    "                \n",
    "            for j in range(1,len(IndexList)):\n",
    "                B=np.vstack((B,k1[IndexList[j]+(Nu_2+1),:]));   \n",
    "\n",
    "            #print \"Rank of new B is:\",self.rankCorankViaNorm2(B)[0];\n",
    "            #print \"The size of B is : \",B.shape[0],B.shape[1];\n",
    "\n",
    "            if self.foncIsFullRank(B)==False:\n",
    "                print \"ERROR: for this Mrep there is no place to find A,B having full rank!\";\n",
    "                invNumber=invNumber+1;\n",
    "                print \"New inversion is being called, with bigger regularity - inversion :\",invNumber;\n",
    "                (mrep,n2,rankMrep)=self.PreInversion(M,sD,uD,1,False,\"jacobian\");\n",
    "                return self.inversion(mrep,n1,n2,rankMrep,False);\n",
    "                #return \"ERROR: for this Mrep there is no place to find A,B having full rank!\";\n",
    "    \n",
    "                      \n",
    "\n",
    "        inv3 = datetime.now();            \n",
    "\n",
    "        #print \"Determinant of the choosen two matrices are :\",linalg.det(A-2*B);\n",
    "        \n",
    "        eigvals, eigvecs = eig(A,B); #it computes eigenvalues and right eigen vectors\n",
    "        #The normalized (unit “length”) eigenvectors, such that the column v[:,i] is the eigenvector corresponding to the eigenvalue w[i].\n",
    "        \n",
    "        inv4 = datetime.now() \n",
    "                \n",
    "        sayac2=[]; #to count which eigen-values are real valued.\n",
    "        EV=[];\n",
    "        \n",
    "        #To verify the eigenvalues and eigenvectors\n",
    "        #for i in range(len(eigvals)):\n",
    "            #print linalg.det(A-eigvals[i]*B);\n",
    "        #Sols=eigvals;                             \n",
    "\n",
    "        \n",
    "        for k in range(len(eigvals)):\n",
    "            if imag(eigvals[k])==0:\n",
    "                Sols.append((eigvals[k]).real);\n",
    "                #print \"is it really an eigen value, i.e. Is it zero? :\", linalg.det(A-eigvals[k]*B)\n",
    "                sayac2.append(k);  \n",
    "                \n",
    "        #EV has all real valued eigenvectors                      \n",
    "        for w in range(len(sayac2)):\n",
    "            EV.append((eigvecs)[:,sayac2[w]]);\n",
    "       \n",
    "        #(B*coord[7])[:,7]  \n",
    "        for w in range(len(sayac2)):\n",
    "            if coRankOfMrep==1:\n",
    "                if imag(EV[w][0])==0:\n",
    "                    Coords.append(((B*eigvecs)[:,sayac2[w]]).item((0, 0)));\n",
    "            else:\n",
    "                Coords.append(((B*eigvecs)[:,sayac2[w]]));\n",
    "                \n",
    "        vList=[];\n",
    "        #print len(len(sayac2));\n",
    "        if (len(sayac2)==1):\n",
    "            #vList.append(((Coords[0].item((0,0))).real));\n",
    "            vList.append(((Coords[0]).real));        \n",
    "        if len(sayac2)>1:\n",
    "            for i in range(len(sayac2)):\n",
    "                vList.append(((Coords[i].real[1]/Coords[i].real[0])).item((0,0)));           \n",
    "            \n",
    "            \n",
    "        #lets compute also the complex ucoordinates, hence all\n",
    "        for i in range(len(eigvals)):\n",
    "            uVecs.append((B*eigvecs)[:,i]);\n",
    "            \n",
    "        for i in range(len(eigvals)):\n",
    "            if len(uCoords)>1:\n",
    "                uCoords.append((uVecs[i][0]/uVecs[i][1]).item((0,0)));      \n",
    "            else:\n",
    "                uCoords.append((uVecs[i][0]).item((0,0)));                              \n",
    "                #uCoords.append((uVecs[i][0]).item((0,0)));  \n",
    "        \n",
    "        inv5 = datetime.now();\n",
    "        \n",
    "        timeOfComputationOfCoKernel = inv2-inv1;\n",
    "        timeOfComputation.append(timeOfComputationOfCoKernel);\n",
    "\n",
    "        timeOfComputationOfAB = inv3 - inv2;\n",
    "        timeOfComputation.append(timeOfComputationOfAB);\n",
    "        \n",
    "        timeOfComputationOfEigen = inv4- inv3;\n",
    "        timeOfComputation.append(timeOfComputationOfEigen);\n",
    "                \n",
    "        timeOfComputationOfReelCoord= inv5-inv4;\n",
    "        timeOfComputation.append(timeOfComputationOfReelCoord);                     \n",
    "        if (writeDetails==True):\n",
    "            #print \"Chosen first matrix is: \",A;\n",
    "            #print \"Chosen second matrix is:\",B; \n",
    "            print \"timeOfComputationOfCoKernel is :\",timeOfComputationOfCoKernel;\n",
    "            print \"timeOfComputationOfAB is :\",timeOfComputationOfAB;\n",
    "            print \"timeOfComputationOfEigen is :\",timeOfComputationOfEigen;\n",
    "            print \"timeOfComputationOfReelCoord is :\",timeOfComputationOfReelCoord;\n",
    "            print \"timeOfTotalComputation is :\",inv5-inv1;\n",
    "            #return vector(eigvals),vector(Coords),len(Coords),p,k1,A,B,eigvecs;\n",
    "            return (Sols),(vList),len(Coords),uCoords,p,k1,A,B,eigvecs,eigvals,inv5-inv1;\n",
    "        elif (writeDetails==False):\n",
    "            #return vector(eigvals),vector(Coords),len(Coords),p,k1,A,B,eigvecs;\n",
    "            return (Sols),(vList),len(Coords),uCoords,p,k1,A,B,eigvecs,eigvals,inv5-inv1;\n",
    "        else:\n",
    "            print \"You can choose either 'True' or 'False' !\"\n",
    "            return \"You can choose either 'True' or 'False' !\"\n",
    "            \n",
    "            \n",
    "\n",
    "    def Plot(self,coord1,coord2,p=[]):\n",
    "        \"\"\"\n",
    "        INPUT : 1)coord1- list, u-coordinates of the orthogonal projections of the point p onto the surface,\n",
    "                2)coord2- list, v-coordinates of the orthogonal projections of the point p onto the surface,\n",
    "                3)p=(float,float,float), point in R^3, p=[] by default. It takes the \n",
    "                point coming from the class properties.\n",
    "\n",
    "        OUTPUT : green point is p,\n",
    "                 yellow points are orthogonal projections of p onto the surface,\n",
    "                 red lines are the lines combining p and its orthogonal projections.\n",
    "                 \n",
    "\n",
    "        \"\"\"   \n",
    "        if p==[]:\n",
    "            p=self._list[4];#the point in the space that we project onto the curve  \n",
    "\n",
    "        S=parametric_plot3d([self._list[1]/self._list[0],self._list[2]/self._list[0],self._list[3]/self._list[0]],(u,-3,3),(v,-3,3));    \n",
    "\n",
    "        length=len(coord1);  \n",
    "        if length>1:\n",
    "            coordx=[];coordy=[];coordz=[];\n",
    "            pointList=[];\n",
    "            for i in range(length):\n",
    "                coordx.append(self._list[1].subs(u=coord1[i],v=coord2[i])/self._list[0].subs(u=coord1[i],v=coord2[i]));\n",
    "                coordy.append(self._list[2].subs(u=coord1[i],v=coord2[i])/self._list[0].subs(u=coord1[i],v=coord2[i]));\n",
    "                coordz.append(self._list[3].subs(u=coord1[i],v=coord2[i])/self._list[0].subs(u=coord1[i],v=coord2[i]));\n",
    "                pointList.append([coordx[i],coordy[i],coordz[i]])\n",
    "                #print pointList;\n",
    "                S=S+line([p,pointList[i]],thickness=5,linestyle=\"solid\",color=\"red\");\n",
    "            S=S+point(pointList,size=(20),color=\"yellow\")+point(p,thickness=5,linestyle=\"solid\",color=\"green\");\n",
    "            return S.show();  \n",
    "        else:\n",
    "                coordx=0;coordy=0;coordz=0;\n",
    "                coordx=self._list[1].subs(u=coord1[0],v=coord2[0])/self._list[0].subs(u=coord1[0],v=coord2[0]);\n",
    "                coordy=self._list[2].subs(u=coord1[0],v=coord2[0])/self._list[0].subs(u=coord1[0],v=coord2[0]);\n",
    "                coordz=self._list[3].subs(u=coord1[0],v=coord2[0])/self._list[0].subs(u=coord1[0],v=coord2[0]);\n",
    "                S=S+point([coordx,coordy,coordz])+point(p,thickness=5,linestyle=\"solid\",color=\"green\")+line([(coordx,coordy,coordz),p],thickness=5,linestyle=\"solid\",color=\"red\");\n",
    "                return S.show();\n",
    "     \n",
    "            \n",
    "            \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 8.6",
   "language": "",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}