| Line | Exclusive | Inclusive | Code |
|---|---|---|---|
| 1 | # This file is a part of Julia. License is MIT: https://julialang.org/license | ||
| 2 | |||
| 3 | ## number-theoretic functions ## | ||
| 4 | |||
| 5 | """ | ||
| 6 | gcd(x,y) | ||
| 7 | |||
| 8 | Greatest common (positive) divisor (or zero if `x` and `y` are both zero). | ||
| 9 | The arguments may be integer and rational numbers. | ||
| 10 | |||
| 11 | !!! compat "Julia 1.4" | ||
| 12 | Rational arguments require Julia 1.4 or later. | ||
| 13 | |||
| 14 | # Examples | ||
| 15 | ```jldoctest | ||
| 16 | julia> gcd(6,9) | ||
| 17 | 3 | ||
| 18 | |||
| 19 | julia> gcd(6,-9) | ||
| 20 | 3 | ||
| 21 | ``` | ||
| 22 | """ | ||
| 23 | function gcd(a::T, b::T) where T<:Integer | ||
| 24 | while b != 0 | ||
| 25 | t = b | ||
| 26 | b = rem(a, b) | ||
| 27 | a = t | ||
| 28 | end | ||
| 29 | checked_abs(a) | ||
| 30 | end | ||
| 31 | |||
| 32 | # binary GCD (aka Stein's) algorithm | ||
| 33 | # about 1.7x (2.1x) faster for random Int64s (Int128s) | ||
| 34 | function gcd(a::T, b::T) where T<:Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128} | ||
| 35 | a == 0 && return abs(b) | ||
| 36 | b == 0 && return abs(a) | ||
| 37 | za = trailing_zeros(a) | ||
| 38 | zb = trailing_zeros(b) | ||
| 39 | k = min(za, zb) | ||
| 40 | u = unsigned(abs(a >> za)) | ||
| 41 | v = unsigned(abs(b >> zb)) | ||
| 42 | while u != v | ||
| 43 | if u > v | ||
| 44 | u, v = v, u | ||
| 45 | end | ||
| 46 | v -= u | ||
| 47 | v >>= trailing_zeros(v) | ||
| 48 | end | ||
| 49 | r = u << k | ||
| 50 | # T(r) would throw InexactError; we want OverflowError instead | ||
| 51 | r > typemax(T) && __throw_gcd_overflow(a, b) | ||
| 52 | r % T | ||
| 53 | end | ||
| 54 | @noinline __throw_gcd_overflow(a, b) = throw(OverflowError("gcd($a, $b) overflows")) | ||
| 55 | |||
| 56 | """ | ||
| 57 | lcm(x,y) | ||
| 58 | |||
| 59 | Least common (non-negative) multiple. | ||
| 60 | The arguments may be integer and rational numbers. | ||
| 61 | |||
| 62 | !!! compat "Julia 1.4" | ||
| 63 | Rational arguments require Julia 1.4 or later. | ||
| 64 | |||
| 65 | # Examples | ||
| 66 | ```jldoctest | ||
| 67 | julia> lcm(2,3) | ||
| 68 | 6 | ||
| 69 | |||
| 70 | julia> lcm(-2,3) | ||
| 71 | 6 | ||
| 72 | ``` | ||
| 73 | """ | ||
| 74 | function lcm(a::T, b::T) where T<:Integer | ||
| 75 | # explicit a==0 test is to handle case of lcm(0,0) correctly | ||
| 76 | if a == 0 | ||
| 77 | return a | ||
| 78 | else | ||
| 79 | return checked_abs(checked_mul(a, div(b, gcd(b,a)))) | ||
| 80 | end | ||
| 81 | end | ||
| 82 | |||
| 83 | gcd(a::Union{Integer,Rational}) = a | ||
| 84 | lcm(a::Union{Integer,Rational}) = a | ||
| 85 | gcd(a::Real, b::Real) = gcd(promote(a,b)...) | ||
| 86 | lcm(a::Real, b::Real) = lcm(promote(a,b)...) | ||
| 87 | gcd(a::Real, b::Real, c::Real...) = gcd(a, gcd(b, c...)) | ||
| 88 | lcm(a::Real, b::Real, c::Real...) = lcm(a, lcm(b, c...)) | ||
| 89 | gcd(a::T, b::T) where T<:Real = throw(MethodError(gcd, (a,b))) | ||
| 90 | lcm(a::T, b::T) where T<:Real = throw(MethodError(lcm, (a,b))) | ||
| 91 | |||
| 92 | gcd(abc::AbstractArray{<:Real}) = reduce(gcd, abc; init=zero(eltype(abc))) | ||
| 93 | lcm(abc::AbstractArray{<:Real}) = reduce(lcm, abc; init=one(eltype(abc))) | ||
| 94 | |||
| 95 | function gcd(abc::AbstractArray{<:Integer}) | ||
| 96 | a = zero(eltype(abc)) | ||
| 97 | for b in abc | ||
| 98 | a = gcd(a,b) | ||
| 99 | if a == 1 | ||
| 100 | return a | ||
| 101 | end | ||
| 102 | end | ||
| 103 | return a | ||
| 104 | end | ||
| 105 | |||
| 106 | # return (gcd(a,b),x,y) such that ax+by == gcd(a,b) | ||
| 107 | """ | ||
| 108 | gcdx(x,y) | ||
| 109 | |||
| 110 | Computes the greatest common (positive) divisor of `x` and `y` and their Bézout | ||
| 111 | coefficients, i.e. the integer coefficients `u` and `v` that satisfy | ||
| 112 | ``ux+vy = d = gcd(x,y)``. ``gcdx(x,y)`` returns ``(d,u,v)``. | ||
| 113 | |||
| 114 | The arguments may be integer and rational numbers. | ||
| 115 | |||
| 116 | !!! compat "Julia 1.4" | ||
| 117 | Rational arguments require Julia 1.4 or later. | ||
| 118 | |||
| 119 | # Examples | ||
| 120 | ```jldoctest | ||
| 121 | julia> gcdx(12, 42) | ||
| 122 | (6, -3, 1) | ||
| 123 | |||
| 124 | julia> gcdx(240, 46) | ||
| 125 | (2, -9, 47) | ||
| 126 | ``` | ||
| 127 | |||
| 128 | !!! note | ||
| 129 | Bézout coefficients are *not* uniquely defined. `gcdx` returns the minimal | ||
| 130 | Bézout coefficients that are computed by the extended Euclidean algorithm. | ||
| 131 | (Ref: D. Knuth, TAoCP, 2/e, p. 325, Algorithm X.) | ||
| 132 | For signed integers, these coefficients `u` and `v` are minimal in | ||
| 133 | the sense that ``|u| < |y/d|`` and ``|v| < |x/d|``. Furthermore, | ||
| 134 | the signs of `u` and `v` are chosen so that `d` is positive. | ||
| 135 | For unsigned integers, the coefficients `u` and `v` might be near | ||
| 136 | their `typemax`, and the identity then holds only via the unsigned | ||
| 137 | integers' modulo arithmetic. | ||
| 138 | """ | ||
| 139 | function gcdx(a::T, b::T) where T<:Integer | ||
| 140 | # a0, b0 = a, b | ||
| 141 | s0, s1 = oneunit(T), zero(T) | ||
| 142 | t0, t1 = s1, s0 | ||
| 143 | # The loop invariant is: s0*a0 + t0*b0 == a | ||
| 144 | while b != 0 | ||
| 145 | q = div(a, b) | ||
| 146 | a, b = b, rem(a, b) | ||
| 147 | s0, s1 = s1, s0 - q*s1 | ||
| 148 | t0, t1 = t1, t0 - q*t1 | ||
| 149 | end | ||
| 150 | a < 0 ? (-a, -s0, -t0) : (a, s0, t0) | ||
| 151 | end | ||
| 152 | gcdx(a::Real, b::Real) = gcdx(promote(a,b)...) | ||
| 153 | gcdx(a::T, b::T) where T<:Real = throw(MethodError(gcdx, (a,b))) | ||
| 154 | |||
| 155 | # multiplicative inverse of n mod m, error if none | ||
| 156 | |||
| 157 | """ | ||
| 158 | invmod(x,m) | ||
| 159 | |||
| 160 | Take the inverse of `x` modulo `m`: `y` such that ``x y = 1 \\pmod m``, | ||
| 161 | with ``div(x,y) = 0``. This is undefined for ``m = 0``, or if | ||
| 162 | ``gcd(x,m) \\neq 1``. | ||
| 163 | |||
| 164 | # Examples | ||
| 165 | ```jldoctest | ||
| 166 | julia> invmod(2,5) | ||
| 167 | 3 | ||
| 168 | |||
| 169 | julia> invmod(2,3) | ||
| 170 | 2 | ||
| 171 | |||
| 172 | julia> invmod(5,6) | ||
| 173 | 5 | ||
| 174 | ``` | ||
| 175 | """ | ||
| 176 | function invmod(n::T, m::T) where T<:Integer | ||
| 177 | g, x, y = gcdx(n, m) | ||
| 178 | g != 1 && throw(DomainError((n, m), "Greatest common divisor is $g.")) | ||
| 179 | m == 0 && throw(DomainError(m, "`m` must not be 0.")) | ||
| 180 | # Note that m might be negative here. | ||
| 181 | # For unsigned T, x might be close to typemax; add m to force a wrap-around. | ||
| 182 | r = mod(x + m, m) | ||
| 183 | # The postcondition is: mod(r * n, m) == mod(T(1), m) && div(r, m) == 0 | ||
| 184 | r | ||
| 185 | end | ||
| 186 | invmod(n::Integer, m::Integer) = invmod(promote(n,m)...) | ||
| 187 | |||
| 188 | # ^ for any x supporting * | ||
| 189 | to_power_type(x) = convert(Base._return_type(*, Tuple{typeof(x), typeof(x)}), x) | ||
| 190 | @noinline throw_domerr_powbysq(::Any, p) = throw(DomainError(p, | ||
| 191 | string("Cannot raise an integer x to a negative power ", p, '.', | ||
| 192 | "\nConvert input to float."))) | ||
| 193 | @noinline throw_domerr_powbysq(::Integer, p) = throw(DomainError(p, | ||
| 194 | string("Cannot raise an integer x to a negative power ", p, '.', | ||
| 195 | "\nMake x or $p a float by adding a zero decimal ", | ||
| 196 | "(e.g., 2.0^$p or 2^$(float(p)) instead of 2^$p), ", | ||
| 197 | "or write 1/x^$(-p), float(x)^$p, x^float($p) or (x//1)^$p"))) | ||
| 198 | @noinline throw_domerr_powbysq(::AbstractMatrix, p) = throw(DomainError(p, | ||
| 199 | string("Cannot raise an integer matrix x to a negative power ", p, '.', | ||
| 200 | "\nMake x a float matrix by adding a zero decimal ", | ||
| 201 | "(e.g., [2.0 1.0;1.0 0.0]^$p instead ", | ||
| 202 | "of [2 1;1 0]^$p), or write float(x)^$p or Rational.(x)^$p"))) | ||
| 203 | function power_by_squaring(x_, p::Integer) | ||
| 204 | x = to_power_type(x_) | ||
| 205 | if p == 1 | ||
| 206 | return copy(x) | ||
| 207 | elseif p == 0 | ||
| 208 | return one(x) | ||
| 209 | elseif p == 2 | ||
| 210 | return x*x | ||
| 211 | elseif p < 0 | ||
| 212 | isone(x) && return copy(x) | ||
| 213 | isone(-x) && return iseven(p) ? one(x) : copy(x) | ||
| 214 | throw_domerr_powbysq(x, p) | ||
| 215 | end | ||
| 216 | t = trailing_zeros(p) + 1 | ||
| 217 | p >>= t | ||
| 218 | while (t -= 1) > 0 | ||
| 219 | x *= x | ||
| 220 | end | ||
| 221 | y = x | ||
| 222 | while p > 0 | ||
| 223 | t = trailing_zeros(p) + 1 | ||
| 224 | p >>= t | ||
| 225 | while (t -= 1) >= 0 | ||
| 226 | 50 (6 %) |
50 (100 %)
samples spent calling
*
x *= x
|
|
| 227 | end | ||
| 228 | 783 (94 %) |
783 (100 %)
samples spent calling
*
y *= x
|
|
| 229 | end | ||
| 230 | return y | ||
| 231 | end | ||
| 232 | power_by_squaring(x::Bool, p::Unsigned) = ((p==0) | x) | ||
| 233 | function power_by_squaring(x::Bool, p::Integer) | ||
| 234 | p < 0 && !x && throw_domerr_powbysq(x, p) | ||
| 235 | return (p==0) | x | ||
| 236 | end | ||
| 237 | |||
| 238 | ^(x::T, p::T) where {T<:Integer} = power_by_squaring(x,p) | ||
| 239 | ^(x::Number, p::Integer) = power_by_squaring(x,p) | ||
| 240 | |||
| 241 | # x^p for any literal integer p is lowered to Base.literal_pow(^, x, Val(p)) | ||
| 242 | # to enable compile-time optimizations specialized to p. | ||
| 243 | # However, we still need a fallback that calls the function ^ which may either | ||
| 244 | # mean Base.^ or something else, depending on context. | ||
| 245 | # We mark these @inline since if the target is marked @inline, | ||
| 246 | # we want to make sure that gets propagated, | ||
| 247 | # even if it is over the inlining threshold. | ||
| 248 | @inline literal_pow(f, x, ::Val{p}) where {p} = f(x,p) | ||
| 249 | |||
| 250 | # Restrict inlining to hardware-supported arithmetic types, which | ||
| 251 | # are fast enough to benefit from inlining. | ||
| 252 | const HWReal = Union{Int8,Int16,Int32,Int64,UInt8,UInt16,UInt32,UInt64,Float32,Float64} | ||
| 253 | const HWNumber = Union{HWReal, Complex{<:HWReal}, Rational{<:HWReal}} | ||
| 254 | |||
| 255 | # Core.Compiler has complicated logic to inline x^2 and x^3 for | ||
| 256 | # numeric types. In terms of Val we can do it much more simply. | ||
| 257 | # (The first argument prevents unexpected behavior if a function ^ | ||
| 258 | # is defined that is not equal to Base.^) | ||
| 259 | @inline literal_pow(::typeof(^), x::HWNumber, ::Val{0}) = one(x) | ||
| 260 | @inline literal_pow(::typeof(^), x::HWNumber, ::Val{1}) = x | ||
| 261 | @inline literal_pow(::typeof(^), x::HWNumber, ::Val{2}) = x*x | ||
| 262 | @inline literal_pow(::typeof(^), x::HWNumber, ::Val{3}) = x*x*x | ||
| 263 | |||
| 264 | # don't use the inv(x) transformation here since float^p is slightly more accurate | ||
| 265 | @inline literal_pow(::typeof(^), x::AbstractFloat, ::Val{p}) where {p} = x^p | ||
| 266 | @inline literal_pow(::typeof(^), x::AbstractFloat, ::Val{-1}) = inv(x) | ||
| 267 | |||
| 268 | # for other types, define x^-n as inv(x)^n so that negative literal powers can | ||
| 269 | # be computed in a type-stable way even for e.g. integers. | ||
| 270 | @inline @generated function literal_pow(f::typeof(^), x, ::Val{p}) where {p} | ||
| 271 | if p < 0 | ||
| 272 | :(literal_pow(^, inv(x), $(Val{-p}()))) | ||
| 273 | else | ||
| 274 | :(f(x,$p)) | ||
| 275 | end | ||
| 276 | end | ||
| 277 | |||
| 278 | # note: it is tempting to add optimized literal_pow(::typeof(^), x, ::Val{n}) | ||
| 279 | # methods here for various n, but this easily leads to method ambiguities | ||
| 280 | # if anyone has defined literal_pow(::typeof(^), x::T, ::Val). | ||
| 281 | |||
| 282 | # b^p mod m | ||
| 283 | |||
| 284 | """ | ||
| 285 | powermod(x::Integer, p::Integer, m) | ||
| 286 | |||
| 287 | Compute ``x^p \\pmod m``. | ||
| 288 | |||
| 289 | # Examples | ||
| 290 | ```jldoctest | ||
| 291 | julia> powermod(2, 6, 5) | ||
| 292 | 4 | ||
| 293 | |||
| 294 | julia> mod(2^6, 5) | ||
| 295 | 4 | ||
| 296 | |||
| 297 | julia> powermod(5, 2, 20) | ||
| 298 | 5 | ||
| 299 | |||
| 300 | julia> powermod(5, 2, 19) | ||
| 301 | 6 | ||
| 302 | |||
| 303 | julia> powermod(5, 3, 19) | ||
| 304 | 11 | ||
| 305 | ``` | ||
| 306 | """ | ||
| 307 | function powermod(x::Integer, p::Integer, m::T) where T<:Integer | ||
| 308 | p < 0 && return powermod(invmod(x, m), -p, m) | ||
| 309 | p == 0 && return mod(one(m),m) | ||
| 310 | (m == 1 || m == -1) && return zero(m) | ||
| 311 | b = oftype(m,mod(x,m)) # this also checks for divide by zero | ||
| 312 | |||
| 313 | t = prevpow(2, p) | ||
| 314 | r::T = 1 | ||
| 315 | while true | ||
| 316 | if p >= t | ||
| 317 | r = mod(widemul(r,b),m) | ||
| 318 | p -= t | ||
| 319 | end | ||
| 320 | t >>>= 1 | ||
| 321 | t <= 0 && break | ||
| 322 | r = mod(widemul(r,r),m) | ||
| 323 | end | ||
| 324 | return r | ||
| 325 | end | ||
| 326 | |||
| 327 | # optimization: promote the modulus m to BigInt only once (cf. widemul in generic powermod above) | ||
| 328 | powermod(x::Integer, p::Integer, m::Union{Int128,UInt128}) = oftype(m, powermod(x, p, big(m))) | ||
| 329 | |||
| 330 | _nextpow2(x::Unsigned) = oneunit(x)<<((sizeof(x)<<3)-leading_zeros(x-oneunit(x))) | ||
| 331 | _nextpow2(x::Integer) = reinterpret(typeof(x),x < 0 ? -_nextpow2(unsigned(-x)) : _nextpow2(unsigned(x))) | ||
| 332 | _prevpow2(x::Unsigned) = one(x) << unsigned((sizeof(x)<<3)-leading_zeros(x)-1) | ||
| 333 | _prevpow2(x::Integer) = reinterpret(typeof(x),x < 0 ? -_prevpow2(unsigned(-x)) : _prevpow2(unsigned(x))) | ||
| 334 | |||
| 335 | """ | ||
| 336 | ispow2(n::Integer) -> Bool | ||
| 337 | |||
| 338 | Test whether `n` is a power of two. | ||
| 339 | |||
| 340 | # Examples | ||
| 341 | ```jldoctest | ||
| 342 | julia> ispow2(4) | ||
| 343 | true | ||
| 344 | |||
| 345 | julia> ispow2(5) | ||
| 346 | false | ||
| 347 | ``` | ||
| 348 | """ | ||
| 349 | ispow2(x::Integer) = x > 0 && count_ones(x) == 1 | ||
| 350 | |||
| 351 | """ | ||
| 352 | nextpow(a, x) | ||
| 353 | |||
| 354 | The smallest `a^n` not less than `x`, where `n` is a non-negative integer. `a` must be | ||
| 355 | greater than 1, and `x` must be greater than 0. | ||
| 356 | |||
| 357 | # Examples | ||
| 358 | ```jldoctest | ||
| 359 | julia> nextpow(2, 7) | ||
| 360 | 8 | ||
| 361 | |||
| 362 | julia> nextpow(2, 9) | ||
| 363 | 16 | ||
| 364 | |||
| 365 | julia> nextpow(5, 20) | ||
| 366 | 25 | ||
| 367 | |||
| 368 | julia> nextpow(4, 16) | ||
| 369 | 16 | ||
| 370 | ``` | ||
| 371 | |||
| 372 | See also [`prevpow`](@ref). | ||
| 373 | """ | ||
| 374 | function nextpow(a::Real, x::Real) | ||
| 375 | x <= 0 && throw(DomainError(x, "`x` must be positive.")) | ||
| 376 | # Special case fast path for x::Integer, a == 2. | ||
| 377 | # This is a very common case. Constant prop will make sure that a call site | ||
| 378 | # specified as `nextpow(2, x)` will get this special case inlined. | ||
| 379 | a == 2 && isa(x, Integer) && return _nextpow2(x) | ||
| 380 | a <= 1 && throw(DomainError(a, "`a` must be greater than 1.")) | ||
| 381 | x <= 1 && return one(a) | ||
| 382 | n = ceil(Integer,log(a, x)) | ||
| 383 | p = a^(n-1) | ||
| 384 | # guard against roundoff error, e.g., with a=5 and x=125 | ||
| 385 | p >= x ? p : a^n | ||
| 386 | end | ||
| 387 | |||
| 388 | """ | ||
| 389 | prevpow(a, x) | ||
| 390 | |||
| 391 | The largest `a^n` not greater than `x`, where `n` is a non-negative integer. | ||
| 392 | `a` must be greater than 1, and `x` must not be less than 1. | ||
| 393 | |||
| 394 | # Examples | ||
| 395 | ```jldoctest | ||
| 396 | julia> prevpow(2, 7) | ||
| 397 | 4 | ||
| 398 | |||
| 399 | julia> prevpow(2, 9) | ||
| 400 | 8 | ||
| 401 | |||
| 402 | julia> prevpow(5, 20) | ||
| 403 | 5 | ||
| 404 | |||
| 405 | julia> prevpow(4, 16) | ||
| 406 | 16 | ||
| 407 | ``` | ||
| 408 | See also [`nextpow`](@ref). | ||
| 409 | """ | ||
| 410 | function prevpow(a::Real, x::Real) | ||
| 411 | x < 1 && throw(DomainError(x, "`x` must be ≥ 1.")) | ||
| 412 | # See comment in nextpos() for a == special case. | ||
| 413 | a == 2 && isa(x, Integer) && return _prevpow2(x) | ||
| 414 | a <= 1 && throw(DomainError(a, "`a` must be greater than 1.")) | ||
| 415 | n = floor(Integer,log(a, x)) | ||
| 416 | p = a^(n+1) | ||
| 417 | p <= x ? p : a^n | ||
| 418 | end | ||
| 419 | |||
| 420 | ## ndigits (number of digits) in base 10 ## | ||
| 421 | |||
| 422 | # decimal digits in an unsigned integer | ||
| 423 | const powers_of_ten = [ | ||
| 424 | 0x0000000000000001, 0x000000000000000a, 0x0000000000000064, 0x00000000000003e8, | ||
| 425 | 0x0000000000002710, 0x00000000000186a0, 0x00000000000f4240, 0x0000000000989680, | ||
| 426 | 0x0000000005f5e100, 0x000000003b9aca00, 0x00000002540be400, 0x000000174876e800, | ||
| 427 | 0x000000e8d4a51000, 0x000009184e72a000, 0x00005af3107a4000, 0x00038d7ea4c68000, | ||
| 428 | 0x002386f26fc10000, 0x016345785d8a0000, 0x0de0b6b3a7640000, 0x8ac7230489e80000, | ||
| 429 | ] | ||
| 430 | function bit_ndigits0z(x::Base.BitUnsigned64) | ||
| 431 | lz = (sizeof(x)<<3)-leading_zeros(x) | ||
| 432 | nd = (1233*lz)>>12+1 | ||
| 433 | nd -= x < powers_of_ten[nd] | ||
| 434 | end | ||
| 435 | function bit_ndigits0z(x::UInt128) | ||
| 436 | n = 0 | ||
| 437 | while x > 0x8ac7230489e80000 | ||
| 438 | x = div(x,0x8ac7230489e80000) | ||
| 439 | n += 19 | ||
| 440 | end | ||
| 441 | return n + ndigits0z(UInt64(x)) | ||
| 442 | end | ||
| 443 | |||
| 444 | ndigits0z(x::BitSigned) = bit_ndigits0z(unsigned(abs(x))) | ||
| 445 | ndigits0z(x::BitUnsigned) = bit_ndigits0z(x) | ||
| 446 | ndigits0z(x::Integer) = ndigits0zpb(x, 10) | ||
| 447 | |||
| 448 | ## ndigits with specified base ## | ||
| 449 | |||
| 450 | # The suffix "nb" stands for "negative base" | ||
| 451 | function ndigits0znb(x::Integer, b::Integer) | ||
| 452 | d = 0 | ||
| 453 | if x isa Unsigned | ||
| 454 | d += (x != 0)::Bool | ||
| 455 | x = -signed(fld(x, -b)) | ||
| 456 | end | ||
| 457 | # precondition: b < -1 && !(typeof(x) <: Unsigned) | ||
| 458 | while x != 0 | ||
| 459 | x = cld(x,b) | ||
| 460 | d += 1 | ||
| 461 | end | ||
| 462 | return d | ||
| 463 | end | ||
| 464 | |||
| 465 | # do first division before conversion with signed here, which can otherwise overflow | ||
| 466 | ndigits0znb(x::Bool, b::Integer) = x % Int | ||
| 467 | |||
| 468 | # The suffix "pb" stands for "positive base" | ||
| 469 | function ndigits0zpb(x::Integer, b::Integer) | ||
| 470 | # precondition: b > 1 | ||
| 471 | x == 0 && return 0 | ||
| 472 | b = Int(b) | ||
| 473 | x = abs(x) | ||
| 474 | if x isa Base.BitInteger | ||
| 475 | x = unsigned(x)::Unsigned | ||
| 476 | b == 2 && return sizeof(x)<<3 - leading_zeros(x) | ||
| 477 | b == 8 && return (sizeof(x)<<3 - leading_zeros(x) + 2) ÷ 3 | ||
| 478 | b == 16 && return sizeof(x)<<1 - leading_zeros(x)>>2 | ||
| 479 | b == 10 && return bit_ndigits0z(x) | ||
| 480 | if ispow2(b) | ||
| 481 | dv, rm = divrem(sizeof(x)<<3 - leading_zeros(x), trailing_zeros(b)) | ||
| 482 | return iszero(rm) ? dv : dv + 1 | ||
| 483 | end | ||
| 484 | end | ||
| 485 | |||
| 486 | d = 0 | ||
| 487 | while x > typemax(Int) | ||
| 488 | x = div(x,b) | ||
| 489 | d += 1 | ||
| 490 | end | ||
| 491 | x = div(x,b) | ||
| 492 | d += 1 | ||
| 493 | |||
| 494 | m = 1 | ||
| 495 | while m <= x | ||
| 496 | m *= b | ||
| 497 | d += 1 | ||
| 498 | end | ||
| 499 | return d | ||
| 500 | end | ||
| 501 | |||
| 502 | ndigits0zpb(x::Bool, b::Integer) = x % Int | ||
| 503 | |||
| 504 | # The suffix "0z" means that the output is 0 on input zero (cf. #16841) | ||
| 505 | """ | ||
| 506 | ndigits0z(n::Integer, b::Integer=10) | ||
| 507 | |||
| 508 | Return 0 if `n == 0`, otherwise compute the number of digits in | ||
| 509 | integer `n` written in base `b` (i.e. equal to `ndigits(n, base=b)` | ||
| 510 | in this case). | ||
| 511 | The base `b` must not be in `[-1, 0, 1]`. | ||
| 512 | |||
| 513 | # Examples | ||
| 514 | ```jldoctest | ||
| 515 | julia> Base.ndigits0z(0, 16) | ||
| 516 | 0 | ||
| 517 | |||
| 518 | julia> Base.ndigits(0, base=16) | ||
| 519 | 1 | ||
| 520 | |||
| 521 | julia> Base.ndigits0z(0) | ||
| 522 | 0 | ||
| 523 | |||
| 524 | julia> Base.ndigits0z(10, 2) | ||
| 525 | 4 | ||
| 526 | |||
| 527 | julia> Base.ndigits0z(10) | ||
| 528 | 2 | ||
| 529 | ``` | ||
| 530 | |||
| 531 | See also [`ndigits`](@ref). | ||
| 532 | """ | ||
| 533 | function ndigits0z(x::Integer, b::Integer) | ||
| 534 | if b < -1 | ||
| 535 | ndigits0znb(x, b) | ||
| 536 | elseif b > 1 | ||
| 537 | ndigits0zpb(x, b) | ||
| 538 | else | ||
| 539 | throw(DomainError(b, "The base must not be in `[-1, 0, 1]`.")) | ||
| 540 | end | ||
| 541 | end | ||
| 542 | |||
| 543 | """ | ||
| 544 | ndigits(n::Integer; base::Integer=10, pad::Integer=1) | ||
| 545 | |||
| 546 | Compute the number of digits in integer `n` written in base `base` | ||
| 547 | (`base` must not be in `[-1, 0, 1]`), optionally padded with zeros | ||
| 548 | to a specified size (the result will never be less than `pad`). | ||
| 549 | |||
| 550 | # Examples | ||
| 551 | ```jldoctest | ||
| 552 | julia> ndigits(12345) | ||
| 553 | 5 | ||
| 554 | |||
| 555 | julia> ndigits(1022, base=16) | ||
| 556 | 3 | ||
| 557 | |||
| 558 | julia> string(1022, base=16) | ||
| 559 | "3fe" | ||
| 560 | |||
| 561 | julia> ndigits(123, pad=5) | ||
| 562 | 5 | ||
| 563 | ``` | ||
| 564 | """ | ||
| 565 | ndigits(x::Integer; base::Integer=10, pad::Integer=1) = max(pad, ndigits0z(x, base)) | ||
| 566 | |||
| 567 | ## integer to string functions ## | ||
| 568 | |||
| 569 | function bin(x::Unsigned, pad::Integer, neg::Bool) | ||
| 570 | i = neg + max(pad,sizeof(x)<<3-leading_zeros(x)) | ||
| 571 | a = StringVector(i) | ||
| 572 | while i > neg | ||
| 573 | @inbounds a[i] = 48+(x&0x1) | ||
| 574 | x >>= 1 | ||
| 575 | i -= 1 | ||
| 576 | end | ||
| 577 | if neg; @inbounds a[1]=0x2d; end | ||
| 578 | String(a) | ||
| 579 | end | ||
| 580 | |||
| 581 | function oct(x::Unsigned, pad::Integer, neg::Bool) | ||
| 582 | i = neg + max(pad,div((sizeof(x)<<3)-leading_zeros(x)+2,3)) | ||
| 583 | a = StringVector(i) | ||
| 584 | while i > neg | ||
| 585 | @inbounds a[i] = 48+(x&0x7) | ||
| 586 | x >>= 3 | ||
| 587 | i -= 1 | ||
| 588 | end | ||
| 589 | if neg; @inbounds a[1]=0x2d; end | ||
| 590 | String(a) | ||
| 591 | end | ||
| 592 | |||
| 593 | function dec(x::Unsigned, pad::Integer, neg::Bool) | ||
| 594 | i = neg + ndigits(x, base=10, pad=pad) | ||
| 595 | a = StringVector(i) | ||
| 596 | while i > neg | ||
| 597 | @inbounds a[i] = 48+rem(x,10) | ||
| 598 | x = oftype(x,div(x,10)) | ||
| 599 | i -= 1 | ||
| 600 | end | ||
| 601 | if neg; @inbounds a[1]=0x2d; end | ||
| 602 | String(a) | ||
| 603 | end | ||
| 604 | |||
| 605 | function hex(x::Unsigned, pad::Integer, neg::Bool) | ||
| 606 | i = neg + max(pad,(sizeof(x)<<1)-(leading_zeros(x)>>2)) | ||
| 607 | a = StringVector(i) | ||
| 608 | while i > neg | ||
| 609 | d = x & 0xf | ||
| 610 | @inbounds a[i] = 48+d+39*(d>9) | ||
| 611 | x >>= 4 | ||
| 612 | i -= 1 | ||
| 613 | end | ||
| 614 | if neg; @inbounds a[1]=0x2d; end | ||
| 615 | String(a) | ||
| 616 | end | ||
| 617 | |||
| 618 | const base36digits = ['0':'9';'a':'z'] | ||
| 619 | const base62digits = ['0':'9';'A':'Z';'a':'z'] | ||
| 620 | |||
| 621 | function _base(b::Integer, x::Integer, pad::Integer, neg::Bool) | ||
| 622 | (x >= 0) | (b < 0) || throw(DomainError(x, "For negative `x`, `b` must be negative.")) | ||
| 623 | 2 <= abs(b) <= 62 || throw(DomainError(b, "base must satisfy 2 ≤ abs(base) ≤ 62")) | ||
| 624 | digits = abs(b) <= 36 ? base36digits : base62digits | ||
| 625 | i = neg + ndigits(x, base=b, pad=pad) | ||
| 626 | a = StringVector(i) | ||
| 627 | @inbounds while i > neg | ||
| 628 | if b > 0 | ||
| 629 | a[i] = digits[1+rem(x,b)] | ||
| 630 | x = div(x,b) | ||
| 631 | else | ||
| 632 | a[i] = digits[1+mod(x,-b)] | ||
| 633 | x = cld(x,b) | ||
| 634 | end | ||
| 635 | i -= 1 | ||
| 636 | end | ||
| 637 | if neg; a[1]='-'; end | ||
| 638 | String(a) | ||
| 639 | end | ||
| 640 | |||
| 641 | split_sign(n::Integer) = unsigned(abs(n)), n < 0 | ||
| 642 | split_sign(n::Unsigned) = n, false | ||
| 643 | |||
| 644 | """ | ||
| 645 | string(n::Integer; base::Integer = 10, pad::Integer = 1) | ||
| 646 | |||
| 647 | Convert an integer `n` to a string in the given `base`, | ||
| 648 | optionally specifying a number of digits to pad to. | ||
| 649 | |||
| 650 | ```jldoctest | ||
| 651 | julia> string(5, base = 13, pad = 4) | ||
| 652 | "0005" | ||
| 653 | |||
| 654 | julia> string(13, base = 5, pad = 4) | ||
| 655 | "0023" | ||
| 656 | ``` | ||
| 657 | """ | ||
| 658 | function string(n::Integer; base::Integer = 10, pad::Integer = 1) | ||
| 659 | if base == 2 | ||
| 660 | (n_positive, neg) = split_sign(n) | ||
| 661 | bin(n_positive, pad, neg) | ||
| 662 | elseif base == 8 | ||
| 663 | (n_positive, neg) = split_sign(n) | ||
| 664 | oct(n_positive, pad, neg) | ||
| 665 | elseif base == 10 | ||
| 666 | (n_positive, neg) = split_sign(n) | ||
| 667 | dec(n_positive, pad, neg) | ||
| 668 | elseif base == 16 | ||
| 669 | (n_positive, neg) = split_sign(n) | ||
| 670 | hex(n_positive, pad, neg) | ||
| 671 | else | ||
| 672 | _base(base, base > 0 ? unsigned(abs(n)) : convert(Signed, n), pad, (base>0) & (n<0)) | ||
| 673 | end | ||
| 674 | end | ||
| 675 | |||
| 676 | string(b::Bool) = b ? "true" : "false" | ||
| 677 | |||
| 678 | """ | ||
| 679 | bitstring(n) | ||
| 680 | |||
| 681 | A string giving the literal bit representation of a number. | ||
| 682 | |||
| 683 | # Examples | ||
| 684 | ```jldoctest | ||
| 685 | julia> bitstring(4) | ||
| 686 | "0000000000000000000000000000000000000000000000000000000000000100" | ||
| 687 | |||
| 688 | julia> bitstring(2.2) | ||
| 689 | "0100000000000001100110011001100110011001100110011001100110011010" | ||
| 690 | ``` | ||
| 691 | """ | ||
| 692 | function bitstring end | ||
| 693 | |||
| 694 | bitstring(x::Union{Bool,Int8,UInt8}) = string(reinterpret(UInt8,x), pad = 8, base = 2) | ||
| 695 | bitstring(x::Union{Int16,UInt16,Float16}) = string(reinterpret(UInt16,x), pad = 16, base = 2) | ||
| 696 | bitstring(x::Union{Char,Int32,UInt32,Float32}) = string(reinterpret(UInt32,x), pad = 32, base = 2) | ||
| 697 | bitstring(x::Union{Int64,UInt64,Float64}) = string(reinterpret(UInt64,x), pad = 64, base = 2) | ||
| 698 | bitstring(x::Union{Int128,UInt128}) = string(reinterpret(UInt128,x), pad = 128, base = 2) | ||
| 699 | |||
| 700 | """ | ||
| 701 | digits([T<:Integer], n::Integer; base::T = 10, pad::Integer = 1) | ||
| 702 | |||
| 703 | Return an array with element type `T` (default `Int`) of the digits of `n` in the given | ||
| 704 | base, optionally padded with zeros to a specified size. More significant digits are at | ||
| 705 | higher indices, such that `n == sum([digits[k]*base^(k-1) for k=1:length(digits)])`. | ||
| 706 | |||
| 707 | # Examples | ||
| 708 | ```jldoctest | ||
| 709 | julia> digits(10, base = 10) | ||
| 710 | 2-element Array{Int64,1}: | ||
| 711 | 0 | ||
| 712 | 1 | ||
| 713 | |||
| 714 | julia> digits(10, base = 2) | ||
| 715 | 4-element Array{Int64,1}: | ||
| 716 | 0 | ||
| 717 | 1 | ||
| 718 | 0 | ||
| 719 | 1 | ||
| 720 | |||
| 721 | julia> digits(10, base = 2, pad = 6) | ||
| 722 | 6-element Array{Int64,1}: | ||
| 723 | 0 | ||
| 724 | 1 | ||
| 725 | 0 | ||
| 726 | 1 | ||
| 727 | 0 | ||
| 728 | 0 | ||
| 729 | ``` | ||
| 730 | """ | ||
| 731 | digits(n::Integer; base::Integer = 10, pad::Integer = 1) = | ||
| 732 | digits(typeof(base), n, base = base, pad = pad) | ||
| 733 | |||
| 734 | function digits(T::Type{<:Integer}, n::Integer; base::Integer = 10, pad::Integer = 1) | ||
| 735 | digits!(zeros(T, ndigits(n, base=base, pad=pad)), n, base=base) | ||
| 736 | end | ||
| 737 | |||
| 738 | """ | ||
| 739 | hastypemax(T::Type) -> Bool | ||
| 740 | |||
| 741 | Return `true` if and only if `typemax(T)` is defined. | ||
| 742 | """ | ||
| 743 | hastypemax(::Base.BitIntegerType) = true | ||
| 744 | hastypemax(::Type{T}) where {T} = applicable(typemax, T) | ||
| 745 | |||
| 746 | """ | ||
| 747 | digits!(array, n::Integer; base::Integer = 10) | ||
| 748 | |||
| 749 | Fills an array of the digits of `n` in the given base. More significant digits are at higher | ||
| 750 | indices. If the array length is insufficient, the least significant digits are filled up to | ||
| 751 | the array length. If the array length is excessive, the excess portion is filled with zeros. | ||
| 752 | |||
| 753 | # Examples | ||
| 754 | ```jldoctest | ||
| 755 | julia> digits!([2,2,2,2], 10, base = 2) | ||
| 756 | 4-element Array{Int64,1}: | ||
| 757 | 0 | ||
| 758 | 1 | ||
| 759 | 0 | ||
| 760 | 1 | ||
| 761 | |||
| 762 | julia> digits!([2,2,2,2,2,2], 10, base = 2) | ||
| 763 | 6-element Array{Int64,1}: | ||
| 764 | 0 | ||
| 765 | 1 | ||
| 766 | 0 | ||
| 767 | 1 | ||
| 768 | 0 | ||
| 769 | 0 | ||
| 770 | ``` | ||
| 771 | """ | ||
| 772 | function digits!(a::AbstractVector{T}, n::Integer; base::Integer = 10) where T<:Integer | ||
| 773 | 2 <= abs(base) || throw(DomainError(base, "base must be ≥ 2 or ≤ -2")) | ||
| 774 | hastypemax(T) && abs(base) - 1 > typemax(T) && | ||
| 775 | throw(ArgumentError("type $T too small for base $base")) | ||
| 776 | isempty(a) && return a | ||
| 777 | |||
| 778 | if base > 0 | ||
| 779 | if ispow2(base) && n >= 0 && n isa Base.BitInteger && base <= typemax(Int) | ||
| 780 | base = Int(base) | ||
| 781 | k = trailing_zeros(base) | ||
| 782 | c = base - 1 | ||
| 783 | for i in eachindex(a) | ||
| 784 | a[i] = (n >> (k * (i - firstindex(a)))) & c | ||
| 785 | end | ||
| 786 | else | ||
| 787 | for i in eachindex(a) | ||
| 788 | n, d = divrem(n, base) | ||
| 789 | a[i] = d | ||
| 790 | end | ||
| 791 | end | ||
| 792 | else | ||
| 793 | # manually peel one loop iteration for type stability | ||
| 794 | n, d = fldmod(n, -base) | ||
| 795 | a[firstindex(a)] = d | ||
| 796 | n = -signed(n) | ||
| 797 | for i in firstindex(a)+1:lastindex(a) | ||
| 798 | n, d = fldmod(n, -base) | ||
| 799 | a[i] = d | ||
| 800 | n = -n | ||
| 801 | end | ||
| 802 | end | ||
| 803 | return a | ||
| 804 | end | ||
| 805 | |||
| 806 | """ | ||
| 807 | isqrt(n::Integer) | ||
| 808 | |||
| 809 | Integer square root: the largest integer `m` such that `m*m <= n`. | ||
| 810 | |||
| 811 | ```jldoctest | ||
| 812 | julia> isqrt(5) | ||
| 813 | 2 | ||
| 814 | ``` | ||
| 815 | """ | ||
| 816 | isqrt(x::Integer) = oftype(x, trunc(sqrt(x))) | ||
| 817 | |||
| 818 | function isqrt(x::Union{Int64,UInt64,Int128,UInt128}) | ||
| 819 | x==0 && return x | ||
| 820 | s = oftype(x, trunc(sqrt(x))) | ||
| 821 | # fix with a Newton iteration, since conversion to float discards | ||
| 822 | # too many bits. | ||
| 823 | s = (s + div(x,s)) >> 1 | ||
| 824 | s*s > x ? s-1 : s | ||
| 825 | end | ||
| 826 | |||
| 827 | """ | ||
| 828 | factorial(n::Integer) | ||
| 829 | |||
| 830 | Factorial of `n`. If `n` is an [`Integer`](@ref), the factorial is computed as an | ||
| 831 | integer (promoted to at least 64 bits). Note that this may overflow if `n` is not small, | ||
| 832 | but you can use `factorial(big(n))` to compute the result exactly in arbitrary precision. | ||
| 833 | |||
| 834 | # Examples | ||
| 835 | ```jldoctest | ||
| 836 | julia> factorial(6) | ||
| 837 | 720 | ||
| 838 | |||
| 839 | julia> factorial(21) | ||
| 840 | ERROR: OverflowError: 21 is too large to look up in the table; consider using `factorial(big(21))` instead | ||
| 841 | Stacktrace: | ||
| 842 | [...] | ||
| 843 | |||
| 844 | julia> factorial(big(21)) | ||
| 845 | 51090942171709440000 | ||
| 846 | ``` | ||
| 847 | |||
| 848 | # See also | ||
| 849 | * [`binomial`](@ref) | ||
| 850 | |||
| 851 | # External links | ||
| 852 | * [Factorial](https://en.wikipedia.org/wiki/Factorial) on Wikipedia. | ||
| 853 | """ | ||
| 854 | function factorial(n::Integer) | ||
| 855 | n < 0 && throw(DomainError(n, "`n` must be nonnegative.")) | ||
| 856 | f::typeof(n*n) = 1 | ||
| 857 | for i::typeof(n*n) = 2:n | ||
| 858 | f *= i | ||
| 859 | end | ||
| 860 | return f | ||
| 861 | end | ||
| 862 | |||
| 863 | """ | ||
| 864 | binomial(n::Integer, k::Integer) | ||
| 865 | |||
| 866 | The _binomial coefficient_ ``\\binom{n}{k}``, being the coefficient of the ``k``th term in | ||
| 867 | the polynomial expansion of ``(1+x)^n``. | ||
| 868 | |||
| 869 | If ``n`` is non-negative, then it is the number of ways to choose `k` out of `n` items: | ||
| 870 | ```math | ||
| 871 | \\binom{n}{k} = \\frac{n!}{k! (n-k)!} | ||
| 872 | ``` | ||
| 873 | where ``n!`` is the [`factorial`](@ref) function. | ||
| 874 | |||
| 875 | If ``n`` is negative, then it is defined in terms of the identity | ||
| 876 | ```math | ||
| 877 | \\binom{n}{k} = (-1)^k \\binom{k-n-1}{k} | ||
| 878 | ``` | ||
| 879 | |||
| 880 | # Examples | ||
| 881 | ```jldoctest | ||
| 882 | julia> binomial(5, 3) | ||
| 883 | 10 | ||
| 884 | |||
| 885 | julia> factorial(5) ÷ (factorial(5-3) * factorial(3)) | ||
| 886 | 10 | ||
| 887 | |||
| 888 | julia> binomial(-5, 3) | ||
| 889 | -35 | ||
| 890 | ``` | ||
| 891 | |||
| 892 | # See also | ||
| 893 | * [`factorial`](@ref) | ||
| 894 | |||
| 895 | # External links | ||
| 896 | * [Binomial coeffient](https://en.wikipedia.org/wiki/Binomial_coefficient) on Wikipedia. | ||
| 897 | """ | ||
| 898 | function binomial(n::T, k::T) where T<:Integer | ||
| 899 | n0, k0 = n, k | ||
| 900 | k < 0 && return zero(T) | ||
| 901 | sgn = one(T) | ||
| 902 | if n < 0 | ||
| 903 | n = -n + k -1 | ||
| 904 | if isodd(k) | ||
| 905 | sgn = -sgn | ||
| 906 | end | ||
| 907 | end | ||
| 908 | k > n && return zero(T) | ||
| 909 | (k == 0 || k == n) && return sgn | ||
| 910 | k == 1 && return sgn*n | ||
| 911 | if k > (n>>1) | ||
| 912 | k = (n - k) | ||
| 913 | end | ||
| 914 | x::T = nn = n - k + 1 | ||
| 915 | nn += 1 | ||
| 916 | rr = 2 | ||
| 917 | while rr <= k | ||
| 918 | xt = div(widemul(x, nn), rr) | ||
| 919 | x = xt % T | ||
| 920 | x == xt || throw(OverflowError("binomial($n0, $k0) overflows")) | ||
| 921 | rr += 1 | ||
| 922 | nn += 1 | ||
| 923 | end | ||
| 924 | convert(T, copysign(x, sgn)) | ||
| 925 | end |