CSharp examples for System:Math Statistics
Returns the k_th number in the Lucas Sequence reduced modulo n
// All rights reserved. using System;//from w ww . j a va 2 s .com public class Main{ //*********************************************************************** // Returns the k_th number in the Lucas Sequence reduced modulo n. // // Uses index doubling to speed up the process. For example, to calculate V(k), // we maintain two numbers in the sequence V(n) and V(n+1). // // To obtain V(2n), we use the identity // V(2n) = (V(n) * V(n)) - (2 * Q^n) // To obtain V(2n+1), we first write it as // V(2n+1) = V((n+1) + n) // and use the identity // V(m+n) = V(m) * V(n) - Q * V(m-n) // Hence, // V((n+1) + n) = V(n+1) * V(n) - Q^n * V((n+1) - n) // = V(n+1) * V(n) - Q^n * V(1) // = V(n+1) * V(n) - Q^n * P // // We use k in its binary expansion and perform index doubling for each // bit position. For each bit position that is set, we perform an // index doubling followed by an index addition. This means that for V(n), // we need to update it to V(2n+1). For V(n+1), we need to update it to // V((2n+1)+1) = V(2*(n+1)) // // This function returns // [0] = U(k) // [1] = V(k) // [2] = Q^n // // Where U(0) = 0 % n, U(1) = 1 % n // V(0) = 2 % n, V(1) = P % n //*********************************************************************** public static BigInteger[] LucasSequence(BigInteger P, BigInteger Q, BigInteger k, BigInteger n) { if (k.dataLength == 1 && k.data[0] == 0) { BigInteger[] result = new BigInteger[3]; result[0] = 0; result[1] = 2 % n; result[2] = 1 % n; return result; } // calculate constant = b^(2k) / m // for Barrett Reduction BigInteger constant = new BigInteger(); int nLen = n.dataLength << 1; constant.data[nLen] = 0x00000001; constant.dataLength = nLen + 1; constant = constant / n; // calculate values of s and t int s = 0; for (int index = 0; index < k.dataLength; index++) { uint mask = 0x01; for (int i = 0; i < 32; i++) { if ((k.data[index] & mask) != 0) { index = k.dataLength; // to break the outer loop break; } mask <<= 1; s++; } } BigInteger t = k >> s; //Console.WriteLine("s = " + s + " t = " + t); return LucasSequenceHelper(P, Q, t, n, constant, s); } //*********************************************************************** // Performs the calculation of the kth term in the Lucas Sequence. // For details of the algorithm, see reference [9]. // // k must be odd. i.e LSB == 1 //*********************************************************************** private static BigInteger[] LucasSequenceHelper(BigInteger P, BigInteger Q, BigInteger k, BigInteger n, BigInteger constant, int s) { BigInteger[] result = new BigInteger[3]; if ((k.data[0] & 0x00000001) == 0) throw (new ArgumentException("Argument k must be odd.")); int numbits = k.bitCount(); uint mask = (uint)0x1 << ((numbits & 0x1F) - 1); // v = v0, v1 = v1, u1 = u1, Q_k = Q^0 BigInteger v = 2 % n, Q_k = 1 % n, v1 = P % n, u1 = Q_k; bool flag = true; for (int i = k.dataLength - 1; i >= 0; i--) // iterate on the binary expansion of k { //Console.WriteLine("round"); while (mask != 0) { if (i == 0 && mask == 0x00000001) // last bit break; if ((k.data[i] & mask) != 0) // bit is set { // index doubling with addition u1 = (u1 * v1) % n; v = ((v * v1) - (P * Q_k)) % n; v1 = n.BarrettReduction(v1 * v1, n, constant); v1 = (v1 - ((Q_k * Q) << 1)) % n; if (flag) flag = false; else Q_k = n.BarrettReduction(Q_k * Q_k, n, constant); Q_k = (Q_k * Q) % n; } else { // index doubling u1 = ((u1 * v) - Q_k) % n; v1 = ((v * v1) - (P * Q_k)) % n; v = n.BarrettReduction(v * v, n, constant); v = (v - (Q_k << 1)) % n; if (flag) { Q_k = Q % n; flag = false; } else Q_k = n.BarrettReduction(Q_k * Q_k, n, constant); } mask >>= 1; } mask = 0x80000000; } // at this point u1 = u(n+1) and v = v(n) // since the last bit always 1, we need to transform u1 to u(2n+1) and v to v(2n+1) u1 = ((u1 * v) - Q_k) % n; v = ((v * v1) - (P * Q_k)) % n; if (flag) flag = false; else Q_k = n.BarrettReduction(Q_k * Q_k, n, constant); Q_k = (Q_k * Q) % n; for (int i = 0; i < s; i++) { // index doubling u1 = (u1 * v) % n; v = ((v * v) - (Q_k << 1)) % n; if (flag) { Q_k = Q % n; flag = false; } else Q_k = n.BarrettReduction(Q_k * Q_k, n, constant); } result[0] = u1; result[1] = v; result[2] = Q_k; return result; } //*********************************************************************** // Fast calculation of modular reduction using Barrett's reduction. // Requires x < b^(2k), where b is the base. In this case, base is // 2^32 (uint). // // Reference [4] //*********************************************************************** private BigInteger BarrettReduction(BigInteger x, BigInteger n, BigInteger constant) { int k = n.dataLength, kPlusOne = k + 1, kMinusOne = k - 1; BigInteger q1 = new BigInteger(); // q1 = x / b^(k-1) for (int i = kMinusOne, j = 0; i < x.dataLength; i++, j++) q1.data[j] = x.data[i]; q1.dataLength = x.dataLength - kMinusOne; if (q1.dataLength <= 0) q1.dataLength = 1; BigInteger q2 = q1 * constant; BigInteger q3 = new BigInteger(); // q3 = q2 / b^(k+1) for (int i = kPlusOne, j = 0; i < q2.dataLength; i++, j++) q3.data[j] = q2.data[i]; q3.dataLength = q2.dataLength - kPlusOne; if (q3.dataLength <= 0) q3.dataLength = 1; // r1 = x mod b^(k+1) // i.e. keep the lowest (k+1) words BigInteger r1 = new BigInteger(); int lengthToCopy = (x.dataLength > kPlusOne) ? kPlusOne : x.dataLength; for (int i = 0; i < lengthToCopy; i++) r1.data[i] = x.data[i]; r1.dataLength = lengthToCopy; // r2 = (q3 * n) mod b^(k+1) // partial multiplication of q3 and n BigInteger r2 = new BigInteger(); for (int i = 0; i < q3.dataLength; i++) { if (q3.data[i] == 0) continue; ulong mcarry = 0; int t = i; for (int j = 0; j < n.dataLength && t < kPlusOne; j++, t++) { // t = i + j ulong val = ((ulong)q3.data[i] * (ulong)n.data[j]) + (ulong)r2.data[t] + mcarry; r2.data[t] = (uint)(val & 0xFFFFFFFF); mcarry = (val >> 32); } if (t < kPlusOne) r2.data[t] = (uint)mcarry; } r2.dataLength = kPlusOne; while (r2.dataLength > 1 && r2.data[r2.dataLength - 1] == 0) r2.dataLength--; r1 -= r2; if ((r1.data[maxLength - 1] & 0x80000000) != 0) // negative { BigInteger val = new BigInteger(); val.data[kPlusOne] = 0x00000001; val.dataLength = kPlusOne + 1; r1 += val; } while (r1 >= n) r1 -= n; return r1; } //*********************************************************************** // Returns the position of the most significant bit in the BigInteger. // // Eg. The result is 0, if the value of BigInteger is 0...0000 0000 // The result is 1, if the value of BigInteger is 0...0000 0001 // The result is 2, if the value of BigInteger is 0...0000 0010 // The result is 2, if the value of BigInteger is 0...0000 0011 // //*********************************************************************** public int bitCount() { while (dataLength > 1 && data[dataLength - 1] == 0) dataLength--; uint value = data[dataLength - 1]; uint mask = 0x80000000; int bits = 32; while (bits > 0 && (value & mask) == 0) { bits--; mask >>= 1; } bits += ((dataLength - 1) << 5); return bits; } }