Compact Position Reporting (CPR) is a way of reducing the number of bits needed to transmit position whilst maintaining high position resolution (~5.1 meters for airborne encoding). Ordinary encoding would require 45 bits, but CPR uses 35 bits thus saving 10 bits. To do so, the world is divided up into a number of zones, both for latitude and for longitude and the position within the zone is encoded. Since transmitting the zone numbers would cost as many bits as that would have been saved a clever trick is used. Two types of encoding are used, each with a different number of zones. From the difference in zone coordinates between the two types of encoding one can determine which zone the aircraft is in. Therefore one needs to receive an even and an odd encoded message within 10 seconds from each other. Once the initial position decoding is successful every subsequent position message will give the updated position.
In CPR decoding, you first determine the latitude of the aircraft based on a pair of odd and even position messages received within 10 seconds from each other.
Latitude is always encoded in the same number of zones between the equator and the poles (15 for even encoding, 14.75 for odd encoding). Due to the roundness of the earth the numbers of zones used for longitude encoding depends on the latitude. Longitude even encoding uses 59 zones at the equator and only 1 at the poles. Odd encoding uses always 1 zone less, except near the pole where 1 zone is used like in even encoding. The number of longitude zones used for even encoding can be found by the $NL()$ function.
The transition latitudes, i.e. where the number of zones for longitude encoding change, can be calculated by the inverse $NL()$ function. This function can be used to create a lookup table of transition latitudes which is more efficient to use than the original $NL()$ function.
Without knowing anything about your implementation it is difficult to assist further, the formula is depicted correctly in your question. Note that input latitude is given in degrees, but the trigonometric function are assumed to work on radians hence the $\frac{\pi}{180}$ factor.
Edit: on further inspection I noticed a mistake in the formula: it should be :
$$
NL = int \left( 2 \pi \left[ arccos \left( 1 - { 1 - cos \left( {\pi \over 2NZ} \right) \over cos^2 \left( { \pi \over 180} \mid lat \mid \right) } \right) \right] ^{-1} \right)
$$
(note that the conversion from degrees to radians is corrected)
For those that really want to understand CPR coding and decoding in detail, you'll have to obtain a copy of RTCA document DO-260B or Eurocae document ED-102A (unfortunately they're expensive, I think about $700) and read Appendix A.1.7 and Appendix T.
