diff options
Diffstat (limited to 'lib/pure/math.nim')
| -rw-r--r-- | lib/pure/math.nim | 292 |
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 |