Ok, I now have some time, so I'll lay out the basics.
First step: how to "fast forward" your construction.
One observation is that your operation is entirely bitwise-linear; that is, there exists a $128\times 128$ matrix that, when multiplied by the input, generates the output. And, that matrix is easy to compute.
Obviously, computing the output using this matrix multiplication is slower than using the obvious code; however "though this be madness, yet there is method in't".
The second observation is that if we multiply two such $128 \times 128$ matricies $A \times B$ together, when we multiply the product by a 128 bit input, this has the effect of performing the operation $B$ first, and then operation $A$. That is, if $A$ was 'performing your operation $a$ times, and $B$ was 'performing your operation $b$ times', then the effect of the matrix $A \times B$ is 'performing your operation $a+b$ times'.
This allows us to perform the operation in faster-than-linear time. In the simplest case, to advance your operation 4 times, we can multiply by the matrix $A \times A \times A \times A = (A \times A) \times (A \times A)$ - this requires just 3 matrix multications, rather than 4. That might not sound like such an impressive speed-up; however consider if we want to compute the matrix that performs your operation 1024 times, all we need to do is take your original matrix, and square it 10 times (where squaring a matrix means multiplying it by itelf). More generally, to compute it $k$ times (for an arbitrary positive integer $k$), we just compute an addition chain that corresponds to $k$; this allows us to compute the matrix in $O(\log k)$ matrix multiplications. For example, to fast forward the operation $2^{128}-1$ steps, instead of calling your fast code $2^{128}-1$ times, we can use the binary method and perform 255 matrix multiplications (and the binary method is not the most efficient; however it's a good place to start, and is only a factor of less than 2 from being optimal).
Now, with that done, we look at how we can verify whether your operation (with specific shifts) actually achieves maximal cycle length.
Now, we know that, for your operation, the all zero vector will always map to the all zero vector, and so that's a cycle of one. Hence, the largest possible cycle for everything else is a cycle of length $2^{128}-1$ (and we know that is achievable with a linear operation, because of the existence of LFSR with a primitive feedback polynomial).
So, the first thing we can check is "what does your operation cycled $2^{128}-1$ look like? If the matrix corresponding to that is anything other than the identity matrix (that is, one that maps anything to itself), then we know that an input must have a cycle length less than $2^{128}-1$ (because it cannot have a larger one).
If that is the identity matrix, then the cycle length must be $(2^{128}-1)/k$ for some integer $k$; the next thing to check is if $k=1$.
For that, we can just check the operation iterated $(2^{128}-1)/p$ times (for all primes $p$ that are divisors of $2^{128}-1$); if $k$ has $p$ as a factor, then that will be the identity matrix. Hence, if it is not the identity matrix for all prime factors $p$, then $k$ must not be divisible by any factor, and so the only possibility is $k=1$, and we are done...
The only question I didn't try to address is "what are the prime factors of $2^{128}-1$ anyways?" Well, $2^{128}-1 = 3 \cdot 5 \cdot 17 \cdot 257 \cdot 641 \cdot 65537 \cdot 274177 \cdot 6700417 \cdot 67280421310721$
Ok, actually implementing this is a bit of work - but it is feasible, and will allow you to perform the test you are looking for...
