How does the squeeze function in the NaCl Poly1305 implementation work?
There are two ideas going on here:

Radix 256 arithmetic with delayed carries.
Here we represent an integer $x$ by $x_0 + 2^8 x_1 + 2^{16} x_2 + \dots + 2^{128} x_{16}$.  In canonical form, the digits $x_i$ lie in $\{0,1,2,\dots,255\}$, but in the Poly1305 computation, we delay the propagation of carries as long as we reliably can in order to save computation, so the $x_i$ might be larger, giving multiple possible representations for each integer $x$ in this form.

Poly1305 reduction $2^{130} x_{\mathrm{hi}} + x_{\mathrm{lo}} \equiv 5 x_{\mathrm{hi}} + x_{\mathrm{lo}} \pmod p$.
We are working modulo $p = 2^{130} - 5$, which has $2^{130} - 5 \equiv 0 \pmod p$ and so $2^{130} \equiv 5 \pmod p$.  The idea is that if we have an integer $x$, we can take its 130 low bits $x_{\mathrm{lo}} = x \mathbin\&amp; (2^{130} - 1)$, and the high bits $x_{\mathrm{hi}} = x \mathbin\gg 130$, so that $$x = (x_{\mathrm{hi}} \ll 130) \mathbin| x_{\mathrm{lo}} = 2^{130} x_{\mathrm{hi}} + x_{\mathrm{lo}},$$ and then do a reduction step by computing $$x \equiv 5x_{\mathrm{hi}} + x_{\mathrm{lo}} \pmod p;$$ i.e., shift/multiply/add.


The two ideas are done simultaneously.  I recommend that you try to work through it yourself from those two ideas, but—spoiler alert!—here's the full gory details if you're still completely stuck:

static void squeeze(unsigned h[17])

On entry, $h$ represents an integer $x$ in radix 256, with digits that may exceed 255 because of delayed carries, as long as they have some bound $B$ so the carries aren't too excessive.
$x = h[0] + 2^8 h[1] + 2^{16} h[2] + \dots + 2^{128} h[16] \\
h[0], h[1], \dots, h[16] &lt; B$
First, we propagate all the delayed carries and bound each digit by the radix up to 128 bits:
u = 0;
for (j = 0; j &lt; 16; j++)


Loop invariants:
$x = h[0] + 2^8 h[1] + \dots + 2^{8j} (h[j] + u) + \dots + 2^{128} h[16] \\
h[0], h[1], \dots, h[j] &lt; 256 \\
h[j + 1], \dots, h[16] &lt; B$


    { u += h[j]; h[j] = u &amp; 255; u &gt;&gt;= 8; }

Now $h[0], \dots, h[15]$ are reduced, leaving us with $h[16]$ and a 128-bit carry.  We'll go two bits further to get a 130-bit carry:
$x = h[0] + 2^8 h[1] + \dots + 2^{128} (h[16] + u) \\
h[0], h[1], \dots, h[15] &lt; 256 \\
h[16] &lt; B$
u += h[16]; h[16] = u &amp; 3;

Now note that $2^{130} \equiv 5 \pmod p$.  So—after we shift $u \gg 2$ below—when we have
$x = h[0] + 2^8 h[1] + \dots + 2^{128} h[16] + 2^{130} u \\
h[0], h[1], \dots, h[15] &lt; 256 \\
h[16] &lt; 4$
we can reduce $x = h[0] + \dots + 2^{130} u \equiv h[0] + \dots + 5 u \pmod p$.
u = 5 * (u &gt;&gt; 2);

Now we can add $5 u$ starting at the first digit, and propagate the carries:
$x = h[0] + 2^8 h[1] + \dots + 2^{128} h[16] + 5 u \\
h[0], h[1], \dots, h[15] &lt; 256 \\
h[16] &lt; 4$
for (j = 0;j &lt; 16;++j)


Loop invariants:
$x = h[0] + 2^8 h[1] + \dots + 2^{8j} (h[j] + u) + \dots + 2^{128} h[16] \\
h[0], h[1], \dots, h[15] &lt; 256 \\
h[16] &lt; 4$


    { u += h[j]; h[j] = u &amp; 255; u &gt;&gt;= 8; }

Now we have a 128-bit carry again, but there's plenty of room because $h[16]$ is only two bits at this point, from the above h[16] = (u &amp; 3).
$x = h[0] + 2^8 h[1] + \dots + 2^{128} (h[16] + u)
h[0], h[1], \dots, h[15] &lt; 256 \\
h[16] &lt; 4$
u += h[16]; h[16] = u;

Finally, $h$ again represents an integer in radix 256, with digits that are bounded by the radix:
$x = h[0] + 2^8 h[1] + \dots + 2^{128} h[16] \\
h[0], h[1], \dots, h[16] &lt; 256$
However, it is not necessarily the least representative modulo $p$.  For that, we need to freeze; squeeze just propagates a batch of delayed carries and does a single reduction step so there's room for more arithmetic.
(Figuring out what $B$ is and setting corresponding bounds on the carry $u$ in each step left as an exercise for the reader.)


Also, if I would want to change the radix here, say from 2^8 bits to 2^16, would that require me to change the above mentioned fragments? I would need to change al the 8's to 16, and the 255's to 2^16-1, but would that be it?

The array $h$ represents the integer $x = h[0] + 2^8 h[1] + 2^{16} h[2] + \dots + 2^{128} h[16]$ on entry, and $x - (2^{130} - 5) k$ for some $k \geq 0$ on exit after the reduction step where we replace $2^{130} u$ by $5 u$.  (Casting out $2^{130} - 5$'s.)
If you want to change the radix, say from $2^8$ to $2^{17}$, then the array $h$ will instead represent $x = h[0] + 2^{17} h[1] + 2^{34} h[2] + \dots + 2^{119} h[7] + 2^{136} h[8]$.  The goal is to compute the same reduction, $x - (2^{130} - 5) k$ for some $k \geq 0$, but finding multiples of $2^{130}$ to replace by multiples of $5$ for the reduction step will happen at a different digit position, $h[i]$, and a different bit position within that digit or carry, $u \gg j$.
It is an easy exercise to find these positions.  For a simpler exercise, you could try writing arithmetic modulo $2^{31} - 1$ in radix $2^8$, in which you can easily exhaustively test a reduction step for correctness.  I recommend that you carry out this exercise rather than asking pseudonymous strangers on the internet to spoon-feed you answers!