{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Formulas\n",
    "\n",
    "$$d^{(j)}(z,t) = exp \\left( -\\sum_{n=1}^\\infty \\frac{z^n}{n} \\sum_{|\\underline i|=n} \\frac {p_{\\underline i}}{1-\\lambda_{\\underline i}^{-2}} \\left(\\prod_{l=0}^{n-1} \\frac {\\| A_j x_{\\sigma^l \\underline i} \\|}{\\|x_{\\sigma^l \\underline i}\\|}\\right)^t \\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$d^{(j)}(z,0) = exp \\left( -\\sum_{n=1}^\\infty \\frac{z^n}{n} \\sum_{|\\underline i|=n} \\frac {p_{\\underline i}}{1-\\lambda_{\\underline i}^{-2}} \\right) = \\color{green}{exp(-S_0)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \\frac{\\partial}{\\partial t} d^{(j)}(z,t)_{|(z,0)} = \\left(-\\sum_{n=1}^\\infty \\frac{z^n}{n} \\sum_{|\\underline i|=n} \\frac {p_{\\underline i}}{1-\\lambda_{\\underline i}^{-2}} \\log \\left(\\prod_{l=0}^{n-1} \\frac {\\| A_j x_{\\sigma^l \\underline i} \\|}{\\|x_{\\sigma^l \\underline i}\\|}\\right) \\right) \\cdot exp \\left( -\\sum_{n=1}^\\infty \\frac{z^n}{n} \\sum_{|\\underline i|=n} \\frac {p_{\\underline i}}{1-\\lambda_{\\underline i}^{-2}} \\right) = \\color{green}{-S_1 \\exp(-S_0)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$x_{i_2, \\dots, i_n, i_1} = A_{i_1}^{-1} x_{i_1, \\dots, i_n}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "def terms(A,p,N):\n",
    "    Ainv = [M.inverse() for M in A]\n",
    "    R.<z> = CBF['z']\n",
    "    k = len(A)\n",
    "\n",
    "    S0,S1 = 0,[0]*k\n",
    "\n",
    "    for n in range(1,N):\n",
    "        cn_S0 = 0\n",
    "        cn_S1 = [0]*k\n",
    "        \n",
    "        for i in itertools.product(range(k), repeat=n):\n",
    "            M = prod(A[s] for s in i)\n",
    "            l, x = perr_frob(M)\n",
    "            scalar = prod(p[s] for s in i)/(1-l^(-2))\n",
    "\n",
    "            P = [1]*k\n",
    "            for l in range(n):\n",
    "                for j in range(k):\n",
    "                    P[j] *= CBF((A[j]*x).norm(1)/x.norm(1))\n",
    "                x = Ainv[i[l]]*x\n",
    "                \n",
    "            # Use Arb polynomials\n",
    "            cn_S0 += scalar\n",
    "            for j in range(k):\n",
    "                cn_S1[j] += scalar*log(P[j])\n",
    "\n",
    "        for j in range(k):\n",
    "            S1[j] += z^n/n*cn_S1[j]\n",
    "        S0 += z^n/n*cn_S0\n",
    "    \n",
    "    an = (-S0)._exp_series(N)\n",
    "    cn = [f1._mul_trunc_(-S1[j],N) for j in range(k)]\n",
    "    \n",
    "    return (an, cn)\n",
    "\n",
    "def perr_frob(M):\n",
    "    ev = M.eigenvectors_right()\n",
    "    val, pos = 0, -1\n",
    "    for i in range(len(ev)):\n",
    "        v = ev[i]\n",
    "        if v[0] > val :\n",
    "            val = v[0]\n",
    "            pos = i\n",
    "    return(ev[pos][0], ev[pos][1][0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$d^{(j)}(z,0) = 1 + \\sum_{n=1}^\\infty a_n^{(j)}(0) z^n  = 1 + \\sum_{n=1}^\\infty a_n(0) z^n$$\n",
    "\n",
    "$$\\frac{d}{dz} d^{(j)}(z,0) = \\sum_{n=1}^\\infty n a_n^{(j)}(0) z^{n-1}$$\n",
    "$$\\frac{d}{dz} d^{(j)}(1,0) = \\sum_{n=1}^\\infty b_n^{(j)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \\frac{\\partial}{\\partial t} d^{(j)}(z,t)_{|(z,0)} = \\sum_{n=1}^\\infty \\frac{\\partial a_n^{(j)}}{\\partial t}(0) z^{n}$$\n",
    "\n",
    "$$ \\frac{\\partial}{\\partial t} d^{(j)}(z,t)_{|(1,0)} = \\sum_{n=1}^\\infty c_n^{(j)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\lambda_n = \\sum_{j=1}^k \\frac{\\sum_{i=1}^n c_i^{(j)}}{\\sum_{i=1}^n b_i^{(j)}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.14331103505173208520461882520715612957215148442230056673747095506825521805 +/- 7.07e-75]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prec = 256\n",
    "CBF = ComplexBallField(prec)\n",
    "\n",
    "A = [\n",
    "    matrix([[2,1],[1,1]]),\n",
    "    matrix([[3,1],[2,1]])\n",
    "]\n",
    "\n",
    "p = [1/2,1/2]\n",
    "\n",
    "an, cn = terms(A,p,5)\n",
    "sum(cn[j].subs(z=1) for j in range(2))/(2*an.differentiate().subs(z=1))\n",
    "# Where does the factor 2 come from ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.1433110351029492458432518536555882994025847950806845414580494534622463138 +/- 5.87e-74]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "an, cn = terms(A,p,10)\n",
    "sum(cn[j].subs(z=1) for j in range(2))/(2*an.differentiate().subs(z=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1.1433110351029492458432518536555882994025"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.2",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}