選択できるのは25トピックまでです。 トピックは、先頭が英数字で、英数字とダッシュ('-')を使用した35文字以内のものにしてください。
 
 
 

178 行
3.4 KiB

  1. (function (root, factory) {
  2. "use strict";
  3. // AMD
  4. if (typeof define === 'function' && define.amd) {
  5. define([], factory);
  6. }
  7. // CommonJS
  8. else if (typeof exports === 'object') {
  9. module.exports = factory();
  10. }
  11. // Browser
  12. else {
  13. root.add = factory();
  14. }
  15. })(this, function () {
  16. "use strict";
  17. // The minimum machine rounding error
  18. var Epsilon = Math.pow(2, -53)
  19. , EpsilonReciprocal = (1 / Epsilon)
  20. /// The smallest positive number that can be represented
  21. , Eta = Math.pow(2, -1074)
  22. // limitB is a constant used in the transform function
  23. , limitB = 0.5 * EpsilonReciprocal * Eta
  24. /**
  25. * S. M. RUMP, T. OGITA AND S. OISHI
  26. * http://www.ti3.tu-harburg.de/paper/rump/RuOgOi07I.pdf
  27. */
  28. // Page 8
  29. // x is result, y is error
  30. // third is so the array is allocated for 4 spaces
  31. // it speeds up transform
  32. function fastTwoSum(a, b) {
  33. var x = a + b
  34. , q = x - a
  35. , y = b - q
  36. return [x, y, null]
  37. }
  38. // Page 12
  39. // p = q + p'
  40. // sigma is a power of 2 greater than or equal to |p|
  41. function extractScalar(sigma, p) {
  42. var q = (sigma + p) - sigma
  43. , pPrime = p - q
  44. return [q, pPrime]
  45. }
  46. // Page 12
  47. function extractVector(sigma, p) {
  48. var tau = 0.0
  49. , extracted
  50. , i = 0
  51. , ii = p.length
  52. , pPrime = new Array(ii)
  53. for(; i<ii; ++i) {
  54. extracted = extractScalar(sigma, p[i])
  55. pPrime[i] = extracted[1]
  56. tau += extracted[0]
  57. }
  58. return [tau, pPrime]
  59. }
  60. // Finds the immediate power of 2 that is larger than p
  61. //// in a fast way
  62. function nextPowerTwo (p) {
  63. var q = EpsilonReciprocal * p
  64. , L = Math.abs((q + p) - q)
  65. if(L === 0)
  66. return Math.abs(p)
  67. return L
  68. }
  69. // Helper, gets the maximum of the absolute values of an array
  70. function maxAbs(arr) {
  71. var i = 0
  72. , ii = arr.length
  73. , best = -1
  74. for(; i<ii; ++i) {
  75. if(Math.abs(arr[i]) > best) {
  76. best = arr[i]
  77. }
  78. }
  79. return best
  80. }
  81. function transform (p) {
  82. var mu = maxAbs(p)
  83. , M
  84. , sigmaPrime
  85. , tPrime
  86. , t
  87. , tau
  88. , sigma
  89. , extracted
  90. , res
  91. // Not part of the original paper, here for optimization
  92. , temp
  93. , bigPow
  94. , limitA
  95. , twoToTheM
  96. if(mu === 0) {
  97. return [0, 0, p, 0]
  98. }
  99. M = nextPowerTwo(p.length + 2)
  100. twoToTheM = Math.pow(2, M)
  101. bigPow = 2 * twoToTheM // equiv to Math.pow(2, 2 * M), faster
  102. sigmaPrime = twoToTheM * nextPowerTwo(mu)
  103. tPrime = 0
  104. do {
  105. t = tPrime
  106. sigma = sigmaPrime
  107. extracted = extractVector(sigma, p)
  108. tau = extracted[0]
  109. tPrime = t + tau
  110. p = extracted[1]
  111. if(tPrime === 0) {
  112. return transform(p)
  113. }
  114. temp = Epsilon * sigma
  115. sigmaPrime = twoToTheM * temp
  116. limitA = bigPow * temp
  117. }
  118. while( Math.abs(tPrime) < limitA && sigma > limitB )
  119. // res already allocated for 4
  120. res = fastTwoSum(t, tau)
  121. res[2] = p
  122. return res
  123. }
  124. function dumbSum(p) {
  125. var i, ii, sum = 0.0
  126. for(i=0, ii=p.length; i<ii; ++i) {
  127. sum += p[i]
  128. }
  129. return sum
  130. }
  131. function accSum(p) {
  132. // Zero length array, or all values are zeros
  133. if(p.length === 0 || maxAbs(p) === 0) {
  134. return 0
  135. }
  136. var tfmd = transform(p)
  137. return tfmd[0] + (tfmd[1] +dumbSum(tfmd[2]))
  138. }
  139. // exports
  140. accSum.dumbSum = dumbSum;
  141. accSum.fastTwoSum = fastTwoSum;
  142. accSum.nextPowerTwo = nextPowerTwo;
  143. return accSum;
  144. });