Merge pull request #2 from nim-lang/devel

Update

1 files changed, 212 insertions, 80 deletions

diff --git a/lib/pure/math.nim b/lib/pure/math.nim
index cbd04a145..8ea8ee203 100644
--- a/lib/pure/math.nim
+++ b/lib/pure/math.nim
@@ -21,6 +21,8 @@ include "system/inclrtl"
 {.push debugger:off .} # the user does not want to trace a part
                        # of the standard library!
 
+import bitops
+
 proc binom*(n, k: int): int {.noSideEffect.} =
   ## Computes the binomial coefficient
   if k <= 0: return 1
@@ -29,11 +31,21 @@ proc binom*(n, k: int): int {.noSideEffect.} =
   for i in countup(2, k):
     result = (result * (n + 1 - i)) div i
 
-proc fac*(n: int): int {.noSideEffect.} =
+proc createFactTable[N: static[int]]: array[N, int] =
+  result[0] = 1
+  for i in 1 ..< N:
+    result[i] = result[i - 1] * i
+
+proc fac*(n: int): int =
   ## Computes the faculty/factorial function.
-  result = 1
-  for i in countup(2, n):
-    result = result * i
+  const factTable =
+    when sizeof(int) == 4:
+      createFactTable[13]()
+    else:
+      createFactTable[21]()
+  assert(n >= 0, $n & " must not be negative.")
+  assert(n < factTable.len, $n & " is too large to look up in the table")
+  factTable[n]
 
 {.push checks:off, line_dir:off, stack_trace:off.}
 
@@ -117,8 +129,14 @@ proc sum*[T](x: openArray[T]): T {.noSideEffect.} =
   ## If `x` is empty, 0 is returned.
   for i in items(x): result = result + i
 
+proc prod*[T](x: openArray[T]): T {.noSideEffect.} =
+  ## Computes the product of the elements in ``x``.
+  ## If ``x`` is empty, 1 is returned.
+  result = 1.T
+  for i in items(x): result = result * i
+
 {.push noSideEffect.}
-when not defined(JS):
+when not defined(JS): # C
   proc sqrt*(x: float32): float32 {.importc: "sqrtf", header: "<math.h>".}
   proc sqrt*(x: float64): float64 {.importc: "sqrt", header: "<math.h>".}
     ## Computes the square root of `x`.
@@ -132,12 +150,33 @@ when not defined(JS):
   proc log10*(x: float32): float32 {.importc: "log10f", header: "<math.h>".}
   proc log10*(x: float64): float64 {.importc: "log10", header: "<math.h>".}
     ## Computes the common logarithm (base 10) of `x`
-  proc log2*[T: float32|float64](x: T): T = return ln(x) / ln(2.0)
+  proc log2*(x: float32): float32 {.importc: "log2f", header: "<math.h>".}
+  proc log2*(x: float64): float64 {.importc: "log2", header: "<math.h>".}
     ## Computes the binary logarithm (base 2) of `x`
   proc exp*(x: float32): float32 {.importc: "expf", header: "<math.h>".}
   proc exp*(x: float64): float64 {.importc: "exp", header: "<math.h>".}
     ## Computes the exponential function of `x` (pow(E, x))
 
+  proc sin*(x: float32): float32 {.importc: "sinf", header: "<math.h>".}
+  proc sin*(x: float64): float64 {.importc: "sin", header: "<math.h>".}
+    ## Computes the sine of `x`
+  proc cos*(x: float32): float32 {.importc: "cosf", header: "<math.h>".}
+  proc cos*(x: float64): float64 {.importc: "cos", header: "<math.h>".}
+    ## Computes the cosine of `x`
+  proc tan*(x: float32): float32 {.importc: "tanf", header: "<math.h>".}
+  proc tan*(x: float64): float64 {.importc: "tan", header: "<math.h>".}
+    ## Computes the tangent of `x`
+
+  proc sinh*(x: float32): float32 {.importc: "sinhf", header: "<math.h>".}
+  proc sinh*(x: float64): float64 {.importc: "sinh", header: "<math.h>".}
+    ## Computes the hyperbolic sine of `x`
+  proc cosh*(x: float32): float32 {.importc: "coshf", header: "<math.h>".}
+  proc cosh*(x: float64): float64 {.importc: "cosh", header: "<math.h>".}
+    ## Computes the hyperbolic cosine of `x`
+  proc tanh*(x: float32): float32 {.importc: "tanhf", header: "<math.h>".}
+  proc tanh*(x: float64): float64 {.importc: "tanh", header: "<math.h>".}
+    ## Computes the hyperbolic tangent of `x`
+
   proc arccos*(x: float32): float32 {.importc: "acosf", header: "<math.h>".}
   proc arccos*(x: float64): float64 {.importc: "acos", header: "<math.h>".}
     ## Computes the arc cosine of `x`
@@ -154,33 +193,80 @@ when not defined(JS):
     ## results even when the resulting angle is near pi/2 or -pi/2
     ## (`x` near 0).
 
-  proc cos*(x: float32): float32 {.importc: "cosf", header: "<math.h>".}
-  proc cos*(x: float64): float64 {.importc: "cos", header: "<math.h>".}
-    ## Computes the cosine of `x`
+  proc arcsinh*(x: float32): float32 {.importc: "asinhf", header: "<math.h>".}
+  proc arcsinh*(x: float64): float64 {.importc: "asinh", header: "<math.h>".}
+    ## Computes the inverse hyperbolic sine of `x`
+  proc arccosh*(x: float32): float32 {.importc: "acoshf", header: "<math.h>".}
+  proc arccosh*(x: float64): float64 {.importc: "acosh", header: "<math.h>".}
+    ## Computes the inverse hyperbolic cosine of `x`
+  proc arctanh*(x: float32): float32 {.importc: "atanhf", header: "<math.h>".}
+  proc arctanh*(x: float64): float64 {.importc: "atanh", header: "<math.h>".}
+    ## Computes the inverse hyperbolic tangent of `x`
+
+else: # JS
+  proc sqrt*(x: float32): float32 {.importc: "Math.sqrt", nodecl.}
+  proc sqrt*(x: float64): float64 {.importc: "Math.sqrt", nodecl.}
 
-  proc cosh*(x: float32): float32 {.importc: "coshf", header: "<math.h>".}
-  proc cosh*(x: float64): float64 {.importc: "cosh", header: "<math.h>".}
-    ## Computes the hyperbolic cosine of `x`
+  proc ln*(x: float32): float32 {.importc: "Math.log", nodecl.}
+  proc ln*(x: float64): float64 {.importc: "Math.log", nodecl.}
+  proc log10*(x: float32): float32 {.importc: "Math.log10", nodecl.}
+  proc log10*(x: float64): float64 {.importc: "Math.log10", nodecl.}
+  proc log2*(x: float32): float32 {.importc: "Math.log2", nodecl.}
+  proc log2*(x: float64): float64 {.importc: "Math.log2", nodecl.}
+  proc exp*(x: float32): float32 {.importc: "Math.exp", nodecl.}
+  proc exp*(x: float64): float64 {.importc: "Math.exp", nodecl.}
 
+  proc sin*[T: float32|float64](x: T): T {.importc: "Math.sin", nodecl.}
+  proc cos*[T: float32|float64](x: T): T {.importc: "Math.cos", nodecl.}
+  proc tan*[T: float32|float64](x: T): T {.importc: "Math.tan", nodecl.}
+
+  proc sinh*[T: float32|float64](x: T): T {.importc: "Math.sinh", nodecl.}
+  proc cosh*[T: float32|float64](x: T): T {.importc: "Math.cosh", nodecl.}
+  proc tanh*[T: float32|float64](x: T): T {.importc: "Math.tanh", nodecl.}
+
+  proc arcsin*[T: float32|float64](x: T): T {.importc: "Math.asin", nodecl.}
+  proc arccos*[T: float32|float64](x: T): T {.importc: "Math.acos", nodecl.}
+  proc arctan*[T: float32|float64](x: T): T {.importc: "Math.atan", nodecl.}
+  proc arctan2*[T: float32|float64](y, x: T): T {.importC: "Math.atan2", nodecl.}
+
+  proc arcsinh*[T: float32|float64](x: T): T {.importc: "Math.asinh", nodecl.}
+  proc arccosh*[T: float32|float64](x: T): T {.importc: "Math.acosh", nodecl.}
+  proc arctanh*[T: float32|float64](x: T): T {.importc: "Math.atanh", nodecl.}
+
+proc cot*[T: float32|float64](x: T): T = 1.0 / tan(x)
+  ## Computes the cotangent of `x`
+proc sec*[T: float32|float64](x: T): T = 1.0 / cos(x)
+  ## Computes the secant of `x`.
+proc csc*[T: float32|float64](x: T): T = 1.0 / sin(x)
+  ## Computes the cosecant of `x`
+
+proc coth*[T: float32|float64](x: T): T = 1.0 / tanh(x)
+  ## Computes the hyperbolic cotangent of `x`
+proc sech*[T: float32|float64](x: T): T = 1.0 / cosh(x)
+  ## Computes the hyperbolic secant of `x`
+proc csch*[T: float32|float64](x: T): T = 1.0 / sinh(x)
+  ## Computes the hyperbolic cosecant of `x`
+
+proc arccot*[T: float32|float64](x: T): T = arctan(1.0 / x)
+  ## Computes the inverse cotangent of `x`
+proc arcsec*[T: float32|float64](x: T): T = arccos(1.0 / x)
+  ## Computes the inverse secant of `x`
+proc arccsc*[T: float32|float64](x: T): T = arcsin(1.0 / x)
+  ## Computes the inverse cosecant of `x`
+
+proc arccoth*[T: float32|float64](x: T): T = arctanh(1.0 / x)
+  ## Computes the inverse hyperbolic cotangent of `x`
+proc arcsech*[T: float32|float64](x: T): T = arccosh(1.0 / x)
+  ## Computes the inverse hyperbolic secant of `x`
+proc arccsch*[T: float32|float64](x: T): T = arcsinh(1.0 / x)
+  ## Computes the inverse hyperbolic cosecant of `x`
+
+when not defined(JS): # C
   proc hypot*(x, y: float32): float32 {.importc: "hypotf", header: "<math.h>".}
   proc hypot*(x, y: float64): float64 {.importc: "hypot", header: "<math.h>".}
     ## Computes the hypotenuse of a right-angle triangle with `x` and
     ## `y` as its base and height. Equivalent to ``sqrt(x*x + y*y)``.
 
-  proc sinh*(x: float32): float32 {.importc: "sinhf", header: "<math.h>".}
-  proc sinh*(x: float64): float64 {.importc: "sinh", header: "<math.h>".}
-    ## Computes the hyperbolic sine of `x`
-  proc sin*(x: float32): float32 {.importc: "sinf", header: "<math.h>".}
-  proc sin*(x: float64): float64 {.importc: "sin", header: "<math.h>".}
-    ## Computes the sine of `x`
-
-  proc tan*(x: float32): float32 {.importc: "tanf", header: "<math.h>".}
-  proc tan*(x: float64): float64 {.importc: "tan", header: "<math.h>".}
-    ## Computes the tangent of `x`
-  proc tanh*(x: float32): float32 {.importc: "tanhf", header: "<math.h>".}
-  proc tanh*(x: float64): float64 {.importc: "tanh", header: "<math.h>".}
-    ## Computes the hyperbolic tangent of `x`
-
   proc pow*(x, y: float32): float32 {.importc: "powf", header: "<math.h>".}
   proc pow*(x, y: float64): float64 {.importc: "pow", header: "<math.h>".}
     ## computes x to power raised of y.
@@ -194,12 +280,18 @@ when not defined(JS):
   proc erfc*(x: float64): float64 {.importc: "erfc", header: "<math.h>".}
     ## The complementary error function
 
+  proc gamma*(x: float32): float32 {.importc: "tgammaf", header: "<math.h>".}
+  proc gamma*(x: float64): float64 {.importc: "tgamma", header: "<math.h>".}
+    ## The gamma function
+  proc tgamma*(x: float32): float32
+    {.deprecated: "use gamma instead", importc: "tgammaf", header: "<math.h>".}
+  proc tgamma*(x: float64): float64
+    {.deprecated: "use gamma instead", importc: "tgamma", header: "<math.h>".}
+    ## The gamma function
+    ## **Deprecated since version 0.19.0**: Use ``gamma`` instead.
   proc lgamma*(x: float32): float32 {.importc: "lgammaf", header: "<math.h>".}
   proc lgamma*(x: float64): float64 {.importc: "lgamma", header: "<math.h>".}
     ## Natural log of the gamma function
-  proc tgamma*(x: float32): float32 {.importc: "tgammaf", header: "<math.h>".}
-  proc tgamma*(x: float64): float64 {.importc: "tgamma", header: "<math.h>".}
-    ## The gamma function
 
   proc floor*(x: float32): float32 {.importc: "floorf", header: "<math.h>".}
   proc floor*(x: float64): float64 {.importc: "floor", header: "<math.h>".}
@@ -283,57 +375,31 @@ when not defined(JS):
       ## .. code-block:: nim
       ##  echo trunc(PI) # 3.0
 
-  proc fmod*(x, y: float32): float32 {.importc: "fmodf", header: "<math.h>".}
-  proc fmod*(x, y: float64): float64 {.importc: "fmod", header: "<math.h>".}
+  proc fmod*(x, y: float32): float32 {.deprecated, importc: "fmodf", header: "<math.h>".}
+  proc fmod*(x, y: float64): float64 {.deprecated, importc: "fmod", header: "<math.h>".}
     ## Computes the remainder of `x` divided by `y`
     ##
     ## .. code-block:: nim
     ##  echo fmod(-2.5, 0.3) ## -0.1
 
-else:
-  proc trunc*(x: float32): float32 {.importc: "Math.trunc", nodecl.}
-  proc trunc*(x: float64): float64 {.importc: "Math.trunc", nodecl.}
+  proc `mod`*(x, y: float32): float32 {.importc: "fmodf", header: "<math.h>".}
+  proc `mod`*(x, y: float64): float64 {.importc: "fmod", header: "<math.h>".}
+    ## Computes the modulo operation for float operators.
+else: # JS
+  proc hypot*[T: float32|float64](x, y: T): T = return sqrt(x*x + y*y)
+  proc pow*(x, y: float32): float32 {.importC: "Math.pow", nodecl.}
+  proc pow*(x, y: float64): float64 {.importc: "Math.pow", nodecl.}
   proc floor*(x: float32): float32 {.importc: "Math.floor", nodecl.}
   proc floor*(x: float64): float64 {.importc: "Math.floor", nodecl.}
   proc ceil*(x: float32): float32 {.importc: "Math.ceil", nodecl.}
   proc ceil*(x: float64): float64 {.importc: "Math.ceil", nodecl.}
-
-  proc sqrt*(x: float32): float32 {.importc: "Math.sqrt", nodecl.}
-  proc sqrt*(x: float64): float64 {.importc: "Math.sqrt", nodecl.}
-  proc ln*(x: float32): float32 {.importc: "Math.log", nodecl.}
-  proc ln*(x: float64): float64 {.importc: "Math.log", nodecl.}
-  proc log10*[T: float32|float64](x: T): T = return ln(x) / ln(10.0)
-  proc log2*[T: float32|float64](x: T): T = return ln(x) / ln(2.0)
-
-  proc exp*(x: float32): float32 {.importc: "Math.exp", nodecl.}
-  proc exp*(x: float64): float64 {.importc: "Math.exp", nodecl.}
   proc round0(x: float): float {.importc: "Math.round", nodecl.}
+  proc trunc*(x: float32): float32 {.importc: "Math.trunc", nodecl.}
+  proc trunc*(x: float64): float64 {.importc: "Math.trunc", nodecl.}
 
-  proc pow*(x, y: float32): float32 {.importC: "Math.pow", nodecl.}
-  proc pow*(x, y: float64): float64 {.importc: "Math.pow", nodecl.}
-
-  proc arccos*(x: float32): float32 {.importc: "Math.acos", nodecl.}
-  proc arccos*(x: float64): float64 {.importc: "Math.acos", nodecl.}
-  proc arcsin*(x: float32): float32 {.importc: "Math.asin", nodecl.}
-  proc arcsin*(x: float64): float64 {.importc: "Math.asin", nodecl.}
-  proc arctan*(x: float32): float32 {.importc: "Math.atan", nodecl.}
-  proc arctan*(x: float64): float64 {.importc: "Math.atan", nodecl.}
-  proc arctan2*(y, x: float32): float32 {.importC: "Math.atan2", nodecl.}
-  proc arctan2*(y, x: float64): float64 {.importc: "Math.atan2", nodecl.}
-
-  proc cos*(x: float32): float32 {.importc: "Math.cos", nodecl.}
-  proc cos*(x: float64): float64 {.importc: "Math.cos", nodecl.}
-  proc cosh*(x: float32): float32 = return (exp(x)+exp(-x))*0.5
-  proc cosh*(x: float64): float64 = return (exp(x)+exp(-x))*0.5
-  proc hypot*[T: float32|float64](x, y: T): T = return sqrt(x*x + y*y)
-  proc sinh*[T: float32|float64](x: T): T = return (exp(x)-exp(-x))*0.5
-  proc sin*(x: float32): float32 {.importc: "Math.sin", nodecl.}
-  proc sin*(x: float64): float64 {.importc: "Math.sin", nodecl.}
-  proc tan*(x: float32): float32 {.importc: "Math.tan", nodecl.}
-  proc tan*(x: float64): float64 {.importc: "Math.tan", nodecl.}
-  proc tanh*[T: float32|float64](x: T): T =
-    var y = exp(2.0*x)
-    return (y-1.0)/(y+1.0)
+  proc `mod`*(x, y: float32): float32 {.importcpp: "# % #".}
+  proc `mod`*(x, y: float64): float64 {.importcpp: "# % #".}
+  ## Computes the modulo operation for float operators.
 
 proc round*[T: float32|float64](x: T, places: int = 0): T =
   ## Round a floating point number.
@@ -350,6 +416,21 @@ proc round*[T: float32|float64](x: T, places: int = 0): T =
     var mult = pow(10.0, places.T)
     result = round0(x*mult)/mult
 
+proc floorDiv*[T: SomeInteger](x, y: T): T =
+  ## Floor division is conceptually defined as ``floor(x / y)``.
+  ## This is different from the ``div`` operator, which is defined
+  ## as ``trunc(x / y)``. That is, ``div`` rounds towards ``0`` and ``floorDiv``
+  ## rounds down.
+  result = x div y
+  let r = x mod y
+  if (r > 0 and y < 0) or (r < 0 and y > 0): result.dec 1
+
+proc floorMod*[T: SomeNumber](x, y: T): T =
+  ## Floor modulus is conceptually defined as ``x - (floorDiv(x, y) * y).
+  ## This proc behaves the same as the ``%`` operator in python.
+  result = x mod y
+  if (result > 0 and y < 0) or (result < 0 and y > 0): result += y
+
 when not defined(JS):
   proc c_frexp*(x: float32, exponent: var int32): float32 {.
     importc: "frexp", header: "<math.h>".}
@@ -414,15 +495,6 @@ proc sgn*[T: SomeNumber](x: T): int {.inline.} =
   ## `NaN`.
   ord(T(0) < x) - ord(x < T(0))
 
-proc `mod`*[T: float32|float64](x, y: T): T =
-  ## Computes the modulo operation for float operators. Equivalent
-  ## to ``x - y * floor(x/y)``. Note that the remainder will always
-  ## have the same sign as the divisor.
-  ##
-  ## .. code-block:: nim
-  ##  echo (4.0 mod -3.1) # -2.2
-  result = if y == 0.0: x else: x - y * (x/y).floor
-
 {.pop.}
 {.pop.}
 
@@ -445,16 +517,42 @@ proc `^`*[T](x: T, y: Natural): T =
     x *= x
 
 proc gcd*[T](x, y: T): T =
-  ## Computes the greatest common divisor of ``x`` and ``y``.
+  ## Computes the greatest common (positive) divisor of ``x`` and ``y``.
   ## Note that for floats, the result cannot always be interpreted as
   ## "greatest decimal `z` such that ``z*N == x and z*M == y``
   ## where N and M are positive integers."
-  var (x,y) = (x,y)
+  var (x, y) = (x, y)
   while y != 0:
     x = x mod y
     swap x, y
   abs x
 
+proc gcd*(x, y: SomeInteger): SomeInteger =
+  ## Computes the greatest common (positive) divisor of ``x`` and ``y``.
+  ## Using binary GCD (aka Stein's) algorithm.
+  when x is SomeSignedInt:
+    var x = abs(x)
+  else:
+    var x = x
+  when y is SomeSignedInt:
+    var y = abs(y)
+  else:
+    var y = y
+
+  if x == 0:
+    return y
+  if y == 0:
+    return x
+
+  let shift = countTrailingZeroBits(x or y)
+  y = y shr countTrailingZeroBits(y)
+  while x != 0:
+    x = x shr countTrailingZeroBits(x)
+    if y > x:
+      swap y, x
+    x -= y
+  y shl shift
+
 proc lcm*[T](x, y: T): T =
   ## Computes the least common multiple of ``x`` and ``y``.
   x div gcd(x, y) * y
@@ -465,6 +563,7 @@ when isMainModule and not defined(JS):
     return sqrt(num)
 
   # check gamma function
+  assert(gamma(5.0) == 24.0) # 4!
   assert($tgamma(5.0) == $24.0) # 4!
   assert(lgamma(1.0) == 0.0) # ln(1.0) == 0.0
   assert(erf(6.0) > erf(5.0))
@@ -474,6 +573,12 @@ when isMainModule:
   # Function for approximate comparison of floats
   proc `==~`(x, y: float): bool = (abs(x-y) < 1e-9)
 
+  block: # prod
+    doAssert prod([1, 2, 3, 4]) == 24
+    doAssert prod([1.5, 3.4]) == 5.1
+    let x: seq[float] = @[]
+    doAssert prod(x) == 1.0
+
   block: # round() tests
     # Round to 0 decimal places
     doAssert round(54.652) ==~ 55.0
@@ -550,3 +655,30 @@ when isMainModule:
     assert sgn(Inf) == 1
     assert sgn(NaN) == 0
 
+  block: # fac() tests
+    try:
+      discard fac(-1)
+    except AssertionError:
+      discard
+
+    doAssert fac(0) == 1
+    doAssert fac(1) == 1
+    doAssert fac(2) == 2
+    doAssert fac(3) == 6
+    doAssert fac(4) == 24
+
+  block: # floorMod/floorDiv
+    doAssert floorDiv(8, 3) == 2
+    doAssert floorMod(8, 3) == 2
+
+    doAssert floorDiv(8, -3) == -3
+    doAssert floorMod(8, -3) == -1
+
+    doAssert floorDiv(-8, 3) == -3
+    doAssert floorMod(-8, 3) == 1
+
+    doAssert floorDiv(-8, -3) == 2
+    doAssert floorMod(-8, -3) == -2
+
+    doAssert floorMod(8.0, -3.0) ==~ -1.0
+    doAssert floorMod(-8.5, 3.0) ==~ 0.5