new.java

  public static class CountedLogSumExp extends CountedCompleter<Double> {
    final int lo;
    final int hi;
    final double[] array;
    final int blockSize;
    double sum;
    double offset;

    CountedLogSumExp[] siblings;

    CountedLogSumExp(final double[] array, final int blockSize, final int maxSiblings) {
      this(null, array, Doubles.max(array),  0, array.length, blockSize, maxSiblings);
    }

    CountedLogSumExp(final CountedCompleter<?> parent,
                     final double[] array, final double offset,
                     final int lo, final int hi,
                     final int blockSize, final int maxSiblings) {
      super(parent);
      this.array = array;
      this.offset = offset;
      this.lo = lo;
      this.hi = hi;
      this.blockSize = blockSize;
      this.siblings = new CountedLogSumExp[maxSiblings];
    }

    @Override
    public void compute() {
      final int l = lo;
      int h = hi;
      int j = 0;
      while (h - l > blockSize) {
        final int mid = (l + h + 1) / 2;
        addToPendingCount(1);
        siblings[j++] = (CountedLogSumExp) new CountedLogSumExp(this,
                                                                array, offset, mid, h, blockSize,
                                                                siblings.length - j).fork();
        h = mid;
      }
      
      if (h > l) {
        double acc = 0.;
        double c = 0.;  // compensation.
        for (int i = l; i < h; i++) {
          final double y = FastMath.exp(array[i] - offset) - c;
          final double t = acc + y;
          c = (t - acc) - y;
          acc = t;
        }

        sum = acc;
      }
      tryComplete();
    }

    @Override
    public void onCompletion(final CountedCompleter<?> caller) {
      if (caller != this) {
        for (final CountedLogSumExp sibling : siblings) {
          if (sibling != null) {
            sum += sibling.sum;
          }
        }
      }
    }

    @Override
    public Double getRawResult() {
      return Math.log(sum) + offset;
    }

    public static double apply(final double[] array, final int pieces) {
      final int blockSize = (array.length + pieces - 1) / pieces;
      final int siblingCount = (int) (Math.log(pieces) / Math.log(2));
      return new CountedLogSumExp(array, blockSize, siblingCount).invoke();
    }
  }

old.java

  public static class CountedLogSumExp extends CountedCompleter<Double> {
    final int lo;
    final int hi;
    final double[] array;
    final int blockSize;
    double sum;
    double offset;

    CountedLogSumExp[] siblings;

    CountedLogSumExp(final double[] array, final int blockSize, final int maxSiblings) {
      this(null, array, 0, array.length, blockSize, maxSiblings);
    }

    CountedLogSumExp(final CountedCompleter<?> parent,
                     final double[] array, final int lo, final int hi,
                     final int blockSize, final int maxSiblings) {
      super(parent);
      this.array = array;
      this.lo = lo;
      this.hi = hi;
      this.blockSize = blockSize;
      this.siblings = new CountedLogSumExp[maxSiblings];
    }

    @Override
    public void compute() {
      final int l = lo;
      int h = hi;
      int j = 0;
      while (h - l > blockSize) {
        final int mid = (l + h + 1) / 2;
        addToPendingCount(1);
        siblings[j++] = (CountedLogSumExp) new CountedLogSumExp(this,
                                                                array, mid, h, blockSize,
                                                                siblings.length - j).fork();
        h = mid;
      }
      if (h > l) {
        double acc = 0.;
        double c = 0.;  // compensation.
        double offset = Double.MIN_VALUE;
        for (int i = l; i < h; i++) {
          if (array[i] > offset) {
            offset = array[i];
          }
        }
        for (int i = l; i < h; i++) {
          final double y = FastMath.exp(array[i] - offset) - c;
          final double t = acc + y;
          c = (t - acc) - y;
          acc = t;
        }
        this.sum = acc;
        this.offset = offset;
      }
      tryComplete();
    }

    @Override
    public void onCompletion(final CountedCompleter<?> caller) {
      if (caller != this) {
        double maxOffset = this.offset;
        for (int i = 0; i < siblings.length; i++) {
          if (siblings[i] != null) {
            if (siblings[i].offset > maxOffset) {
              maxOffset = siblings[i].offset;
            }
          }
        }
        sum = Math.exp(Math.log(sum) + offset - maxOffset);
        offset = maxOffset;
        for (final CountedLogSumExp sibling : siblings) {
          if (sibling != null) {
            sum += Math.exp(Math.log(sibling.sum) + sibling.offset - offset);
          }
        }
      }
    }

    @Override
    public Double getRawResult() {
      return Math.log(sum) + offset;
    }

    public static double apply(final double[] array, final int pieces) {
      final int blockSize = (array.length + pieces - 1) / pieces;
      final int siblingCount = new Double(Math.floor(Math.log(1. / pieces) / Math.log(0.5))).intValue();
      return new CountedLogSumExp(array, blockSize, siblingCount).invoke();
    }
  }