Java, custom big integer class: 32.9 (120000000 / 365000)
The main class is pretty straightforward:
import java.util.*;
public class PPCG37270 {
public static void main(String[] args) {
long start = System.nanoTime();
int n = 12000000;
if (args.length == 1) n = Integer.parseInt(args[0]);
boolean[] sieve = new boolean[n + 1];
int[] remaining = new int[n + 1];
int[] count = new int[n + 1];
for (int p = 2; p <= n; p++) {
if (sieve[p]) continue;
long p2 = p * (long)p;
if (p2 > n) continue;
for (int i = (int)p2; i <= n; i += p) sieve[i] = true;
}
for (int i = 2; i <= n; i++) remaining[i] = i;
for (int p = 2; p <= n; p++) {
if (sieve[p]) continue;
for (int i = p; i <= n; i += p) {
while (remaining[i] % p == 0) {
remaining[i] /= p;
count[p]++;
if (i <= n/2) count[p] -= 2;
}
}
}
count[2] -= count[5];
count[5] = 0;
List<BigInt> partialProd = new ArrayList<BigInt>();
long accum = 1;
for (int i = 2; i <= n; i++) {
for (int j = count[i]; j > 0; j--) {
long tmp = accum * i;
if (tmp < 1000000000L) accum = tmp;
else {
partialProd.add(new BigInt((int)accum));
accum = i;
}
}
}
partialProd.add(new BigInt((int)accum));
System.out.println(prod(partialProd).digsum());
System.out.println((System.nanoTime() - start) / 1000000 + "ms");
}
private static BigInt prod(List<BigInt> vals) {
while (vals.size() > 1) {
int n = vals.size();
List<BigInt> next = new ArrayList<BigInt>();
for (int i = 0; i < n; i += 2) {
if (i == n - 1) next.add(vals.get(i));
else next.add(vals.get(i).mul(vals.get(i+1)));
}
vals = next;
}
return vals.get(0);
}
}
It relies on a big integer class which is optimised for multiplication and toString()
, both of which are significant bottlenecks in an implementation with java.math.BigInteger
.
/**
* A big integer class which is optimised for conversion to decimal.
* For use in simple applications where BigInteger.toString() is a bottleneck.
*/
public class BigInt {
// The base of the representation.
private static final int B = 1000000000;
// The number of decimal digits per digit of the representation.
private static final int LOG10_B = 9;
public static final BigInt ZERO = new BigInt(0);
public static final BigInt ONE = new BigInt(1);
// We use sign-magnitude representation.
private final boolean negative;
// Least significant digit is at val[off]; most significant is at val[off + len - 1]
// Unless len == 1 we guarantee that val[off + len - 1] is non-zero.
private final int[] val;
private final int off;
private final int len;
// Toom-style multiplication parameters from
// Zuras, D. (1994). More on squaring and multiplying large integers. IEEE Transactions on Computers, 43(8), 899-908.
private static final int[][][] Q = new int[][][]{
{},
{},
{{1, -1}},
{{4, 2, 1}, {1, 1, 1}, {1, 2, 4}},
{{8, 4, 2, 1}, {-8, 4, -2, 1}, {1, 1, 1, 1}, {1, -2, 4, -8}, {1, 2, 4, 8}}
};
private static final int[][][] R = new int[][][]{
{},
{},
{{1, -1, 1}},
{{-21, 2, -12, 1, -6}, {7, -1, 10, -1, 7}, {-6, 1, -12, 2, -21}},
{{-180, 6, 2, -80, 1, 3, -180}, {-510, 4, 4, 0, -1, -1, 120}, {1530, -27, -7, 680, -7, -27, 1530}, {120, -1, -1, 0, 4, 4, -510}, {-180, 3, 1, -80, 2, 6, -180}}
};
private static final int[][] S = new int[][]{
{},
{},
{1, 1, 1},
{1, 6, 2, 6, 1},
{1, 180, 120, 360, 120, 180, 1}
};
/**
* Constructs a big version of an integer value.
* @param x The value to represent.
*/
public BigInt(int x) {
this(Integer.toString(x));
}
/**
* Constructs a big version of a long value.
* @param x The value to represent.
*/
public BigInt(long x) {
this(Long.toString(x));
}
/**
* Parses a decimal representation of an integer.
* @param str The value to represent.
*/
public BigInt(String str) {
this(str.charAt(0) == '-', split(str));
}
/**
* Constructs a sign-magnitude representation taking the entire span of the array as the range of interest.
* @param neg Is the value negative?
* @param val The base-B digits, least significant first.
*/
private BigInt(boolean neg, int[] val) {
this(neg, val, 0, val.length);
}
/**
* Constructs a sign-magnitude representation taking a range of an array as the magnitude.
* @param neg Is the value negative?
* @param val The base-B digits, least significant at offset off, most significant at off + val - 1.
* @param off The offset within the array.
* @param len The number of base-B digits.
*/
private BigInt(boolean neg, int[] val, int off, int len) {
// Bounds checks
if (val == null) throw new IllegalArgumentException("val");
if (off < 0 || off >= val.length) throw new IllegalArgumentException("off");
if (len < 1 || off + len > val.length) throw new IllegalArgumentException("len");
this.negative = neg;
this.val = val;
this.off = off;
// Enforce the invariant that this.len is 1 or val[off + len - 1] is non-zero.
while (len > 1 && val[off + len - 1] == 0) len--;
this.len = len;
// Sanity check
for (int i = 0; i < len; i++) {
if (val[off + i] < 0) throw new IllegalArgumentException("val contains negative digits");
}
}
/**
* Splits a string into base-B digits.
* @param str The string to parse.
* @return An array which can be passed to the (boolean, int[]) constructor.
*/
private static int[] split(String str) {
if (str.charAt(0) == '-') str = str.substring(1);
int[] arr = new int[(str.length() + LOG10_B - 1) / LOG10_B];
int i, off;
// Each element of arr represents LOG10_B characters except (probably) the last one.
for (i = 0, off = str.length() - LOG10_B; off > 0; off -= LOG10_B) {
arr[i++] = Integer.parseInt(str.substring(off, off + LOG10_B));
}
arr[i] = Integer.parseInt(str.substring(0, off + LOG10_B));
return arr;
}
public boolean isZero() {
return len == 1 && val[off] == 0;
}
public BigInt negate() {
return new BigInt(!negative, val, off, len);
}
public BigInt add(BigInt that) {
// If the signs differ, then since we use sign-magnitude representation we want to do a subtraction.
boolean isSubtraction = negative ^ that.negative;
BigInt left, right;
if (len < that.len) {
left = that;
right = this;
}
else {
left = this;
right = that;
// For addition I just care about the lengths of the arrays.
// For subtraction I want the largest absolute value on the left.
if (isSubtraction && len == that.len) {
int cmp = compareAbsolute(that);
if (cmp == 0) return ZERO; // Cheap special case
if (cmp < 0) {
left = that;
right = this;
}
}
}
if (right.isZero()) return left;
BigInt result;
if (!isSubtraction) {
int[] sum = new int[left.len + 1];
// A copy here rather than using left.val in the main loops and copying remaining values
// at the end gives a small performance boost, probably due to cache locality.
System.arraycopy(left.val, left.off, sum, 0, left.len);
int carry = 0, k = 0;
for (; k < right.len; k++) {
int a = sum[k] + right.val[right.off + k] + carry;
sum[k] = a % B;
carry = a / B;
}
for (; carry > 0 && k < left.len; k++) {
int a = sum[k] + carry;
sum[k] = a % B;
carry = a / B;
}
sum[left.len] = carry;
result = new BigInt(negative, sum);
}
else {
int[] diff = new int[left.len];
System.arraycopy(left.val, left.off, diff, 0, left.len);
int carry = 0, k = 0;
for (; k < right.len; k++) {
int a = diff[k] - right.val[right.off + k] + carry;
// Why did anyone ever think that rounding positive and negative divisions differently made sense?
if (a < 0) {
diff[k] = a + B;
carry = -1;
}
else {
diff[k] = a % B;
carry = a / B;
}
}
for (; carry != 0 && k < left.len; k++) {
int a = diff[k] + carry;
if (a < 0) {
diff[k] = a + B;
carry = -1;
}
else {
diff[k] = a % B;
carry = a / B;
}
}
result = new BigInt(left.negative, diff, 0, k > left.len ? k : left.len);
}
return result;
}
private int compareAbsolute(BigInt that) {
if (len > that.len) return 1;
if (len < that.len) return -1;
for (int i = len - 1; i >= 0; i--) {
if (val[off + i] > that.val[that.off + i]) return 1;
if (val[off + i] < that.val[that.off + i]) return -1;
}
return 0;
}
public BigInt mul(BigInt that) {
if (isZero() || that.isZero()) return ZERO;
if (len == 1) return that.mulSmall(negative ? -val[off] : val[off]);
if (that.len == 1) return mulSmall(that.negative ? -that.val[that.off] : that.val[that.off]);
int shorter = len < that.len ? len : that.len;
BigInt result;
// Cutoffs have been hand-tuned.
if (shorter > 300) result = mulToom(3, that);
else if (shorter > 28) result = mulToom(2, that);
else result = mulNaive(that);
return result;
}
BigInt mulSmall(int m) {
if (m == 0) return ZERO;
if (m == 1) return this;
if (m == -1) return negate();
// We want to do the magnitude calculation with a positive multiplicand.
boolean neg = negative;
if (m < 0) {
neg = !neg;
m = -m;
}
int[] pr = new int[len + 1];
int carry = 0;
for (int i = 0; i < len; i++) {
long t = val[off + i] * (long)m + carry;
pr[i] = (int)(t % B);
carry = (int)(t / B);
}
pr[len] = carry;
return new BigInt(neg, pr);
}
// NB This truncates.
BigInt divSmall(int d) {
if (d == 0) throw new ArithmeticException();
if (d == 1) return this;
if (d == -1) return negate();
// We want to do the magnitude calculation with a positive divisor.
boolean neg = negative;
if (d < 0) {
neg = !neg;
d = -d;
}
int[] div = new int[len];
int rem = 0;
for (int i = len - 1; i >= 0; i--) {
long t = val[off + i] + rem * (long)B;
div[i] = (int)(t / d);
rem = (int)(t % d);
}
return new BigInt(neg, div);
}
BigInt mulNaive(BigInt that) {
int[] rv = new int[len + that.len];
// Naive multiplication
for (int i = 0; i < len; i++) {
for (int j = 0; j < that.len; j++) {
int k = i + j;
long c = val[off + i] * (long)that.val[that.off + j];
while (c > 0) {
c += rv[k];
rv[k] = (int)(c % B);
c /= B;
k++;
}
}
}
return new BigInt(this.negative ^ that.negative, rv);
}
private BigInt mulToom(int k, BigInt that) {
// We split each number into k parts of m base-B digits each.
// m = ceil(longer / k)
int m = ((len > that.len ? len : that.len) + k - 1) / k;
// Perform the splitting and evaluation steps of Toom-Cook.
BigInt[] f1 = this.toomFwd(k, m);
BigInt[] f2 = that.toomFwd(k, m);
// Pointwise multiplication.
for (int i = 0; i < f1.length; i++) f1[i] = f1[i].mul(f2[i]);
// Inverse (or interpolation) and recomposition.
return toomBk(k, m, f1, negative ^ that.negative, val[off], that.val[that.off]);
}
// Splits a number into k parts of m base-B digits each and does the polynomial evaluation.
private BigInt[] toomFwd(int k, int m) {
// Split.
BigInt[] a = new BigInt[k];
for (int i = 0; i < k; i++) {
int o = i * m;
if (o >= len) a[i] = ZERO;
else {
int l = m;
if (o + l > len) l = len - o;
// Ignore signs for now.
a[i] = new BigInt(false, val, off + o, l);
}
}
// Evaluate
return transform(Q[k], a);
}
private BigInt toomBk(int k, int m, BigInt[] f, boolean neg, int lsd1, int lsd2) {
// Inverse (or interpolation).
BigInt[] b = transform(R[k], f);
// Recomposition: add at suitable offsets, dividing by the normalisation factors
BigInt prod = ZERO;
int[] s = S[k];
for (int i = 0; i < b.length; i++) {
int[] shifted = new int[i * m + b[i].len];
System.arraycopy(b[i].val, b[i].off, shifted, i * m, b[i].len);
prod = prod.add(new BigInt(neg ^ b[i].negative, shifted).divSmall(s[i]));
}
// Handle the remainders.
// In the worst case the absolute value of the sum of the remainders is s.length, so pretty small.
// It should be easy enough to work out whether to go up or down.
int lsd = (int)((lsd1 * (long)lsd2) % B);
int err = lsd - prod.val[prod.off];
if (err > B / 2) err -= B / 2;
if (err < -B / 2) err += B / 2;
return prod.add(new BigInt(err));
}
/**
* Multiplies a matrix of small integers and a vector of big ones.
* The matrix has a implicit leading row [1 0 ... 0] and an implicit trailing row [0 ... 0 1].
* @param m The matrix.
* @param v The vector.
* @return m v
*/
private BigInt[] transform(int[][] m, BigInt[] v) {
BigInt[] b = new BigInt[m.length + 2];
b[0] = v[0];
for (int i = 0; i < m.length; i++) {
BigInt s = ZERO;
for (int j = 0; j < m[i].length; j++) s = s.add(v[j].mulSmall(m[i][j]));
b[i + 1] = s;
}
b[b.length - 1] = v[v.length - 1];
return b;
}
/**
* Sums the digits of this integer.
* @return The sum of the digits of this integer.
*/
public long digsum() {
long rv = 0;
for (int i = 0; i < len; i++) {
int x = val[off + i];
while (x > 0) {
rv += x % 10;
x /= 10;
}
}
return rv;
}
}
The big bottleneck is naïve multiplication (60%), followed by the other multiplication (37%) and the sieving (3%). The digsum()
call is insignificant.
Performance measured with OpenJDK 7 (64 bit).
2Are we supposed to submit solutions in python or in a language of choice? – None – 2014-09-03T19:31:47.480
It's possible to write a routine that does this on a modern computer and takes
n
well into the millions, while I doubt the Python function would handle anything greater thann = 1e5
without choking. – COTO – 2014-09-03T20:05:45.910@Alessandro You can use any language of your choice. The only restriction was that you can't use builtin functions to compute the coefficients. – None – 2014-09-03T20:27:28.877
@COTO That sounds very interesting. I look forward to seeing it! – None – 2014-09-03T20:29:30.627
It was just an observation, @Lembik. I'm unfortunately splitting my time between work and one of the existing challenges at this point. I'll upvote the question to get it a bit of exposure, however. ;) – COTO – 2014-09-03T20:34:43.007
What if my highest n is n=1? (What I want to say is that someone might profit from having a much lower/higher n than someone else in regard of the time complexity of each algorithm) I'd suggest that the algorithm has to work up to a certain (fixed) n and that the timing is calculate according to that. – flawr – 2014-09-03T21:02:36.480
I want to point of Lucas' Theorem and this generalization to higher powers , though unfortunately 10 is not prime.
– xnor – 2014-09-03T21:31:33.467@COTO Well into the millions is an understatement. On my old Core Quad 6600 I can get a result for n = 100,000,000 in 32 seconds in C with mpz_bin_uiui(). I actually takes twice as long to stream the digits to the terminal... – Comintern – 2014-09-04T02:37:13.827
@Comintern Part of the challenge is not being able to use a built-in function. – qwr – 2014-09-04T02:43:00.313
@qwr - That's why I didn't post it as an answer. ;-) Mainly wanted to give a frame of reference for performance. – Comintern – 2014-09-04T02:44:35.433
2Are factorial functions allowed? I assumed not since they could be "quickly transformed into binomial coefficients" (the whole thing is just one factorial divided by another factorial squared), but since an answer is using one now, clarity would be nice. – Geobits – 2014-09-04T03:00:08.387
@Geobits You are right that using a built in factorial function is cheating. – None – 2014-09-04T07:05:59.520
1@Comintern: I have successfully replicated that point of reference with 287mil in 1 min, or 169mil in 35 secs! :) – justhalf – 2014-09-05T04:15:12.873
Just realized that the 287mil in 1 min includes the time to convert the result into string. The calculation itself takes only a third of the time (like Comintern noticed in his/her comment): around 18s. – justhalf – 2014-09-08T06:59:29.127