Unbalanced numeric Feistel network implementation
Does the implementation follow accurately the described network generalization?

No. The Feistel network analyzed in the paper assumes independent round functions, but the code does not give that.
The issue is where the code reads random.seed(input + key + round). The corresponding pseudocode either has the same issue, or is at at least ambiguous with
$$\text{return }\mathit{hash}(\mathit{data}+\mathit{key}+\mathit{round})$$
when perhaps it means
$$\text{return }\mathit{hash}(\widehat{\mathit{data}}\mathbin\|\widehat{\mathit{key}}\mathbin\|\widehat{\mathit{round}})$$
where $\widehat x$ is a representation of integer $x$ per (say) big-endian binary, over a fixed-width bitstring chosen large enough for the maximum possible value.
If we use the same hash at each round, it is desirable that the input of $\mathit{hash}$ is unique for any triple $(\mathit{data},\mathit{key},\mathit{round})$. We can do this with just numbers, as
$$\text{return }\mathit{hash}(\mathit{data}+(\mathit{key}\cdot\mathit{nbrounds}+\mathit{round})\cdot\mathit{second})$$
Problem with the current proposition is: there is very little difference in the round functions for the different rounds. For a fixed $\mathit{key}$, if we note $F$ the function $x\mapsto\mathit{hash}(\mathit{key}+x)$, then the Feistel round function $F_i$ at round $i$ is $x\mapsto F_i(x)=F(i+x)$, which is very far from independent functions.
This is a good context for slide attack when $\mathit{first}=\mathit{second}$ and these are small (that is, $2^\mathit{first}$ is enumerable). At least a related-key attack works: no mater how many rounds, we can query two black boxes that implement the cipher and tell if they have their keys offset by one or not. And because the round function's result is added, which is associative with the $i+x$ operation at the input of $F$, I think it extends to a distinguishing attack in the single random key setup; perhaps even a way to tabulate the $F$ function and recover the functional equivalent of the key in sizably less than $4^\mathit{first}$ queries.


Is there any way of reversing the modulus addition ⊞ on the encryption function other than by using modulus subtraction ⊟ on the decryption function?

I see nothing simpler that works.


What would be the security implications of dynamically computing $\mathit{first}$ and $\mathit{second}$ as the $\mathit{first}\cdot\mathit{second}$ finite field is strictly equal to the input $\mathit{data}$?

As is, with $\mathit{data}$ the input block, that's not functional. This would mean the maximum value of the ciphertext is $\mathit{data}$. Thus, by induction, either the cipher is identity (which it won't be in all likelyhood), or it is not injective thus is not decipherable.
Therefore, we can't do without at least another input, which is the $\mathit{top}$ of the range $[0,\mathit{top})$ for $\mathit{data}$. And then we have a serious problem when $\mathit{top}$ is prime.
In Format Preserving Encryption (the academic name of what the question studies), one standard solution to this is cycle walking, which essentially increases $\mathit{top}$ a little until finding $\mathit{top'}=\mathit{first}\cdot\mathit{second}$ with $\min(\mathit{first},\mathit{second})$ above some reasonableness thresold, then fixes encryption by encrypting again when the ciphertext is in range $[\mathit{top},\mathit{top'})$. That's compensated at decryption by decrypting again when the decrypted plaintext is in range $[\mathit{top},\mathit{top'})$.

I have not touched side channel attacks, which are heavily language-dependent anyway.