sumpw.js (4678B)
1 /** 2 * @license Apache-2.0 3 * 4 * Copyright (c) 2020 The Stdlib Authors. 5 * 6 * Licensed under the Apache License, Version 2.0 (the "License"); 7 * you may not use this file except in compliance with the License. 8 * You may obtain a copy of the License at 9 * 10 * http://www.apache.org/licenses/LICENSE-2.0 11 * 12 * Unless required by applicable law or agreed to in writing, software 13 * distributed under the License is distributed on an "AS IS" BASIS, 14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 15 * See the License for the specific language governing permissions and 16 * limitations under the License. 17 */ 18 19 'use strict'; 20 21 // MODULES // 22 23 var isnan = require( '@stdlib/math/base/assert/is-nan' ); 24 var floor = require( '@stdlib/math/base/special/floor' ); 25 26 27 // VARIABLES // 28 29 // Blocksize for pairwise summation (NOTE: decreasing the blocksize decreases rounding error as more pairs are summed, but also decreases performance. Because the inner loop is unrolled eight times, the blocksize is effectively `16`.): 30 var BLOCKSIZE = 128; 31 32 33 // MAIN // 34 35 /** 36 * Computes the sum of a double-precision floating-point strided array elements, ignoring `NaN` values and using pairwise summation. 37 * 38 * ## Method 39 * 40 * - This implementation uses pairwise summation, which accrues rounding error `O(log2 N)` instead of `O(N)`. The recursion depth is also `O(log2 N)`. 41 * 42 * ## References 43 * 44 * - Higham, Nicholas J. 1993. "The Accuracy of Floating Point Summation." _SIAM Journal on Scientific Computing_ 14 (4): 783–99. doi:[10.1137/0914050](https://doi.org/10.1137/0914050). 45 * 46 * @private 47 * @param {PositiveInteger} N - number of indexed elements 48 * @param {Float64Array} x - input array 49 * @param {integer} strideX - `x` stride length 50 * @param {NonNegativeInteger} offsetX - `x` starting index 51 * @param {Float64Array} out - two-element output array whose first element is the accumulated sum and whose second element is the accumulated number of summed values 52 * @param {integer} strideOut - `out` stride length 53 * @param {NonNegativeInteger} offsetOut - `out` starting index 54 * @returns {Float64Array} output array 55 * 56 * @example 57 * var Float64Array = require( '@stdlib/array/float64' ); 58 * var floor = require( '@stdlib/math/base/special/floor' ); 59 * 60 * var x = new Float64Array( [ 2.0, 1.0, 2.0, -2.0, -2.0, 2.0, 3.0, 4.0, NaN, NaN ] ); 61 * var N = floor( x.length / 2 ); 62 * 63 * var out = new Float64Array( [ 0.0, 0 ] ); 64 * var v = sumpw( N, x, 2, 1, out, 1, 0 ); 65 * // returns <Float64Array>[ 5.0, 4 ] 66 */ 67 function sumpw( N, x, strideX, offsetX, out, strideOut, offsetOut ) { 68 var ix; 69 var io; 70 var s0; 71 var s1; 72 var s2; 73 var s3; 74 var s4; 75 var s5; 76 var s6; 77 var s7; 78 var M; 79 var s; 80 var n; 81 var v; 82 var i; 83 84 if ( N <= 0 ) { 85 return out; 86 } 87 ix = offsetX; 88 io = offsetOut; 89 if ( N === 1 || strideX === 0 ) { 90 if ( isnan( x[ ix ] ) ) { 91 return out; 92 } 93 out[ io ] += x[ ix ]; 94 out[ io+strideOut ] += 1; 95 return out; 96 } 97 if ( N < 8 ) { 98 // Use simple summation... 99 s = 0.0; 100 n = 0; 101 for ( i = 0; i < N; i++ ) { 102 v = x[ ix ]; 103 if ( v === v ) { 104 s += v; 105 n += 1; 106 } 107 ix += strideX; 108 } 109 out[ io ] += s; 110 out[ io+strideOut ] += n; 111 return out; 112 } 113 if ( N <= BLOCKSIZE ) { 114 // Sum a block with 8 accumulators (by loop unrolling, we lower the effective blocksize to 16)... 115 s0 = 0.0; 116 s1 = 0.0; 117 s2 = 0.0; 118 s3 = 0.0; 119 s4 = 0.0; 120 s5 = 0.0; 121 s6 = 0.0; 122 s7 = 0.0; 123 n = 0; 124 125 M = N % 8; 126 for ( i = 0; i < N-M; i += 8 ) { 127 v = x[ ix ]; 128 if ( v === v ) { 129 s0 += v; 130 n += 1; 131 } 132 ix += strideX; 133 v = x[ ix ]; 134 if ( v === v ) { 135 s1 += v; 136 n += 1; 137 } 138 ix += strideX; 139 v = x[ ix ]; 140 if ( v === v ) { 141 s2 += v; 142 n += 1; 143 } 144 ix += strideX; 145 v = x[ ix ]; 146 if ( v === v ) { 147 s3 += v; 148 n += 1; 149 } 150 ix += strideX; 151 v = x[ ix ]; 152 if ( v === v ) { 153 s4 += v; 154 n += 1; 155 } 156 ix += strideX; 157 v = x[ ix ]; 158 if ( v === v ) { 159 s5 += v; 160 n += 1; 161 } 162 ix += strideX; 163 v = x[ ix ]; 164 if ( v === v ) { 165 s6 += v; 166 n += 1; 167 } 168 ix += strideX; 169 v = x[ ix ]; 170 if ( v === v ) { 171 s7 += v; 172 n += 1; 173 } 174 ix += strideX; 175 } 176 // Pairwise sum the accumulators: 177 s = ((s0+s1) + (s2+s3)) + ((s4+s5) + (s6+s7)); 178 179 // Clean-up loop... 180 for ( i; i < N; i++ ) { 181 v = x[ ix ]; 182 if ( v === v ) { 183 s += v; 184 n += 1; 185 } 186 ix += strideX; 187 } 188 out[ io ] += s; 189 out[ io+strideOut ] += n; 190 return out; 191 } 192 // Recurse by dividing by two, but avoiding non-multiples of unroll factor... 193 n = floor( N/2 ); 194 n -= n % 8; 195 sumpw( n, x, strideX, ix, out, strideOut, offsetOut ); 196 sumpw( N-n, x, strideX, ix+(n*strideX), out, strideOut, offsetOut ); 197 return out; 198 } 199 200 201 // EXPORTS // 202 203 module.exports = sumpw;