There is a simple geometric way to think about the fact that $P_n^{\text{ab}} \cong \mathbb{Z}^{\binom{n}{2}}$: Draw the braid with strands going from left to right across the page. Pick any two of the strands, $p$ and $q$. Let $c$ be a crossing of $p$ and $q$. Any time that $p$ and $q$ cross, either $p$ is going up the page and $q$ down the page or vice versa: Set $\alpha_{c} = \pm 1$ accordingly. And either $p$ is coming above the page and $q$ below the page or vice versa, set $\beta_{c} = \pm 1$ accordingly. Note that the positive braid generator $\sigma$ will have $\alpha_c \beta_c = 1$ and the negative braid generator will have $\alpha_c \beta_c = -1$. Set
$$k_{pq} = \sum_{\text{$p$ and $q$ cross at $c$}} \alpha_c \beta_c.$$
Since our word is in the pure braid group, $k_{pq}$ is even; set $\ell_{pq} = \tfrac{1}{2} k_{pq}$. Then $\ell_{pq}$ is a homomorphism $P_n \to \mathbb{Z}$, and combinining this for all $1 \leq p < q \leq n$ gives the abelianization map $P_n \to \mathbb{Z}^{\binom{n}{2}}$.
I don't know what a reference for this is; I checked it by remembering that $P_n$ is the fundamental group of $\mathcal{A}:=\{ (z_1, z_2, \ldots, z_n) \in \mathbb{C}^n : z_p \neq z_q \ \text{for}\ p \neq q \}$ and remembering that $H^1$ of a hyperplane arrangement is given by the maps to $\mathbb{C}^{\ast}$ given by the defining hyperplanes, so in this case by the map $\mathcal{A} \to (\mathbb{C}^{\ast})^{\binom{n}{2}}$ by $(z_1, z_2, \ldots, z_n) \to (z_1-z_2, z_1-z_3, \ldots, z_{n-1}-z_n)$.
Here is how I would implement that on a computer. Let your word be $\sigma_{i_1}^{\epsilon_1} \sigma_{i_2}^{\epsilon_2} \cdots \sigma_{i_N}^{\epsilon_N}$, with $\epsilon_j = \pm 1$. Let $s_i$ be the element $(i\ i+1)$ in $S_n$ and let $w_j = s_{i_1} s_{i_2} \cdots s_{i_j}$, which you can store compactly as a list $(w_j(1), w_j(2), \ldots, w_j(n))$. Let $a_j = w_j(i_j)$ and $b_j = w_j(i_j+1)$, so $w_{j+1}(i_j) = b_j$ and $w_{j+1}(i_j) = a_j$. Then
$$\ell_{pq}(w) = \frac{1}{2} \sum_{\{ a_j,b_j \} = \{ p,q \}} \epsilon_j = \sum_{(a_j, b_j) = (p,q)} \epsilon_j.$$
Note that the difference between the two formulas is that I have curly braces, meaning unordered pairs, in the subscript of the first summation and round braces, meaning ordered pairs, on the right.