Sum of powers.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"text-align: right\">Paul Novaes<br>June 2018</div> \n",
    "\n",
    "# Sum of Powers\n",
    "\n",
    "We show how to find closed formulas for $$1^k+2^k +\\cdots + n^k$$\n",
    "using a simple program and elementary linear algebra.\n",
    "\n",
    "We then compare our results to those in Bernouilli's  [Ars Conjectandi](https://en.wikipedia.org/wiki/Ars_Conjectandi) and find a typo in [its formula](https://books.google.com/books?id=kD4PAAAAQAAJ&pg=PA97#v=onepage&q&f=false) for $\\int n^9$.\n",
    "\n",
    "![Wikipedia](https://upload.wikimedia.org/wikipedia/commons/7/74/JakobBernoulliSummaePotestatum.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finding a closed formula\n",
    "\n",
    "The goal is to find close formulas for the sum of powers, such as $1+2 +\\cdots + n = {{n(n+1)}\\over 2}$. Of course we could look them up, but it would be more fun to find them by ourselves, rather than to deal with Bernouilli numbers and what not. \n",
    "\n",
    "Let's define $$S_k(n) = 1^k+2^k +\\cdots + n^k$$\n",
    "\n",
    "We will assume that $S_k(n)$ is a polynomial of degree $k+1$. Note that this same assuption had to be made by Bernouilli too, since his formulas were not proven until Jacobi.\n",
    "\n",
    "Our goal is to find a polynomial $$P_k(n) = a_{k+1}n^{k+1} + a_kn^k + \\cdots a_0n^0$$ such that, for all $n \\ge 0$, $P_k(n) = S_k(n)$.\n",
    "\n",
    "Since a polynomial of degree $d$ is uniquely determined by its value on $d+1$ points, we just need to find $P_k$ such that $P_k(n) = S_k(n)$ for $n = 0, 1, \\ldots, k+1$. Or, in matrix form:\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{pmatrix} 0^{k+1} & 0^k & \\cdots & 0^0 \\\\ 1^{k+1} & 1^k & \\cdots & 1^0 \\\\ \\vdots & \\vdots & & \\vdots \\\\ (k+1)^{k+1} & (k+1)^k & \\cdots & (k+1)^0 \\end{pmatrix}\n",
    "\\begin{pmatrix} a_{k+1} \\\\ a_k \\\\ \\vdots \\\\ a_0 \\end{pmatrix}\n",
    "= \\begin{pmatrix} 0 \\\\ 1^k \\\\ \\vdots \\\\ 1^k + 2^k + \\cdots + (k+1)^k\\end{pmatrix} \n",
    "\\end{equation}\n",
    "\n",
    "So we have reduced our problem to solving a system of linear equations of the form $AX = B$, where $A$ and $B$ are known and $X$ (the coefficients of the polynomial) is unknown.\n",
    "\n",
    "Note that in the formula above, we follow Python's convention: $0^0 = 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using sympy to solve the system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "\n",
    "def MatrixA(k):\n",
    "    A = [[0 for x in range(k + 2)] for y in range(k + 2)]\n",
    "    for i in range(k + 2):\n",
    "        for j in range(k + 2):\n",
    "            A[i][j] = i ** (k - j + 1)\n",
    "    return Matrix(A)\n",
    "\n",
    "def MatrixB(k):\n",
    "    B = [0 for x in range(k + 2)]\n",
    "    for i in range(1, k + 2):\n",
    "        B[i] = B[i - 1] + i ** k \n",
    "    return Matrix(B)\n",
    "\n",
    "def formula_sum_of_powers(k):\n",
    "    A = MatrixA(k)\n",
    "    B = MatrixB(k)\n",
    "    X = A.LUsolve(B)\n",
    "    return X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Verifying the formula in simple cases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([\n",
       "[1/2],\n",
       "[1/2],\n",
       "[  0]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula_sum_of_powers(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is, we have the familiar formula $$1+2+\\cdots+n = {{n^2}\\over 2} + {n\\over 2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([\n",
       "[1/3],\n",
       "[1/2],\n",
       "[1/6],\n",
       "[  0]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula_sum_of_powers(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is, $$1^2+2^2+\\cdots+n^2 = {{n^3}\\over3} + {{n^2}\\over 2} + {n\\over 6}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finding a formula for $1^{100} + 2^{100} + \\cdots + n^{100}$\n",
    "\n",
    "Using our program, it is a breeze to show that $1^{100} + 2^{100} + \\cdots + n^{100}$ can be expressed as\n",
    "\n",
    "$${1\\over 101}n^{101} + {1\\over 2}n^{100} + {25\\over 3}n^{99} + \\cdots -{94598037819122125295227433069493721872702841533066936133385696204311395415197247711\\over 33330}n$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([\n",
       "[                                                                                     1/101],\n",
       "[                                                                                       1/2],\n",
       "[                                                                                      25/3],\n",
       "[                                                                                         0],\n",
       "[                                                                                   -2695/2],\n",
       "[                                                                                         0],\n",
       "[                                                                                    298760],\n",
       "[                                                                                         0],\n",
       "[                                                                                 -66698170],\n",
       "[                                                                                         0],\n",
       "[                                                                             43232541100/3],\n",
       "[                                                                                         0],\n",
       "[                                                                            -2987368590010],\n",
       "[                                                                                         0],\n",
       "[                                                                           592545208316600],\n",
       "[                                                                                         0],\n",
       "[                                                                   -1909010939318560231/17],\n",
       "[                                                                                         0],\n",
       "[                                                                    60927624576260699950/3],\n",
       "[                                                                                         0],\n",
       "[                                                                -10503770178403996919537/3],\n",
       "[                                                                                         0],\n",
       "[                                                                  574696979476592789539800],\n",
       "[                                                                                         0],\n",
       "[                                                               -89701724868851868404757210],\n",
       "[                                                                                         0],\n",
       "[                                                             13296745470530926863904296852],\n",
       "[                                                                                         0],\n",
       "[                                                          -1869298239768618416234153813290],\n",
       "[                                                                                         0],\n",
       "[                                                         248870955751990847260270884407400],\n",
       "[                                                                                         0],\n",
       "[                                                 -1065245686771269279784908613651828005/34],\n",
       "[                                                                                         0],\n",
       "[                                                     3723652407297582727619274890591931075],\n",
       "[                                                                                         0],\n",
       "[                                             -10844299000116828980379757772973769420469/26],\n",
       "[                                                                                         0],\n",
       "[                                                 43950288418050613210495571589828389262800],\n",
       "[                                                                                         0],\n",
       "[                                              -4348447505694585428839166185138223415249684],\n",
       "[                                                                                         0],\n",
       "[                                             403139711179170251736670257248480111641926600],\n",
       "[                                                                                         0],\n",
       "[                                          -34944260063316672143127900206016265956506516820],\n",
       "[                                                                                         0],\n",
       "[                                         2825393845314372316887963516463774302220183874480],\n",
       "[                                                                                         0],\n",
       "[                                -46975128963737486419489164499794297560673231041202090/221],\n",
       "[                                                                                         0],\n",
       "[                                     14838677083702274079364344314958348779821717660137900],\n",
       "[                                                                                         0],\n",
       "[                                   -958463607702055700662952442255558159708363131670845130],\n",
       "[                                                                                         0],\n",
       "[                                  57102248319760803358389861613269821873306861123018474800],\n",
       "[                                                                                         0],\n",
       "[                               -3127153223428501512572222954783017354698250212191091945572],\n",
       "[                                                                                         0],\n",
       "[                              156839198684220220595062954768424384099662857573873761392200],\n",
       "[                                                                                         0],\n",
       "[                        -93273006623793637434656479802977293641893710056414115793830132/13],\n",
       "[                                                                                         0],\n",
       "[                           298055222117767447988694875776788702575308931452828672542296400],\n",
       "[                                                                                         0],\n",
       "[                    -380420681562789081339436627697748498619486609696130138245054547645/34],\n",
       "[                                                                                         0],\n",
       "[                        377511069257143967197314136886615170865786408764196373626268649635],\n",
       "[                                                                                         0],\n",
       "[                   -22758671683254934243234770245768111655371809025564559292966948184145/2],\n",
       "[                                                                                         0],\n",
       "[                     304383493005866429515905139920856508495517057169598168611307541714600],\n",
       "[                                                                                         0],\n",
       "[               -93215398532963113031284566771666289746833192047887884379325411297276490/13],\n",
       "[                                                                                         0],\n",
       "[                  147482537396120605270641214092320919218664984298368495230915199366025300],\n",
       "[                                                                                         0],\n",
       "[             -13112574861745373977119865065563790341936971501269597493806710307275562234/5],\n",
       "[                                                                                         0],\n",
       "[                39857448167724145521712375634202721882854746999642509687138048831663608600],\n",
       "[                                                                                         0],\n",
       "[          -8684587426344073606489504709883050170678418405635672010429984157277666084023/17],\n",
       "[                                                                                         0],\n",
       "[           16304634906294224258925539055711850026937478314146500355475130833028845746350/3],\n",
       "[                                                                                         0],\n",
       "[        -612067818251686839120746668057583204003010059302601007211600982133792394509295/13],\n",
       "[                                                                                         0],\n",
       "[            324388433491984944852419075920450620931865462636226130063940688368734730163320],\n",
       "[                                                                                         0],\n",
       "[          -1725539552167813151807602972961047761861449443882732089672204732211086637590530],\n",
       "[                                                                                         0],\n",
       "[      225010984573490896048358130480328881791747873028668541553377868767760647319540900/33],\n",
       "[                                                                                         0],\n",
       "[       -56995948223784756021697384698478769187481110881793715780290371409596995575985270/3],\n",
       "[                                                                                         0],\n",
       "[          34649381621107623485288357725441359210991177541965976048926944361209977534643400],\n",
       "[                                                                                         0],\n",
       "[  -16293234618989521508515025064456465992824384487957638029599182473343901462949018943/442],\n",
       "[                                                                                         0],\n",
       "[          18674771685049011296614057325260991542019103825338734064995339343748060441015725],\n",
       "[                                                                                         0],\n",
       "[-94598037819122125295227433069493721872702841533066936133385696204311395415197247711/33330],\n",
       "[                                                                                         0]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula_sum_of_powers(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A typo in _Ars Conjectandi_\n",
    "\n",
    "In _Ars Conjectandi_, it says, for $\\int n^9$ (Bernouilli's notation for $\\sum_{i=1}^n{i^9}$):\n",
    "\n",
    "$$\\int n^9 = {1\\over 10}n^{10} + {1\\over 2}n^9 +{3\\over 4}n^8 -{7\\over 10}n^6 +{1\\over 2}n^4 -{1\\over 12}nn$$\n",
    "\n",
    "But the actual closed formula is\n",
    "\n",
    "$$\\int n^9 = {1\\over 10}n^{10} + {1\\over 2}n^9 +{3\\over 4}n^8 -{7\\over 10}n^6 +{1\\over 2}n^4 -{3\\over 20}nn$$\n",
    "\n",
    "Indeed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([\n",
       "[ 1/10],\n",
       "[  1/2],\n",
       "[  3/4],\n",
       "[    0],\n",
       "[-7/10],\n",
       "[    0],\n",
       "[  1/2],\n",
       "[    0],\n",
       "[-3/20],\n",
       "[    0],\n",
       "[    0]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula_sum_of_powers(9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bernouilli numbers\n",
    "\n",
    "Finally, let's note that the coefficient of $n$ in the polynomial for $S_k(n)$ is the k-th [Bernouilli number](https://en.wikipedia.org/wiki/Bernoulli_number). Therefore, we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bernouilli_number(k):\n",
    "    return formula_sum_of_powers(k)[k]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which can be used to find the first 20 Bernouilli numbers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "B_1 = 1/2\n",
      "B_2 = 1/6\n",
      "B_3 = 0\n",
      "B_4 = -1/30\n",
      "B_5 = 0\n",
      "B_6 = 1/42\n",
      "B_7 = 0\n",
      "B_8 = -1/30\n",
      "B_9 = 0\n",
      "B_10 = 5/66\n",
      "B_11 = 0\n",
      "B_12 = -691/2730\n",
      "B_13 = 0\n",
      "B_14 = 7/6\n",
      "B_15 = 0\n",
      "B_16 = -3617/510\n",
      "B_17 = 0\n",
      "B_18 = 43867/798\n",
      "B_19 = 0\n",
      "B_20 = -174611/330\n"
     ]
    }
   ],
   "source": [
    "for i in range(1, 21):\n",
    "    print(\"B_\", i, \" = \", bernouilli_number(i), sep=\"\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}