handwritten asm.js code for logsumexp

logsumexp-asm.js

// MIT License. Copyright (c) 2016 Hajime Senuma.
/**
 * asm.js module for logsumexp
 */
function asmModule(stdlib, foreign, heap) {
  'use asm';

  var exp = stdlib.Math.exp;
  var log = stdlib.Math.log;

  var f32heap = new stdlib.Float32Array(heap);

  /**
   * Returns the largest number of one or more 32-bit floats.
   *
   * @param {number} valueP - byte offset to the 32-bit float sequence.
   * @param {number} len - length of the float sequence.
   *   If this value is less than 1, the behavior is undefined.
   * @returns {number} - 64-bit float (not 32-bit)
   */
  function maxFloat32(valueP, len) {
    // Type annotations
    valueP = valueP | 0;
    len = len | 0;
  
    // Local variables
    var p = 0;
    var end = 0;
    var v = 0.0;
    var result = 0.0;
  
    // Main
    p = valueP;
    end = (valueP + (len << 2)) | 0;
    result = +f32heap[p >> 2];
    p = (p + 4) | 0;
  
    for (; (p | 0) < (end | 0); p = (p + 4) | 0) {
      v = +f32heap[p >> 2];

      if (v >= result) {
        result = v;
      }
    }
  
    return +result;
  }

  /**
   * Returns the logsumexp of one or more 32-bit floats.
   *
   * @param {number} valueP - byte offset to the 32-bit float sequence.
   * @param {number} len - length of the float sequence.
   *   If this value is less than 1, the behavior is undefined.
   * @returns {number} - 64-bit float (not 32-bit)
   */
  function logsumexp(valueP, len) {
    // Type annotations
    valueP = valueP | 0;
    len = len | 0;

    // Local variables
    var p = 0;
    var end = 0;
    var v = 0.0;
    var maxValue = 0.0;
    var result = 0.0;

    // Main
    maxValue = +maxFloat32(valueP, len);
    p = valueP;
    end = (valueP + (len << 2)) | 0;
  
    for (; (p | 0) < (end | 0); p = (p + 4) | 0) {
      v = +f32heap[p >> 2];

      // exp(-20) = 2.06e-9, machine epsilon for float32 = 5.96e-08
      if (v - maxValue > -16.0) {
        result = +(result + +exp(v - maxValue));
      }
    }

    return +(maxValue + +log(result));
  } 

  return {
    maxFloat32: maxFloat32,
    logsumexp: logsumexp
  };
}