{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Factorization into q-integers\n",
    "\n",
    "In many cases we have a nice product formula and its q-analog.\n",
    "For example, the number of standard Young tableaux of shape $\\lambda\\vdash n$ is\n",
    "$$\n",
    "f^\\lambda = \\frac{n!}{\\prod_{u\\in \\lambda} h(u)},\n",
    "$$\n",
    "where $h(u)$ is the hook length of the cell $u$ in the Young diagram $\\lambda$. Its $q$-analog is\n",
    "$$\n",
    "\\sum_{T\\in SYT(\\lambda)} q^{maj(T)} = \\frac{[n]_q!}{\\prod_{u\\in \\lambda} [h(u)]_q},\n",
    "$$\n",
    "\n",
    "\n",
    "In this Sage Notebook we factor a rational function in $q$ into $q$-integers $[k]_q = q^k-1$. (Note that we use a rather different convention for $q$-integers.)\n",
    "\n",
    "In the code we use cyclotomic polynomials.\n",
    "Recall that the cyclotomic polynomials $\\Phi_n(x)$ satisfy\n",
    "$$\n",
    "x^N - 1 = \\prod_{d|N} \\Phi_d(x).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "max_cyclotomic_index = 100   \n",
    "var('q','x','y','z')\n",
    "var('a','b','c','d','e','k','m')\n",
    "\n",
    "def find_cyclotomic_index(p,x,N=max_cyclotomic_index):\n",
    "    r\"\"\"\n",
    "    Return n if p is the nth cyclotomic polynomial for some 1<=n<=N and 0 otherwise.\n",
    "        \n",
    "    EXAMPLES::\n",
    "\n",
    "        sage: find_cyclotomic_index(1+q^4,q)\n",
    "        8\n",
    "        sage: find_cyclotomic_index(1-q,q)\n",
    "        0\n",
    "        sage: find_cyclotomic_index(q-1,q)\n",
    "        1\n",
    "\n",
    "    \"\"\"\n",
    "    for i in range(1,N+1):\n",
    "        if p == cyclotomic_polynomial(i,x):\n",
    "            return i\n",
    "    return 0\n",
    "\n",
    "def find_qint(f,q):\n",
    "    \"\"\"\n",
    "    Return x if f = q^x-1 and 0 otherwise\n",
    "    \"\"\"\n",
    "    if f == q-1:\n",
    "        return 1\n",
    "    if (f+1).operator() == (x^y).operator():\n",
    "        A = (f+1).operands()\n",
    "        if A[0] == q:\n",
    "            return A[1]\n",
    "    return 0\n",
    "\n",
    "def qint_expression(f,q,N=max_cyclotomic_index):\n",
    "    r\"\"\"\n",
    "    Return [const,qnum,qden] where if F = C*NUM/DEN, qnum and qden are the lists of factors in NUM and DEN which are products of factors of the form (1-q^k).\n",
    "        \n",
    "    EXAMPLES::\n",
    "\n",
    "        sage: qint_expression(-q^2*(1-q)*(1-q^6)^2/(1-q^4)/(1-q^18),q)\n",
    "        [q^2, [6, 6, 1], [18, 4]]\n",
    "        sage: qint_expression(-q^2*(1+q)^2/(1-q^2+q^4),q)\n",
    "        [-q^2, [6, 4, 2], [12]]\n",
    "        sage: qint_expression(-q^2*(1+q)^2*(1-q^6)/(1-q^4),q)\n",
    "        [-q^2, [6, 2, 2], [4]]\n",
    "\n",
    "    \"\"\"\n",
    "    const,qnum,qden,num,den = 1,[],[],[],[]\n",
    "    for A in f.factor_list():\n",
    "        k = find_qint(A[0],q)\n",
    "        if k != 0:\n",
    "            if A[1] > 0:\n",
    "                qnum += [k for i in range(A[1])]\n",
    "            else:\n",
    "                qden += [k for i in range(-A[1])]\n",
    "        else:\n",
    "            n = find_cyclotomic_index(A[0],q,N)\n",
    "            if n == 0:\n",
    "                const *= A[0]^A[1]\n",
    "            elif A[1] > 0:\n",
    "                num += [n for i in range(A[1])]\n",
    "            else:\n",
    "                den += [n for i in range(-A[1])]\n",
    "    while len(num)+len(den) > 0:\n",
    "        num.sort()\n",
    "        den.sort()\n",
    "        if num == []:\n",
    "            k = den.pop()\n",
    "            qden.append(k)\n",
    "            num += divisors(k)[:-1]\n",
    "        elif den == []:\n",
    "            k = num.pop()\n",
    "            qnum.append(k)\n",
    "            den += divisors(k)[:-1]\n",
    "        elif num[-1] == den[-1]:\n",
    "            num.pop()\n",
    "            den.pop()\n",
    "        elif num[-1] < den[-1]:\n",
    "            k = den.pop()\n",
    "            qden.append(k)\n",
    "            num += divisors(k)[:-1]\n",
    "        elif num[-1] > den[-1]:\n",
    "            k = num.pop()\n",
    "            qnum.append(k)\n",
    "            den += divisors(k)[:-1]\n",
    "    while 1 in qden and 1 in qnum:\n",
    "        qden.remove(1)\n",
    "        qnum.remove(1)\n",
    "    qden.sort(reverse=True)\n",
    "    qnum.sort(reverse=True)\n",
    "    return [const,qnum,qden]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(q^6 + q^5 + q^4 + q^3 + q^2 + q + 1)*(q^4 + q^3 + q^2 + q + 1)*(q^4 + 1)*q^5"
      ]
     },
     "execution_count": 121,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var('q')\n",
    "A=add(q^(T.standard_major_index()) for T in StandardTableaux([4,3,1]))\n",
    "factor(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above example, it is not easy to see directly the q-integer factors. More obvious example would be $[n]_q!$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(q^6 + q^5 + q^4 + q^3 + q^2 + q + 1)*(q^6 + q^3 + 1)*(q^4 + q^3 + q^2 + q + 1)*(q^4 + 1)*(q^2 + q + 1)^3*(q^2 - q + 1)*(q^2 + 1)^2*(q + 1)^4*(q - 1)^9"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = factor(mul(q^k-1 for k in range(1,10))); B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we used \"qint_expression\" method we can see the q-integer factors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[q^5, [8, 7, 5], [4, 1, 1]]"
      ]
     },
     "execution_count": 122,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qint_expression(A,q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above output means that\n",
    "$$\n",
    "A = q^5 \\frac{[8]_q[7]_q[5]_q}{[4]_q[1]_q[1]_q}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, [9, 8, 7, 6, 5, 4, 3, 2, 1], []]"
      ]
     },
     "execution_count": 124,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qint_expression(B,q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 8.7",
   "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
}