implementing Lucas Primality Test

Im trying to implement the Lucas Primality test as described in FIPS 186-3 C.3.3. However, ive hit a problem:

im just testing my code with the value 33, im getting $D = 5$ as the first value for $\left(\frac{D}{33}\right) = -1$. with that, at step 6.2 for starting value $i = 4$ im getting $V_{temp} = {{{1 * 1 + 5 * 1 * 1} \over 2} \mod 33} = 3$. then, in the if/else block, there is a $"V_i ="$ step. since at $i = 4, k_i = 0$, the step taken is $V_i = V_{temp}$, how do i stuff 2 bits into a $V_i$? am i supposed to add? bitwise or? do i mod the value by 2 before i set the bit?

// integer is a custom arbitrary precision integer type i wrote
// it's supposed works like a normal signed int
// integer::bits() returns the length of the value's binary string
// integer::operator[] returns the bit at whatever given index, with index 0 being the least significant digit and bits() - 1 being the most significant digit

bool Lucas_FIPS186(integer C){
if (perfect_square(C))
return 0;
integer D = -3;
bool flip = false;
integer j = 1;
while (j != -1){
D = abs(D) + 2;
if (flip)
D = -D;
j = jacobi(D, C);
if (j == 0)
return 0;
flip ^= 1;
}
integer K = C + 1;
integer U = 1;
integer V = 1;
for(unsigned int i = K.bits() - 1; i > 0; i--){
integer Utemp = (U * V) % C;
integer Vtemp = (((V * V) + (D * U * U)) / 2) % C;
if (K[i - 1] == 1){
U = (((Utemp + Vtemp) / 2) % C);
V = (((Vtemp + D * Utemp) / 2) % C);
}
else{
U = Utemp;
V = Vtemp;
}
}
return !(U == 0);
}


EDIT: The code has been updated. 33 is still a problem. I also tested a very large known composite, and it returned prime

-
Maybe this would go over better on the programming website. –  Gerry Myerson Jun 26 '12 at 4:02

There's no reason to treat U and V as bit sequences. each $U_i$ and $V_i$ is a residue $\bmod C$ that should be stored as an integer.

[Edit in response to the edit in the question:]

You're taking the division by $2$ too literally. Note the explanation following the algorithm in the file you linked to:

Steps 6.2, 6.3.1 and 6.3.2 contain expressions of the form $A/2 \bmod C$, where $A$ is an integer, and $C$ is an odd integer. If $A/2$ is not an integer (i.e., $A$ is odd), then $A/2 \bmod C$ may be calculated as $(A+C)/2 \bmod C$. Alternatively, $A/2 \bmod C = A\cdot(C+1)/2 \bmod C$, for any integer $A$, without regard to $A$ being odd or even.

-
5, -7, 9, -11, 13, -15, 17, -19, 21 seems right to me, unless i missed something. also, +D means "force the sign to be positive", so if it were only D + 2, i would get -7 -> -5. –  calccrypto Jun 26 '12 at 5:50
@calccrypto: Sorry, I didn't realize this was your custom type with a custom + operator; I removed that part of my answer. I would recommend not to define operators such that their action disagrees with what one would expect from similar types; D.abs () would be much clearer. –  joriki Jun 26 '12 at 6:31
sorry. doesnt +int cause a normal int to become positive? I also updated my post. it should be clearer now, with potential fixes –  calccrypto Jun 26 '12 at 6:51
@calccrypto: No, $+x$ is just $x$ for a normal int, and for the same reason: It would be bad style to have a programming language gratuitously give a different meaning to $+x$ than mathematics. I updated the answer in response to your update in the post. –  joriki Jun 26 '12 at 7:03
how did i miss that? –  calccrypto Jun 26 '12 at 7:09