poly + POMAIN BRAINY POLRT.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:9665b14192a155e7b8f34879f4aad866060189672b6103818cebb5cd58f74f55"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "import sys\n",
      "# sys.path.append(\"/Users/becqg/Documents/code/Robinson\")\n",
      "sys.path.append(\"C:\\Users\\Guillaume\\Documents\\projets\\code\\Robinson\")\n",
      "import ER\n",
      "import imp"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2> POLYEV</h2>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imp.reload(ER)\n",
      "LA = 3\n",
      "A = numpy.array([3., 2., 1.])\n",
      "Z = complex(1, 1.0)\n",
      "AZ = numpy.empty((LA))\n",
      "AZ = ER.POLYEV(3, A, Z, AZ) \n",
      "print(AZ)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2>POMAIN</h2>\n",
      "p.162"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imp.reload(ER)\n",
      "N = 2\n",
      "LA = 3\n",
      "A = numpy.array([1, 0, 0, 1, -3, 14, 21, 1, 2, 5, -6, 14])\n",
      "LADJ = (N - 1) * (LA - 1) + 1\n",
      "LDET = N * (LA - 1) + 1\n",
      "ADJ = empty(N * N * LADJ)\n",
      "P = empty(N * LDET)\n",
      "DET = empty(LDET)\n",
      "S = empty(N * N * LDET)\n",
      "(ADJ, P, DET, S) = ER.POMAIN(N, LA, A, ADJ, P, DET, S)\n",
      "print(\"ADJ\", ADJ)\n",
      "print(\"P\", P)\n",
      "print(\"DET\", DET)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "('ADJ', array([  1.-0.j,  -0.+0.j,  -0.+0.j,   1.-0.j,   1.-0.j, -14.+0.j,\n",
        "       -21.+0.j,  -3.+0.j,  14.-0.j,  -5.+0.j,   6.-0.j,   2.-0.j]))\n",
        "('P', array([   2.+0.j,   -1.+0.j,   -2.+0.j,    2.+0.j,   16.+0.j,  281.+0.j,\n",
        "         61.+0.j,   61.+0.j,  -58.+0.j,  -58.+0.j]))\n",
        "('DET', array([   1.-0.j,   -2.+0.j, -281.+0.j,  -61.+0.j,   58.-0.j]))\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2> BRAINY </h2>\n",
      "p. 154"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "NRA = 2\n",
      "NCA = 2\n",
      "LA = 3\n",
      "A = numpy.array([0, 1, 0, -1, 2, 0, 1, 1, 1, 1, 0, 0], \"float\")\n",
      "NRB = 2\n",
      "NCB = 2\n",
      "LB = 3\n",
      "B = numpy.array([0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0], \"float\")\n",
      "C = ()\n",
      "C = ER.BRAINY(NRA, NCA, LA, A, NRB, NCB, LB, B, C)\n",
      "print(C)\n",
      "C = ER.BRAINY(NRB, NCB, LB, B, NRA, NCA, LA, A, C)\n",
      "print(C)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0. -1.  0.  0.  1.  0.  0.  0.  1.  2.  3.  2.  2.  0.  3.  1.  1.  1.\n",
        "  1.  1.]\n",
        "[ 0.  0.  0.  0.  1.  3. -1.  0.  1.  3.  0.  2.  3.  2.  2.  0.  2.  0.\n",
        "  0.  0.]\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2> POLRT </h2>\n",
      "p. 35"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imp.reload(ER)\n",
      "XCOF = numpy.array([36., 0, -13, 0, 1])\n",
      "M = 4\n",
      "COF = numpy.empty(M + 1)\n",
      "(ROOTR, ROOTI, IER) = ([], [], [])\n",
      "(COF, ROOTR, ROOTI, IER) = ER.POLRT(XCOF, COF, M, ROOTR, ROOTI, IER)\n",
      "print(COF, ROOTR, ROOTI, IER)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(array([  1.,   0.,  -4.,  12.,  36.]), array([-3.,  3., -2.,  2.]), array([ 0.,  0.,  0.,  0.]), 0)\n"
       ]
      }
     ],
     "prompt_number": 58
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imp.reload(ER)\n",
      "XCOF = numpy.array([83, -278., 720, -238, 120])\n",
      "M = 4\n",
      "COF = numpy.empty(M + 1)\n",
      "(ROOTR, ROOTI, IER) = ([], [], [])\n",
      "(COF, ROOTR, ROOTI, IER) = ER.POLRT(XCOF, COF, M, ROOTR, ROOTI, IER)\n",
      "print(COF, ROOTR, ROOTI, IER)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(array([  1.20000000e+02,  -5.68434189e-14,   6.27718661e+02,\n",
        "        -2.78000000e+02,   8.30000000e+01]), array([ 0.20146262,  0.20146262,  0.79020404,  0.79020404]), array([-0.30271712,  0.30271712, -2.14629132,  2.14629132]), 0)\n"
       ]
      }
     ],
     "prompt_number": 57
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "___\n",
      "test pas \u00e0 pas sur le premier exemple pour d\u00e9bugguer \n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imp.reload(ER)\n",
      "\n",
      "XCOF = [36., 0, -13, 0, 1]\n",
      "COF = [0., 0, 0, 0, 0]\n",
      "IFIT = 0\n",
      "N = M = 4\n",
      "IER = 0\n",
      "IER = 1\n",
      "NX = N\n",
      "NXX = N + 1\n",
      "N2 = 1\n",
      "KJ1 = N + 1\n",
      "print(\"COF: \", COF)\n",
      "for L in range(1, KJ1+1): \n",
      "    MT = KJ1 - L + 1\n",
      "    COF[MT-1] = XCOF[L-1]\n",
      "print(\"COF: \", COF)\n",
      "ROOTI = zeros(4)\n",
      "ROOTR = zeros(4)\n",
      "\n",
      "\n",
      "print(\"__________________________________________________\")\n",
      "XO = 0.00500101\n",
      "YO = 0.01000101\n",
      "IN = 0\n",
      "X = XO\n",
      "XO = -10.0 * YO\n",
      "YO = -10.0 * X\n",
      "X = XO\n",
      "Y = YO\n",
      "ICT = 0\n",
      "IN = IN + 1\n",
      "while (ICT < 500): \n",
      "    UX = 0.0\n",
      "    UY = 0.0\n",
      "    V = 0.0\n",
      "    YT = 0.0\n",
      "    XT = 1.0\n",
      "    U = COF[N + 1 - 1]\n",
      "    for I in range(1, N + 1): \n",
      "        L = N - I + 1\n",
      "        XT2 = X * XT - Y * YT\n",
      "        YT2 = X * YT + Y * XT\n",
      "        U = U + COF[L - 1] * XT2\n",
      "        V = V + COF[L - 1] * YT2\n",
      "        FI = I\n",
      "        UX = UX + FI * XT * COF[L - 1]\n",
      "        UY = UY - FI * YT * COF[L - 1]\n",
      "        XT = XT2\n",
      "        YT = YT2\n",
      "    SUMSQ = UX * UX + UY * UY\n",
      "#    print(SUMSQ)\n",
      "#    print(X, Y)\n",
      "    if (SUMSQ == 0): \n",
      "        break\n",
      "    DX = (V * UY - U * UX) / SUMSQ\n",
      "    X = X + DX\n",
      "    DY = -(U * UY + V * UX) / SUMSQ\n",
      "    Y = Y + DY\n",
      "    if (abs(DY) + abs(DX) < 1e-5): \n",
      "        break\n",
      "    ICT = ICT + 1\n",
      "print(\"COF, XCOF: \", COF, XCOF)\n",
      "for L in range(1, NXX + 1): \n",
      "    MT = KJ1 - L + 1\n",
      "    print(\"L, MT: \", L, MT)\n",
      "    TEMP = XCOF[MT - 1]\n",
      "    XCOF[MT - 1] = COF[L - 1]\n",
      "    COF[L - 1] = TEMP\n",
      "print(\"COF, XCOF: \", COF, XCOF)\n",
      "ITEMP = N\n",
      "N = NX\n",
      "NX = ITEMP\n",
      "print(IFIT)\n",
      "IFIT = 1\n",
      "XPR = X\n",
      "YPR = Y\n",
      "while ICT < 500: \n",
      "    UX = 0.0\n",
      "    UY = 0.0\n",
      "    V = 0.0\n",
      "    YT = 0.0\n",
      "    XT = 1.0\n",
      "    U = COF[N + 1 - 1]\n",
      "    for I in range(1, N + 1): \n",
      "        L = N - I + 1\n",
      "        XT2 = X * XT - Y * YT\n",
      "        YT2 = X * YT + Y * XT\n",
      "        U = U + COF[L - 1] * XT2\n",
      "        V = V + COF[L - 1] * YT2\n",
      "        FI = I\n",
      "        UX = UX + FI * XT * COF[L - 1]\n",
      "        UY = UY - FI * YT * COF[L - 1]\n",
      "        XT = XT2\n",
      "        YT = YT2\n",
      "    SUMSQ = UX * UX + UY * UY\n",
      "    if (SUMSQ == 0): \n",
      "        break\n",
      "    DX = (V * UY - U * UX) / SUMSQ\n",
      "    X = X + DX\n",
      "    DY = -(U * UY + V * UX) / SUMSQ\n",
      "    Y = Y + DY\n",
      "    if (abs(DY) + abs(DX) < 1e-5): \n",
      "        break\n",
      "    ICT = ICT + 1\n",
      "print(\"X, Y\", X, Y)\n",
      "print(\"IFIT\", IFIT)\n",
      "print(\"COF, XCOF: \", COF, XCOF)\n",
      "for L in range(1, NXX + 1): \n",
      "    MT = KJ1 - L + 1\n",
      "    print(\"L, MT: \", L, MT)\n",
      "    TEMP = XCOF[MT - 1]\n",
      "    XCOF[MT - 1] = COF[L - 1]\n",
      "    COF[L - 1] = TEMP\n",
      "print(\"COF, XCOF: \", COF, XCOF)\n",
      "ITEMP = N\n",
      "N = NX\n",
      "NX = ITEMP\n",
      "IFIT = 0\n",
      "print(abs(Y / X) < 1e-4)\n",
      "Y = 0.0\n",
      "SUMSQ = 0.0\n",
      "ALPHA = X\n",
      "N = N - 1\n",
      "COF[2 - 1] = COF[2 - 1] + ALPHA * COF[1 - 1]\n",
      "for L in range(2, N + 1): \n",
      "    COF[L + 1 - 1] = COF[L + 1 - 1] + ALPHA * COF[L - 1] - SUMSQ * COF[L - 1]\n",
      "print(COF)\n",
      "ROOTI[N2 - 1] = Y\n",
      "ROOTR[N2 - 1] = X\n",
      "N2 = N2 + 1\n",
      "print(SUMSQ, N)\n",
      "print(\"ROOT: \", ROOTR, ROOTI)\n",
      "print(\"__________________________________________________\")\n",
      "print(\"__________________________________________________\")\n",
      "XO = 0.00500101\n",
      "YO = 0.01000101\n",
      "IN = 0\n",
      "X = XO\n",
      "XO = -10.0 * YO\n",
      "YO = -10.0 * X\n",
      "X = XO\n",
      "Y = YO\n",
      "ICT = 0\n",
      "IN = IN + 1\n",
      "while (ICT < 500): \n",
      "    UX = 0.0\n",
      "    UY = 0.0\n",
      "    V = 0.0\n",
      "    YT = 0.0\n",
      "    XT = 1.0\n",
      "    U = COF[N + 1 - 1]\n",
      "    for I in range(1, N + 1): \n",
      "        L = N - I + 1\n",
      "        XT2 = X * XT - Y * YT\n",
      "        YT2 = X * YT + Y * XT\n",
      "        U = U + COF[L - 1] * XT2\n",
      "        V = V + COF[L - 1] * YT2\n",
      "        FI = I\n",
      "        UX = UX + FI * XT * COF[L - 1]\n",
      "        UY = UY - FI * YT * COF[L - 1]\n",
      "        XT = XT2\n",
      "        YT = YT2\n",
      "    SUMSQ = UX * UX + UY * UY\n",
      "#    print(SUMSQ)\n",
      "#    print(X, Y)\n",
      "    if (SUMSQ == 0): \n",
      "        break\n",
      "    DX = (V * UY - U * UX) / SUMSQ\n",
      "    X = X + DX\n",
      "    DY = -(U * UY + V * UX) / SUMSQ\n",
      "    Y = Y + DY\n",
      "    if (abs(DY) + abs(DX) < 1e-5): \n",
      "        break\n",
      "    ICT = ICT + 1\n",
      "print(\"COF, XCOF: \", COF, XCOF)\n",
      "for L in range(1, NXX + 1): \n",
      "    MT = KJ1 - L + 1\n",
      "    print(\"L, MT: \", L, MT)\n",
      "    TEMP = XCOF[MT - 1]\n",
      "    XCOF[MT - 1] = COF[L - 1]\n",
      "    COF[L - 1] = TEMP\n",
      "print(\"COF, XCOF: \", COF, XCOF)\n",
      "ITEMP = N\n",
      "N = NX\n",
      "NX = ITEMP\n",
      "print(IFIT)\n",
      "IFIT = 1\n",
      "XPR = X\n",
      "YPR = Y\n",
      "while ICT < 500: \n",
      "    UX = 0.0\n",
      "    UY = 0.0\n",
      "    V = 0.0\n",
      "    YT = 0.0\n",
      "    XT = 1.0\n",
      "    U = COF[N + 1 - 1]\n",
      "    for I in range(1, N + 1): \n",
      "        L = N - I + 1\n",
      "        XT2 = X * XT - Y * YT\n",
      "        YT2 = X * YT + Y * XT\n",
      "        U = U + COF[L - 1] * XT2\n",
      "        V = V + COF[L - 1] * YT2\n",
      "        FI = I\n",
      "        UX = UX + FI * XT * COF[L - 1]\n",
      "        UY = UY - FI * YT * COF[L - 1]\n",
      "        XT = XT2\n",
      "        YT = YT2\n",
      "    SUMSQ = UX * UX + UY * UY\n",
      "    if (SUMSQ == 0): \n",
      "        break\n",
      "    DX = (V * UY - U * UX) / SUMSQ\n",
      "    X = X + DX\n",
      "    DY = -(U * UY + V * UX) / SUMSQ\n",
      "    Y = Y + DY\n",
      "    if (abs(DY) + abs(DX) < 1e-5): \n",
      "        break\n",
      "    ICT = ICT + 1\n",
      "print(\"X, Y\", X, Y)\n",
      "print(\"IFIT\", IFIT)\n",
      "print(\"COF, XCOF: \", COF, XCOF)\n",
      "for L in range(1, NXX + 1): \n",
      "    MT = KJ1 - L + 1\n",
      "    print(\"L, MT: \", L, MT)\n",
      "    TEMP = XCOF[MT - 1]\n",
      "    XCOF[MT - 1] = COF[L - 1]\n",
      "    COF[L - 1] = TEMP\n",
      "print(\"COF, XCOF: \", COF, XCOF)\n",
      "ITEMP = N\n",
      "N = NX\n",
      "NX = ITEMP\n",
      "IFIT = 0\n",
      "print(abs(Y / X) < 1e-4)\n",
      "Y = 0.0\n",
      "SUMSQ = 0.0\n",
      "ALPHA = X\n",
      "N = N - 1\n",
      "COF[2 - 1] = COF[2 - 1] + ALPHA * COF[1 - 1]\n",
      "for L in range(2, N + 1): \n",
      "    COF[L + 1 - 1] = COF[L + 1 - 1] + ALPHA * COF[L - 1] - SUMSQ * COF[L - 1]\n",
      "print(COF)\n",
      "ROOTI[N2 - 1] = Y\n",
      "ROOTR[N2 - 1] = X\n",
      "N2 = N2 + 1\n",
      "print(SUMSQ, N)\n",
      "print(\"ROOT: \", ROOTR, ROOTI)\n",
      "print(\"__________________________________________________\")\n",
      "print(\"__________________________________________________\")\n",
      "XO = 0.00500101\n",
      "YO = 0.01000101\n",
      "IN = 0\n",
      "X = XO\n",
      "XO = -10.0 * YO\n",
      "YO = -10.0 * X\n",
      "X = XO\n",
      "Y = YO\n",
      "ICT = 0\n",
      "IN = IN + 1\n",
      "while (ICT < 500): \n",
      "    UX = 0.0\n",
      "    UY = 0.0\n",
      "    V = 0.0\n",
      "    YT = 0.0\n",
      "    XT = 1.0\n",
      "    U = COF[N + 1 - 1]\n",
      "    for I in range(1, N + 1): \n",
      "        L = N - I + 1\n",
      "        XT2 = X * XT - Y * YT\n",
      "        YT2 = X * YT + Y * XT\n",
      "        U = U + COF[L - 1] * XT2\n",
      "        V = V + COF[L - 1] * YT2\n",
      "        FI = I\n",
      "        UX = UX + FI * XT * COF[L - 1]\n",
      "        UY = UY - FI * YT * COF[L - 1]\n",
      "        XT = XT2\n",
      "        YT = YT2\n",
      "    SUMSQ = UX * UX + UY * UY\n",
      "#    print(SUMSQ)\n",
      "#    print(X, Y)\n",
      "    if (SUMSQ == 0): \n",
      "        break\n",
      "    DX = (V * UY - U * UX) / SUMSQ\n",
      "    X = X + DX\n",
      "    DY = -(U * UY + V * UX) / SUMSQ\n",
      "    Y = Y + DY\n",
      "    if (abs(DY) + abs(DX) < 1e-5): \n",
      "        break\n",
      "    ICT = ICT + 1\n",
      "print(\"COF, XCOF: \", COF, XCOF)\n",
      "for L in range(1, NXX + 1): \n",
      "    MT = KJ1 - L + 1\n",
      "    print(\"L, MT: \", L, MT)\n",
      "    TEMP = XCOF[MT - 1]\n",
      "    XCOF[MT - 1] = COF[L - 1]\n",
      "    COF[L - 1] = TEMP\n",
      "print(\"COF, XCOF: \", COF, XCOF)\n",
      "ITEMP = N\n",
      "N = NX\n",
      "NX = ITEMP\n",
      "print(IFIT)\n",
      "IFIT = 1\n",
      "XPR = X\n",
      "YPR = Y\n",
      "while ICT < 500: \n",
      "    UX = 0.0\n",
      "    UY = 0.0\n",
      "    V = 0.0\n",
      "    YT = 0.0\n",
      "    XT = 1.0\n",
      "    U = COF[N + 1 - 1]\n",
      "    for I in range(1, N + 1): \n",
      "        L = N - I + 1\n",
      "        XT2 = X * XT - Y * YT\n",
      "        YT2 = X * YT + Y * XT\n",
      "        U = U + COF[L - 1] * XT2\n",
      "        V = V + COF[L - 1] * YT2\n",
      "        FI = I\n",
      "        UX = UX + FI * XT * COF[L - 1]\n",
      "        UY = UY - FI * YT * COF[L - 1]\n",
      "        XT = XT2\n",
      "        YT = YT2\n",
      "    SUMSQ = UX * UX + UY * UY\n",
      "    if (SUMSQ == 0): \n",
      "        break\n",
      "    DX = (V * UY - U * UX) / SUMSQ\n",
      "    X = X + DX\n",
      "    DY = -(U * UY + V * UX) / SUMSQ\n",
      "    Y = Y + DY\n",
      "    if (abs(DY) + abs(DX) < 1e-5): \n",
      "        break\n",
      "    ICT = ICT + 1\n",
      "print(\"X, Y\", X, Y)\n",
      "print(\"IFIT\", IFIT)\n",
      "print(\"COF, XCOF: \", COF, XCOF)\n",
      "for L in range(1, NXX + 1): \n",
      "    MT = KJ1 - L + 1\n",
      "    print(\"L, MT: \", L, MT)\n",
      "    TEMP = XCOF[MT - 1]\n",
      "    XCOF[MT - 1] = COF[L - 1]\n",
      "    COF[L - 1] = TEMP\n",
      "print(\"COF, XCOF: \", COF, XCOF)\n",
      "ITEMP = N\n",
      "N = NX\n",
      "NX = ITEMP\n",
      "IFIT = 0\n",
      "print(abs(Y / X) < 1e-4)\n",
      "Y = 0.0\n",
      "SUMSQ = 0.0\n",
      "ALPHA = X\n",
      "N = N - 1\n",
      "COF[2 - 1] = COF[2 - 1] + ALPHA * COF[1 - 1]\n",
      "for L in range(2, N + 1): \n",
      "    COF[L + 1 - 1] = COF[L + 1 - 1] + ALPHA * COF[L - 1] - SUMSQ * COF[L - 1]\n",
      "print(COF)\n",
      "ROOTI[N2 - 1] = Y\n",
      "ROOTR[N2 - 1] = X\n",
      "N2 = N2 + 1\n",
      "print(SUMSQ, N)\n",
      "print(\"ROOT: \", ROOTR, ROOTI)\n",
      "print(\"__________________________________________________\")\n",
      "print(\"__________________________________________________\")\n",
      "XO = 0.00500101\n",
      "YO = 0.01000101\n",
      "IN = 0\n",
      "X = XO\n",
      "XO = -10.0 * YO\n",
      "YO = -10.0 * X\n",
      "X = XO\n",
      "Y = YO\n",
      "ICT = 0\n",
      "IN = IN + 1\n",
      "while (ICT < 500): \n",
      "    UX = 0.0\n",
      "    UY = 0.0\n",
      "    V = 0.0\n",
      "    YT = 0.0\n",
      "    XT = 1.0\n",
      "    U = COF[N + 1 - 1]\n",
      "    for I in range(1, N + 1): \n",
      "        L = N - I + 1\n",
      "        XT2 = X * XT - Y * YT\n",
      "        YT2 = X * YT + Y * XT\n",
      "        U = U + COF[L - 1] * XT2\n",
      "        V = V + COF[L - 1] * YT2\n",
      "        FI = I\n",
      "        UX = UX + FI * XT * COF[L - 1]\n",
      "        UY = UY - FI * YT * COF[L - 1]\n",
      "        XT = XT2\n",
      "        YT = YT2\n",
      "    SUMSQ = UX * UX + UY * UY\n",
      "#    print(SUMSQ)\n",
      "#    print(X, Y)\n",
      "    if (SUMSQ == 0): \n",
      "        break\n",
      "    DX = (V * UY - U * UX) / SUMSQ\n",
      "    X = X + DX\n",
      "    DY = -(U * UY + V * UX) / SUMSQ\n",
      "    Y = Y + DY\n",
      "    if (abs(DY) + abs(DX) < 1e-5): \n",
      "        break\n",
      "    ICT = ICT + 1\n",
      "print(\"COF, XCOF: \", COF, XCOF)\n",
      "for L in range(1, NXX + 1): \n",
      "    MT = KJ1 - L + 1\n",
      "    print(\"L, MT: \", L, MT)\n",
      "    TEMP = XCOF[MT - 1]\n",
      "    XCOF[MT - 1] = COF[L - 1]\n",
      "    COF[L - 1] = TEMP\n",
      "print(\"COF, XCOF: \", COF, XCOF)\n",
      "ITEMP = N\n",
      "N = NX\n",
      "NX = ITEMP\n",
      "print(IFIT)\n",
      "IFIT = 1\n",
      "XPR = X\n",
      "YPR = Y\n",
      "while ICT < 500: \n",
      "    UX = 0.0\n",
      "    UY = 0.0\n",
      "    V = 0.0\n",
      "    YT = 0.0\n",
      "    XT = 1.0\n",
      "    U = COF[N + 1 - 1]\n",
      "    print(U)\n",
      "    for I in range(1, N + 1): \n",
      "        L = N - I + 1\n",
      "        XT2 = X * XT - Y * YT\n",
      "        YT2 = X * YT + Y * XT\n",
      "        U = U + COF[L - 1] * XT2\n",
      "        V = V + COF[L - 1] * YT2\n",
      "        FI = I\n",
      "        UX = UX + FI * XT * COF[L - 1]\n",
      "        UY = UY - FI * YT * COF[L - 1]\n",
      "        XT = XT2\n",
      "        YT = YT2\n",
      "    SUMSQ = UX * UX + UY * UY\n",
      "    if (SUMSQ == 0): \n",
      "        break\n",
      "    DX = (V * UY - U * UX) / SUMSQ\n",
      "    X = X + DX\n",
      "    DY = -(U * UY + V * UX) / SUMSQ\n",
      "    Y = Y + DY\n",
      "    if (abs(DY) + abs(DX) < 1e-5): \n",
      "        break\n",
      "    ICT = ICT + 1\n",
      "print(\"X, Y\", X, Y)\n",
      "print(\"IFIT\", IFIT)\n",
      "print(\"COF, XCOF: \", COF, XCOF)\n",
      "for L in range(1, NXX + 1): \n",
      "    MT = KJ1 - L + 1\n",
      "    print(\"L, MT: \", L, MT)\n",
      "    TEMP = XCOF[MT - 1]\n",
      "    XCOF[MT - 1] = COF[L - 1]\n",
      "    COF[L - 1] = TEMP\n",
      "print(\"COF, XCOF: \", COF, XCOF)\n",
      "ITEMP = N\n",
      "N = NX\n",
      "NX = ITEMP\n",
      "IFIT = 0\n",
      "print(abs(Y / X) < 1e-4)\n",
      "Y = 0.0\n",
      "SUMSQ = 0.0\n",
      "ALPHA = X\n",
      "N = N - 1\n",
      "COF[2 - 1] = COF[2 - 1] + ALPHA * COF[1 - 1]\n",
      "for L in range(2, N + 1): \n",
      "    COF[L + 1 - 1] = COF[L + 1 - 1] + ALPHA * COF[L - 1] - SUMSQ * COF[L - 1]\n",
      "print(COF)\n",
      "ROOTI[N2 - 1] = Y\n",
      "ROOTR[N2 - 1] = X\n",
      "N2 = N2 + 1\n",
      "print(SUMSQ, N)\n",
      "print(\"ROOT: \", ROOTR, ROOTI)\n",
      "print(\"__________________________________________________\")\n",
      "\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "('COF: ', [0.0, 0, 0, 0, 0])\n",
        "('COF: ', [1, 0, -13, 0, 36.0])\n",
        "__________________________________________________\n",
        "('COF, XCOF: ', [1, 0, -13, 0, 36.0], [36.0, 0, -13, 0, 1])\n",
        "('L, MT: ', 1, 5)\n",
        "('L, MT: ', 2, 4)\n",
        "('L, MT: ', 3, 3)\n",
        "('L, MT: ', 4, 2)\n",
        "('L, MT: ', 5, 1)\n",
        "('COF, XCOF: ', [1, 0, -13, 0, 36.0], [36.0, 0, -13, 0, 1])\n",
        "0\n",
        "('X, Y', -2.9999999999999996, 3.527785968148365e-26)\n",
        "('IFIT', 1)\n",
        "('COF, XCOF: ', [1, 0, -13, 0, 36.0], [36.0, 0, -13, 0, 1])\n",
        "('L, MT: ', 1, 5)\n",
        "('L, MT: ', 2, 4)\n",
        "('L, MT: ', 3, 3)\n",
        "('L, MT: ', 4, 2)\n",
        "('L, MT: ', 5, 1)\n",
        "('COF, XCOF: ', [1, 0, -13, 0, 36.0], [36.0, 0, -13, 0, 1])\n",
        "True\n",
        "[1, -2.9999999999999996, -4.000000000000002, 12.000000000000004, 36.0]\n",
        "(0.0, 3)\n",
        "('ROOT: ', array([-3.,  0.,  0.,  0.]), array([ 0.,  0.,  0.,  0.]))\n",
        "__________________________________________________\n",
        "__________________________________________________\n",
        "('COF, XCOF: ', [1, -2.9999999999999996, -4.000000000000002, 12.000000000000004, 36.0], [36.0, 0, -13, 0, 1])\n",
        "('L, MT: ', 1, 5)\n",
        "('L, MT: ', 2, 4)\n",
        "('L, MT: ', 3, 3)\n",
        "('L, MT: ', 4, 2)\n",
        "('L, MT: ', 5, 1)\n",
        "('COF, XCOF: ', [1, 0, -13, 0, 36.0], [36.0, 12.000000000000004, -4.000000000000002, -2.9999999999999996, 1])\n",
        "0\n",
        "('X, Y', 2.9999999999999996, -1.5331986524126435e-21)\n",
        "('IFIT', 1)\n",
        "('COF, XCOF: ', [1, 0, -13, 0, 36.0], [36.0, 12.000000000000004, -4.000000000000002, -2.9999999999999996, 1])\n",
        "('L, MT: ', 1, 5)\n",
        "('L, MT: ', 2, 4)\n",
        "('L, MT: ', 3, 3)\n",
        "('L, MT: ', 4, 2)\n",
        "('L, MT: ', 5, 1)\n",
        "('COF, XCOF: ', [1, -2.9999999999999996, -4.000000000000002, 12.000000000000004, 36.0], [36.0, 0, -13, 0, 1])\n",
        "True\n",
        "[1, 0.0, -4.000000000000002, 12.000000000000004, 36.0]\n",
        "(0.0, 2)\n",
        "('ROOT: ', array([-3.,  3.,  0.,  0.]), array([ 0.,  0.,  0.,  0.]))\n",
        "__________________________________________________\n",
        "__________________________________________________\n",
        "('COF, XCOF: ', [1, 0.0, -4.000000000000002, 12.000000000000004, 36.0], [36.0, 0, -13, 0, 1])\n",
        "('L, MT: ', 1, 5)\n",
        "('L, MT: ', 2, 4)\n",
        "('L, MT: ', 3, 3)\n",
        "('L, MT: ', 4, 2)\n",
        "('L, MT: ', 5, 1)\n",
        "('COF, XCOF: ', [1, 0, -13, 0, 36.0], [36.0, 12.000000000000004, -4.000000000000002, 0.0, 1])\n",
        "0\n",
        "('X, Y', -2.0, -4.70197740328915e-38)\n",
        "('IFIT', 1)\n",
        "('COF, XCOF: ', [1, 0, -13, 0, 36.0], [36.0, 12.000000000000004, -4.000000000000002, 0.0, 1])\n",
        "('L, MT: ', 1, 5)\n",
        "('L, MT: ', 2, 4)\n",
        "('L, MT: ', 3, 3)\n",
        "('L, MT: ', 4, 2)\n",
        "('L, MT: ', 5, 1)\n",
        "('COF, XCOF: ', [1, 0.0, -4.000000000000002, 12.000000000000004, 36.0], [36.0, 0, -13, 0, 1])\n",
        "True\n",
        "[1, -2.0, -4.000000000000002, 12.000000000000004, 36.0]\n",
        "(0.0, 1)\n",
        "('ROOT: ', array([-3.,  3., -2.,  0.]), array([ 0.,  0.,  0.,  0.]))\n",
        "__________________________________________________\n",
        "__________________________________________________\n",
        "('COF, XCOF: ', [1, -2.0, -4.000000000000002, 12.000000000000004, 36.0], [36.0, 0, -13, 0, 1])\n",
        "('L, MT: ', 1, 5)\n",
        "('L, MT: ', 2, 4)\n",
        "('L, MT: ', 3, 3)\n",
        "('L, MT: ', 4, 2)\n",
        "('L, MT: ', 5, 1)\n",
        "('COF, XCOF: ', [1, 0, -13, 0, 36.0], [36.0, 12.000000000000004, -4.000000000000002, -2.0, 1])\n",
        "0\n",
        "36.0\n",
        "('X, Y', 2.0, 0.0)\n",
        "('IFIT', 1)\n",
        "('COF, XCOF: ', [1, 0, -13, 0, 36.0], [36.0, 12.000000000000004, -4.000000000000002, -2.0, 1])\n",
        "('L, MT: ', 1, 5)\n",
        "('L, MT: ', 2, 4)\n",
        "('L, MT: ', 3, 3)\n",
        "('L, MT: ', 4, 2)\n",
        "('L, MT: ', 5, 1)\n",
        "('COF, XCOF: ', [1, -2.0, -4.000000000000002, 12.000000000000004, 36.0], [36.0, 0, -13, 0, 1])\n",
        "True\n",
        "[1, 0.0, -4.000000000000002, 12.000000000000004, 36.0]\n",
        "(0.0, 0)\n",
        "('ROOT: ', array([-3.,  3., -2.,  2.]), array([ 0.,  0.,  0.,  0.]))\n",
        "__________________________________________________\n"
       ]
      }
     ],
     "prompt_number": 32
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "numpy.polynomial.polynomial.polyroots((36, 0, -13, 0, 1))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 14,
       "text": [
        "array([-3., -2.,  2.,  3.])"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "numpy.polynomial.polynomial.polyroots((83, -278, 720, -238, 120))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 15,
       "text": [
        "array([ 0.20146262-0.30271712j,  0.20146262+0.30271712j,\n",
        "        0.79020404-2.14629132j,  0.79020404+2.14629132j])"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imp.reload(ER)\n",
      "COF = numpy.array([1., 2])\n",
      "X = 1\n",
      "NX = 1\n",
      "NXX = 1\n",
      "Y = 1\n",
      "SUMSQ = 0\n",
      "ALPHA = 0\n",
      "M = N = 1\n",
      "(U, X, Y, SUMSQ, N, DX, DY) = ER.POLRT_eval_pol_and_deriv(COF, X, Y, SUMSQ, N)\n",
      "print(U, X, Y, SUMSQ, N, DX, DY) \n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(1, -2.0, 0.0, 1.0, 1, -3.0, -1.0)\n"
       ]
      }
     ],
     "prompt_number": 52
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}