How does the squeeze function in the NaCl Poly1305 implementation work? The NaCl ref implementation of the Poly1305 algorithm uses the following reduction function (which is called squeeze()):

static void squeeze(unsigned int h[17])
{
    unsigned int j;
    unsigned int u;
    u = 0;
    for (j = 0;j &lt; 16;++j) { u += h[j]; h[j] = u &amp; 255; u &gt;&gt;= 8; }
    u += h[16]; h[16] = u &amp; 3;
    u = 5 * (u &gt;&gt; 2);
    for (j = 0;j &lt; 16;++j) { u += h[j]; h[j] = u &amp; 255; u &gt;&gt;= 8; }
    u += h[16]; h[16] = u;
}


I understand globally what is done in the first for loop, namely that each integer in h is reduced to 256 bits, with the carry being carried over to the next integer. I only don't grasp what happens from the h[16] = u &amp; 3; on. Where does the u &amp; 3; come from? and the u = 5 * (u &gt;&gt; 2);? 

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?

Does anyone know how this squeeze function works, and provide an explanation, as well as an answer to my second question?

The implementation itself provides no documentation, and the answer on How does NaCl Poly1305 implementation do modular multiplication? does not delve into the squeeze function. See this stackexchange question for the full code.
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!