1  pragma solidity >=0.6.0;
  2  
  3  import "./SafeMath.sol";
  4  import "./SignedSafeMath.sol";
  5  import "./FixedPointMathLib.sol";
  6  
  7  library MoreMath {
  8  
  9      using SafeMath for uint;
 10      using SignedSafeMath for int;
 11  
 12      using FixedPointMathLib for int256;
 13      using FixedPointMathLib for uint256;
 14  
 15      uint256 internal constant HALF_WAD = 0.5 ether;
 16      uint256 internal constant PI = 3_141592653589793238;
 17      int256 internal constant SQRT_2PI = 2_506628274631000502;
 18      int256 internal constant SIGN = -1;
 19      int256 internal constant SCALAR = 1e18;
 20      int256 internal constant HALF_SCALAR = 1e9;
 21      int256 internal constant SCALAR_SQRD = 1e36;
 22      int256 internal constant HALF = 5e17;
 23      int256 internal constant ONE = 1e18;
 24      int256 internal constant TWO = 2e18;
 25      int256 internal constant NEGATIVE_TWO = -2e18;
 26      int256 internal constant SQRT2 = 1_414213562373095048; // √2 with 18 decimals of precision.
 27      int256 internal constant ERFC_A = 1_265512230000000000;
 28      int256 internal constant ERFC_B = 1_000023680000000000;
 29      int256 internal constant ERFC_C = 374091960000000000; // 1e-1
 30      int256 internal constant ERFC_D = 96784180000000000; // 1e-2
 31      int256 internal constant ERFC_E = -186288060000000000; // 1e-1
 32      int256 internal constant ERFC_F = 278868070000000000; // 1e-1
 33      int256 internal constant ERFC_G = -1_135203980000000000;
 34      int256 internal constant ERFC_H = 1_488515870000000000;
 35      int256 internal constant ERFC_I = -822152230000000000; // 1e-1
 36      int256 internal constant ERFC_J = 170872770000000000; // 1e-1
 37      int256 internal constant IERFC_A = -707110000000000000; // 1e-1
 38      int256 internal constant IERFC_B = 2_307530000000000000;
 39      int256 internal constant IERFC_C = 270610000000000000; // 1e-1
 40      int256 internal constant IERFC_D = 992290000000000000; // 1e-1
 41      int256 internal constant IERFC_E = 44810000000000000; // 1e-2
 42      int256 internal constant IERFC_F = 1_128379167095512570;
 43  
 44      //see: https://ethereum.stackexchange.com/questions/8086/logarithm-math-operation-in-solidity
 45  
 46      /// @dev log2(e) as a signed 59.18-decimal fixed-point number.
 47      int256 internal constant LOG2_E = 1_442695040888963407;
 48  
 49      /// @dev Half the SCALE number.
 50      int256 internal constant HALF_SCALE = 5e17;
 51  
 52      int256 internal constant SCALE = 1e18;
 53  
 54      int256 internal constant OLD_PI = 3_141592653589793238;
 55  
 56      /// @notice Finds the zero-based index of the first one in the binary representation of x.
 57      /// @dev See the note on msb in the "Find First Set" Wikipedia article https://en.wikipedia.org/wiki/Find_first_set
 58      /// @param x The uint256 number for which to find the index of the most significant bit.
 59      /// @return msb The index of the most significant bit as an uint256.
 60      function mostSignificantBit(uint256 x) internal pure returns (uint256 msb) {
 61          if (x >= 2**128) {
 62              x >>= 128;
 63              msb += 128;
 64          }
 65          if (x >= 2**64) {
 66              x >>= 64;
 67              msb += 64;
 68          }
 69          if (x >= 2**32) {
 70              x >>= 32;
 71              msb += 32;
 72          }
 73          if (x >= 2**16) {
 74              x >>= 16;
 75              msb += 16;
 76          }
 77          if (x >= 2**8) {
 78              x >>= 8;
 79              msb += 8;
 80          }
 81          if (x >= 2**4) {
 82              x >>= 4;
 83              msb += 4;
 84          }
 85          if (x >= 2**2) {
 86              x >>= 2;
 87              msb += 2;
 88          }
 89          if (x >= 2**1) {
 90              // No need to shift x any more.
 91              msb += 1;
 92          }
 93      }
 94      /// @notice Calculates the binary logarithm of x.
 95      ///
 96      /// @dev Based on the iterative approximation algorithm.
 97      /// https://en.wikipedia.org/wiki/Binary_logarithm#Iterative_approximation
 98      ///
 99      /// Requirements:
100      /// - x must be greater than zero.
101      ///
102      /// Caveats:
103      /// - The results are nor perfectly accurate to the last digit, due to the lossy precision of the iterative approximation.
104      ///
105      /// @param x The signed 59.18-decimal fixed-point number for which to calculate the binary logarithm.
106      /// @return result The binary logarithm as a signed 59.18-decimal fixed-point number.
107      function old_log2(int256 x) internal pure returns (int256 result) {
108          require(x > 0);
109          // This works because log2(x) = -log2(1/x).
110          int256 sign;
111          if (x >= SCALE) {
112              sign = 1;
113          } else {
114              sign = -1;
115              // Do the fixed-point inversion inline to save gas. The numerator is SCALE * SCALE.
116              assembly {
117                  x := div(1000000000000000000000000000000000000, x)
118              }
119          }
120  
121          // Calculate the integer part of the logarithm and add it to the result and finally calculate y = x * 2^(-n).
122          uint256 n = mostSignificantBit(uint256(x / SCALE));
123  
124          // The integer part of the logarithm as a signed 59.18-decimal fixed-point number. The operation can't overflow
125          // because n is maximum 255, SCALE is 1e18 and sign is either 1 or -1.
126          result = int256(n) * SCALE;
127  
128          // This is y = x * 2^(-n).
129          int256 y = x >> n;
130  
131          // If y = 1, the fractional part is zero.
132          if (y == SCALE) {
133              return result * sign;
134          }
135  
136          // Calculate the fractional part via the iterative approximation.
137          // The "delta >>= 1" part is equivalent to "delta /= 2", but shifting bits is faster.
138          for (int256 delta = int256(HALF_SCALE); delta > 0; delta >>= 1) {
139              y = (y * y) / SCALE;
140  
141              // Is y^2 > 2 and so in the range [2,4)?
142              if (y >= 2 * SCALE) {
143                  // Add the 2^(-m) factor to the logarithm.
144                  result += delta;
145  
146                  // Corresponds to z/2 on Wikipedia.
147                  y >>= 1;
148              }
149          }
150          result *= sign;
151      }
152  
153      /// @notice Calculates the natural logarithm of x.
154      ///
155      /// @dev Based on the insight that ln(x) = log2(x) / log2(e).
156      ///
157      /// Requirements:
158      /// - All from "log2".
159      ///
160      /// Caveats:
161      /// - All from "log2".
162      /// - This doesn't return exactly 1 for 2718281828459045235, for that we would need more fine-grained precision.
163      ///
164      /// @param x The signed 59.18-decimal fixed-point number for which to calculate the natural logarithm.
165      /// @return result The natural logarithm as a signed 59.18-decimal fixed-point number.
166      function ln(int256 x) internal pure returns (int256 result) {
167          // Do the fixed-point multiplication inline to save gas. This is overflow-safe because the maximum value that log2(x)
168          // can return is 195205294292027477728.
169          result = (old_log2(x) * SCALE) / LOG2_E;
170      }
171  
172      /*
173       * notice Approximation of the Complimentary Error Function.
174       * Related to the Error Function: `erfc(x) = 1 - erf(x)`.
175       * Both cumulative distribution and error functions are integrals
176       * which cannot be expressed in elementary terms. They are called special functions.
177       * The error and complimentary error functions have numerical approximations
178       * which is what is used in this library to compute the cumulative distribution function.
179       *
180       *  dev This is a special function with its own identities.
181       * Identity: `erfc(-x) = 2 - erfc(x)`.
182       * Special Values:
183       * erfc(-infinity)  =   2
184       * erfc(0)          =   1
185       * erfc(infinity)   =   0
186       *
187       * custom:epsilon Fractional error less than 1.2e-7.
188       * custom:source Numerical Recipes in C 2e p221.
189       * custom:source https://mathworld.wolfram.com/Erfc.html.
190       */
191      function erfc(int256 input) internal pure returns (int256 output) {
192          uint256 z = abs(input);
193          int256 t;
194          int256 step;
195          int256 k;
196          assembly {
197              let quo := sdiv(mul(z, ONE), TWO) // 1 / (1 + z / 2).
198              let den := add(ONE, quo)
199              t := sdiv(SCALAR_SQRD, den)
200  
201              function muli(pxn, pxd) -> res {
202                  res := sdiv(mul(pxn, pxd), ONE)
203              }
204  
205              {
206                  step := add(
207                      ERFC_F,
208                      muli(
209                          t,
210                          add(
211                              ERFC_G,
212                              muli(
213                                  t,
214                                  add(
215                                      ERFC_H,
216                                      muli(t, add(ERFC_I, muli(t, ERFC_J)))
217                                  )
218                              )
219                          )
220                      )
221                  )
222              }
223              {
224                  step := muli(
225                      t,
226                      add(
227                          ERFC_B,
228                          muli(
229                              t,
230                              add(
231                                  ERFC_C,
232                                  muli(
233                                      t,
234                                      add(
235                                          ERFC_D,
236                                          muli(t, add(ERFC_E, muli(t, step)))
237                                      )
238                                  )
239                              )
240                          )
241                      )
242                  )
243              }
244  
245              k := add(sub(mul(SIGN, muli(z, z)), ERFC_A), step)
246          }
247  
248          int256 expWad = FixedPointMathLib.expWad(k);
249          int256 r;
250          assembly {
251              r := sdiv(mul(t, expWad), ONE)
252              switch iszero(slt(input, 0))
253              case 0 {
254                  output := sub(TWO, r)
255              }
256              case 1 {
257                  output := r
258              }
259          }
260      }
261  
262      /**
263       * notice Approximation of the Cumulative Distribution Function.
264       *
265       * dev Equal to `D(x) = 0.5[ 1 + erf((x - µ) / σ√2)]`.
266       * Only computes cdf of a distribution with µ = 0 and σ = 1.
267       *
268       * custom:error Maximum error of 1.2e-7.
269       * custom:source https://mathworld.wolfram.com/NormalDistribution.html.
270       */
271      function cdf(int256 x) internal pure returns (int256 z) {
272          int256 negated;
273          assembly {
274              let res := sdiv(mul(x, ONE), SQRT2)
275              negated := add(not(res), 1)
276          }
277  
278          int256 _erfc = erfc(negated);
279          assembly {
280              z := sdiv(mul(ONE, _erfc), TWO)
281          }
282      }
283  
284      /*
285       * notice Approximation of the Probability Density Function.
286       *
287       * dev Equal to `Z(x) = (1 / σ√2π)e^( (-(x - µ)^2) / 2σ^2 )`.
288       * Only computes pdf of a distribution with µ = 0 and σ = 1.
289       *
290       * custom:error Maximum error of 1.2e-7.
291       * custom:source https://mathworld.wolfram.com/ProbabilityDensityFunction.html.
292       */
293      function pdf(int256 x) internal pure returns (int256 z) {
294          int256 e;
295          assembly {
296              e := sdiv(mul(add(not(x), 1), x), TWO) // (-x * x) / 2.
297          }
298          e = FixedPointMathLib.expWad(e);
299  
300          assembly {
301              z := sdiv(mul(e, ONE), SQRT_2PI)
302          }
303      }
304  
305  
306      // rounds "v" considering a base "b"
307      function round(uint v, uint b) internal pure returns (uint) {
308  
309          return v.div(b).add((v % b) >= b.div(2) ? 1 : 0);
310      }
311  
312      // calculates {[(n/d)^e]*f}
313      function powAndMultiply(uint n, uint d, uint e, uint f) internal pure returns (uint) {
314          
315          if (e == 0) {
316              return 1;
317          } else if (e == 1) {
318              return f.mul(n).div(d);
319          } else {
320              uint p = powAndMultiply(n, d, e.div(2), f);
321              p = p.mul(p).div(f);
322              if (e.mod(2) == 1) {
323                  p = p.mul(n).div(d);
324              }
325              return p;
326          }
327      }
328  
329      // calculates (n^e)
330      function pow(uint n, uint e) internal pure returns (uint) {
331          
332          if (e == 0) {
333              return 1;
334          } else if (e == 1) {
335              return n;
336          } else {
337              uint p = pow(n, e.div(2));
338              p = p.mul(p);
339              if (e.mod(2) == 1) {
340                  p = p.mul(n);
341              }
342              return p;
343          }
344      }
345  
346      // calculates {n^(e/b)}
347      function powDecimal(uint n, uint e, uint b) internal pure returns (uint v) {
348          
349          if (e == 0) {
350              return b;
351          }
352  
353          if (e > b) {
354              return n.mul(powDecimal(n, e.sub(b), b)).div(b);
355          }
356  
357          v = b;
358          uint f = b;
359          uint aux = 0;
360          uint rootN = n;
361          uint rootB = sqrt(b);
362          while (f > 1) {
363              f = f.div(2);
364              rootN = sqrt(rootN).mul(rootB);
365              if (aux.add(f) < e) {
366                  aux = aux.add(f);
367                  v = v.mul(rootN).div(b);
368              }
369          }
370      }
371      
372      // calculates ceil(n/d)
373      function divCeil(uint n, uint d) internal pure returns (uint v) {
374          
375          v = n.div(d);
376          if (n.mod(d) > 0) {
377              v = v.add(1);
378          }
379      }
380      
381      // calculates the square root of "x" and multiplies it by "f"
382      function sqrtAndMultiply(uint x, uint f) internal pure returns (uint y) {
383      
384          y = sqrt(x.mul(1e18)).mul(f).div(1e9);
385      }
386      
387      // calculates the square root of "x"
388      function sqrt(uint x) internal pure returns (uint y) {
389      
390          uint z = (x.div(2)).add(1);
391          y = x;
392          while (z < y) {
393              y = z;
394              z = (x.div(z).add(z)).div(2);
395          }
396      }
397  
398      // calculates the standard deviation
399      function std(int[] memory array) internal pure returns (uint _std) {
400  
401          int avg = sum(array).div(int(array.length));
402          uint x2 = 0;
403          for (uint i = 0; i < array.length; i++) {
404              int p = array[i].sub(avg);
405              x2 = x2.add(uint(p.mul(p)));
406          }
407          _std = sqrt(x2 / array.length);
408      }
409  
410      function sum(int[] memory array) internal pure returns (int _sum) {
411  
412          for (uint i = 0; i < array.length; i++) {
413              _sum = _sum.add(array[i]);
414          }
415      }
416  
417      function abs(int a) internal pure returns (uint) {
418  
419          return uint(a < 0 ? -a : a);
420      }
421      
422      function max(int a, int b) internal pure returns (int) {
423          
424          return a > b ? a : b;
425      }
426      
427      function max(uint a, uint b) internal pure returns (uint) {
428          
429          return a > b ? a : b;
430      }
431      
432      function min(int a, int b) internal pure returns (int) {
433          
434          return a < b ? a : b;
435      }
436      
437      function min(uint a, uint b) internal pure returns (uint) {
438          
439          return a < b ? a : b;
440      }
441  
442      function toString(uint v) internal pure returns (string memory str) {
443  
444          str = toString(v, true);
445      }
446      
447      function toString(uint v, bool scientific) internal pure returns (string memory str) {
448  
449          if (v == 0) {
450              return "0";
451          }
452  
453          uint maxlength = 100;
454          bytes memory reversed = new bytes(maxlength);
455          uint i = 0;
456          
457          while (v != 0) {
458              uint remainder = v % 10;
459              v = v / 10;
460              reversed[i++] = byte(uint8(48 + remainder));
461          }
462  
463          uint zeros = 0;
464          if (scientific) {
465              for (uint k = 0; k < i; k++) {
466                  if (reversed[k] == '0') {
467                      zeros++;
468                  } else {
469                      break;
470                  }
471              }
472          }
473  
474          uint len = i - (zeros > 2 ? zeros : 0);
475          bytes memory s = new bytes(len);
476          for (uint j = 0; j < len; j++) {
477              s[j] = reversed[i - j - 1];
478          }
479  
480          str = string(s);
481  
482          if (scientific && zeros > 2) {
483              str = string(abi.encodePacked(s, "e", toString(zeros, false)));
484          }
485      }
486  }