Correlators_ext_complex.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import division, print_function\n",
    "\n",
    "import numpy, math, cmath\n",
    "np = numpy\n",
    "import mpmath, scipy.special\n",
    "import operator\n",
    "from cmath import sin, sinh, exp\n",
    "from scipy.interpolate import InterpolatedUnivariateSpline\n",
    "\n",
    "\n",
    "from Auxiliary_classes import complex_quadrature, numerical\n",
    "from CFT import Charge, Dimension\n",
    "from Blocks import Block, BlockNum "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Upsilon functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Upsilon function $\\Upsilon_b(P)$ is defined as \n",
    "$$\n",
    "\\Upsilon_b(P) = \\exp \\int_0^\\infty \\frac{dt}{t} \\left[\\frac{\\sin^2 Pt}{\\sinh bt\\sinh\\frac{t}{b}} - P^2 e^{-2t}\\right]\n",
    "$$\n",
    "which is valid if $|\\Im P |\\leq \\frac{Q}{2}$ with $Q = b + \\frac{1}{b}$. For momentums outside this strip, we should use shift equations such as\n",
    "$$\n",
    "\\frac{\\Upsilon_b(P-ib)}{\\Upsilon_b(P)} = b^{-b^2-2ibP}\\gamma(\\tfrac{bQ}{2} + ibP)\n",
    "$$\n",
    "In the case $b=1$ we have \n",
    "$$\n",
    "\\Upsilon_1(iP) = G(1+P)G(1-P)\n",
    "$$\n",
    "where we use Barnes' $G$-function.\n",
    "\n",
    "There is an alternative representation as the infinite product\n",
    "$$\n",
    "\\Upsilon_b(P) = e^{\\lambda_b P^2} \\prod_{m,n=0}^\\infty F\\left(\\frac{P^2}{\\left(\\frac{Q}{2}+mb+nb^{-1}\\right)^2}\\right) \\qquad \\mathrm{where} \\qquad F(x) = (1+x)\\exp(-x)\n",
    "$$\n",
    "where we define\n",
    "$$\n",
    "\\lambda_b = \\int_0^\\infty \\frac{dt}{t} \\left[\\frac{t^2}{\\sinh bt\\sinh\\frac{t}{b}} - e^{-2t}\\right]\n",
    "$$\n",
    "The factor $e^{\\lambda_b P^2}$ plays no role in most calculations. We neglect it by default, while providing\n",
    "a static method for computing it. We have $\\lambda_1 = 1 + \\gamma$ where $\\gamma$ is the Euler-Mascheroni constant.\n",
    "\n",
    "The product representation can be improved by approximating the neglected factors with (the exponential of) an integral, according to the formula\n",
    "$$\n",
    "\\Upsilon_b(P) \\sim \\prod_{m=0}^M\\prod_{n=0}^N \\left(1+\\frac{P^2}{y_{m,n}^2}\\right) \n",
    "\\times\n",
    "e^{ - \\sum_{m=0}^M\\sum_{n=0}^N \\frac{P^2}{y_{m,n}^2} - Pb\\sum_{m=0}^M \\varphi'\\left(\\frac{y_{m,N+\\frac12}}{P}\\right) - Pb^{-1}\\sum_{n=0}^N\\varphi'\\left(\\frac{y_{M+\\frac12,n}}{P}\\right) + P^2 \\varphi\\left(\\frac{y_{M+\\frac12,N+\\frac12}}{P}\\right) }\n",
    "$$\n",
    "where we define $y_{m,n} = \\frac{Q}{2}+mb+nb^{-1}$ and \n",
    "$$\n",
    "\\varphi(y) = \\tfrac32 + \\tfrac12 (y^2-1)\\log(1+\\tfrac{1}{y^2}) + iy\\log\\tfrac{y+i}{y-i}\n",
    "\\\\\n",
    " \\varphi'(y) = y\\log(1+\\tfrac{1}{y^2}) +\\tfrac{1}{y} +i\\log\\tfrac{y+i}{y-i}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Upsilons:\n",
    "    \"\"\" Computing products and ratios of Upsilon functions. \"\"\"\n",
    "    \n",
    "    def __init__(self, b = 1, Plist = [[1, True]]):\n",
    "        \"\"\" We need the parameter b, and lists of Ps for numerator \n",
    "        (True) or denominator (False) Upsilon functions. \"\"\" \n",
    "        self.b = b\n",
    "        self.Plist = Plist\n",
    "        \n",
    "    @staticmethod\n",
    "    def phi(y):\n",
    "        return 3/2 + (y**2 - 1)/2 * cmath.log(1 + 1/y**2) + 1j*y*cmath.log((y+1j)/(y-1j))\n",
    "    \n",
    "    @staticmethod\n",
    "    def phip(y):\n",
    "        return y * cmath.log(1 + 1/y**2) + 1 / y + 1j*cmath.log((y+1j)/(y-1j))\n",
    "    \n",
    "    def product(self, L = 20, improved = True, normalized = False):\n",
    "        \"\"\" Computes the product, where m and n are bounded by integer cutoffs M and N, chosen so that \n",
    "        we compute O(L^2) factors and so that Mb = Nb^{-1}.\n",
    "        If improved, includes the approximation of the neglected factors. \n",
    "        If normalized, agrees with the integral representation.\n",
    "        \"\"\"\n",
    "            \n",
    "        b = self.b\n",
    "        binv = 1 / b\n",
    "        q = ( b + binv ) / 2      \n",
    "        \n",
    "        M = int(L/abs(b)) + 1\n",
    "        N = int(L*abs(b)) + 1\n",
    "        \n",
    "        Plist = [pair for pair in self.Plist if pair[0] != 0]  \n",
    "        \"\"\" We remove the zeros from the list, as they give trivial contributions and nevertheless\n",
    "        cause errors if improved. \"\"\"\n",
    "        Plists = [[pair[0]**2, pair[1]] for pair in Plist]\n",
    "        \n",
    "        prod = 1\n",
    "        exponent = 0\n",
    "        \n",
    "        for m in range(M + 1):\n",
    "            for n in range(N + 1):\n",
    "                y = q + m*b + n*binv\n",
    "                x = 1 / y**2\n",
    "                exponent += x\n",
    "                for pair in Plists:\n",
    "                    factor = 1 + pair[0] * x\n",
    "                    prod = prod * factor if pair[1] else prod / factor\n",
    "\n",
    "        expFactor = 0\n",
    "        for pair in Plists:\n",
    "            expFactor = expFactor + pair[0] if pair[1] else expFactor - pair[0]\n",
    "        exponent *= expFactor\n",
    "        \n",
    "        if improved:\n",
    "            \n",
    "            addterm = 0\n",
    "            for m in range(M + 1):\n",
    "                y = q + m*b + (N+1/2)*binv\n",
    "                for pair in Plist:\n",
    "                    P = pair[0]\n",
    "                    term = P * Upsilons.phip(y / P)\n",
    "                    addterm = addterm + term if pair[1] else addterm - term\n",
    "            addterm *= b \n",
    "            exponent += addterm\n",
    "            \n",
    "            addterm = 0\n",
    "            for n in range(N + 1):\n",
    "                y = q + (M+1/2)*b + n*binv\n",
    "                for pair in Plist:\n",
    "                    P = pair[0]\n",
    "                    term = P * Upsilons.phip(y / P)\n",
    "                    addterm = addterm + term if pair[1] else addterm - term\n",
    "            addterm *= binv\n",
    "            exponent += addterm\n",
    "            \n",
    "            y = q + (M+1/2)*b + (N+1/2)*binv\n",
    "            for pair in Plist:\n",
    "                term = Upsilons.phi(y / pair[0])\n",
    "                exponent = exponent - term if pair[1] else exponent + term\n",
    "                \n",
    "        if normalized:\n",
    "            shift = 0\n",
    "            for pair in Plists:\n",
    "                shift = shift + pair[0] if pair[1] else shift - pair[0]\n",
    "            exponent -= shift * Upsilons.lambda_b(b)\n",
    "                \n",
    "        return prod * cmath.exp(- exponent)\n",
    "    \n",
    "    @staticmethod\n",
    "    def lambda_b(b, cutoff = 10):\n",
    "        \"\"\" Computes the number lambda_b which controls \n",
    "        the discrepancy between the product and integral representations. \"\"\"\n",
    "        def integrand(t):\n",
    "            return (t**2 / cmath.sinh(b*t) / cmath.sinh(t/b) - cmath.exp(-2*t))/t  \n",
    "        return complex_quadrature(integrand, 0, cutoff)[0]    \n",
    "    \n",
    "    def integral(self, cutoff = 100):\n",
    "        \"\"\" Directly computes the integral over t. Uses recursion relations if need be. \"\"\"\n",
    "        result = 1\n",
    "        for pair in self.Plist:\n",
    "            factor = Upsilons.upsilon(pair[0], self.b)\n",
    "            result = result * factor if pair[1] else result / factor\n",
    "        return result      \n",
    "\n",
    "    @staticmethod \n",
    "    def upsilon(P, b, cutoff = 100):\n",
    "        \"\"\" Computes the Upsilon function, using recursion relations if need be. \"\"\"\n",
    "        bs = [b, -b, 1/b, -1/b]\n",
    "        for oneb in bs:\n",
    "            if oneb.real > b.real:\n",
    "                b = oneb\n",
    "        Q = b + 1/b\n",
    "        q = Q.real/2\n",
    "        alpha = Q/2 + 1j*P\n",
    "        \n",
    "        def integrand(t):\n",
    "            return ((cmath.sin(P*t))**2 / cmath.sinh(b*t) / cmath.sinh(t/b) - P**2*cmath.exp(-2*t))/t  \n",
    "        \n",
    "        if P.imag >= q:\n",
    "            return numerical(Upsilons.upsilon(P - 1j*b, b, cutoff) * mpmath.gamma(1-b*alpha) \n",
    "                           / mpmath.gamma(b*alpha) * b**(2*b*alpha - 1))\n",
    "        elif P.imag <= -q:\n",
    "            return numerical(Upsilons.upsilon(P + 1j*b, b, cutoff) * mpmath.gamma(b*(alpha-b))\n",
    "                           / mpmath.gamma(1-b*(alpha-b)) * b**(1-2*b*(alpha - b)))\n",
    "        else:    \n",
    "            return cmath.exp(complex_quadrature(integrand, 0, cutoff)[0]) \n",
    "        \n",
    "\n",
    "from mpmath import gamma\n",
    "class Upsilonprime:\n",
    "    \"\"\" calcuates 1 / derivatives of the upsilon function\n",
    "    at alpha = - m b - n / b, as function of alpha\"\"\"\n",
    "    def __init__(self, b, mmax = 10, nmax = 10):\n",
    "        self.b=b\n",
    "        Q = b + 1 / b\n",
    "        def smallgamma(z):\n",
    "            try:\n",
    "                res =  gamma(z) / gamma(1 - z)\n",
    "                return numerical(res)\n",
    "            except ValueError:\n",
    "                return 0.j\n",
    "        upsb = Upsilons(b, [[(b - Q/2) / 1.j, True]]).integral()\n",
    "        #print( upsb )\n",
    "        table = np.zeros((mmax, nmax), dtype = complex)\n",
    "        table[0,0] = 1.  /  upsb\n",
    "        for j in range(0, mmax):\n",
    "            x = -(j * b)\n",
    "            if j > 0:\n",
    "                table[j,0] = numerical(table[j - 1, 0] * (b ** (1 - 2 * b * x)) * smallgamma(b * x))\n",
    "            for k in range(1, nmax):\n",
    "                x = - k / b - j * b\n",
    "                table[j,k] = numerical(table[j, k - 1] * (b ** (2 * x / b - 1)) * smallgamma(x / b))\n",
    "        self.table = table\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Liouville theory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reflection coefficient"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reflection coefficient associated to a field is \n",
    "$$\n",
    "R = - \\frac{\\Gamma(b(2\\alpha-Q)) \\Gamma(b^{-1}(2\\alpha-Q))}{\\Gamma(b(Q-2\\alpha)) \\Gamma(b^{-1}(Q-2\\alpha))}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Reflection:\n",
    "    \"\"\" Liouville reflection coefficient. \"\"\"\n",
    "    \n",
    "    def __init__(self, field):\n",
    "        \n",
    "        alpha = field.get('alpha')\n",
    "        b = field.charge.get('b')\n",
    "        x = 2*alpha - b - 1/b\n",
    "        self.value = - ( scipy.special.gamma(b*x) * scipy.special.gamma(x/b) \n",
    "                    / scipy.special.gamma(-b*x) / scipy.special.gamma(-x/b) )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Three-point correlation function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Liouville three-point function is given by the DOZZ formula,\n",
    "$$\n",
    "C(P_1, P_2, P_3) = \\frac{\\prod_{i=1}^3 \\Upsilon_b(2P_i -i\\frac{Q}{2})}{\\prod_{\\pm\\pm}\\Upsilon_b(P_1\\pm P_2 \\pm P_3)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ThreePoint:\n",
    "    \"\"\" Liouville three-point function. \"\"\"\n",
    "    \n",
    "    def __init__(self, Ps, b, smart = -1):\n",
    "        \"\"\" We encode the three-point function as an Upsilons object.\n",
    "            when smart = 0,1,2, only the Ps[smart] dependent factors will be calcuated. \"\"\"\n",
    "        Q = b + 1/b\n",
    "        Plist = ([[2*P - 1j * Q/2, True] for P in ([Ps[smart]] if smart in [0,1,2] else Ps)] \n",
    "                 + [[P, False] for P in Z3.Psums(Ps)])\n",
    "        self.upsilons = Upsilons(b, Plist)\n",
    "class ThreePointInt:\n",
    "    \"\"\" Liouville three-point function, without the two external Upsilon(2alpha) 's \"\"\"\n",
    "    \n",
    "    def __init__(self, Ps, b):\n",
    "        \"\"\" We encode the three-point function as an Upsilons object.\n",
    "            when smart = 0,1,2, only the Ps[smart] dependent factors will be calcuated. \"\"\"\n",
    "        Q = b + 1/b\n",
    "        Plist = ([[2*Ps[2] - 1j * Q/2, True]] \n",
    "                 + [[P, False] for P in Z3.Psums(Ps)])\n",
    "        self.upsilons = Upsilons(b, Plist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Four-point correlation function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We consider four fields, parametrized by real momentums $P_1,P_2,P_3,P_4$. We want to compute the corresponding four-point function $Z$ in Liouville theory, by using its decomposition in some channel,\n",
    "$$\n",
    "Z = \\prod_{i=1}^4 \\frac{\\Upsilon_b(2P_i-i\\frac{Q}{2})}{\\Upsilon_b(2P_i)} \\int_0^\\infty dP\\ K(P) |B(P)|^2\n",
    "$$\n",
    "where $K(P)$ is a combination of structure constants, $B(P)$ is a conformal block, where complex conjugation only acts on the position $x$. In the case of the $s$-channel, we have\n",
    "$$\n",
    "K_s(P) = P^2 \\exp \\int_0^\\infty \\frac{dt}{t}\\frac{4\\varphi(t)}{\\sinh bt\\sinh b^{-1}t }\n",
    "$$\n",
    "where we used\n",
    "$$\n",
    "\\varphi(t) = 4S(1-S)\\sinh^2 \\tfrac12(b-\\tfrac{1}{b})t-2S^2\n",
    "+S\\left(2\\sum_{i=1}^4S_i-4S_1S_2-4S_3S_4\\right)-(S_1-S_2)^2-(S_3-S_4)^2\n",
    "$$\n",
    "where $S_i=\\sin^2P_it$ and $S=\\sin^2Pt$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class FourPoint:\n",
    "    \"\"\" Liouville four-point function. <\\prod_i V_{Q/2 + ip_i}(z_i)>, z_i = (0, x, 1, \\infty)\"\"\"\n",
    "    \n",
    "       \n",
    "    def __init__(self, t_channel = False, Ps = [i + .5 for i in range(4)], b = math.sqrt(2) - 1j,\n",
    "                 Nmax = 4, Pcutoff = 10, Kcutoff = 40,\n",
    "                 reduced = False, spline = 0, normalized = False, useK2 = 0):\n",
    "        \"\"\" We need the values of b, four momentums P, the order Nmax of the blocks' expansion, \n",
    "        two integration cutoffs, and the indication whether to use reduced blocks. \n",
    "        We also give the number of points in the interpolating spline\n",
    "        for K(P): if this number is zero, then no interpolation is used. The interpolating\n",
    "        spline dispenses us from computing K(P) again whenever we change N or q. For precision computations\n",
    "        the spline should be about 2000. For graphical display, much lower values such as 100 may be enough.\n",
    "        There is an option to include the normalization prefactor, which is absent by default.\n",
    "        useK2: 0 - use K(P), -1: use K2(P) with prefactor, 2:use K2(P) without prefactor.\n",
    "        the K(P) trick works only if |imag(P)| < Q / 4, so will switch to K2 automatically.\n",
    "        \"\"\"\n",
    "        useK2 = 0\n",
    "        for P in Ps:\n",
    "            if abs(P.imag) > (b + 1 / b).real / 4:\n",
    "                useK2 = 2\n",
    "                break \n",
    "        print( \"use K2\", useK2 )\n",
    "        self.useK2 = useK2\n",
    "        self.b = b\n",
    "        self.t_channel = t_channel\n",
    "        self.Kcutoff = Kcutoff\n",
    "        self.Pcutoff = Pcutoff\n",
    "        self.charge = Charge('b', b)        \n",
    "        self.Ps = Ps\n",
    "        self.computeBlock(Nmax)\n",
    "        self.Nmax = Nmax\n",
    "        self.reduced = reduced\n",
    "        self.prefactor = 1.\n",
    "        self.normalized = normalized\n",
    "        self.upsprime = Upsilonprime(b, 5, 5)\n",
    "        if normalized and useK2 in [0, 2]:\n",
    "            if useK2 == 0:\n",
    "                Plist = [[2*P - 1j * (b+1/b)/2, True] for P in Ps] + [[2*P, False] for P in Ps]\n",
    "            if useK2 == 2:\n",
    "                Plist = [[2*P - 1j * (b+1/b)/2, True] for P in Ps]\n",
    "            self.upsilons = Upsilons(b, Plist)\n",
    "            self.prefactor = self.upsilons.integral()\n",
    "            if useK2 == 0:\n",
    "                self.prefactor *= (self.upsprime.table[0,0] ** (-2)) * 4\n",
    "        self.spline = spline\n",
    "        self.useSpline = spline > 0            \n",
    "        if self.useSpline:\n",
    "            self.Kspline = FourPoint.complex_spline([i*Pcutoff/spline for i in range(spline+1)], \n",
    "                                                    [self.K, self.K2, self.K21][self.useK2])\n",
    "        self.polelist = self.getcrossedpoles()\n",
    "        self.CCdiscrete = self.getCCdiscrete(self.polelist) \n",
    "    def computeBlock(self, N):\n",
    "        \"\"\" Computing the relevant Block. This method is here so that we can change N without\n",
    "        reconstructing the object.\n",
    "        \"\"\"       \n",
    "        self.block = Block([Dimension('P', self.Ps[i], self.charge) for i in range(4)], N, self.t_channel)\n",
    "\n",
    "    @staticmethod\n",
    "    def complex_spline(Pvals, K):\n",
    "        \"\"\" Computing spline functions for real and imaginary parts of a function K\n",
    "        on argument values Pvals. \"\"\"\n",
    "        Kvals = [K(Pval) for Pval in Pvals]\n",
    "        Kreals = [[Kval.real for Kval in Kvals], [Kval.imag for Kval in Kvals]]\n",
    "        return [InterpolatedUnivariateSpline(Pvals, Kreals[0]), \n",
    "                InterpolatedUnivariateSpline(Pvals, Kreals[1])]\n",
    "    def getallpoles(self, nmax = 10, mmax = 10):\n",
    "        \"\"\"We look for all poles in Ps of the structure constant \"\"\"\n",
    "        Ps = self.Ps\n",
    "        if self.t_channel:\n",
    "            Ps = [Ps[1], Ps[2], Ps[3], Ps[0]]\n",
    "        b = self.b\n",
    "        Q = b + 1. / b\n",
    "        R = np.array([[b * k + j / b for k in range(nmax)] for j in range(mmax)])\n",
    "        Poles = np.zeros((16, nmax, mmax), dtype = complex)\n",
    "        for series in range(16):\n",
    "            PA, PB = (Ps[0], Ps[1]) if series < 8 else (Ps[2], Ps[3])\n",
    "            signs = [1,1,1]\n",
    "            s1 = series\n",
    "            for k in range(3):\n",
    "                signs[k] = -1 if s1 % 2 == 1 else 1\n",
    "                s1 = s1 // 2\n",
    "            Poles[series, :, :] = (R + Q / 2) * 1.j * signs[0] + PA * signs[1] + PB * signs[2] \n",
    "        return Poles\n",
    "    \n",
    "    def getcrossedpoles(self):\n",
    "        \"\"\"return the poles that have crossed the spectrum. Return a list of (p, k, j, series), \n",
    "           where p is the p-position of the pole, k, j are the coordinate k b + j / b in the lattice\n",
    "           and series gives the mupsilon that produces the pole. \"\"\"\n",
    "        res = []\n",
    "        Ps = self.Ps\n",
    "        if self.t_channel:\n",
    "            Ps = [Ps[1], Ps[2], Ps[3], Ps[0]]\n",
    "        b = self.b\n",
    "        Q = b + 1. / b\n",
    "        for series in range(16):\n",
    "            #the lowest bit is \\pm P_int, the next two are \\pm P_1, P_2, /comp \n",
    "            PA, PB = (Ps[0], Ps[1]) if series < 8 else (Ps[2], Ps[3])\n",
    "            signs = [1,1,1]\n",
    "            s1 = series\n",
    "            for k in range(3):\n",
    "                signs[k] = -1 if s1 % 2 == 1 else 1\n",
    "                s1 = s1 // 2\n",
    "            leading = (Q / 2) * 1.j * signs[0] + PA * signs[1] + PB * signs[2]\n",
    "            borne = -leading.imag * signs[0] \n",
    "            if borne < 0: \n",
    "                continue\n",
    "            else: \n",
    "                print( \"crossing in series\", series, 'signs:', signs )\n",
    "            for j in range(int(borne / b.real) + 1):\n",
    "                jb = borne - j * b\n",
    "                for k in  range(int(jb.imag * b.real) + 1):\n",
    "                    pos = (j * b + k / b + Q / 2) * 1.j * signs[0] + signs[1] * PA + signs[2] * PB\n",
    "                    res.append((pos, j, k, series))\n",
    "        return res\n",
    "    def getCCdiscrete(self, polelist):\n",
    "        \"\"\"calculate the product of structure constants at the crossed poles.\"\"\"\n",
    "        Ps = self.Ps\n",
    "        b = self.b\n",
    "        Q = b + 1. / b\n",
    "        if self.t_channel:\n",
    "            Ps = [Ps[1], Ps[2], Ps[3], Ps[0]]\n",
    "        res = []\n",
    "        for (pos, j, k, series) in polelist:\n",
    "            Plists = [[pos * 2 + Q / 2.j, True], [- pos * 2 + Q / 2.j, True]]\n",
    "            signs = [[1, 1],[-1, 1],[1, -1],[-1, -1]]\n",
    "            for s in range(8):\n",
    "                if s ==  series // 2:\n",
    "                    continue\n",
    "                (Pa, Pb) = (Ps[0], Ps[1]) if s < 4 else (Ps[2], Ps[3])\n",
    "                newp = -Pa * signs[s % 4][0] - Pb * signs[s % 4][1] + pos\n",
    "                Plists.append([newp, False])\n",
    "            prod = Upsilons(b, Plists).integral()\n",
    "            res.append(2 * np.pi * self.upsprime.table[j,k] * prod)\n",
    "        return res \n",
    "    def K(self, P, *args):\n",
    "        \"\"\" We compute the combination of structure constants. \"\"\"\n",
    "        Ps = self.Ps\n",
    "        if self.t_channel:\n",
    "            Ps = [Ps[1], Ps[2], Ps[3], Ps[0]]\n",
    "            \n",
    "        def integrand(t):\n",
    "            Ss = [sin(Ps[i]*t)**2 for i in range(4)]\n",
    "            S = sin(P*t)**2\n",
    "            sh = sinh((self.b - 1/self.b)/2*t)**2\n",
    "            phi = (4*S*(1-S)*sh - 2*S**2 + S*(2*(Ss[0]+Ss[1]+Ss[2]+Ss[3]) \n",
    "                                        - 4*Ss[0]*Ss[1] - 4*Ss[2]*Ss[3]) - (Ss[0]-Ss[1])**2 - (Ss[2]-Ss[3])**2)\n",
    "            return 4*phi/t/sinh(self.b*t)/sinh(t/self.b)\n",
    "        \n",
    "        return P**2*exp(complex_quadrature(integrand, 0, self.Kcutoff)[0])\n",
    "    def K2(self, P):\n",
    "        Ps = self.Ps\n",
    "        if self.t_channel:\n",
    "            Ps = [Ps[1], Ps[2], Ps[3], Ps[0]]\n",
    "        return (ThreePoint([Ps[0], Ps[1], P], self.b).upsilons.integral() * \n",
    "                ThreePoint([Ps[2], Ps[3], -P], self.b).upsilons.integral())\n",
    "    def K21(self, P):\n",
    "        Ps = self.Ps\n",
    "        if self.t_channel:\n",
    "            Ps = [Ps[1], Ps[2], Ps[3], Ps[0]]\n",
    "        return (ThreePointInt([Ps[0], Ps[1], P], self.b).upsilons.integral() * \n",
    "                ThreePointInt([Ps[2], Ps[3], -P], self.b).upsilons.integral())\n",
    "\n",
    "    def value(self, q = .05):\n",
    "        \"\"\" Computes the numerical value of the four-point function. The integrand decreases quickly at large\n",
    "        P so the cutoff can be quite small. Anyway, computing K(P >= 12) is problematic. (Warning about number\n",
    "        of subdivisions from integration over t.) However this value of Pcutoff might have to depend on the \n",
    "        parameters b and Ps. We take an asymmetric range (0, Pcutoff) because the integrand is an even \n",
    "        function of P.   \n",
    "        Discrete fusion contritbutions added. \n",
    "        \"\"\"\n",
    "        \n",
    "        blockNum = BlockNum(self.block, q)\n",
    "        blockNum_c = BlockNum(self.block, q.conjugate())\n",
    "        if self.useSpline:\n",
    "            kfunc = lambda P: self.Kspline[0](P) + self.Kspline[1](P)*1j\n",
    "        else: \n",
    "            kfunc = [self.K, self.K2, self.K21][self.useK2]\n",
    "        integrand = lambda P: numerical(kfunc(P) * blockNum.value(P**2, self.reduced) \n",
    "                                        * blockNum_c.value(P**2, self.reduced))\n",
    "        \"\"\" We apply the method 'numerical' because 'mpmath' can produce objects of the \n",
    "        'mpc' and 'mpf' classes which are not recognized as complex numbers and floats \n",
    "        respectively by the 'quad' integration method.\n",
    "        \"\"\"\n",
    "        res = complex_quadrature(integrand, 0, self.Pcutoff)[0]\n",
    "          \n",
    "        for i in range(len(self.polelist)):\n",
    "            pos = self.polelist[i][0]\n",
    "            res += (.5 *  self.CCdiscrete[i] * blockNum.value(pos**2, self.reduced) \n",
    "                       * blockNum_c.value(pos**2, self.reduced))\n",
    "            # .5 because we integrate over P = 0 to +infty, so we take the half of the residue    \n",
    "        if self.normalized:\n",
    "            res *= self.prefactor\n",
    "        \n",
    "        return res\n",
    "    def value_discrete(self, q = .05):\n",
    "        \"\"\"return the list of discrete contributions to 4p function.\"\"\"\n",
    "        blockNum = BlockNum(self.block, q)\n",
    "        blockNum_c = BlockNum(self.block, q.conjugate())\n",
    "        res = [] \n",
    "        for i in range(len(self.polelist)):\n",
    "            pos = self.polelist[i][0]\n",
    "            contribution = (self.CCdiscrete[i] * blockNum.value(pos**2, self.reduced)\n",
    "                            * blockNum_c.value(pos**2, self.reduced))\n",
    "            if self.normalized:\n",
    "                contribution *= self.prefactor\n",
    "            res.append(contribution)\n",
    "        return res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generalized minimal models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Three-point function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The three-point structure constant of a generalized minimal model is given by\n",
    "$$\n",
    "C_{\\{(r_i,s_i)\\}} = \\frac{(-1)^{r_0+s_0+1}b^{2(r_0-s_0)}}{P_{r_0}(b^2)P_{s_0}(b^{-2})Q_{r_0,s_0}(b)} \n",
    "\\prod_{i=1}^3 \\frac{P_{r_i-1}(b^2)P_{s_i-1}(b^{-2}) Q_{r_i-1,s_i-1}(b)}{P_{r_0-r_i}(b^2)P_{s_0-s_i}(b^{-2})Q_{r_0-r_i,s_0-s_i}(b)}\n",
    "$$\n",
    "where we introduced the integers\n",
    "$$\n",
    "r_0 = \\frac{r_1+r_2+r_3-1}{2} \\quad , \\quad s_0 = \\frac{s_1+s_2+s_3-1}{2}\n",
    "$$\n",
    "and the functions \n",
    "$$\n",
    "P_r(x) = \\prod_{i=1}^r \\gamma(1+ix) \\quad ,\\quad Q_{r,s}(b) = \\prod_{i=1}^r \\prod_{j=1}^s(ib+jb^{-1})^2\n",
    "$$\n",
    "where $\\gamma(x) =\\frac{\\Gamma(x)}{\\Gamma(1-x)}$.\n",
    "\n",
    "The conditions for the three fields to be intertwined by fusion are \n",
    "$$\n",
    "r_0, s_0 \\in \\mathbb{N} \\quad , \\quad \\left\\{\\begin{array}{l} r_i \\leq r_0 \\\\ s_i \\leq s_0 \\end{array} \\right.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Parray:\n",
    "    \"\"\" An array of values of P_r(x) for various values of r. \"\"\"\n",
    "    \n",
    "    def __init__(self, x):\n",
    "        \"\"\" We need only specify x. \"\"\" \n",
    "        \n",
    "        self.x = x\n",
    "        self.values = [1]\n",
    "        \n",
    "    def get(self, r):\n",
    "        \"\"\" Given an integer r, lengthens the lists of values if appropriate, and returns\n",
    "        the desired value. \"\"\"\n",
    "        \n",
    "        value = self.values[-1]\n",
    "        \n",
    "        for i in range(len(self.values), r+1): \n",
    "            value *= scipy.special.gamma(1+i*self.x) / scipy.special.gamma(-i*self.x)\n",
    "            self.values.append(value)\n",
    "            \n",
    "        return self.values[r] \n",
    "    \n",
    "class Qarray:\n",
    "    \"\"\" An array of values of Q_{r,s}(x) for various values of r, s. \"\"\"\n",
    "    \n",
    "    def __init__(self, x):\n",
    "        \"\"\" We need only specify x. \"\"\"\n",
    "        \n",
    "        self.x = x\n",
    "        self.values = [[1]]\n",
    "        \n",
    "    def get(self, r, s):\n",
    "        \"\"\" Given r and s, lengthens the list of values as appropriate, and returns the \n",
    "        desired value. \"\"\"\n",
    "        \n",
    "        for i in range(len(self.values)):\n",
    "            line = self.values[i]\n",
    "            value = line[-1]\n",
    "            for j in range(len(line), s+1):\n",
    "                for k in range(1, i+1):\n",
    "                    value *= (k*self.x + j/self.x)**2\n",
    "                line.append(value)\n",
    "            \n",
    "        line = list(self.values[-1])    \n",
    "        \n",
    "        for i in range(len(self.values), r+1):\n",
    "            for j in range(len(line)):\n",
    "                for k in range(1, j+1):\n",
    "                    line[j] *= (i*self.x + k/self.x)**2\n",
    "            self.values.append(list(line))\n",
    "        \n",
    "        return self.values[r][s]\n",
    "    \n",
    "class GMM3:\n",
    "    \"\"\" Three degenerate fields, and the corresponding three-point structure constant. \n",
    "    We also determine whether the fields are intertwined by fusion. That r0 and s0 are \n",
    "    integer is not explicitly tested. \n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, fields, minimal = False, forced = True):\n",
    "        \"\"\" We need an array of three degenerate fields. We should also indicate \n",
    "        whether we really want to compute the three-point function if the fields \n",
    "        are not intertwined by fusion. And we have the option to normalize as in \n",
    "        minimal models.\n",
    "        \"\"\"\n",
    "        \n",
    "        pairs = [field.get('degenerate') for field in fields]\n",
    "        r0 = 0\n",
    "        s0 = 0\n",
    "        for pair in pairs:\n",
    "            (r, s) = pair\n",
    "            r0 += r\n",
    "            s0 += s\n",
    "        r0 = (r0 - 1)//2\n",
    "        s0 = (s0 - 1)//2\n",
    "        \n",
    "        I = True\n",
    "        for pair in pairs:\n",
    "            (r, s) = pair\n",
    "            I = I and r <= r0 and s <= s0\n",
    "        self.canFuse = I            \n",
    "\n",
    "        if I or forced:\n",
    "            \n",
    "            charge = fields[0].charge\n",
    "            b = charge.get('b')\n",
    "            Ps = Parray(b**2)\n",
    "            Pinvs = Parray(1/b**2)\n",
    "            Qs = Qarray(b)\n",
    "            \n",
    "            C = (-1)**(r0 + s0 + 1) * b**(2*(r0-s0)) / Ps.get(r0) / Pinvs.get(s0) / Qs.get(r0, s0)\n",
    "            for pair in pairs:\n",
    "                (r, s) = pair\n",
    "                C *= Ps.get(r-1) * Pinvs.get(s-1) * Qs.get(r-1, s-1)\n",
    "                C /= Ps.get(r0-r) * Pinvs.get(s0-s) * Qs.get(r0-r, s0-s)\n",
    "            if minimal:\n",
    "                C *= cmath.sqrt( Reflection(Dimension('degenerate', [1, 1], charge)).value )\n",
    "                for field in fields:\n",
    "                    C /= cmath.sqrt( Reflection(field).value ) \n",
    "                                \n",
    "            self.value = numerical(C)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Four-point function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GMM4:\n",
    "    \"\"\" Computing the four-point function. \"\"\"\n",
    "    \n",
    "    def __init__(self, fields, t_channel = False,  Nmax = 4, reduced = False):\n",
    "        \"\"\" We need an array of four degenerate fields. \"\"\"\n",
    "        \n",
    "        self.t_channel = t_channel\n",
    "        self.original_fields = fields\n",
    "        self.Nmax = Nmax\n",
    "        self.computeBlock(Nmax)\n",
    "        self.fields = [fields[1], fields[2], fields[3], fields[0]] if t_channel else fields\n",
    "        self.reduced = reduced\n",
    "        \n",
    "        prop = []\n",
    "        for field in self.fields[0].fuse(self.fields[1]):\n",
    "            if GMM3([field, self.fields[2], self.fields[3]], False, False).canFuse:\n",
    "                prop.append(field)\n",
    "        self.propagating = prop\n",
    "\n",
    "    def computeBlock(self, N):\n",
    "        \"\"\" Computing the relevant Block. This method is here so that we can change N without\n",
    "        reconstructing the object.\n",
    "        \"\"\"       \n",
    "        self.block = Block(self.original_fields, N, self.t_channel) \n",
    "        \n",
    "    def value(self, q = .05):\n",
    "        \"\"\" The value of the four-point function. \"\"\"\n",
    "        total = 0\n",
    "        blockNum = BlockNum(self.block, q)\n",
    "        blockNum_c = BlockNum(self.block, q.conjugate())\n",
    "        \n",
    "        for field in self.propagating:\n",
    "            \n",
    "            factor = ( GMM3([field, self.fields[0], self.fields[1]]).value \n",
    "                      * GMM3([field, self.fields[2], self.fields[3]]).value / Reflection(field).value )\n",
    "            delta = field.get('delta')\n",
    "            total += factor * blockNum.value(delta, self.reduced) * blockNum_c.value(delta, self.reduced)\n",
    "        \n",
    "        return numerical(total)\n",
    "    \n",
    "    def getFields(self):\n",
    "        \"\"\" The list of propagating fields. \"\"\"\n",
    "        return [field.get('degenerate') for field in self.propagating]\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conformal field theory with $c<1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Three-point structure constants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute the three-point structure constant\n",
    "$$\n",
    "C(P_1,P_2,P_3) = \n",
    "\\frac{\\prod_{\\pm\\pm}\\Upsilon_\\beta(i(P_1\\pm P_2 \\pm P_3))}{\\prod_{i=1}^3\\Upsilon_\\beta(i(q'+2P_i))}\n",
    "$$\n",
    "where we defined\n",
    "$$\n",
    "\\beta = ib \\quad , \\quad q' = \\frac12(\\beta^{-1}-\\beta)\n",
    "$$\n",
    "In the statistical physics application, the normalized three point structure constants have the extra normalization factor $A(\\beta)\\prod_{i=1}^3 B(\\beta, P_i) $ where \n",
    "$$\n",
    "A(\\beta)=\\beta^{-1-\\beta^2+\\beta^{-2}}\\frac{\\sqrt{\\gamma(\\beta^{-2}-1)\\gamma(\\beta^2)}}{\\Upsilon_{\\beta}(iq')}\n",
    "$$\n",
    "$$\n",
    "B(\\beta, P) = \\sqrt{\\frac{\\Upsilon_\\beta(i(q'+2P))}{\\Upsilon_\\beta(i(q'-2P))}} = \\beta^{4\\beta^{-1}P-\\beta^{-2}-1}\\sqrt{\\gamma(4\\beta^{-1}P) \\gamma(4\\beta^{-1}P -\\beta^{-2})}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Z3:\n",
    "    \"\"\" Computing the three-point function. \"\"\"\n",
    "    \n",
    "    def __init__(self, Ps, beta): \n",
    "        \"\"\" We need a list of three fields. \"\"\"\n",
    "        \n",
    "        self.beta = beta\n",
    "        self.Ps = Ps\n",
    "        self.qp = (1/beta - beta)/ 2\n",
    "        Plist = [[-1j * P, True] for P in Z3.Psums(Ps)] + [[-1j * (self.qp + 2*P), False] for P in Ps]\n",
    "        self.upsilons = Upsilons(beta, Plist)\n",
    "\n",
    "    @staticmethod\n",
    "    def Psums(Ps):\n",
    "        \"\"\" We compute four combinations of momenums. \"\"\"\n",
    "        return [Ps[0] + Ps[1] + Ps[2], Ps[0] - Ps[1] + Ps[2], Ps[0] + Ps[1] - Ps[2], - Ps[0] + Ps[1] + Ps[2]]\n",
    "        \n",
    "    def upsilons(self):\n",
    "        \"\"\" We return an Upsilon object, which can then be evaluated in various ways. \"\"\"\n",
    "        return self.upsilons\n",
    "    \n",
    "    @staticmethod\n",
    "    def overall_factor(beta):\n",
    "        \"\"\" We compute the factor A(beta). \"\"\"\n",
    "        result = 1   \n",
    "        result *= Upsilons(beta, [[1j*(1/beta - beta)/2, False]]).integral()\n",
    "        result *= beta**(beta**(-2) - beta**2 -1)\n",
    "        result *= cmath.sqrt(mpmath.gamma(beta**2) / mpmath.gamma(1 - beta**2) * mpmath.gamma(beta**(-2) - 1)\n",
    "                             / mpmath.gamma(2 - beta**(-2)))        \n",
    "        return result\n",
    "    \n",
    "    @staticmethod\n",
    "    def momentum_factor(beta, P):\n",
    "        \"\"\" We compute the momentum-dependent normalization. \"\"\"\n",
    "        (p1, p2) = (4*P/beta, 4*P/beta - 1/beta**2)\n",
    "        result = beta**(p1 - 1/beta**2 - 1)\n",
    "        result *= cmath.sqrt(mpmath.gamma(p1) / mpmath.gamma(1 - p1) * mpmath.gamma(p2) / mpmath.gamma(1 - p2))\n",
    "        return result\n",
    "    \n",
    "    def normalized_value(self):\n",
    "        \"\"\" Computes the value of the three-point function, as normalized for statistical physics\n",
    "        applications.\n",
    "        \"\"\"\n",
    "        result = self.upsilons.integral()\n",
    "        result *= Z3.overall_factor(self.beta)\n",
    "        for P in self.Ps:\n",
    "            result *= Z3.momentum_factor(self.beta, P)\n",
    "        return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Four-point functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute four-point functions with $c<1$ i.e. $\\beta\\in \\mathbb{R}$. In the $s$-channel, this is \n",
    "$$\n",
    "Z = \\int_{\\mathcal{C_\\epsilon}} dP\\ C(P_1, P_2, P) C(-P, P_3, P_4)\\ |B_s(P)|^2\n",
    "$$\n",
    "for a contour $\\mathcal{C}_\\epsilon = \\mathbb{R} + i\\epsilon$.\n",
    "\n",
    "By default we neglect $P_i$-dependent prefactors, and use the product formula for Upsilon functions. \n",
    "To obtain the statistical physics normalization, we should multiply with the normalization factor\n",
    "$$\n",
    "\\frac{1}{\\prod_{i=1}^4 \\Upsilon_\\beta(i(q'+2P_i))} \\cdot A(\\beta)^2 \\prod_{i=1}^4 B(\\beta, P_i) \\cdot e^{\\lambda_b(6q'^2 + 8q'\\sum_{i=1}^4 P_i)}\n",
    "$$\n",
    "where the first factor contains the neglected $P_i$-dependent prefactors, the second factor brings us to the statistical physics normalization, and the third factor accounts for our use of products instead of integrals for the Upsilon functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Z4:\n",
    "    \"\"\" Computing the four-point function. \"\"\"\n",
    "    \n",
    "    def __init__(self, t_channel = False, Ps = [i + .5 for i in range(4)], beta = math.sqrt(2) - 1, \n",
    "                 Nmax = 4, Pcutoff = 10, Ucutoff = 20, reduced = False, spline = 0, epsilon = .1, \n",
    "                 normalized = False):\n",
    "        \"\"\" We need the values of beta, four momentums P, the order Nmax of the blocks' expansion, \n",
    "        cutoffs for the P integration and the Upsilon product formula, and the indication whether to \n",
    "        use reduced blocks. \n",
    "        We also give the number of points in the interpolating spline: if this number is zero, then no \n",
    "        interpolation is used. For precision computations the spline should be about 2000. \n",
    "        For graphical display, much lower values such as 100 may be enough.\n",
    "        We also give a parameter epsilon for the position of the integration line.\n",
    "        We indicate whether the result should be normalized for statistical physics applications.\n",
    "        \"\"\"  \n",
    "        \n",
    "        self.beta = beta\n",
    "        self.qp = (1/beta - beta)/2\n",
    "        self.t_channel = t_channel\n",
    "        self.Ucutoff = Ucutoff\n",
    "        self.Pcutoff = Pcutoff\n",
    "        self.charge = Charge('beta', beta)        \n",
    "        self.Ps = Ps\n",
    "        self.Nmax = Nmax\n",
    "        self.computeBlock(Nmax)\n",
    "        self.reduced = reduced\n",
    "        self.epsilon = epsilon\n",
    "        \n",
    "        self.useSpline = spline > 0\n",
    "        if self.useSpline:\n",
    "            Pvals = [i*Pcutoff/spline for i in range(-spline, spline+1)]\n",
    "            self.Kspline = FourPoint.complex_spline(Pvals, self.K)   \n",
    "            \n",
    "        self.normalized = normalized\n",
    "        if normalized:\n",
    "            upsilons = [[-1j * (self.qp + 2*P), False] for P in Ps]\n",
    "            factor = Upsilons(beta, upsilons).product(Ucutoff)\n",
    "            factor *= Z3.overall_factor(beta)**2\n",
    "            for P in Ps:\n",
    "                factor *= Z3.momentum_factor(beta, P)\n",
    "            factor *= cmath.exp(Upsilons.lambda_b(beta) * (6 * self.qp**2 + 8 * self.qp * sum(Ps)))\n",
    "            self.factor = factor                \n",
    "\n",
    "    def computeBlock(self, N):\n",
    "        \"\"\" Computing the relevant Block. This method is here so that we can change N without\n",
    "        reconstructing the object.\n",
    "        \"\"\"       \n",
    "        self.block = Block([Dimension('P', self.Ps[i], self.charge) for i in range(4)], N, self.t_channel)\n",
    "        \n",
    "    def K(self, p0):\n",
    "        \"\"\" We compute the prefactor. We need to shift p at this point, so that we can compute the spline. \"\"\"\n",
    "        Ps = self.Ps\n",
    "        p = self.shift(p0)\n",
    "        if self.t_channel:\n",
    "            Ps = [Ps[1], Ps[2], Ps[3], Ps[0]]\n",
    "        Psums = Z3.Psums([Ps[0], Ps[1], p]) + Z3.Psums([Ps[2], Ps[3], -p])\n",
    "        Plist = [[-1j * P, True] for P in Psums] + [[-1j * (self.qp + 2*P), False] for P in [p, -p]]\n",
    "        upsilons = Upsilons(self.beta, Plist)\n",
    "        # return upsilons.integral()\n",
    "        return upsilons.product(self.Ucutoff)\n",
    "    \n",
    "    def shift(self, P):\n",
    "        \"\"\" Shifts P away from the real line. \"\"\"\n",
    "        return P + 1j * self.epsilon\n",
    "    \n",
    "    def value(self, q = .05):\n",
    "        \"\"\" Computes the numerical value of the four-point function. The integrand decreases quickly at large\n",
    "        P so the cutoff can be quite small. \n",
    "        \"\"\"\n",
    "        \n",
    "        blockNum = BlockNum(self.block, q)\n",
    "        blockNum_c = BlockNum(self.block, q.conjugate())\n",
    "        \n",
    "        def integrand(P):\n",
    "            if self.useSpline:\n",
    "                K = self.Kspline[0](P) + self.Kspline[1](P)*1j\n",
    "            else:\n",
    "                K = self.K(P)\n",
    "            Delta = self.shift(P)**2\n",
    "            return numerical(K * blockNum.value(Delta, self.reduced) * blockNum_c.value(Delta, self.reduced))\n",
    "            \"\"\" We apply the method 'numerical' because 'mpmath' can produce objects of the \n",
    "            'mpc' and 'mpf' classes which are not recognized as complex numbers and floats \n",
    "            respectively by the 'quad' integration method.\n",
    "            \"\"\"\n",
    "        result = complex_quadrature(integrand, -self.Pcutoff, self.Pcutoff)[0]      \n",
    "        if self.normalized:\n",
    "            result *= self.factor\n",
    "        return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rational central charge"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us have \n",
    "$$\n",
    "\\beta^2 = \\frac{q}{p}\n",
    "$$\n",
    "We consider the three-point structure constants \n",
    "$$\n",
    "\\tilde{C}(P_1,P_2,P_3) = C(P_1, P_2, P_3) \\sigma\n",
    "$$\n",
    "where $C$ is the usual, smooth three-point structure constant, and\n",
    "$$\n",
    "\\sigma = 1_E =\n",
    "\\left\\{\n",
    "\\begin{array}{l} \n",
    "1 \\quad \\text{if} \\quad \n",
    "s(P_1+P_2+P_3)s(P_1+P_2-P_3)s(P_1-P_2+P_3)s(-P_1+P_2+P_3) > 0\n",
    "\\\\\n",
    "0 \\quad \\text{otherwise}\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "with \n",
    "$$\n",
    "s(x) = \\sin \\pi \\left(\\sqrt{pq} x + \\tfrac12(p-q)\\right)\n",
    "$$\n",
    "The periodic function $\\sigma$ is the characteristic function of the set $E$, which is a union of intervals. In a four-point function we have a factor $\\sigma\\sigma' = 1_{E\\cap E'}$, and we call the corresponding integral the intersection integral. We also consider the xor integral with the factor $1_{E\\sqcup E'}$, which would appear if we considered three-point function of the type $C(P_1, P_2, P_3)(1-\\nu + \\nu \\sigma)$ for some number $\\nu$. \n",
    "\n",
    "Let us examine how we compute $\\sigma$ as a function of $P_1$ in the case $p=q=1$. We call $P_-, P_+$ the two smallest numbers of the set of four numbers $\\{\\pm P_2 \\pm P_3 \\operatorname{mod} 1\\} \\subset (0,1)$. Then, as a function of $P_1$, $\\sigma$ obeys\n",
    "$$\n",
    "P_1\\in (0, 1) \\ \\text{and} \\ \\sigma(P_1) = 1 \\quad \\Leftrightarrow \\quad P_1 \\in (P_-, P_+) \\cup (1-P_+, 1-P_-)\n",
    "$$\n",
    "Then let us consider a product of two three-point functions, with intervals $(P_-, P_+)$ and $(P'_-, P'_+)$. Let us order these numbers so that $\\{P_-, P_+, P'_-, P'_+\\} = \\{a, b, c, d\\}$ with $a<b<c<d$. Then\n",
    "$$\n",
    "P_1\\in (0, 1) \\ \\text{and} \\ \\sigma(P_1) \\neq \\sigma'(P_1)\\quad \\Leftrightarrow \\quad P_1 \\in (a, b) \\cup (c, d) \\cup (1-d, 1-c) \\cup (1-b, 1-a)\n",
    "$$\n",
    "and \n",
    "$$\n",
    "P_1 \\in (0, 1) \\ \\text{and} \\ \\sigma(P_1)=\\sigma'(P_1) = 1 \\quad \\Leftrightarrow \\quad P_1 \\in \n",
    "\\left\\{\\begin{array}{ll} \n",
    "\\emptyset & \\text{if}\\ \\  \\operatorname{sign}(P_\\pm - P'_\\pm) = \\text{constant}\n",
    "\\\\ \n",
    "(b, c) \\cup (1-c, 1-b) & \\text{else} \\end{array}\\right.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Four-point function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class R4:\n",
    "    \"\"\" Computing the four-point function. \"\"\"\n",
    "    \n",
    "    def __init__(self, t_channel = False, Ps = [1/4.37, 1/11.73, 1/3.81, 1/12.46], \n",
    "                 charge = Charge('minimal', [1, 2]), eta = .00001274, Nmax = 4, \n",
    "                 Pcutoff = 2.25, Ucutoff = 20, reduced = False, epsilon = .1, barnes = False, disjoint = False):\n",
    "        \"\"\" We need a minimal central charge, four momentums P, the order Nmax of the blocks' expansion, \n",
    "        cutoffs for the P integration and the Upsilon product formula, and the indication whether to \n",
    "        use reduced blocks. \n",
    "        We also give a parameter epsilon for the position of the integration line, and\n",
    "        the deviation eta of the central charge from its rational value. \n",
    "        If c = 1 we can use Barnes' G-function for computing structure constants, but this is slower.\n",
    "        We indicate whether we want to compute the integral on an intersection (default) or disjoint union.\n",
    "        \"\"\"  \n",
    "        \n",
    "        self.beta = charge.get('beta')\n",
    "        self.qp = (1/self.beta - self.beta)/ 2\n",
    "        self.charge = Charge('beta', self.beta + 1j * eta)\n",
    "        (self.p, self.q) = charge.get('minimal')\n",
    "        self.even = self.p - self.q % 2 == 0\n",
    "        self.barnes = self.p == self.q and barnes\n",
    "        self.factor = 1 / math.sqrt(self.p * self.q)\n",
    "        self.t_channel = t_channel\n",
    "        self.Ucutoff = Ucutoff\n",
    "        self.Pcutoff = Pcutoff       \n",
    "        self.Ps = Ps\n",
    "        self.Nmax = Nmax\n",
    "        self.computeBlock(Nmax)\n",
    "        self.reduced = reduced\n",
    "        self.epsilon = epsilon\n",
    "        self.disjoint = disjoint\n",
    "        \n",
    "    def computeBlock(self, N):\n",
    "        \"\"\" Computing the relevant Block. This method is here so that we can change N without\n",
    "        reconstructing the object.\n",
    "        \"\"\"       \n",
    "        self.block = Block([Dimension('P', self.Ps[i], self.charge) for i in range(4)], N, self.t_channel)\n",
    "    \n",
    "    @staticmethod \n",
    "    def upsilon(P):\n",
    "        return mpmath.barnesg(1 + P) * mpmath.barnesg(1 - P)\n",
    "    \n",
    "    def Ksmooth(self, P):\n",
    "        \"\"\" We compute the analytic factor of the three-point function. \"\"\"\n",
    "        Ps = self.Ps\n",
    "        if self.t_channel:\n",
    "            Ps = [Ps[1], Ps[2], Ps[3], Ps[0]]\n",
    "        Psums = Z3.Psums([Ps[0], Ps[1], P]) + Z3.Psums([Ps[2], Ps[3], -P])\n",
    "        if self.barnes:\n",
    "            result = 1 \n",
    "            for P1 in Psums:\n",
    "                result *= R4.upsilon(P1)\n",
    "            result /= R4.upsilon(2*P)**2\n",
    "            return result\n",
    "        Plist = [[-1j * P1, True] for P1 in Psums] + [[-1j * (self.qp + 2*P1), False] for P1 in [P, -P]]\n",
    "        upsilons = Upsilons(self.beta, Plist)  \n",
    "        return upsilons.product(self.Ucutoff)\n",
    "    \n",
    "    @staticmethod\n",
    "    def interval(P2, P3, factor = 1):\n",
    "        \"\"\" We compute the first interval where the non-analytic factor does not vanish. \"\"\"\n",
    "        (P2, P3) = (P2.real, P3.real)\n",
    "        return list(numpy.sort(map(lambda x: x % factor, [P2 + P3, P2 - P3, P3 - P2, - P2 - P3]))[0:2])\n",
    "    \n",
    "    @staticmethod\n",
    "    def intersections(interval1, interval2, display = False):\n",
    "        \"\"\" Computes the xor and intersection of the two intervals, as a list of four bounds,\n",
    "        plus two boolean which are True if the xor and intersections are respectively empty. \n",
    "        \"\"\"\n",
    "        diffs = map(lambda x: x > 0, [i1 - i2 for i1 in interval1 for i2 in interval2])\n",
    "        bounds = interval1 + interval2\n",
    "        bounds = numpy.sort(bounds)\n",
    "        if display:\n",
    "            R4.display([interval1, interval2, bounds], sort = False)        \n",
    "        empty_intersection = diffs[0] == diffs[1] and diffs[0] == diffs[2] and diffs[0] == diffs[3]\n",
    "        empty_xor = bounds[0] == bounds[1] and bounds[2] == bounds[3]\n",
    "        return [bounds, empty_xor, empty_intersection]\n",
    "    \n",
    "    @staticmethod\n",
    "    def extend(interval_list, period, interval):\n",
    "        \"\"\" Extends interval_list into interval modulo period and reflection. \"\"\"\n",
    "        reflected = [[period - inter[1], period - inter[0]] for inter in interval_list]\n",
    "        shifts = range(int(interval[0] / period), int(interval[1] / period))\n",
    "        return [[inter[0] + i*period, inter[1] + i*period] \n",
    "                for inter in interval_list + reflected for i in shifts]\n",
    "    \n",
    "    @staticmethod\n",
    "    def display(interval_list, sort = True):\n",
    "        \"\"\" Orders the interval list, truncates the bounds, and prints. \"\"\"\n",
    "        if sort: \n",
    "            interval_list.sort(key = lambda interval: interval[0])\n",
    "        print( [[round(bound, 3) for bound in interval] for interval in interval_list] )\n",
    "    \n",
    "    def value(self, q = .02):\n",
    "        \"\"\" Computes the numerical value of the four-point function. The non-smooth terms are computed\n",
    "        as sums of integrals over the intervals where they do not vanish.\n",
    "        \"\"\"\n",
    "        \n",
    "        blockNum = BlockNum(self.block, q)\n",
    "        blockNum_c = BlockNum(self.block, q.conjugate())\n",
    "        \n",
    "        Ps = self.Ps\n",
    "        period = self.factor\n",
    "        if self.t_channel:\n",
    "            Ps = [Ps[1], Ps[2], Ps[3], Ps[0]]\n",
    "        if not self.even:\n",
    "            Ps = [Ps[0] + period / 2, Ps[1], Ps[2], Ps[3] + period / 2]\n",
    "        inters = R4.intersections(R4.interval(Ps[0], Ps[1], period), R4.interval(Ps[2], Ps[3], period))  \n",
    "        bounds = inters[0]\n",
    "        interval = [0, self.Pcutoff]\n",
    "        xor = [] if inters[1] else R4.extend([[bounds[0], bounds[1]], [bounds[2], bounds[3]]], period, interval)\n",
    "        intersection = [] if inters[2] else R4.extend([[bounds[1], bounds[2]]], period, interval)\n",
    "        \n",
    "        def integrand(P):\n",
    "            Delta = P**2 \n",
    "            return numerical(self.Ksmooth(P) * blockNum.value(Delta, self.reduced) \n",
    "                           * blockNum_c.value(Delta, self.reduced))\n",
    "        \n",
    "        intervals = xor if self.disjoint else intersection\n",
    "        # self.xor_int = 2*sum([complex_quadrature(integrand, inter[0], inter[1])[0] for inter in xor])\n",
    "        # self.intersection_int = 2*sum([complex_quadrature(integrand, inter[0], inter[1])[0] for inter in intersection])\n",
    "        # self.smooth_int = complex_quadrature(lambda P: integrand(P + 1j * self.epsilon), -self.Pcutoff, self.Pcutoff)[0]\n",
    "        return 2*sum([complex_quadrature(integrand, inter[0], inter[1])[0] for inter in intervals])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" Test of crossing symmetry for Liouville for c \\notin (-\\infty,1)\"\"\"\n",
    "b = .31 + 0.32j\n",
    "Q = b + 1. / b\n",
    "Ps = [0.2, 0.123, 0.3,0.17]\n",
    "crossratio =  0.32\n",
    "sch = FourPoint(b = b, Ps = Ps, Nmax = 8, normalized = True)\n",
    "tch = FourPoint(b = b, Ps = Ps, t_channel=True,  Nmax = 8, normalized = True)\n",
    "print( \"crossing sym.\", sch.value(mpmath.qfrom(m=crossratio)), tch.value(mpmath.qfrom(m= crossratio)) )\n",
    "#ls = sch.value_discrete() \n",
    "#print( \"discrete contributions\", ls )\n",
    "#polelist = sch.getcrossedpoles()\n",
    "#sch.getCCdiscrete(polelist)\n",
    "#tch = FourPoint(b = b, Ps = Ps, t_channel=True,  Nmax = 8, normalized = True)\n",
    "#sch.value(mpmath.qfrom(m=.5)), tch.value(mpmath.qfrom(m=.5))\n",
    "#ls = sch.value_discrete() \n",
    "#print( \"discrete contributions\", ls )\n",
    "#import matplotlib.pyplot as plt\n",
    "#poles = sch.getallpoles(3,3)\n",
    "#print( poles )\n",
    "#allpoles = poles.flatten()\n",
    "#for j in range(16):\n",
    "#    plt.plot(poles[j,:,:].flatten().real, poles[j,:,:].flatten().imag, \"+-\")\n",
    "#plt.plot([-1, 2], [0,0], color = \"black\")\n",
    "#plt.plot([-1, 8], [Q / 2,0], color = \"black\")\n",
    "#plt.show()\n",
    "#sch.getcrossedpoles()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Edit Metadata",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}