esnext.math.sum-precise.js 4.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151
  1. 'use strict';
  2. // based on Shewchuk's algorithm for exactly floating point addition
  3. // adapted from https://github.com/tc39/proposal-math-sum/blob/3513d58323a1ae25560e8700aa5294500c6c9287/polyfill/polyfill.mjs
  4. var $ = require('../internals/export');
  5. var uncurryThis = require('../internals/function-uncurry-this');
  6. var iterate = require('../internals/iterate');
  7. var $RangeError = RangeError;
  8. var $TypeError = TypeError;
  9. var $Infinity = Infinity;
  10. var $NaN = NaN;
  11. var abs = Math.abs;
  12. var pow = Math.pow;
  13. var push = uncurryThis([].push);
  14. var POW_2_1023 = pow(2, 1023);
  15. var MAX_SAFE_INTEGER = pow(2, 53) - 1; // 2 ** 53 - 1 === 9007199254740992
  16. var MAX_DOUBLE = Number.MAX_VALUE; // 2 ** 1024 - 2 ** (1023 - 52) === 1.79769313486231570815e+308
  17. var MAX_ULP = pow(2, 971); // 2 ** (1023 - 52) === 1.99584030953471981166e+292
  18. var NOT_A_NUMBER = {};
  19. var MINUS_INFINITY = {};
  20. var PLUS_INFINITY = {};
  21. var MINUS_ZERO = {};
  22. var FINITE = {};
  23. // prerequisite: abs(x) >= abs(y)
  24. var twosum = function (x, y) {
  25. var hi = x + y;
  26. var lo = y - (hi - x);
  27. return { hi: hi, lo: lo };
  28. };
  29. // `Math.sumPrecise` method
  30. // https://github.com/tc39/proposal-math-sum
  31. $({ target: 'Math', stat: true }, {
  32. // eslint-disable-next-line max-statements -- ok
  33. sumPrecise: function sumPrecise(items) {
  34. var numbers = [];
  35. var count = 0;
  36. var state = MINUS_ZERO;
  37. iterate(items, function (n) {
  38. if (++count >= MAX_SAFE_INTEGER) throw new $RangeError('Maximum allowed index exceeded');
  39. if (typeof n != 'number') throw new $TypeError('Value is not a number');
  40. if (state !== NOT_A_NUMBER) {
  41. // eslint-disable-next-line no-self-compare -- NaN check
  42. if (n !== n) state = NOT_A_NUMBER;
  43. else if (n === $Infinity) state = state === MINUS_INFINITY ? NOT_A_NUMBER : PLUS_INFINITY;
  44. else if (n === -$Infinity) state = state === PLUS_INFINITY ? NOT_A_NUMBER : MINUS_INFINITY;
  45. else if ((n !== 0 || (1 / n) === $Infinity) && (state === MINUS_ZERO || state === FINITE)) {
  46. state = FINITE;
  47. push(numbers, n);
  48. }
  49. }
  50. });
  51. switch (state) {
  52. case NOT_A_NUMBER: return $NaN;
  53. case MINUS_INFINITY: return -$Infinity;
  54. case PLUS_INFINITY: return $Infinity;
  55. case MINUS_ZERO: return -0;
  56. }
  57. var partials = [];
  58. var overflow = 0; // conceptually 2 ** 1024 times this value; the final partial is biased by this amount
  59. var x, y, sum, hi, lo, tmp;
  60. for (var i = 0; i < numbers.length; i++) {
  61. x = numbers[i];
  62. var actuallyUsedPartials = 0;
  63. for (var j = 0; j < partials.length; j++) {
  64. y = partials[j];
  65. if (abs(x) < abs(y)) {
  66. tmp = x;
  67. x = y;
  68. y = tmp;
  69. }
  70. sum = twosum(x, y);
  71. hi = sum.hi;
  72. lo = sum.lo;
  73. if (abs(hi) === $Infinity) {
  74. var sign = hi === $Infinity ? 1 : -1;
  75. overflow += sign;
  76. x = (x - (sign * POW_2_1023)) - (sign * POW_2_1023);
  77. if (abs(x) < abs(y)) {
  78. tmp = x;
  79. x = y;
  80. y = tmp;
  81. }
  82. sum = twosum(x, y);
  83. hi = sum.hi;
  84. lo = sum.lo;
  85. }
  86. if (lo !== 0) partials[actuallyUsedPartials++] = lo;
  87. x = hi;
  88. }
  89. partials.length = actuallyUsedPartials;
  90. if (x !== 0) push(partials, x);
  91. }
  92. // compute the exact sum of partials, stopping once we lose precision
  93. var n = partials.length - 1;
  94. hi = 0;
  95. lo = 0;
  96. if (overflow !== 0) {
  97. var next = n >= 0 ? partials[n] : 0;
  98. n--;
  99. if (abs(overflow) > 1 || (overflow > 0 && next > 0) || (overflow < 0 && next < 0)) {
  100. return overflow > 0 ? $Infinity : -$Infinity;
  101. }
  102. // here we actually have to do the arithmetic
  103. // drop a factor of 2 so we can do it without overflow
  104. // assert(abs(overflow) === 1)
  105. sum = twosum(overflow * POW_2_1023, next / 2);
  106. hi = sum.hi;
  107. lo = sum.lo;
  108. lo *= 2;
  109. if (abs(2 * hi) === $Infinity) {
  110. // rounding to the maximum value
  111. if (hi > 0) {
  112. return (hi === POW_2_1023 && lo === -(MAX_ULP / 2) && n >= 0 && partials[n] < 0) ? MAX_DOUBLE : $Infinity;
  113. } return (hi === -POW_2_1023 && lo === (MAX_ULP / 2) && n >= 0 && partials[n] > 0) ? -MAX_DOUBLE : -$Infinity;
  114. }
  115. if (lo !== 0) {
  116. partials[++n] = lo;
  117. lo = 0;
  118. }
  119. hi *= 2;
  120. }
  121. while (n >= 0) {
  122. sum = twosum(hi, partials[n--]);
  123. hi = sum.hi;
  124. lo = sum.lo;
  125. if (lo !== 0) break;
  126. }
  127. if (n >= 0 && ((lo < 0 && partials[n] < 0) || (lo > 0 && partials[n] > 0))) {
  128. y = lo * 2;
  129. x = hi + y;
  130. if (y === x - hi) hi = x;
  131. }
  132. return hi;
  133. }
  134. });