Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

Add new implementation of C# BigDecimal and BigInteger, also implemen…

…t sqrt correctly
  • Loading branch information...
commit ccbf538b4227e6f2db9f1b7dd1d6222fc40f0a4c 1 parent b6dec38
@olabini olabini authored
View
26 src/ikc/main/Ioke.Lang/Decimal.cs
@@ -12,7 +12,7 @@ public class Decimal : IokeData {
}
public Decimal(string textRepresentation) {
- this.value = new BigDecimal(textRepresentation);
+ this.value = new BigDecimal(textRepresentation).stripTrailingZeros();
}
public Decimal(BigDecimal value) {
@@ -32,12 +32,18 @@ public class Decimal : IokeData {
}
public string AsNativeString() {
- string s = value.ToString(MathContext.PLAIN);
- if(s[s.Length-1] == '0' && s.IndexOf('.') != -1 && s.IndexOf('e') == -1 && s.IndexOf('E') == -1) {
- int end = s.Length-1;
- while(s[end] == '0' && s[end-1] != '.') end--;
- if(s[end-1] == '.') end++;
- return s.Substring(0, end);
+ string s = value.toPlainString();
+ if(s.IndexOf('.') != -1) {
+ if(s[s.Length-1] == '0' && s.IndexOf('e') == -1 && s.IndexOf('E') == -1) {
+ int end = s.Length-1;
+ while(s[end] == '0' && s[end-1] != '.') end--;
+ if(s[end-1] == '.') end++;
+ return s.Substring(0, end);
+ }
+ } else {
+ if(s.IndexOf('e') == -1 && s.IndexOf('E') == -1) {
+ return s + ".0";
+ }
}
return s;
}
@@ -93,6 +99,12 @@ public class Decimal : IokeData {
return runtime.NewText(on.ToString());
})));
+ obj.RegisterMethod(runtime.NewNativeMethod("returns the square root of the receiver. this should return the same result as calling ** with 0.5",
+ new TypeCheckingNativeMethod.WithNoArguments("sqrt", obj,
+ (method, on, args, keywords, context, message) => {
+ return runtime.NewDecimal(new BigSquareRoot().Get(((Decimal)IokeObject.dataOf(on)).value));
+ })));
+
obj.RegisterMethod(obj.runtime.NewNativeMethod("Returns a text inspection of the object",
new TypeCheckingNativeMethod.WithNoArguments("inspect", obj,
(method, on, args, keywords, context, message) => {
View
21 src/ikc/main/Ioke.Lang/Number.cs
@@ -152,6 +152,27 @@ public class Number : IokeData {
return context.runtime.NewNumber(Number.GetValue(on).GetHashCode());
})));
+ number.RegisterMethod(runtime.NewNativeMethod("returns the square root of the receiver. this should return the same result as calling ** with 0.5",
+ new NativeMethod.WithNoArguments("sqrt", (method, context, message, on, outer) => {
+ outer.ArgumentsDefinition.CheckArgumentCount(context, message, on);
+
+ RatNum value = Number.GetValue(on);
+ if(value is IntFraction) {
+ IntNum num = value.numerator();
+ IntNum den = value.denominator();
+ BigDecimal nums = new BigSquareRoot().Get(num.AsBigDecimal());
+ BigDecimal dens = new BigSquareRoot().Get(den.AsBigDecimal());
+ try {
+ num = IntNum.valueOf(nums.toBigIntegerExact().ToString());
+ den = IntNum.valueOf(dens.toBigIntegerExact().ToString());
+ return context.runtime.NewNumber(new IntFraction(num, den));
+ } catch(ArithmeticException e) {
+ // Ignore and fall through
+ }
+ }
+ return context.runtime.NewDecimal(new BigSquareRoot().Get(value.AsBigDecimal()));
+ })));
+
number.RegisterMethod(runtime.NewNativeMethod("returns true if the left hand side number is equal to the right hand side number.",
new TypeCheckingNativeMethod("==", TypeCheckingArgumentsDefinition.builder()
.ReceiverMustMimic(runtime.Number)
View
3,563 src/ikc/main/Ioke.Math/BigDecimal.cs
1,502 additions, 2,061 deletions not shown
View
854 src/ikc/main/Ioke.Math/BigInteger.cs
@@ -0,0 +1,854 @@
+namespace Ioke.Math {
+ using System;
+ using System.Text;
+
+ public class BigInteger {
+ internal int[] digits;
+ internal int numberLength;
+ internal int sign;
+
+ public static readonly BigInteger ZERO = new BigInteger(0, 0);
+ public static readonly BigInteger ONE = new BigInteger(1, 1);
+ public static readonly BigInteger TEN = new BigInteger(1, 10);
+ internal static readonly BigInteger MINUS_ONE = new BigInteger(-1, 1);
+ internal const int EQUALS = 0;
+ internal const int GREATER = 1;
+ internal const int LESS = -1;
+ static readonly BigInteger[] SMALL_VALUES = { ZERO, ONE, new BigInteger(1, 2),
+ new BigInteger(1, 3), new BigInteger(1, 4), new BigInteger(1, 5),
+ new BigInteger(1, 6), new BigInteger(1, 7), new BigInteger(1, 8),
+ new BigInteger(1, 9), TEN };
+
+ static readonly BigInteger[] TWO_POWS;
+
+ static BigInteger() {
+ TWO_POWS = new BigInteger[32];
+ for(int i = 0; i < TWO_POWS.Length; i++) {
+ TWO_POWS[i] = BigInteger.valueOf(1L<<i);
+ }
+ }
+
+ private int firstNonzeroDigit = -2;
+ private int hashCode = 0;
+
+ public BigInteger(int numBits, Random rnd) {
+ if (numBits < 0) {
+ throw new ArgumentException("numBits must be non-negative");
+ }
+ if (numBits == 0) {
+ sign = 0;
+ numberLength = 1;
+ digits = new int[] { 0 };
+ } else {
+ sign = 1;
+ numberLength = (numBits + 31) >> 5;
+ digits = new int[numberLength];
+ for (int i = 0; i < numberLength; i++) {
+ digits[i] = rnd.Next();
+ }
+ // Using only the necessary bits
+ digits[numberLength - 1] = (int)(((uint)digits[numberLength - 1]) >> (-numBits) & 31);
+ cutOffLeadingZeroes();
+ }
+ }
+
+ public BigInteger(int bitLength, int certainty, Random rnd) {
+ if (bitLength < 2) {
+ throw new ArithmeticException("bitLength < 2");
+ }
+ BigInteger me = Primality.consBigInteger(bitLength, certainty, rnd);
+ sign = me.sign;
+ numberLength = me.numberLength;
+ digits = me.digits;
+ }
+
+ public BigInteger(String val) : this(val, 10) {
+ }
+
+ public BigInteger(String val, int radix) {
+ if (val == null) {
+ throw new NullReferenceException();
+ }
+ if ((radix < 0) || (radix > 26)) {
+ throw new System.FormatException("Radix out of range");
+ }
+ if (val.Length == 0) {
+ throw new System.FormatException("Zero length BigInteger");
+ }
+ setFromString(this, val, radix);
+ }
+
+ public BigInteger(int signum, byte[] magnitude) {
+ if (magnitude == null) {
+ throw new NullReferenceException();
+ }
+ if ((signum < -1) || (signum > 1)) {
+ // math.13=Invalid signum value
+ throw new System.FormatException("Invalid signum value");
+ }
+ if (signum == 0) {
+ foreach (byte element in magnitude) {
+ if (element != 0) {
+ // math.14=signum-magnitude mismatch
+ throw new System.FormatException("signum-magnitude mismatch");
+ }
+ }
+ }
+ if (magnitude.Length == 0) {
+ sign = 0;
+ numberLength = 1;
+ digits = new int[] { 0 };
+ } else {
+ sign = signum;
+ putBytesPositiveToIntegers(magnitude);
+ cutOffLeadingZeroes();
+ }
+ }
+
+ public BigInteger(byte[] val) {
+ if (val.Length == 0) {
+ // math.12=Zero length BigInteger
+ throw new System.FormatException("Zero length BigInteger");
+ }
+ if (val[0] < 0) {
+ sign = -1;
+ putBytesNegativeToIntegers(val);
+ } else {
+ sign = 1;
+ putBytesPositiveToIntegers(val);
+ }
+ cutOffLeadingZeroes();
+ }
+
+ internal BigInteger(int sign, int value) {
+ this.sign = sign;
+ numberLength = 1;
+ digits = new int[] { value };
+ }
+
+ internal BigInteger(int sign, int numberLength, int[] digits) {
+ this.sign = sign;
+ this.numberLength = numberLength;
+ this.digits = digits;
+ }
+
+ internal BigInteger(int sign, long val) {
+ // PRE: (val >= 0) && (sign >= -1) && (sign <= 1)
+ this.sign = sign;
+ if (((ulong)val & 0xFFFFFFFF00000000L) == 0) {
+ // It fits in one 'int'
+ numberLength = 1;
+ digits = new int[] { (int) val };
+ } else {
+ numberLength = 2;
+ digits = new int[] { (int) val, (int) (val >> 32) };
+ }
+ }
+
+ internal BigInteger(int signum, int[] digits) {
+ if (digits.Length == 0) {
+ sign = 0;
+ numberLength = 1;
+ this.digits = new int[] { 0 };
+ } else {
+ sign = signum;
+ numberLength = digits.Length;
+ this.digits = digits;
+ cutOffLeadingZeroes();
+ }
+ }
+
+ public static BigInteger valueOf(long val) {
+ if (val < 0) {
+ if (val != -1) {
+ return new BigInteger(-1, -val);
+ }
+ return MINUS_ONE;
+ } else if (val <= 10) {
+ return SMALL_VALUES[(int) val];
+ } else {// (val > 10)
+ return new BigInteger(1, val);
+ }
+ }
+
+ public byte[] toByteArray() {
+ if (this.sign == 0) {
+ return new byte[] { 0 };
+ }
+ BigInteger temp = this;
+ int bitLen = bitLength();
+ int iThis = getFirstNonzeroDigit();
+ int bytesLen = (bitLen >> 3) + 1;
+ /*
+ * Puts the little-endian int array representing the magnitude of this
+ * BigInteger into the big-endian byte array.
+ */
+ byte[] bytes = new byte[bytesLen];
+ int firstByteNumber = 0;
+ int highBytes;
+ int digitIndex = 0;
+ int bytesInInteger = 4;
+ int digit;
+ int hB;
+
+ if (bytesLen - (numberLength << 2) == 1) {
+ bytes[0] = (byte) ((sign < 0) ? -1 : 0);
+ highBytes = 4;
+ firstByteNumber++;
+ } else {
+ hB = bytesLen & 3;
+ highBytes = (hB == 0) ? 4 : hB;
+ }
+
+ digitIndex = iThis;
+ bytesLen -= iThis << 2;
+
+ if (sign < 0) {
+ digit = -temp.digits[digitIndex];
+ digitIndex++;
+ if (digitIndex == numberLength) {
+ bytesInInteger = highBytes;
+ }
+ for (int i = 0; i < bytesInInteger; i++, digit >>= 8) {
+ bytes[--bytesLen] = (byte) digit;
+ }
+ while (bytesLen > firstByteNumber) {
+ digit = ~temp.digits[digitIndex];
+ digitIndex++;
+ if (digitIndex == numberLength) {
+ bytesInInteger = highBytes;
+ }
+ for (int i = 0; i < bytesInInteger; i++, digit >>= 8) {
+ bytes[--bytesLen] = (byte) digit;
+ }
+ }
+ } else {
+ while (bytesLen > firstByteNumber) {
+ digit = temp.digits[digitIndex];
+ digitIndex++;
+ if (digitIndex == numberLength) {
+ bytesInInteger = highBytes;
+ }
+ for (int i = 0; i < bytesInInteger; i++, digit >>= 8) {
+ bytes[--bytesLen] = (byte) digit;
+ }
+ }
+ }
+ return bytes;
+ }
+
+ private static void setFromString(BigInteger bi, String val, int radix) {
+ int sign;
+ int[] digits;
+ int numberLength;
+ int stringLength = val.Length;
+ int startChar;
+ int endChar = stringLength;
+
+ if (val[0] == '-') {
+ sign = -1;
+ startChar = 1;
+ stringLength--;
+ } else {
+ sign = 1;
+ startChar = 0;
+ }
+ /*
+ * We use the following algorithm: split a string into portions of n
+ * characters and convert each portion to an integer according to the
+ * radix. Then convert an exp(radix, n) based number to binary using the
+ * multiplication method. See D. Knuth, The Art of Computer Programming,
+ * vol. 2.
+ */
+
+ int charsPerInt = Conversion.digitFitInInt[radix];
+ int bigRadixDigitsLength = stringLength / charsPerInt;
+ int topChars = stringLength % charsPerInt;
+
+ if (topChars != 0) {
+ bigRadixDigitsLength++;
+ }
+ digits = new int[bigRadixDigitsLength];
+ // Get the maximal power of radix that fits in int
+ int bigRadix = Conversion.bigRadices[radix - 2];
+ // Parse an input string and accumulate the BigInteger's magnitude
+ int digitIndex = 0; // index of digits array
+ int substrEnd = startChar + ((topChars == 0) ? charsPerInt : topChars);
+ int newDigit;
+
+ for (int substrStart = startChar; substrStart < endChar; substrStart = substrEnd, substrEnd = substrStart
+ + charsPerInt) {
+ int bigRadixDigit = Convert.ToInt32(val.Substring(substrStart,
+ substrEnd - substrStart), radix);
+ newDigit = Multiplication.multiplyByInt(digits, digitIndex,
+ bigRadix);
+ newDigit += Elementary
+ .inplaceAdd(digits, digitIndex, bigRadixDigit);
+ digits[digitIndex++] = newDigit;
+ }
+ numberLength = digitIndex;
+ bi.sign = sign;
+ bi.numberLength = numberLength;
+ bi.digits = digits;
+ bi.cutOffLeadingZeroes();
+ }
+
+ public BigInteger abs() {
+ return ((sign < 0) ? new BigInteger(1, numberLength, digits) : this);
+ }
+
+ public BigInteger negate() {
+ return ((sign == 0) ? this
+ : new BigInteger(-sign, numberLength, digits));
+ }
+
+ public BigInteger add(BigInteger val) {
+ return Elementary.add(this, val);
+ }
+
+ public BigInteger subtract(BigInteger val) {
+ return Elementary.subtract(this, val);
+ }
+
+ public int signum() {
+ return sign;
+ }
+
+ public BigInteger shiftRight(int n) {
+ if ((n == 0) || (sign == 0)) {
+ return this;
+ }
+ return ((n > 0) ? BitLevel.shiftRight(this, n) : BitLevel.shiftLeft(
+ this, -n));
+ }
+
+ public BigInteger shiftLeft(int n) {
+ if ((n == 0) || (sign == 0)) {
+ return this;
+ }
+ return ((n > 0) ? BitLevel.shiftLeft(this, n) : BitLevel.shiftRight(
+ this, -n));
+ }
+
+ internal BigInteger shiftLeftOneBit() {
+ return (sign == 0) ? this : BitLevel.shiftLeftOneBit(this);
+ }
+
+ public int bitLength() {
+ return BitLevel.bitLength(this);
+ }
+
+ public bool testBit(int n) {
+ if (n == 0) {
+ return ((digits[0] & 1) != 0);
+ }
+ if (n < 0) {
+ throw new ArithmeticException("Negative bit address");
+ }
+ int intCount = n >> 5;
+ if (intCount >= numberLength) {
+ return (sign < 0);
+ }
+ int digit = digits[intCount];
+ n = (1 << (n & 31)); // int with 1 set to the needed position
+ if (sign < 0) {
+ int firstNonZeroDigit = getFirstNonzeroDigit();
+ if (intCount < firstNonZeroDigit) {
+ return false;
+ } else if (firstNonZeroDigit == intCount) {
+ digit = -digit;
+ } else {
+ digit = ~digit;
+ }
+ }
+ return ((digit & n) != 0);
+ }
+
+ public BigInteger setBit(int n) {
+ if (!testBit(n)) {
+ return BitLevel.flipBit(this, n);
+ }
+ return this;
+ }
+
+ public BigInteger clearBit(int n) {
+ if (testBit(n)) {
+ return BitLevel.flipBit(this, n);
+ }
+ return this;
+ }
+
+ public BigInteger flipBit(int n) {
+ if (n < 0) {
+ throw new ArithmeticException("Negative bit address");
+ }
+ return BitLevel.flipBit(this, n);
+ }
+
+ public int getLowestSetBit() {
+ if (sign == 0) {
+ return -1;
+ }
+ // (sign != 0) implies that exists some non zero digit
+ int i = getFirstNonzeroDigit();
+ return ((i << 5) + BigDecimal.numberOfTrailingZeros(digits[i]));
+ }
+
+ public int bitCount() {
+ return BitLevel.bitCount(this);
+ }
+
+ public BigInteger not() {
+ return Logical.not(this);
+ }
+
+ public BigInteger and(BigInteger val) {
+ return Logical.and(this, val);
+ }
+
+ public BigInteger or(BigInteger val) {
+ return Logical.or(this, val);
+ }
+
+ public BigInteger xor(BigInteger val) {
+ return Logical.xor(this, val);
+ }
+
+ public BigInteger andNot(BigInteger val) {
+ return Logical.andNot(this, val);
+ }
+
+ public int intValue() {
+ return (sign * digits[0]);
+ }
+
+ public long longValue() {
+ long value = (numberLength > 1) ? (((long) digits[1]) << 32)
+ | (digits[0] & 0xFFFFFFFFL) : (digits[0] & 0xFFFFFFFFL);
+ return (sign * value);
+ }
+
+ public float floatValue() {
+ return (float) doubleValue();
+ }
+
+ public double doubleValue() {
+ return Conversion.bigInteger2Double(this);
+ }
+
+ public int compareTo(BigInteger val) {
+ if (sign > val.sign) {
+ return GREATER;
+ }
+ if (sign < val.sign) {
+ return LESS;
+ }
+ if (numberLength > val.numberLength) {
+ return sign;
+ }
+ if (numberLength < val.numberLength) {
+ return -val.sign;
+ }
+ // Equal sign and equal numberLength
+ return (sign * Elementary.compareArrays(digits, val.digits,
+ numberLength));
+ }
+
+ public BigInteger min(BigInteger val) {
+ return ((this.compareTo(val) == LESS) ? this : val);
+ }
+
+ public BigInteger max(BigInteger val) {
+ return ((this.compareTo(val) == GREATER) ? this : val);
+ }
+
+ public override int GetHashCode() {
+ if (hashCode != 0) {
+ return hashCode;
+ }
+ for (int i = 0; i < digits.Length; i++) {
+ hashCode = (hashCode * 33 + (int)(digits[i] & 0xffffffff));
+ }
+ hashCode = hashCode * sign;
+ return hashCode;
+ }
+
+ public override bool Equals(object x) {
+ if (this == x) {
+ return true;
+ }
+ if (x is BigInteger) {
+ BigInteger x1 = (BigInteger) x;
+ return sign == x1.sign && numberLength == x1.numberLength
+ && equalsArrays(x1.digits);
+ }
+ return false;
+ }
+
+ bool equalsArrays(int[] b) {
+ int i;
+ for (i = numberLength - 1; (i >= 0) && (digits[i] == b[i]); i--) {
+ // Empty
+ }
+ return i < 0;
+ }
+
+ public override string ToString(){
+ return Conversion.toDecimalScaledString(this, 0);
+ }
+
+ public String toString(int radix) {
+ return Conversion.bigInteger2String(this, radix);
+ }
+
+ public BigInteger gcd(BigInteger val) {
+ BigInteger val1 = this.abs();
+ BigInteger val2 = val.abs();
+ // To avoid a possible division by zero
+ if (val1.signum() == 0) {
+ return val2;
+ } else if (val2.signum() == 0) {
+ return val1;
+ }
+
+ // Optimization for small operands
+ // (op2.bitLength() < 64) and (op1.bitLength() < 64)
+ if (((val1.numberLength == 1) || ((val1.numberLength == 2) && (val1.digits[1] > 0)))
+ && (val2.numberLength == 1 || (val2.numberLength == 2 && val2.digits[1] > 0))) {
+ return BigInteger.valueOf(Division.gcdBinary(val1.longValue(), val2
+ .longValue()));
+ }
+
+ return Division.gcdBinary(val1.copy(), val2.copy());
+
+ }
+
+ public BigInteger multiply(BigInteger val) {
+ // This let us to throw NullReferenceException when val == null
+ if (val.sign == 0) {
+ return ZERO;
+ }
+ if (sign == 0) {
+ return ZERO;
+ }
+ return Multiplication.multiply(this, val);
+ }
+
+ public BigInteger pow(int exp) {
+ if (exp < 0) {
+ // math.16=Negative exponent
+ throw new ArithmeticException("Negative exponent");
+ }
+ if (exp == 0) {
+ return ONE;
+ } else if (exp == 1 || Equals(ONE) || Equals(ZERO)) {
+ return this;
+ }
+
+ // if even take out 2^x factor which we can
+ // calculate by shifting.
+ if (!testBit(0)) {
+ int x = 1;
+ while (!testBit(x)) {
+ x++;
+ }
+ return getPowerOfTwo(x*exp).multiply(this.shiftRight(x).pow(exp));
+ }
+ return Multiplication.pow(this, exp);
+ }
+
+ public BigInteger[] divideAndRemainder(BigInteger divisor) {
+ int divisorSign = divisor.sign;
+ if (divisorSign == 0) {
+ throw new ArithmeticException("BigInteger divide by zero");
+ }
+ int divisorLen = divisor.numberLength;
+ int[] divisorDigits = divisor.digits;
+ if (divisorLen == 1) {
+ return Division.divideAndRemainderByInteger(this, divisorDigits[0],
+ divisorSign);
+ }
+
+
+ // res[0] is a quotient and res[1] is a remainder:
+ int[] thisDigits = digits;
+ int thisLen = numberLength;
+ int cmp = (thisLen != divisorLen) ? ((thisLen > divisorLen) ? 1 : -1)
+ : Elementary.compareArrays(thisDigits, divisorDigits, thisLen);
+ if (cmp < 0) {
+ return new BigInteger[] { ZERO, this };
+ }
+ int thisSign = sign;
+ int quotientLength = thisLen - divisorLen + 1;
+ int remainderLength = divisorLen;
+ int quotientSign = ((thisSign == divisorSign) ? 1 : -1);
+ int[] quotientDigits = new int[quotientLength];
+ int[] remainderDigits = Division.divide(quotientDigits, quotientLength,
+ thisDigits, thisLen, divisorDigits, divisorLen);
+ BigInteger result0 = new BigInteger(quotientSign, quotientLength,
+ quotientDigits);
+ BigInteger result1 = new BigInteger(thisSign, remainderLength,
+ remainderDigits);
+ result0.cutOffLeadingZeroes();
+ result1.cutOffLeadingZeroes();
+ return new BigInteger[] { result0, result1 };
+ }
+
+ public BigInteger divide(BigInteger divisor) {
+ if (divisor.sign == 0) {
+ throw new ArithmeticException("BigInteger divide by zero");
+ }
+ int divisorSign = divisor.sign;
+ if (divisor.isOne()) {
+ return ((divisor.sign > 0) ? this : this.negate());
+ }
+ int thisSign = sign;
+ int thisLen = numberLength;
+ int divisorLen = divisor.numberLength;
+ if (thisLen + divisorLen == 2) {
+ long val = (digits[0] & 0xFFFFFFFFL)
+ / (divisor.digits[0] & 0xFFFFFFFFL);
+ if (thisSign != divisorSign) {
+ val = -val;
+ }
+ return valueOf(val);
+ }
+ int cmp = ((thisLen != divisorLen) ? ((thisLen > divisorLen) ? 1 : -1)
+ : Elementary.compareArrays(digits, divisor.digits, thisLen));
+ if (cmp == EQUALS) {
+ return ((thisSign == divisorSign) ? ONE : MINUS_ONE);
+ }
+ if (cmp == LESS) {
+ return ZERO;
+ }
+ int resLength = thisLen - divisorLen + 1;
+ int[] resDigits = new int[resLength];
+ int resSign = ((thisSign == divisorSign) ? 1 : -1);
+ if (divisorLen == 1) {
+ Division.divideArrayByInt(resDigits, digits, thisLen,
+ divisor.digits[0]);
+ } else {
+ Division.divide(resDigits, resLength, digits, thisLen,
+ divisor.digits, divisorLen);
+ }
+ BigInteger result = new BigInteger(resSign, resLength, resDigits);
+ result.cutOffLeadingZeroes();
+ return result;
+ }
+
+ public BigInteger remainder(BigInteger divisor) {
+ if (divisor.sign == 0) {
+ throw new ArithmeticException("BigInteger divide by zero");
+ }
+ int thisLen = numberLength;
+ int divisorLen = divisor.numberLength;
+ if (((thisLen != divisorLen) ? ((thisLen > divisorLen) ? 1 : -1)
+ : Elementary.compareArrays(digits, divisor.digits, thisLen)) == LESS) {
+ return this;
+ }
+ int resLength = divisorLen;
+ int[] resDigits = new int[resLength];
+ if (resLength == 1) {
+ resDigits[0] = Division.remainderArrayByInt(digits, thisLen,
+ divisor.digits[0]);
+ } else {
+ int qLen = thisLen - divisorLen + 1;
+ resDigits = Division.divide(null, qLen, digits, thisLen,
+ divisor.digits, divisorLen);
+ }
+ BigInteger result = new BigInteger(sign, resLength, resDigits);
+ result.cutOffLeadingZeroes();
+ return result;
+ }
+
+ public BigInteger modInverse(BigInteger m) {
+ if (m.sign <= 0) {
+ throw new ArithmeticException("BigInteger: modulus not positive");
+ }
+ // If both are even, no inverse exists
+ if (!(testBit(0) || m.testBit(0))) {
+ throw new ArithmeticException("BigInteger not invertible.");
+ }
+ if (m.isOne()) {
+ return ZERO;
+ }
+
+ // From now on: (m > 1)
+ BigInteger res = Division.modInverseMontgomery(abs().mod(m), m);
+ if (res.sign == 0) {
+ throw new ArithmeticException("BigInteger not invertible.");
+ }
+
+ res = ((sign < 0) ? m.subtract(res) : res);
+ return res;
+
+ }
+
+ public BigInteger modPow(BigInteger exponent, BigInteger m) {
+ if (m.sign <= 0) {
+ throw new ArithmeticException("BigInteger: modulus not positive");
+ }
+ BigInteger _base = this;
+
+ if (m.isOne() | (exponent.sign > 0 & _base.sign == 0)) {
+ return BigInteger.ZERO;
+ }
+ if (_base.sign == 0 && exponent.sign == 0) {
+ return BigInteger.ONE;
+ }
+ if (exponent.sign < 0) {
+ _base = modInverse(m);
+ exponent = exponent.negate();
+ }
+ // From now on: (m > 0) and (exponent >= 0)
+ BigInteger res = (m.testBit(0)) ? Division.oddModPow(_base.abs(),
+ exponent, m) : Division.evenModPow(_base.abs(), exponent, m);
+ if ((_base.sign < 0) && exponent.testBit(0)) {
+ // -b^e mod m == ((-1 mod m) * (b^e mod m)) mod m
+ res = m.subtract(BigInteger.ONE).multiply(res).mod(m);
+ }
+ // else exponent is even, so base^exp is positive
+ return res;
+ }
+
+ public BigInteger mod(BigInteger m) {
+ if (m.sign <= 0) {
+ throw new ArithmeticException("BigInteger: modulus not positive");
+ }
+ BigInteger rem = remainder(m);
+ return ((rem.sign < 0) ? rem.add(m) : rem);
+ }
+
+ public bool isProbablePrime(int certainty) {
+ return Primality.isProbablePrime(abs(), certainty);
+ }
+
+ public BigInteger nextProbablePrime() {
+ if (sign < 0) {
+ throw new ArithmeticException("start < 0: " + this);
+ }
+ return Primality.nextProbablePrime(this);
+ }
+
+ public static BigInteger probablePrime(int bitLength, Random rnd) {
+ return new BigInteger(bitLength, 100, rnd);
+ }
+
+ internal void cutOffLeadingZeroes() {
+ while ((numberLength > 0) && (digits[--numberLength] == 0)) {
+ // Empty
+ }
+ if (digits[numberLength++] == 0) {
+ sign = 0;
+ }
+ }
+
+ internal bool isOne() {
+ return ((numberLength == 1) && (digits[0] == 1));
+ }
+
+ private void putBytesPositiveToIntegers(byte[] byteValues) {
+ int bytesLen = byteValues.Length;
+ int highBytes = bytesLen & 3;
+ numberLength = (bytesLen >> 2) + ((highBytes == 0) ? 0 : 1);
+ digits = new int[numberLength];
+ int i = 0;
+ // Put bytes to the int array starting from the end of the byte array
+ while (bytesLen > highBytes) {
+ digits[i++] = (byteValues[--bytesLen] & 0xFF)
+ | (byteValues[--bytesLen] & 0xFF) << 8
+ | (byteValues[--bytesLen] & 0xFF) << 16
+ | (byteValues[--bytesLen] & 0xFF) << 24;
+ }
+ // Put the first bytes in the highest element of the int array
+ for (int j = 0; j < bytesLen; j++) {
+ digits[i] = (digits[i] << 8) | (byteValues[j] & 0xFF);
+ }
+ }
+
+ private void putBytesNegativeToIntegers(byte[] byteValues) {
+ int bytesLen = byteValues.Length;
+ int highBytes = bytesLen & 3;
+ numberLength = (bytesLen >> 2) + ((highBytes == 0) ? 0 : 1);
+ digits = new int[numberLength];
+ int i = 0;
+ // Setting the sign
+ digits[numberLength - 1] = -1;
+ // Put bytes to the int array starting from the end of the byte array
+ while (bytesLen > highBytes) {
+ digits[i] = (byteValues[--bytesLen] & 0xFF)
+ | (byteValues[--bytesLen] & 0xFF) << 8
+ | (byteValues[--bytesLen] & 0xFF) << 16
+ | (byteValues[--bytesLen] & 0xFF) << 24;
+ if (digits[i] != 0) {
+ digits[i] = -digits[i];
+ firstNonzeroDigit = i;
+ i++;
+ while (bytesLen > highBytes) {
+ digits[i] = (byteValues[--bytesLen] & 0xFF)
+ | (byteValues[--bytesLen] & 0xFF) << 8
+ | (byteValues[--bytesLen] & 0xFF) << 16
+ | (byteValues[--bytesLen] & 0xFF) << 24;
+ digits[i] = ~digits[i];
+ i++;
+ }
+ break;
+ }
+ i++;
+ }
+ if (highBytes != 0) {
+ // Put the first bytes in the highest element of the int array
+ if (firstNonzeroDigit != -2) {
+ for (int j = 0; j < bytesLen; j++) {
+ digits[i] = (digits[i] << 8) | (byteValues[j] & 0xFF);
+ }
+ digits[i] = ~digits[i];
+ } else {
+ for (int j = 0; j < bytesLen; j++) {
+ digits[i] = (digits[i] << 8) | (byteValues[j] & 0xFF);
+ }
+ digits[i] = -digits[i];
+ }
+ }
+ }
+
+ internal int getFirstNonzeroDigit() {
+ if (firstNonzeroDigit == -2) {
+ int i;
+ if (this.sign == 0) {
+ i = -1;
+ } else {
+ for (i = 0; digits[i] == 0; i++) {
+ // Empty
+ }
+ }
+ firstNonzeroDigit = i;
+ }
+ return firstNonzeroDigit;
+ }
+
+ /*
+ * Returns a copy of the current instance to achieve immutability
+ */
+ internal BigInteger copy() {
+ int[] copyDigits = new int[numberLength];
+ Array.Copy(digits, copyDigits, numberLength);
+ return new BigInteger(sign, numberLength, copyDigits);
+ }
+
+ internal void unCache() {
+ firstNonzeroDigit = -2;
+ }
+
+ internal static BigInteger getPowerOfTwo(int exp) {
+ if(exp < TWO_POWS.Length) {
+ return TWO_POWS[exp];
+ }
+ int intCount = exp >> 5;
+ int bitN = exp & 31;
+ int[] resDigits = new int[intCount+1];
+ resDigits[intCount] = 1 << bitN;
+ return new BigInteger(1, intCount+1, resDigits);
+ }
+ }
+}
View
60 src/ikc/main/Ioke.Math/BigSquareRoot.cs
@@ -0,0 +1,60 @@
+
+namespace Ioke.Math {
+ using System;
+
+ public class BigSquareRoot {
+ private static readonly BigDecimal ZERO = BigDecimal.ZERO;
+ private static readonly BigDecimal ONE = BigDecimal.ONE;
+ private static readonly BigDecimal TWO = new BigDecimal ("2");
+ public const int DEFAULT_MAX_ITERATIONS = 50;
+ public const int DEFAULT_SCALE = 10;
+
+ private BigDecimal error;
+ private readonly int scale;
+ private readonly int maxIterations;
+
+ public BigSquareRoot() : this(DEFAULT_MAX_ITERATIONS, DEFAULT_SCALE) {}
+
+ public BigSquareRoot(int maxIterations, int scale) {
+ this.maxIterations = maxIterations;
+ this.scale = scale;
+ }
+
+ public BigDecimal Get(BigDecimal n) {
+ if(n.CompareTo(ZERO) <= 0) {
+ throw new System.ArgumentException();
+ }
+
+ BigDecimal initialGuess = GetInitialApproximation(n);
+ BigDecimal lastGuess = ZERO;
+ BigDecimal guess = new BigDecimal(initialGuess.ToString());
+
+ int iterations = 0;
+ bool more = true;
+ while(more) {
+ lastGuess = guess;
+ guess = n.divide(guess, scale, BigDecimal.ROUND_HALF_UP);
+ guess = guess.add(lastGuess);
+ guess = guess.divide(TWO, scale, BigDecimal.ROUND_HALF_UP);
+ error = n.subtract(guess.multiply(guess));
+ if(++iterations >= maxIterations) {
+ more = false;
+ } else if(lastGuess.Equals(guess)) {
+ more = error.abs().CompareTo(ONE) >= 0;
+ }
+ }
+ return guess;
+ }
+
+ private static BigDecimal GetInitialApproximation(BigDecimal n) {
+ BigInteger integerPart = n.toBigInteger();
+ int length = integerPart.ToString().Length;
+ if((length % 2) == 0) {
+ length--;
+ }
+ length /= 2;
+ BigDecimal guess = ONE.movePointRight(length);
+ return guess;
+ }
+ }
+}
View
257 src/ikc/main/Ioke.Math/BitLevel.cs
@@ -0,0 +1,257 @@
+
+namespace Ioke.Math {
+ using System;
+ using System.Text;
+ class BitLevel {
+ private BitLevel() {}
+
+ internal static int bitLength(BigInteger val) {
+ if (val.sign == 0) {
+ return 0;
+ }
+ int bLength = (val.numberLength << 5);
+ int highDigit = val.digits[val.numberLength - 1];
+
+ if (val.sign < 0) {
+ int i = val.getFirstNonzeroDigit();
+ // We reduce the problem to the positive case.
+ if (i == val.numberLength - 1) {
+ highDigit--;
+ }
+ }
+ // Subtracting all sign bits
+ bLength -= BigDecimal.numberOfLeadingZeros(highDigit);
+ return bLength;
+ }
+
+ internal static int bitCount(BigInteger val) {
+ int bCount = 0;
+
+ if (val.sign == 0) {
+ return 0;
+ }
+
+ int i = val.getFirstNonzeroDigit();;
+ if (val.sign > 0) {
+ for ( ; i < val.numberLength; i++) {
+ bCount += BigDecimal.bitCount(val.digits[i]);
+ }
+ } else {// (sign < 0)
+ // this digit absorbs the carry
+ bCount += BigDecimal.bitCount(-val.digits[i]);
+ for (i++; i < val.numberLength; i++) {
+ bCount += BigDecimal.bitCount(~val.digits[i]);
+ }
+ // We take the complement sum:
+ bCount = (val.numberLength << 5) - bCount;
+ }
+ return bCount;
+ }
+
+ internal static bool testBit(BigInteger val, int n) {
+ // PRE: 0 <= n < val.bitLength()
+ return ((val.digits[n >> 5] & (1 << (n & 31))) != 0);
+ }
+
+ internal static bool nonZeroDroppedBits(int numberOfBits, int[] digits) {
+ int intCount = numberOfBits >> 5;
+ int bitCount = numberOfBits & 31;
+ int i;
+
+ for (i = 0; (i < intCount) && (digits[i] == 0); i++) {
+ ;
+ }
+ return ((i != intCount) || (digits[i] << (32 - bitCount) != 0));
+ }
+
+ internal static BigInteger shiftLeft(BigInteger source, int count) {
+ int intCount = count >> 5;
+ count &= 31; // %= 32
+ int resLength = source.numberLength + intCount
+ + ( ( count == 0 ) ? 0 : 1 );
+ int[] resDigits = new int[resLength];
+
+ shiftLeft(resDigits, source.digits, intCount, count);
+ BigInteger result = new BigInteger(source.sign, resLength, resDigits);
+ result.cutOffLeadingZeroes();
+ return result;
+ }
+
+ internal static void inplaceShiftLeft(BigInteger val, int count) {
+ int intCount = count >> 5; // count of integers
+ val.numberLength += intCount
+ + ( BigDecimal
+ .numberOfLeadingZeros(val.digits[val.numberLength - 1])
+ - ( count & 31 ) >= 0 ? 0 : 1 );
+ shiftLeft(val.digits, val.digits, intCount, count & 31);
+ val.cutOffLeadingZeroes();
+ val.unCache();
+ }
+
+ internal static void shiftLeft(int[] result, int[] source, int intCount, int count) {
+ if (count == 0) {
+ Array.Copy(source, 0, result, intCount, result.Length
+ - intCount);
+ } else {
+ int rightShiftCount = 32 - count;
+
+ result[result.Length - 1] = 0;
+ for (int i = result.Length - 1; i > intCount; i--) {
+ result[i] |= (int)(((uint)source[i - intCount - 1]) >> rightShiftCount);
+ result[i - 1] = source[i - intCount - 1] << count;
+ }
+ }
+
+ for (int i = 0; i < intCount; i++) {
+ result[i] = 0;
+ }
+ }
+
+ internal static void shiftLeftOneBit(int[] result, int[] source, int srcLen) {
+ int carry = 0;
+ for (int i = 0; i < srcLen; i++) {
+ int val = source[i];
+ result[i] = (val << 1) | carry;
+ carry = (int)(((uint)val) >> 31);
+ }
+ if (carry != 0) {
+ result[srcLen] = carry;
+ }
+ }
+
+ internal static BigInteger shiftLeftOneBit(BigInteger source) {
+ int srcLen = source.numberLength;
+ int resLen = srcLen + 1;
+ int[] resDigits = new int[resLen];
+ shiftLeftOneBit(resDigits, source.digits, srcLen);
+ BigInteger result = new BigInteger(source.sign, resLen, resDigits);
+ result.cutOffLeadingZeroes();
+ return result;
+ }
+
+ internal static BigInteger shiftRight(BigInteger source, int count) {
+ int intCount = count >> 5; // count of integers
+ count &= 31; // count of remaining bits
+ if (intCount >= source.numberLength) {
+ return ((source.sign < 0) ? BigInteger.MINUS_ONE : BigInteger.ZERO);
+ }
+ int i;
+ int resLength = source.numberLength - intCount;
+ int[] resDigits = new int[resLength + 1];
+
+ shiftRight(resDigits, resLength, source.digits, intCount, count);
+ if (source.sign < 0) {
+ // Checking if the dropped bits are zeros (the remainder equals to
+ // 0)
+ for (i = 0; (i < intCount) && (source.digits[i] == 0); i++) {
+ ;
+ }
+ // If the remainder is not zero, add 1 to the result
+ if ((i < intCount)
+ || ((count > 0) && ((source.digits[i] << (32 - count)) != 0))) {
+ for (i = 0; (i < resLength) && (resDigits[i] == -1); i++) {
+ resDigits[i] = 0;
+ }
+ if (i == resLength) {
+ resLength++;
+ }
+ resDigits[i]++;
+ }
+ }
+ BigInteger result = new BigInteger(source.sign, resLength, resDigits);
+ result.cutOffLeadingZeroes();
+ return result;
+ }
+
+ internal static void inplaceShiftRight(BigInteger val, int count) {
+ int sign = val.signum();
+ if (count == 0 || val.signum() == 0)
+ return;
+ int intCount = count >> 5; // count of integers
+ val.numberLength -= intCount;
+ if (!shiftRight(val.digits, val.numberLength, val.digits, intCount,
+ count & 31)
+ && sign < 0) {
+ // remainder not zero: add one to the result
+ int i;
+ for (i = 0; ( i < val.numberLength ) && ( val.digits[i] == -1 ); i++) {
+ val.digits[i] = 0;
+ }
+ if (i == val.numberLength) {
+ val.numberLength++;
+ }
+ val.digits[i]++;
+ }
+ val.cutOffLeadingZeroes();
+ val.unCache();
+ }
+
+ internal static bool shiftRight(int[] result, int resultLen, int[] source,
+ int intCount, int count) {
+ int i;
+ bool allZero = true;
+ for (i = 0; i < intCount; i++)
+ allZero &= source[i] == 0;
+ if (count == 0) {
+ Array.Copy(source, intCount, result, 0, resultLen);
+ i = resultLen;
+ } else {
+ int leftShiftCount = 32 - count;
+
+ allZero &= ( source[i] << leftShiftCount ) == 0;
+ for (i = 0; i < resultLen - 1; i++) {
+ result[i] = ( (int)(((uint)source[i + intCount]) >> count) )
+ | ( source[i + intCount + 1] << leftShiftCount );
+ }
+ result[i] = (int)(((uint)( source[i + intCount]) >> count) );
+ i++;
+ }
+
+ return allZero;
+ }
+
+ internal static BigInteger flipBit(BigInteger val, int n){
+ int resSign = (val.sign == 0) ? 1 : val.sign;
+ int intCount = n >> 5;
+ int bitN = n & 31;
+ int resLength = Math.Max(intCount + 1, val.numberLength) + 1;
+ int[] resDigits = new int[resLength];
+ int i;
+
+ int bitNumber = 1 << bitN;
+ Array.Copy(val.digits, 0, resDigits, 0, val.numberLength);
+
+ if (val.sign < 0) {
+ if (intCount >= val.numberLength) {
+ resDigits[intCount] = bitNumber;
+ } else {
+ //val.sign<0 y intCount < val.numberLength
+ int firstNonZeroDigit = val.getFirstNonzeroDigit();
+ if (intCount > firstNonZeroDigit) {
+ resDigits[intCount] ^= bitNumber;
+ } else if (intCount < firstNonZeroDigit) {
+ resDigits[intCount] = -bitNumber;
+ for (i=intCount + 1; i < firstNonZeroDigit; i++) {
+ resDigits[i]=-1;
+ }
+ resDigits[i] = resDigits[i]--;
+ } else {
+ i = intCount;
+ resDigits[i] = -((-resDigits[intCount]) ^ bitNumber);
+ if (resDigits[i] == 0) {
+ for (i++; resDigits[i] == -1 ; i++) {
+ resDigits[i] = 0;
+ }
+ resDigits[i]++;
+ }
+ }
+ }
+ } else {//case where val is positive
+ resDigits[intCount] ^= bitNumber;
+ }
+ BigInteger result = new BigInteger(resSign, resLength, resDigits);
+ result.cutOffLeadingZeroes();
+ return result;
+ }
+ }
+}
View
393 src/ikc/main/Ioke.Math/Conversion.cs
@@ -0,0 +1,393 @@
+namespace Ioke.Math {
+ using System;
+ using System.Text;
+
+ class Conversion {
+ private Conversion() {}
+ internal static readonly int[] digitFitInInt = { -1, -1, 31, 19, 15, 13, 11,
+ 11, 10, 9, 9, 8, 8, 8, 8, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 6,
+ 6, 6, 6, 6, 6, 6, 6, 5 };
+
+ internal static readonly int[] bigRadices = { -2147483648, 1162261467,
+ 1073741824, 1220703125, 362797056, 1977326743, 1073741824,
+ 387420489, 1000000000, 214358881, 429981696, 815730721, 1475789056,
+ 170859375, 268435456, 410338673, 612220032, 893871739, 1280000000,
+ 1801088541, 113379904, 148035889, 191102976, 244140625, 308915776,
+ 387420489, 481890304, 594823321, 729000000, 887503681, 1073741824,
+ 1291467969, 1544804416, 1838265625, 60466176 };
+
+
+ internal static String bigInteger2String(BigInteger val, int radix) {
+ int sign = val.sign;
+ int numberLength = val.numberLength;
+ int[] digits = val.digits;
+
+ if (sign == 0) {
+ return "0"; //$NON-NLS-1$
+ }
+ if (numberLength == 1) {
+ int highDigit = digits[numberLength - 1];
+ long v = highDigit & 0xFFFFFFFFL;
+ if (sign < 0) {
+ v = -v;
+ }
+ return Convert.ToString(v, radix);
+ }
+ if ((radix == 10) || (radix < 0)
+ || (radix > 26)) {
+ return val.ToString();
+ }
+ double bitsForRadixDigit;
+ bitsForRadixDigit = Math.Log(radix) / Math.Log(2);
+ int resLengthInChars = (int) (val.abs().bitLength() / bitsForRadixDigit + ((sign < 0) ? 1
+ : 0)) + 1;
+
+ char[] result = new char[resLengthInChars];
+ int currentChar = resLengthInChars;
+ int resDigit;
+ if (radix != 16) {
+ int[] temp = new int[numberLength];
+ Array.Copy(digits, 0, temp, 0, numberLength);
+ int tempLen = numberLength;
+ int charsPerInt = digitFitInInt[radix];
+ int i;
+ // get the maximal power of radix that fits in int
+ int bigRadix = bigRadices[radix - 2];
+ while (true) {
+ // divide the array of digits by bigRadix and convert remainders
+ // to characters collecting them in the char array
+ resDigit = Division.divideArrayByInt(temp, temp, tempLen,
+ bigRadix);
+ int previous = currentChar;
+ do {
+ result[--currentChar] = Convert.ToString(resDigit % radix, radix)[0];
+ } while (((resDigit /= radix) != 0) && (currentChar != 0));
+ int delta = charsPerInt - previous + currentChar;
+ for (i = 0; i < delta && currentChar > 0; i++) {
+ result[--currentChar] = '0';
+ }
+ for (i = tempLen - 1; (i > 0) && (temp[i] == 0); i--) {
+ ;
+ }
+ tempLen = i + 1;
+ if ((tempLen == 1) && (temp[0] == 0)) { // the quotient is 0
+ break;
+ }
+ }
+ } else {
+ // radix == 16
+ for (int i = 0; i < numberLength; i++) {
+ for (int j = 0; (j < 8) && (currentChar > 0); j++) {
+ resDigit = digits[i] >> (j << 2) & 0xf;
+ result[--currentChar] = Convert.ToString(resDigit % radix, 16)[0];
+ }
+ }
+ }
+ while (result[currentChar] == '0') {
+ currentChar++;
+ }
+ if (sign == -1) {
+ result[--currentChar] = '-';
+ }
+ return new String(result, currentChar, resLengthInChars - currentChar);
+ }
+
+ internal static String toDecimalScaledString(BigInteger val, int scale) {
+ int sign = val.sign;
+ int numberLength = val.numberLength;
+ int[] digits = val.digits;
+ int resLengthInChars;
+ int currentChar;
+ char[] result;
+
+ if (sign == 0) {
+ switch (scale) {
+ case 0:
+ return "0"; //$NON-NLS-1$
+ case 1:
+ return "0.0"; //$NON-NLS-1$
+ case 2:
+ return "0.00"; //$NON-NLS-1$
+ case 3:
+ return "0.000"; //$NON-NLS-1$
+ case 4:
+ return "0.0000"; //$NON-NLS-1$
+ case 5:
+ return "0.00000"; //$NON-NLS-1$
+ case 6:
+ return "0.000000"; //$NON-NLS-1$
+ default:
+ StringBuilder result1x = new StringBuilder();
+ if (scale < 0) {
+ result1x.Append("0E+"); //$NON-NLS-1$
+ } else {
+ result1x.Append("0E"); //$NON-NLS-1$
+ }
+ result1x.Append(-scale);
+ return result1x.ToString();
+ }
+ }
+ // one 32-bit unsigned value may contains 10 decimal digits
+ resLengthInChars = numberLength * 10 + 1 + 7;
+ // Explanation why +1+7:
+ // +1 - one char for sign if needed.
+ // +7 - For "special case 2" (see below) we have 7 free chars for
+ // inserting necessary scaled digits.
+ result = new char[resLengthInChars + 1];
+ // allocated [resLengthInChars+1] characters.
+ // a free latest character may be used for "special case 1" (see
+ // below)
+ currentChar = resLengthInChars;
+ if (numberLength == 1) {
+ int highDigit = digits[0];
+ if (highDigit < 0) {
+ long v = highDigit & 0xFFFFFFFFL;
+ do {
+ long prev = v;
+ v /= 10;
+ result[--currentChar] = (char) (0x0030 + ((int) (prev - v * 10)));
+ } while (v != 0);
+ } else {
+ int v = highDigit;
+ do {
+ int prev = v;
+ v /= 10;
+ result[--currentChar] = (char) (0x0030 + (prev - v * 10));
+ } while (v != 0);
+ }
+ } else {
+ int[] temp = new int[numberLength];
+ int tempLen = numberLength;
+ Array.Copy(digits, temp, tempLen);
+ while (true) {
+ // divide the array of digits by bigRadix and convert
+ // remainders
+ // to characters collecting them in the char array
+ long result11 = 0;
+ for (int i1 = tempLen - 1; i1 >= 0; i1--) {
+ long temp1 = (result11 << 32)
+ + (temp[i1] & 0xFFFFFFFFL);
+ long res = divideLongByBillion(temp1);
+ temp[i1] = (int) res;
+ result11 = (int) (res >> 32);
+ }
+ int resDigit = (int) result11;
+ int previous = currentChar;
+ do {
+ result[--currentChar] = (char) (0x0030 + (resDigit % 10));
+ } while (((resDigit /= 10) != 0) && (currentChar != 0));
+ int delta = 9 - previous + currentChar;
+ for (int i = 0; (i < delta) && (currentChar > 0); i++) {
+ result[--currentChar] = '0';
+ }
+ int j = tempLen - 1;
+ bool bigloop = false;
+ for (; temp[j] == 0; j--) {
+ if (j == 0) { // means temp[0] == 0
+ bigloop = true;
+ break;
+ }
+ }
+ if(bigloop) break;
+ tempLen = j + 1;
+ }
+ while (result[currentChar] == '0') {
+ currentChar++;
+ }
+ }
+ bool negNumber = (sign < 0);
+ int exponent = resLengthInChars - currentChar - scale - 1;
+ if (scale == 0) {
+ if (negNumber) {
+ result[--currentChar] = '-';
+ }
+ return new String(result, currentChar, resLengthInChars
+ - currentChar);
+ }
+ if ((scale > 0) && (exponent >= -6)) {
+ if (exponent >= 0) {
+ // special case 1
+ int insertPoint = currentChar + exponent;
+ for (int j = resLengthInChars - 1; j >= insertPoint; j--) {
+ result[j + 1] = result[j];
+ }
+ result[++insertPoint] = '.';
+ if (negNumber) {
+ result[--currentChar] = '-';
+ }
+ return new String(result, currentChar, resLengthInChars
+ - currentChar + 1);
+ }
+ // special case 2
+ for (int j = 2; j < -exponent + 1; j++) {
+ result[--currentChar] = '0';
+ }
+ result[--currentChar] = '.';
+ result[--currentChar] = '0';
+ if (negNumber) {
+ result[--currentChar] = '-';
+ }
+ return new String(result, currentChar, resLengthInChars
+ - currentChar);
+ }
+ int startPoint = currentChar + 1;
+ int endPoint = resLengthInChars;
+ StringBuilder result1 = new StringBuilder(16 + endPoint - startPoint);
+ if (negNumber) {
+ result1.Append('-');
+ }
+ if (endPoint - startPoint >= 1) {
+ result1.Append(result[currentChar]);
+ result1.Append('.');
+ result1.Append(result, currentChar + 1, resLengthInChars
+ - currentChar - 1);
+ } else {
+ result1.Append(result, currentChar, resLengthInChars
+ - currentChar);
+ }
+ result1.Append('E');
+ if (exponent > 0) {
+ result1.Append('+');
+ }
+ result1.Append(exponent.ToString());
+ return result1.ToString();
+ }
+
+ /* can process only 32-bit numbers */
+ internal static String toDecimalScaledString(long value, int scale) {
+ int resLengthInChars;
+ int currentChar;
+ char[] result;
+ bool negNumber = value < 0;
+ if(negNumber) {
+ value = -value;
+ }
+ if (value == 0) {
+ switch (scale) {
+ case 0: return "0"; //$NON-NLS-1$
+ case 1: return "0.0"; //$NON-NLS-1$
+ case 2: return "0.00"; //$NON-NLS-1$
+ case 3: return "0.000"; //$NON-NLS-1$
+ case 4: return "0.0000"; //$NON-NLS-1$
+ case 5: return "0.00000"; //$NON-NLS-1$
+ case 6: return "0.000000"; //$NON-NLS-1$
+ default:
+ StringBuilder result1x = new StringBuilder();
+ if (scale < 0) {
+ result1x.Append("0E+");
+ } else {
+ result1x.Append("0E");
+ }
+ result1x.Append( (scale == Int32.MinValue) ? "2147483648" : (-scale).ToString());
+ return result1x.ToString();
+ }
+ }
+ // one 32-bit unsigned value may contains 10 decimal digits
+ resLengthInChars = 18;
+ // Explanation why +1+7:
+ // +1 - one char for sign if needed.
+ // +7 - For "special case 2" (see below) we have 7 free chars for
+ // inserting necessary scaled digits.
+ result = new char[resLengthInChars+1];
+ // Allocated [resLengthInChars+1] characters.
+ // a free latest character may be used for "special case 1" (see below)
+ currentChar = resLengthInChars;
+ long v = value;
+ do {
+ long prev = v;
+ v /= 10;
+ result[--currentChar] = (char) (0x0030 + (prev - v * 10));
+ } while (v != 0);
+
+ long exponent = (long)resLengthInChars - (long)currentChar - scale - 1L;
+ if (scale == 0) {
+ if (negNumber) {
+ result[--currentChar] = '-';
+ }
+ return new String(result, currentChar, resLengthInChars - currentChar);
+ }
+ if (scale > 0 && exponent >= -6) {
+ if (exponent >= 0) {
+ // special case 1
+ int insertPoint = currentChar + (int) exponent ;
+ for(int j=resLengthInChars-1; j>=insertPoint; j--) {
+ result[j+1] = result[j];
+ }
+ result[++insertPoint]='.';
+ if (negNumber) {
+ result[--currentChar] = '-';
+ }
+ return new String(result, currentChar, resLengthInChars - currentChar + 1);
+ }
+ // special case 2
+ for (int j = 2; j < -exponent + 1; j++) {
+ result[--currentChar] = '0';
+ }
+ result[--currentChar] = '.';
+ result[--currentChar] = '0';
+ if (negNumber) {
+ result[--currentChar] = '-';
+ }
+ return new String(result, currentChar, resLengthInChars - currentChar);
+ }
+ int startPoint = currentChar + 1;
+ int endPoint = resLengthInChars;
+ StringBuilder result1 = new StringBuilder(16+endPoint-startPoint);
+ if (negNumber) {
+ result1.Append('-');
+ }
+ if (endPoint - startPoint >= 1) {
+ result1.Append(result[currentChar]);
+ result1.Append('.');
+ result1.Append(result,currentChar+1,resLengthInChars - currentChar-1);
+ } else {
+ result1.Append(result,currentChar,resLengthInChars - currentChar);
+ }
+ result1.Append('E');
+ if (exponent > 0) {
+ result1.Append('+');
+ }
+ result1.Append(exponent.ToString());
+ return result1.ToString();
+ }
+
+ internal static long divideLongByBillion(long a) {
+ long quot;
+ long rem;
+
+ if (a >= 0) {
+ long bLong = 1000000000L;
+ quot = (a / bLong);
+ rem = (a % bLong);
+ } else {
+ /*
+ * Make the dividend positive shifting it right by 1 bit then get
+ * the quotient an remainder and correct them properly
+ */
+ long aPos = (long)(((ulong)a) >> 1);
+ long bPos = (long)(((ulong)1000000000L) >> 1);
+ quot = aPos / bPos;
+ rem = aPos % bPos;
+ // double the remainder and add 1 if 'a' is odd
+ rem = (rem << 1) + (a & 1);
+ }
+ return ((rem << 32) | (quot & 0xFFFFFFFFL));
+ }
+
+ /** @see BigInteger#doubleValue() */
+ internal static double bigInteger2Double(BigInteger val) {
+ // val.bitLength() < 64
+ if ((val.numberLength < 2)
+ || ((val.numberLength == 2) && (val.digits[1] > 0))) {
+ return val.longValue();
+ }
+ // val.bitLength() >= 33 * 32 > 1024
+ if (val.numberLength > 32) {
+ return ((val.sign > 0) ? double.PositiveInfinity
+ : double.NegativeInfinity);
+ }
+
+ return Convert.ToDouble(val.ToString());
+ }
+ }
+}
View
783 src/ikc/main/Ioke.Math/Division.cs
@@ -0,0 +1,783 @@
+namespace Ioke.Math {
+ using System;
+
+ class Division {
+ internal static int[] divide(int[] quot, int quotLength, int[] a, int aLength, int[] b, int bLength) {
+ int[] normA = new int[aLength + 1]; // the normalized dividend
+ // an extra byte is needed for correct shift
+ int[] normB = new int[bLength + 1]; // the normalized divisor;
+ int normBLength = bLength;
+ /*
+ * Step D1: normalize a and b and put the results to a1 and b1 the
+ * normalized divisor's first digit must be >= 2^31
+ */
+ int divisorShift = BigDecimal.numberOfLeadingZeros(b[bLength - 1]);
+ if (divisorShift != 0) {
+ BitLevel.shiftLeft(normB, b, 0, divisorShift);
+ BitLevel.shiftLeft(normA, a, 0, divisorShift);
+ } else {
+ Array.Copy(a, normA, aLength);
+ Array.Copy(b, normB, bLength);
+ }
+ int firstDivisorDigit = normB[normBLength - 1];
+ // Step D2: set the quotient index
+ int i = quotLength - 1;
+ int j = aLength;
+
+ while (i >= 0) {
+ // Step D3: calculate a guess digit guessDigit
+ int guessDigit = 0;
+ if (normA[j] == firstDivisorDigit) {
+ // set guessDigit to the largest unsigned int value
+ guessDigit = -1;
+ } else {
+ long product = (((normA[j] & 0xffffffffL) << 32) + (normA[j - 1] & 0xffffffffL));
+ long res = Division.divideLongByInt(product, firstDivisorDigit);
+ guessDigit = (int) res; // the quotient of divideLongByInt
+ int rem = (int) (res >> 32); // the remainder of
+ // divideLongByInt
+ // decrease guessDigit by 1 while leftHand > rightHand
+ if (guessDigit != 0) {
+ long leftHand = 0;
+ long rightHand = 0;
+ bool rOverflowed = false;
+ guessDigit++; // to have the proper value in the loop
+ // below
+ do {
+ guessDigit--;
+ if (rOverflowed) {
+ break;
+ }
+ // leftHand always fits in an unsigned long
+ leftHand = (guessDigit & 0xffffffffL)
+ * (normB[normBLength - 2] & 0xffffffffL);
+ /*
+ * rightHand can overflow; in this case the loop
+ * condition will be true in the next step of the loop
+ */
+ rightHand = ((long) rem << 32)
+ + (normA[j - 2] & 0xffffffffL);
+ long longR = (rem & 0xffffffffL)
+ + (firstDivisorDigit & 0xffffffffL);
+ /*
+ * checks that longR does not fit in an unsigned int;
+ * this ensures that rightHand will overflow unsigned
+ * long in the next step
+ */
+ if (BigDecimal.numberOfLeadingZeros((int) ((long)(((ulong)longR) >> 32))) < 32) {
+ rOverflowed = true;
+ } else {
+ rem = (int) longR;
+ }
+ } while (((long)((ulong)leftHand ^ 0x8000000000000000L) > (long)((ulong)rightHand ^ 0x8000000000000000L)));
+ }
+ }
+ // Step D4: multiply normB by guessDigit and subtract the production
+ // from normA.
+ if (guessDigit != 0) {
+ int borrow = Division.multiplyAndSubtract(normA, j
+ - normBLength, normB, normBLength,
+ guessDigit);
+ // Step D5: check the borrow
+ if (borrow != 0) {
+ // Step D6: compensating addition
+ guessDigit--;
+ long carry = 0;
+ for (int k = 0; k < normBLength; k++) {
+ carry += (normA[j - normBLength + k] & 0xffffffffL)
+ + (normB[k] & 0xffffffffL);
+ normA[j - normBLength + k] = (int) carry;
+ carry = (long)(((ulong)carry) >> 32);
+ }
+ }
+ }
+ if (quot != null) {
+ quot[i] = guessDigit;
+ }
+ // Step D7
+ j--;
+ i--;
+ }
+ /*
+ * Step D8: we got the remainder in normA. Denormalize it id needed
+ */
+ if (divisorShift != 0) {
+ // reuse normB
+ BitLevel.shiftRight(normB, normBLength, normA, 0, divisorShift);
+ return normB;
+ }
+ Array.Copy(normA, normB, bLength);
+ return normA;
+ }
+
+ internal static int divideArrayByInt(int[] dest, int[] src, int srcLength, int divisor) {
+ long rem = 0;
+ long bLong = divisor & 0xffffffffL;
+
+ for (int i = srcLength - 1; i >= 0; i--) {
+ long temp = (rem << 32) | (src[i] & 0xffffffffL);
+ long quot;
+ if (temp >= 0) {
+ quot = (temp / bLong);
+ rem = (temp % bLong);
+ } else {
+ /*
+ * make the dividend positive shifting it right by 1 bit then
+ * get the quotient an remainder and correct them properly
+ */
+ long aPos = (long)(((ulong)temp) >> 1);
+ long bPos = (int)(((uint)divisor) >> 1);
+ quot = aPos / bPos;
+ rem = aPos % bPos;
+ // double the remainder and add 1 if a is odd
+ rem = (rem << 1) + (temp & 1);
+ if ((divisor & 1) != 0) {
+ // the divisor is odd
+ if (quot <= rem) {
+ rem -= quot;
+ } else {
+ if (quot - rem <= bLong) {
+ rem += bLong - quot;
+ quot -= 1;
+ } else {
+ rem += (bLong << 1) - quot;
+ quot -= 2;
+ }
+ }
+ }
+ }
+ dest[i] = (int) (quot & 0xffffffffL);
+ }
+ return (int) rem;
+ }
+
+ internal static int remainderArrayByInt(int[] src, int srcLength,
+ int divisor) {
+
+ long result = 0;
+
+ for (int i = srcLength - 1; i >= 0; i--) {
+ long temp = (result << 32) + (src[i] & 0xffffffffL);
+ long res = divideLongByInt(temp, divisor);
+ result = (int) (res >> 32);
+ }
+ return (int) result;
+ }
+
+ internal static int remainder(BigInteger dividend, int divisor) {
+ return remainderArrayByInt(dividend.digits, dividend.numberLength,
+ divisor);
+ }
+
+ internal static long divideLongByInt(long a, int b) {
+ long quot;
+ long rem;
+ long bLong = b & 0xffffffffL;
+
+ if (a >= 0) {
+ quot = (a / bLong);
+ rem = (a % bLong);
+ } else {
+ /*
+ * Make the dividend positive shifting it right by 1 bit then get
+ * the quotient an remainder and correct them properly
+ */
+ long aPos = (long)(((ulong)a) >> 1);
+ long bPos = (int)(((uint)b) >> 1);
+ quot = aPos / bPos;
+ rem = aPos % bPos;
+ // double the remainder and add 1 if a is odd
+ rem = (rem << 1) + (a & 1);
+ if ((b & 1) != 0) { // the divisor is odd
+ if (quot <= rem) {
+ rem -= quot;
+ } else {
+ if (quot - rem <= bLong) {
+ rem += bLong - quot;
+ quot -= 1;
+ } else {
+ rem += (bLong << 1) - quot;
+ quot -= 2;
+ }
+ }
+ }
+ }
+ return (rem << 32) | (quot & 0xffffffffL);
+ }
+
+ internal static BigInteger[] divideAndRemainderByInteger(BigInteger val,
+ int divisor, int divisorSign) {
+ // res[0] is a quotient and res[1] is a remainder:
+ int[] valDigits = val.digits;
+ int valLen = val.numberLength;
+ int valSign = val.sign;
+ if (valLen == 1) {
+ long a = (valDigits[0] & 0xffffffffL);
+ long b = (divisor & 0xffffffffL);
+ long quo = a / b;
+ long rem = a % b;
+ if (valSign != divisorSign) {
+ quo = -quo;
+ }
+ if (valSign < 0) {
+ rem = -rem;
+ }
+ return new BigInteger[] { BigInteger.valueOf(quo),
+ BigInteger.valueOf(rem) };
+ }
+ int quotientLength = valLen;
+ int quotientSign = ((valSign == divisorSign) ? 1 : -1);
+ int[] quotientDigits = new int[quotientLength];
+ int[] remainderDigits;
+ remainderDigits = new int[] { Division.divideArrayByInt(
+ quotientDigits, valDigits, valLen, divisor) };
+ BigInteger result0 = new BigInteger(quotientSign, quotientLength,
+ quotientDigits);
+ BigInteger result1 = new BigInteger(valSign, 1, remainderDigits);
+ result0.cutOffLeadingZeroes();
+ result1.cutOffLeadingZeroes();
+ return new BigInteger[] { result0, result1 };
+ }
+
+ internal static int multiplyAndSubtract(int[] a, int start, int[] b, int bLen, int c) {
+ long carry0 = 0;
+ long carry1 = 0;
+
+ for (int i = 0; i < bLen; i++) {
+ carry0 = Multiplication.unsignedMultAddAdd(b[i], c, (int)carry0, 0);
+ carry1 = (a[start+i] & 0xffffffffL) - (carry0 & 0xffffffffL) + carry1;
+ a[start+i] = (int)carry1;
+ carry1 >>= 32; // -1 or 0
+ carry0 = (long)(((ulong)carry0) >> 32);
+ }
+
+ carry1 = (a[start + bLen] & 0xffffffffL) - carry0 + carry1;
+ a[start + bLen] = (int)carry1;
+ return (int)(carry1 >> 32); // -1 or 0
+ }
+
+ internal static BigInteger gcdBinary(BigInteger op1, BigInteger op2) {
+ // PRE: (op1 > 0) and (op2 > 0)
+
+ /*
+ * Divide both number the maximal possible times by 2 without rounding
+ * gcd(2*a, 2*b) = 2 * gcd(a,b)
+ */
+ int lsb1 = op1.getLowestSetBit();
+ int lsb2 = op2.getLowestSetBit();
+ int pow2Count = Math.Min(lsb1, lsb2);
+
+ BitLevel.inplaceShiftRight(op1, lsb1);
+ BitLevel.inplaceShiftRight(op2, lsb2);
+
+ BigInteger swap;
+ // I want op2 > op1
+ if (op1.compareTo(op2) == BigInteger.GREATER) {
+ swap = op1;
+ op1 = op2;
+ op2 = swap;
+ }
+
+ do { // INV: op2 >= op1 && both are odd unless op1 = 0
+
+ // Optimization for small operands
+ // (op2.bitLength() < 64) implies by INV (op1.bitLength() < 64)
+ if (( op2.numberLength == 1 )
+ || ( ( op2.numberLength == 2 ) && ( op2.digits[1] > 0 ) )) {
+ op2 = BigInteger.valueOf(Division.gcdBinary(op1.longValue(),
+ op2.longValue()));
+ break;
+ }
+
+ // Implements one step of the Euclidean algorithm
+ // To reduce one operand if it's much smaller than the other one
+ if (op2.numberLength > op1.numberLength * 1.2) {
+ op2 = op2.remainder(op1);
+ if (op2.signum() != 0) {
+ BitLevel.inplaceShiftRight(op2, op2.getLowestSetBit());
+ }
+ } else {
+
+ // Use Knuth's algorithm of successive subtract and shifting
+ do {
+ Elementary.inplaceSubtract(op2, op1); // both are odd
+ BitLevel.inplaceShiftRight(op2, op2.getLowestSetBit()); // op2 is even
+ } while (op2.compareTo(op1) >= BigInteger.EQUALS);
+ }
+ // now op1 >= op2
+ swap = op2;
+ op2 = op1;
+ op1 = swap;
+ } while (op1.sign != 0);
+ return op2.shiftLeft(pow2Count);
+ }
+
+ internal static long gcdBinary(long op1, long op2) {
+ // PRE: (op1 > 0) and (op2 > 0)
+ int lsb1 = BigDecimal.numberOfTrailingZeros(op1);
+ int lsb2 = BigDecimal.numberOfTrailingZeros(op2);
+ int pow2Count = Math.Min(lsb1, lsb2);
+
+ if (lsb1 != 0) {
+ op1 = (long)(((ulong)op1) >> lsb1);
+ }
+ if (lsb2 != 0) {
+ op2 = (long)(((ulong)op2) >> lsb2);
+ }
+ do {
+ if (op1 >= op2) {
+ op1 -= op2;
+ op1 = (long)(((ulong)op1) >> BigDecimal.numberOfTrailingZeros(op1));
+ } else {
+ op2 -= op1;
+ op2 = (long)(((ulong)op2) >> BigDecimal.numberOfTrailingZeros(op2));
+ }
+ } while (op1 != 0);
+ return ( op2 << pow2Count );
+ }
+
+ internal static BigInteger modInverseMontgomery(BigInteger a, BigInteger p) {
+
+ if (a.sign == 0){
+ // ZERO hasn't inverse
+ throw new ArithmeticException("BigInteger not invertible");
+ }
+
+
+ if (!p.testBit(0)){
+ // montgomery inverse require even modulo
+ return modInverseLorencz(a, p);
+ }
+
+ int m = p.numberLength * 32;
+ // PRE: a \in [1, p - 1]
+ BigInteger u, v, r, s;
+ u = p.copy(); // make copy to use inplace method
+ v = a.copy();
+ int max = Math.Max(v.numberLength, u.numberLength);
+ r = new BigInteger(1, 1, new int[max + 1]);
+ s = new BigInteger(1, 1, new int[max + 1]);
+ s.digits[0] = 1;
+ // s == 1 && v == 0
+
+ int k = 0;
+
+ int lsbu = u.getLowestSetBit();
+ int lsbv = v.getLowestSetBit();
+ int toShift;
+
+ if (lsbu > lsbv) {
+ BitLevel.inplaceShiftRight(u, lsbu);
+ BitLevel.inplaceShiftRight(v, lsbv);
+ BitLevel.inplaceShiftLeft(r, lsbv);
+ k += lsbu - lsbv;
+ } else {
+ BitLevel.inplaceShiftRight(u, lsbu);
+ BitLevel.inplaceShiftRight(v, lsbv);
+ BitLevel.inplaceShiftLeft(s, lsbu);
+ k += lsbv - lsbu;
+ }
+
+ r.sign = 1;
+ while (v.signum() > 0) {
+ // INV v >= 0, u >= 0, v odd, u odd (except last iteration when v is even (0))
+
+ while (u.compareTo(v) > BigInteger.EQUALS) {
+ Elementary.inplaceSubtract(u, v);
+ toShift = u.getLowestSetBit();
+ BitLevel.inplaceShiftRight(u, toShift);
+ Elementary.inplaceAdd(r, s);
+ BitLevel.inplaceShiftLeft(s, toShift);
+ k += toShift;
+ }
+
+ while (u.compareTo(v) <= BigInteger.EQUALS) {
+ Elementary.inplaceSubtract(v, u);
+ if (v.signum() == 0)
+ break;
+ toShift = v.getLowestSetBit();
+ BitLevel.inplaceShiftRight(v, toShift);
+ Elementary.inplaceAdd(s, r);
+ BitLevel.inplaceShiftLeft(r, toShift);
+ k += toShift;
+ }
+ }
+ if (!u.isOne()){
+ // in u is stored the gcd
+ throw new ArithmeticException("BigInteger not invertible.");
+ }
+ if (r.compareTo(p) >= BigInteger.EQUALS) {
+ Elementary.inplaceSubtract(r, p);
+ }
+
+ r = p.subtract(r);
+
+ // Have pair: ((BigInteger)r, (Integer)k) where r == a^(-1) * 2^k mod (module)
+ int n1 = calcN(p);
+ if (k > m) {
+ r = monPro(r, BigInteger.ONE, p, n1);
+ k = k - m;
+ }
+
+ r = monPro(r, BigInteger.getPowerOfTwo(m - k), p, n1);
+ return r;
+ }
+
+ private static int calcN(BigInteger a) {
+ long m0 = a.digits[0] & 0xFFFFFFFFL;
+ long n2 = 1L; // this is a'[0]
+ long powerOfTwo = 2L;
+ do {
+ if (((m0 * n2) & powerOfTwo) != 0) {
+ n2 |= powerOfTwo;
+ }
+ powerOfTwo <<= 1;
+ } while (powerOfTwo < 0x100000000L);
+ n2 = -n2;
+ return (int)(n2 & 0xFFFFFFFFL);
+ }
+
+ private static bool isPowerOfTwo(BigInteger bi, int exp) {
+ bool result = false;
+ result = ( exp >> 5 == bi.numberLength - 1 )
+ && ( bi.digits[bi.numberLength - 1] == 1 << ( exp & 31 ) );
+ if (result) {
+ for (int i = 0; result && i < bi.numberLength - 1; i++) {
+ result = bi.digits[i] == 0;
+ }
+ }
+ return result;
+ }
+
+ private static int howManyIterations(BigInteger bi, int n) {
+ int i = n - 1;
+ if (bi.sign > 0) {
+ while (!bi.testBit(i))
+ i--;
+ return n - 1 - i;
+ } else {
+ while (bi.testBit(i))
+ i--;
+ return n - 1 - Math.Max(i, bi.getLowestSetBit());
+ }
+
+ }
+
+ internal static BigInteger modInverseLorencz(BigInteger a, BigInteger modulo) {
+ int max = Math.Max(a.numberLength, modulo.numberLength);
+ int[] uDigits = new int[max + 1]; // enough place to make all the inplace operation
+ int[] vDigits = new int[max + 1];
+ Array.Copy(modulo.digits, uDigits, modulo.numberLength);
+ Array.Copy(a.digits, vDigits, a.numberLength);
+ BigInteger u = new BigInteger(modulo.sign, modulo.numberLength,
+ uDigits);
+ BigInteger v = new BigInteger(a.sign, a.numberLength, vDigits);
+
+ BigInteger r = new BigInteger(0, 1, new int[max + 1]); // BigInteger.ZERO;
+ BigInteger s = new BigInteger(1, 1, new int[max + 1]);
+ s.digits[0] = 1;
+ // r == 0 && s == 1, but with enough place
+
+ int coefU = 0, coefV = 0;
+ int n = modulo.bitLength();
+ int k;
+ while (!isPowerOfTwo(u, coefU) && !isPowerOfTwo(v, coefV)) {
+
+ // modification of original algorithm: I calculate how many times the algorithm will enter in the same branch of if
+ k = howManyIterations(u, n);
+
+ if (k != 0) {
+ BitLevel.inplaceShiftLeft(u, k);
+ if (coefU >= coefV) {
+ BitLevel.inplaceShiftLeft(r, k);
+ } else {
+ BitLevel.inplaceShiftRight(s, Math.Min(coefV - coefU, k));
+ if (k - ( coefV - coefU ) > 0) {
+ BitLevel.inplaceShiftLeft(r, k - coefV + coefU);
+ }
+ }
+ coefU += k;
+ }
+
+ k = howManyIterations(v, n);
+ if (k != 0) {
+ BitLevel.inplaceShiftLeft(v, k);
+ if (coefV >= coefU) {
+ BitLevel.inplaceShiftLeft(s, k);
+ } else {
+ BitLevel.inplaceShiftRight(r, Math.Min(coefU - coefV, k));
+ if (k - ( coefU - coefV ) > 0) {
+ BitLevel.inplaceShiftLeft(s, k - coefU + coefV);
+ }
+ }
+ coefV += k;
+
+ }
+
+ if (u.signum() == v.signum()) {
+ if (coefU <= coefV) {
+ Elementary.completeInPlaceSubtract(u, v);
+ Elementary.completeInPlaceSubtract(r, s);
+ } else {
+ Elementary.completeInPlaceSubtract(v, u);
+ Elementary.completeInPlaceSubtract(s, r);
+ }
+ } else {
+ if (coefU <= coefV) {
+ Elementary.completeInPlaceAdd(u, v);
+ Elementary.completeInPlaceAdd(r, s);
+ } else {
+ Elementary.completeInPlaceAdd(v, u);
+ Elementary.completeInPlaceAdd(s, r);
+ }
+ }
+ if (v.signum() == 0 || u.signum() == 0){
+ throw new ArithmeticException("BigInteger not invertible");
+ }
+ }
+
+ if (isPowerOfTwo(v, coefV)) {
+ r = s;
+ if (v.signum() != u.signum())
+ u = u.negate();
+ }
+ if (u.testBit(n)) {
+ if (r.signum() < 0) {
+ r = r.negate();
+ } else {
+ r = modulo.subtract(r);
+ }
+ }
+ if (r.signum() < 0) {
+ r = r.add(modulo);
+ }
+
+ return r;
+ }
+
+ internal static BigInteger squareAndMultiply(BigInteger x2, BigInteger a2, BigInteger exponent,BigInteger modulus, int n2 ){
+ BigInteger res = x2;
+ for (int i = exponent.bitLength() - 1; i >= 0; i--) {
+ res = monPro(res,res,modulus, n2);
+ if (BitLevel.testBit(exponent, i)) {
+ res = monPro(res, a2, modulus, n2);
+ }
+ }
+ return res;
+ }
+
+ internal static BigInteger slidingWindow(BigInteger x2, BigInteger a2, BigInteger exponent,BigInteger modulus, int n2){
+ // fill odd low pows of a2
+ BigInteger[] pows = new BigInteger[8];
+ BigInteger res = x2;
+ int lowexp;
+ BigInteger x3;
+ int acc3;
+ pows[0] = a2;
+
+ x3 = monPro(a2,a2,modulus,n2);
+ for (int i = 1; i <= 7; i++){
+ pows[i] = monPro(pows[i-1],x3,modulus,n2) ;
+ }
+
+ for (int i = exponent.bitLength()-1; i>=0;i--){
+ if( BitLevel.testBit(exponent,i) ) {
+ lowexp = 1;
+ acc3 = i;
+
+ for(int j = Math.Max(i-3,0);j <= i-1 ;j++) {
+ if (BitLevel.testBit(exponent,j)) {
+ if (j<acc3) {
+ acc3 = j;
+ lowexp = (lowexp << (i-j))^1;
+ } else {
+ lowexp = lowexp^(1<<(j-acc3));
+ }
+ }
+ }
+
+ for(int j = acc3; j <= i; j++) {
+ res = monPro(res,res,modulus,n2);
+ }
+ res = monPro(pows[(lowexp-1)>>1], res, modulus,n2);
+ i = acc3 ;
+ }else{
+ res = monPro(res, res, modulus, n2) ;
+ }
+ }
+ return res;
+ }
+
+ internal static BigInteger oddModPow(BigInteger _base, BigInteger exponent,
+ BigInteger modulus) {
+ // PRE: (base > 0), (exponent > 0), (modulus > 0) and (odd modulus)
+ int k = (modulus.numberLength << 5); // r = 2^k
+ // n-residue of base [base * r (mod modulus)]
+ BigInteger a2 = _base.shiftLeft(k).mod(modulus);
+ // n-residue of base [1 * r (mod modulus)]
+ BigInteger x2 = BigInteger.getPowerOfTwo(k).mod(modulus);
+ BigInteger res;
+ // Compute (modulus[0]^(-1)) (mod 2^32) for odd modulus
+
+ int n2 = calcN(modulus);
+ if( modulus.numberLength == 1 ){
+ res = squareAndMultiply(x2,a2, exponent, modulus,n2);
+ } else {
+ res = slidingWindow(x2, a2, exponent, modulus, n2);
+ }
+
+ return monPro(res, BigInteger.ONE, modulus, n2);
+ }
+
+ internal static BigInteger evenModPow(BigInteger _base, BigInteger exponent,
+ BigInteger modulus) {
+ // PRE: (base > 0), (exponent > 0), (modulus > 0) and (modulus even)
+ // STEP 1: Obtain the factorization 'modulus'= q * 2^j.
+ int j = modulus.getLowestSetBit();
+ BigInteger q = modulus.shiftRight(j);
+
+ // STEP 2: Compute x1 := base^exponent (mod q).
+ BigInteger x1 = oddModPow(_base, exponent, q);
+
+ // STEP 3: Compute x2 := base^exponent (mod 2^j).
+ BigInteger x2 = pow2ModPow(_base, exponent, j);
+
+ // STEP 4: Compute q^(-1) (mod 2^j) and y := (x2-x1) * q^(-1) (mod 2^j)
+ BigInteger qInv = modPow2Inverse(q, j);
+ BigInteger y = (x2.subtract(x1)).multiply(qInv);
+ inplaceModPow2(y, j);
+ if (y.sign < 0) {
+ y = y.add(BigInteger.getPowerOfTwo(j));
+ }
+ // STEP 5: Compute and return: x1 + q * y
+ return x1.add(q.multiply(y));
+ }
+
+ internal static BigInteger pow2ModPow(BigInteger _base, BigInteger exponent, int j) {
+ // PRE: (base > 0), (exponent > 0) and (j > 0)
+ BigInteger res = BigInteger.ONE;
+ BigInteger e = exponent.copy();
+ BigInteger baseMod2toN = _base.copy();
+ BigInteger res2;
+ /*
+ * If 'base' is odd then it's coprime with 2^j and phi(2^j) = 2^(j-1);
+ * so we can reduce reduce the exponent (mod 2^(j-1)).
+ */
+ if (_base.testBit(0)) {
+ inplaceModPow2(e, j - 1);
+ }
+ inplaceModPow2(baseMod2toN, j);
+
+ for (int i = e.bitLength() - 1; i >= 0; i--) {
+ res2 = res.copy();
+ inplaceModPow2(res2, j);
+ res = res.multiply(res2);
+ if (BitLevel.testBit(e, i)) {
+ res = res.multiply(baseMod2toN);
+ inplaceModPow2(res, j);
+ }
+ }
+ inplaceModPow2(res, j);
+ return res;
+ }
+
+ private static void monReduction(int[] res, BigInteger modulus, int n2) {
+
+ /* res + m*modulus_digits */
+ int[] modulus_digits = modulus.digits;
+ int modulusLen = modulus.numberLength;
+ long outerCarry = 0;
+
+ for (int i = 0; i < modulusLen; i++){
+ long innerCarry = 0;
+ int m = (int) Multiplication.unsignedMultAddAdd(res[i],n2,0,0);
+ for(int j = 0; j < modulusLen; j++){
+ innerCarry = Multiplication.unsignedMultAddAdd(m, modulus_digits[j], res[i+j], (int)innerCarry);
+ res[i+j] = (int) innerCarry;
+ innerCarry = (long)(((ulong)innerCarry) >> 32);
+ }
+
+ outerCarry += (res[i+modulusLen] & 0xFFFFFFFFL) + innerCarry;
+ res[i+modulusLen] = (int) outerCarry;
+ outerCarry = (long)(((ulong)outerCarry) >> 32);
+ }
+
+ res[modulusLen << 1] = (int) outerCarry;
+
+ /* res / r */
+ for(int j = 0; j < modulusLen+1; j++){
+ res[j] = res[j+modulusLen];
+ }
+ }
+
+ internal static BigInteger monPro(BigInteger a, BigInteger b, BigInteger modulus, int n2) {
+ int modulusLen = modulus.numberLength;
+ int[] res = new int[(modulusLen << 1) + 1];
+ Multiplication.multArraysPAP(a.digits, Math.Min(modulusLen, a.numberLength),
+ b.digits, Math.Min(modulusLen, b.numberLength), res);
+ monReduction(res,modulus,n2);
+ return finalSubtraction(res, modulus);
+
+ }
+
+ internal static BigInteger finalSubtraction(int[] res, BigInteger modulus){
+
+ // skipping leading zeros
+ int modulusLen = modulus.numberLength;
+ bool doSub = res[modulusLen]!=0;
+ if(!doSub) {
+ int[] modulusDigits = modulus.digits;
+ doSub = true;
+ for(int i = modulusLen - 1; i >= 0; i--) {
+ if(res[i] != modulusDigits[i]) {
+ doSub = (res[i] != 0) && ((res[i] & 0xFFFFFFFFL) > (modulusDigits[i] & 0xFFFFFFFFL));
+ break;
+ }
+ }
+ }
+
+ BigInteger result = new BigInteger(1, modulusLen+1, res);
+
+ // if (res >= modulusDigits) compute (res - modulusDigits)
+ if (doSub) {
+ Elementary.inplaceSubtract(result, modulus);
+ }
+
+ result.cutOffLeadingZeroes();
+ return result;
+ }
+
+ internal static BigInteger modPow2Inverse(BigInteger x, int n) {
+ // PRE: (x > 0), (x is odd), and (n > 0)
+ BigInteger y = new BigInteger(1, new int[1 << n]);
+ y.numberLength = 1;
+ y.digits[0] = 1;
+ y.sign = 1;
+
+ for (int i = 1; i < n; i++) {
+ if (BitLevel.testBit(x.multiply(y), i)) {
+ // Adding 2^i to y (setting the i-th bit)
+ y.digits[i >> 5] |= (1 << (i & 31));
+ }
+ }
+ return y;
+ }
+
+ internal static void inplaceModPow2(BigInteger x, int n) {
+ // PRE: (x > 0) and (n >= 0)
+ int fd = n >> 5;
+ int leadingZeros;
+
+ if ((x.numberLength < fd) || (x.bitLength() <= n)) {
+ return;
+ }
+ leadingZeros = 32 - (n & 31);
+ x.numberLength = fd + 1;
+ unchecked {
+ x.digits[fd] &= (leadingZeros < 32) ? ((int)(((uint)-1) >> leadingZeros)) : 0;
+ }
+ x.cutOffLeadingZeroes();
+ }
+
+ }
+}
View
325 src/ikc/main/Ioke.Math/Elementary.cs
@@ -0,0 +1,325 @@
+namespace Ioke.Math {
+ using System;
+ using System.Text;
+
+ class Elementary {
+ private Elementary() {
+ }
+
+ internal static int compareArrays(int[] a, int[] b, int size) {
+ int i;
+ for (i = size - 1; (i >= 0) && (a[i] == b[i]); i--) {
+ ;
+ }
+ return ((i < 0) ? BigInteger.EQUALS
+ : (a[i] & 0xFFFFFFFFL) < (b[i] & 0xFFFFFFFFL) ? BigInteger.LESS
+ : BigInteger.GREATER);
+ }
+
+ internal static BigInteger add(BigInteger op1, BigInteger op2) {
+ int[] resDigits;
+ int resSign;
+ int op1Sign = op1.sign;
+ int op2Sign = op2.sign;
+
+ if (op1Sign == 0) {
+ return op2;
+ }
+ if (op2Sign == 0) {
+ return op1;
+ }