Report index
StatProfilerHTML.jl report
Generated on Mon, 01 Apr 2024 21:01:18
File source code
Line Exclusive Inclusive Code
1 # This file is a part of Julia. License is MIT: https://julialang.org/license
2
3 module Math
4
5 export sin, cos, sincos, tan, sinh, cosh, tanh, asin, acos, atan,
6 asinh, acosh, atanh, sec, csc, cot, asec, acsc, acot,
7 sech, csch, coth, asech, acsch, acoth,
8 sinpi, cospi, sincospi, tanpi, sinc, cosc,
9 cosd, cotd, cscd, secd, sind, tand, sincosd,
10 acosd, acotd, acscd, asecd, asind, atand,
11 rad2deg, deg2rad,
12 log, log2, log10, log1p, exponent, exp, exp2, exp10, expm1,
13 cbrt, sqrt, fourthroot, significand,
14 hypot, max, min, minmax, ldexp, frexp,
15 clamp, clamp!, modf, ^, mod2pi, rem2pi,
16 @evalpoly, evalpoly
17
18 import .Base: log, exp, sin, cos, tan, sinh, cosh, tanh, asin,
19 acos, atan, asinh, acosh, atanh, sqrt, log2, log10,
20 max, min, minmax, ^, exp2, muladd, rem,
21 exp10, expm1, log1p, @constprop, @assume_effects
22
23 using .Base: sign_mask, exponent_mask, exponent_one,
24 exponent_half, uinttype, significand_mask,
25 significand_bits, exponent_bits, exponent_bias,
26 exponent_max, exponent_raw_max
27
28 using Core.Intrinsics: sqrt_llvm
29
30 using .Base: IEEEFloat
31
32 @noinline function throw_complex_domainerror(f::Symbol, x)
33 throw(DomainError(x,
34 LazyString(f," was called with a negative real argument but will only return a complex result if called with a complex argument. Try ", f,"(Complex(x)).")))
35 end
36 @noinline function throw_complex_domainerror_neg1(f::Symbol, x)
37 throw(DomainError(x,
38 LazyString(f," was called with a real argument < -1 but will only return a complex result if called with a complex argument. Try ", f,"(Complex(x)).")))
39 end
40 @noinline function throw_exp_domainerror(x)
41 throw(DomainError(x, LazyString(
42 "Exponentiation yielding a complex result requires a ",
43 "complex argument.\nReplace x^y with (x+0im)^y, ",
44 "Complex(x)^y, or similar.")))
45 end
46
47 # non-type specific math functions
48
49 function two_mul(x::T, y::T) where {T<:Number}
50 xy = x*y
51 xy, fma(x, y, -xy)
52 end
53
54 @assume_effects :consistent @inline function two_mul(x::Float64, y::Float64)
55 if Core.Intrinsics.have_fma(Float64)
56 xy = x*y
57 return xy, fma(x, y, -xy)
58 end
59 return Base.twomul(x,y)
60 end
61
62 @assume_effects :consistent @inline function two_mul(x::T, y::T) where T<: Union{Float16, Float32}
63 if Core.Intrinsics.have_fma(T)
64 xy = x*y
65 return xy, fma(x, y, -xy)
66 end
67 xy = widen(x)*y
68 Txy = T(xy)
69 return Txy, T(xy-Txy)
70 end
71
72 """
73 clamp(x, lo, hi)
74
75 Return `x` if `lo <= x <= hi`. If `x > hi`, return `hi`. If `x < lo`, return `lo`. Arguments
76 are promoted to a common type.
77
78 See also [`clamp!`](@ref), [`min`](@ref), [`max`](@ref).
79
80 !!! compat "Julia 1.3"
81 `missing` as the first argument requires at least Julia 1.3.
82
83 # Examples
84 ```jldoctest
85 julia> clamp.([pi, 1.0, big(10)], 2.0, 9.0)
86 3-element Vector{BigFloat}:
87 3.141592653589793238462643383279502884197169399375105820974944592307816406286198
88 2.0
89 9.0
90
91 julia> clamp.([11, 8, 5], 10, 6) # an example where lo > hi
92 3-element Vector{Int64}:
93 6
94 6
95 10
96 ```
97 """
98 clamp(x::X, lo::L, hi::H) where {X,L,H} =
99 ifelse(x > hi, convert(promote_type(X,L,H), hi),
100 ifelse(x < lo,
101 convert(promote_type(X,L,H), lo),
102 convert(promote_type(X,L,H), x)))
103
104 """
105 clamp(x, T)::T
106
107 Clamp `x` between `typemin(T)` and `typemax(T)` and convert the result to type `T`.
108
109 See also [`trunc`](@ref).
110
111 # Examples
112 ```jldoctest
113 julia> clamp(200, Int8)
114 127
115
116 julia> clamp(-200, Int8)
117 -128
118
119 julia> trunc(Int, 4pi^2)
120 39
121 ```
122 """
123 clamp(x, ::Type{T}) where {T<:Integer} = clamp(x, typemin(T), typemax(T)) % T
124
125
126 """
127 clamp!(array::AbstractArray, lo, hi)
128
129 Restrict values in `array` to the specified range, in-place.
130 See also [`clamp`](@ref).
131
132 !!! compat "Julia 1.3"
133 `missing` entries in `array` require at least Julia 1.3.
134
135 # Examples
136 ```jldoctest
137 julia> row = collect(-4:4)';
138
139 julia> clamp!(row, 0, Inf)
140 1×9 adjoint(::Vector{Int64}) with eltype Int64:
141 0 0 0 0 0 1 2 3 4
142
143 julia> clamp.((-4:4)', 0, Inf)
144 1×9 Matrix{Float64}:
145 0.0 0.0 0.0 0.0 0.0 1.0 2.0 3.0 4.0
146 ```
147 """
148 function clamp!(x::AbstractArray, lo, hi)
149 @inbounds for i in eachindex(x)
150 x[i] = clamp(x[i], lo, hi)
151 end
152 x
153 end
154
155 """
156 clamp(x::Integer, r::AbstractUnitRange)
157
158 Clamp `x` to lie within range `r`.
159
160 !!! compat "Julia 1.6"
161 This method requires at least Julia 1.6.
162 """
163 clamp(x::Integer, r::AbstractUnitRange{<:Integer}) = clamp(x, first(r), last(r))
164
165 """
166 evalpoly(x, p)
167
168 Evaluate the polynomial ``\\sum_k x^{k-1} p[k]`` for the coefficients `p[1]`, `p[2]`, ...;
169 that is, the coefficients are given in ascending order by power of `x`.
170 Loops are unrolled at compile time if the number of coefficients is statically known, i.e.
171 when `p` is a `Tuple`.
172 This function generates efficient code using Horner's method if `x` is real, or using
173 a Goertzel-like [^DK62] algorithm if `x` is complex.
174
175 [^DK62]: Donald Knuth, Art of Computer Programming, Volume 2: Seminumerical Algorithms, Sec. 4.6.4.
176
177 !!! compat "Julia 1.4"
178 This function requires Julia 1.4 or later.
179
180 # Example
181 ```jldoctest
182 julia> evalpoly(2, (1, 2, 3))
183 17
184 ```
185 """
186 1 (2 %)
1 (2 %) samples spent in evalpoly
1 (100 %) (incl.) when called from expm1b_kernel line 78
1 (100 %) samples spent calling macro expansion
function evalpoly(x, p::Tuple)
187 1 (2 %)
1 (2 %) samples spent in macro expansion
1 (100 %) (incl.) when called from evalpoly line 186
1 (100 %) samples spent calling muladd
if @generated
188 N = length(p.parameters::Core.SimpleVector)
189 ex = :(p[end])
190 for i in N-1:-1:1
191 ex = :(muladd(x, $ex, p[$i]))
192 end
193 ex
194 else
195 _evalpoly(x, p)
196 end
197 end
198
199 evalpoly(x, p::AbstractVector) = _evalpoly(x, p)
200
201 function _evalpoly(x, p)
202 Base.require_one_based_indexing(p)
203 N = length(p)
204 ex = p[end]
205 for i in N-1:-1:1
206 ex = muladd(x, ex, p[i])
207 end
208 ex
209 end
210
211 function evalpoly(z::Complex, p::Tuple)
212 if @generated
213 N = length(p.parameters)
214 a = :(p[end])
215 b = :(p[end-1])
216 as = []
217 for i in N-2:-1:1
218 ai = Symbol("a", i)
219 push!(as, :($ai = $a))
220 a = :(muladd(r, $ai, $b))
221 b = :(muladd(-s, $ai, p[$i]))
222 end
223 ai = :a0
224 push!(as, :($ai = $a))
225 C = Expr(:block,
226 :(x = real(z)),
227 :(y = imag(z)),
228 :(r = x + x),
229 :(s = muladd(x, x, y*y)),
230 as...,
231 :(muladd($ai, z, $b)))
232 else
233 _evalpoly(z, p)
234 end
235 end
236 evalpoly(z::Complex, p::Tuple{<:Any}) = p[1]
237
238
239 evalpoly(z::Complex, p::AbstractVector) = _evalpoly(z, p)
240
241 function _evalpoly(z::Complex, p)
242 Base.require_one_based_indexing(p)
243 length(p) == 1 && return p[1]
244 N = length(p)
245 a = p[end]
246 b = p[end-1]
247
248 x = real(z)
249 y = imag(z)
250 r = 2x
251 s = muladd(x, x, y*y)
252 for i in N-2:-1:1
253 ai = a
254 a = muladd(r, ai, b)
255 b = muladd(-s, ai, p[i])
256 end
257 ai = a
258 muladd(ai, z, b)
259 end
260
261 """
262 @horner(x, p...)
263
264 Evaluate `p[1] + x * (p[2] + x * (....))`, i.e. a polynomial via Horner's rule.
265
266 See also [`@evalpoly`](@ref), [`evalpoly`](@ref).
267 """
268 macro horner(x, p...)
269 xesc, pesc = esc(x), esc.(p)
270 :(invoke(evalpoly, Tuple{Any, Tuple}, $xesc, ($(pesc...),)))
271 end
272
273 # Evaluate p[1] + z*p[2] + z^2*p[3] + ... + z^(n-1)*p[n]. This uses
274 # Horner's method if z is real, but for complex z it uses a more
275 # efficient algorithm described in Knuth, TAOCP vol. 2, section 4.6.4,
276 # equation (3).
277
278 """
279 @evalpoly(z, c...)
280
281 Evaluate the polynomial ``\\sum_k z^{k-1} c[k]`` for the coefficients `c[1]`, `c[2]`, ...;
282 that is, the coefficients are given in ascending order by power of `z`. This macro expands
283 to efficient inline code that uses either Horner's method or, for complex `z`, a more
284 efficient Goertzel-like algorithm.
285
286 See also [`evalpoly`](@ref).
287
288 # Examples
289 ```jldoctest
290 julia> @evalpoly(3, 1, 0, 1)
291 10
292
293 julia> @evalpoly(2, 1, 0, 1)
294 5
295
296 julia> @evalpoly(2, 1, 1, 1)
297 7
298 ```
299 """
300 macro evalpoly(z, p...)
301 zesc, pesc = esc(z), esc.(p)
302 :(evalpoly($zesc, ($(pesc...),)))
303 end
304
305 # polynomial evaluation using compensated summation.
306 # much more accurate, especially when lo can be combined with other rounding errors
307 Base.@assume_effects :terminates_locally @inline function exthorner(x, p::Tuple)
308 hi, lo = p[end], zero(x)
309 for i in length(p)-1:-1:1
310 pi = getfield(p, i) # needed to prove consistency
311 prod, err = two_mul(hi,x)
312 hi = pi+prod
313 lo = fma(lo, x, prod - (hi - pi) + err)
314 end
315 return hi, lo
316 end
317
318 """
319 rad2deg(x)
320
321 Convert `x` from radians to degrees.
322
323 See also [`deg2rad`](@ref).
324
325 # Examples
326 ```jldoctest
327 julia> rad2deg(pi)
328 180.0
329 ```
330 """
331 rad2deg(z::AbstractFloat) = z * (180 / oftype(z, pi))
332
333 """
334 deg2rad(x)
335
336 Convert `x` from degrees to radians.
337
338 See also [`rad2deg`](@ref), [`sind`](@ref), [`pi`](@ref).
339
340 # Examples
341 ```jldoctest
342 julia> deg2rad(90)
343 1.5707963267948966
344 ```
345 """
346 deg2rad(z::AbstractFloat) = z * (oftype(z, pi) / 180)
347 rad2deg(z::Real) = rad2deg(float(z))
348 deg2rad(z::Real) = deg2rad(float(z))
349 rad2deg(z::Number) = (z/pi)*180
350 deg2rad(z::Number) = (z*pi)/180
351
352 log(b::T, x::T) where {T<:Number} = log(x)/log(b)
353
354 """
355 log(b,x)
356
357 Compute the base `b` logarithm of `x`. Throws [`DomainError`](@ref) for negative
358 [`Real`](@ref) arguments.
359
360 # Examples
361 ```jldoctest; filter = r"Stacktrace:(\\n \\[[0-9]+\\].*)*"
362 julia> log(4,8)
363 1.5
364
365 julia> log(4,2)
366 0.5
367
368 julia> log(-2, 3)
369 ERROR: DomainError with -2.0:
370 log was called with a negative real argument but will only return a complex result if called with a complex argument. Try log(Complex(x)).
371 Stacktrace:
372 [1] throw_complex_domainerror(::Symbol, ::Float64) at ./math.jl:31
373 [...]
374
375 julia> log(2, -3)
376 ERROR: DomainError with -3.0:
377 log was called with a negative real argument but will only return a complex result if called with a complex argument. Try log(Complex(x)).
378 Stacktrace:
379 [1] throw_complex_domainerror(::Symbol, ::Float64) at ./math.jl:31
380 [...]
381 ```
382
383 !!! note
384 If `b` is a power of 2 or 10, [`log2`](@ref) or [`log10`](@ref) should be used, as these will
385 typically be faster and more accurate. For example,
386
387 ```jldoctest
388 julia> log(100,1000000)
389 2.9999999999999996
390
391 julia> log10(1000000)/2
392 3.0
393 ```
394 """
395 log(b::Number, x::Number) = log(promote(b,x)...)
396
397 # type specific math functions
398
399 const libm = Base.libm_name
400
401 # functions with no domain error
402 """
403 sinh(x)
404
405 Compute hyperbolic sine of `x`.
406 """
407 sinh(x::Number)
408
409 """
410 cosh(x)
411
412 Compute hyperbolic cosine of `x`.
413 """
414 cosh(x::Number)
415
416 """
417 tanh(x)
418
419 Compute hyperbolic tangent of `x`.
420
421 See also [`tan`](@ref), [`atanh`](@ref).
422
423 # Examples
424
425 ```jldoctest
426 julia> tanh.(-3:3f0) # Here 3f0 isa Float32
427 7-element Vector{Float32}:
428 -0.9950548
429 -0.9640276
430 -0.7615942
431 0.0
432 0.7615942
433 0.9640276
434 0.9950548
435
436 julia> tan.(im .* (1:3))
437 3-element Vector{ComplexF64}:
438 0.0 + 0.7615941559557649im
439 0.0 + 0.9640275800758169im
440 0.0 + 0.9950547536867306im
441 ```
442 """
443 tanh(x::Number)
444
445 """
446 atan(y)
447 atan(y, x)
448
449 Compute the inverse tangent of `y` or `y/x`, respectively.
450
451 For one argument, this is the angle in radians between the positive *x*-axis and the point
452 (1, *y*), returning a value in the interval ``[-\\pi/2, \\pi/2]``.
453
454 For two arguments, this is the angle in radians between the positive *x*-axis and the
455 point (*x*, *y*), returning a value in the interval ``[-\\pi, \\pi]``. This corresponds to a
456 standard [`atan2`](https://en.wikipedia.org/wiki/Atan2) function. Note that by convention
457 `atan(0.0,x)` is defined as ``\\pi`` and `atan(-0.0,x)` is defined as ``-\\pi`` when `x < 0`.
458
459 See also [`atand`](@ref) for degrees.
460
461 # Examples
462
463 ```jldoctest
464 julia> rad2deg(atan(-1/√3))
465 -30.000000000000004
466
467 julia> rad2deg(atan(-1, √3))
468 -30.000000000000004
469
470 julia> rad2deg(atan(1, -√3))
471 150.0
472 ```
473 """
474 atan(x::Number)
475
476 """
477 asinh(x)
478
479 Compute the inverse hyperbolic sine of `x`.
480 """
481 asinh(x::Number)
482
483
484 # utility for converting NaN return to DomainError
485 # the branch in nan_dom_err prevents its callers from inlining, so be sure to force it
486 # until the heuristics can be improved
487 @inline nan_dom_err(out, x) = isnan(out) & !isnan(x) ? throw(DomainError(x, "NaN result for non-NaN input.")) : out
488
489 # functions that return NaN on non-NaN argument for domain error
490 """
491 sin(x)
492
493 Compute sine of `x`, where `x` is in radians.
494
495 See also [`sind`](@ref), [`sinpi`](@ref), [`sincos`](@ref), [`cis`](@ref), [`asin`](@ref).
496
497 # Examples
498 ```jldoctest
499 julia> round.(sin.(range(0, 2pi, length=9)'), digits=3)
500 1×9 Matrix{Float64}:
501 0.0 0.707 1.0 0.707 0.0 -0.707 -1.0 -0.707 -0.0
502
503 julia> sind(45)
504 0.7071067811865476
505
506 julia> sinpi(1/4)
507 0.7071067811865475
508
509 julia> round.(sincos(pi/6), digits=3)
510 (0.5, 0.866)
511
512 julia> round(cis(pi/6), digits=3)
513 0.866 + 0.5im
514
515 julia> round(exp(im*pi/6), digits=3)
516 0.866 + 0.5im
517 ```
518 """
519 sin(x::Number)
520
521 """
522 cos(x)
523
524 Compute cosine of `x`, where `x` is in radians.
525
526 See also [`cosd`](@ref), [`cospi`](@ref), [`sincos`](@ref), [`cis`](@ref).
527 """
528 cos(x::Number)
529
530 """
531 tan(x)
532
533 Compute tangent of `x`, where `x` is in radians.
534 """
535 tan(x::Number)
536
537 """
538 asin(x)
539
540 Compute the inverse sine of `x`, where the output is in radians.
541
542 See also [`asind`](@ref) for output in degrees.
543
544 # Examples
545 ```jldoctest
546 julia> asin.((0, 1/2, 1))
547 (0.0, 0.5235987755982989, 1.5707963267948966)
548
549 julia> asind.((0, 1/2, 1))
550 (0.0, 30.000000000000004, 90.0)
551 ```
552 """
553 asin(x::Number)
554
555 """
556 acos(x)
557
558 Compute the inverse cosine of `x`, where the output is in radians
559 """
560 acos(x::Number)
561
562 """
563 acosh(x)
564
565 Compute the inverse hyperbolic cosine of `x`.
566 """
567 acosh(x::Number)
568
569 """
570 atanh(x)
571
572 Compute the inverse hyperbolic tangent of `x`.
573 """
574 atanh(x::Number)
575
576 """
577 log(x)
578
579 Compute the natural logarithm of `x`. Throws [`DomainError`](@ref) for negative
580 [`Real`](@ref) arguments. Use complex negative arguments to obtain complex results.
581
582 See also [`ℯ`](@ref), [`log1p`](@ref), [`log2`](@ref), [`log10`](@ref).
583
584 # Examples
585 ```jldoctest; filter = r"Stacktrace:(\\n \\[[0-9]+\\].*)*"
586 julia> log(2)
587 0.6931471805599453
588
589 julia> log(-3)
590 ERROR: DomainError with -3.0:
591 log was called with a negative real argument but will only return a complex result if called with a complex argument. Try log(Complex(x)).
592 Stacktrace:
593 [1] throw_complex_domainerror(::Symbol, ::Float64) at ./math.jl:31
594 [...]
595
596 julia> log.(exp.(-1:1))
597 3-element Vector{Float64}:
598 -1.0
599 0.0
600 1.0
601 ```
602 """
603 log(x::Number)
604
605 """
606 log2(x)
607
608 Compute the logarithm of `x` to base 2. Throws [`DomainError`](@ref) for negative
609 [`Real`](@ref) arguments.
610
611 See also: [`exp2`](@ref), [`ldexp`](@ref), [`ispow2`](@ref).
612
613 # Examples
614 ```jldoctest; filter = r"Stacktrace:(\\n \\[[0-9]+\\].*)*"
615 julia> log2(4)
616 2.0
617
618 julia> log2(10)
619 3.321928094887362
620
621 julia> log2(-2)
622 ERROR: DomainError with -2.0:
623 log2 was called with a negative real argument but will only return a complex result if called with a complex argument. Try log2(Complex(x)).
624 Stacktrace:
625 [1] throw_complex_domainerror(f::Symbol, x::Float64) at ./math.jl:31
626 [...]
627
628 julia> log2.(2.0 .^ (-1:1))
629 3-element Vector{Float64}:
630 -1.0
631 0.0
632 1.0
633 ```
634 """
635 log2(x)
636
637 """
638 log10(x)
639
640 Compute the logarithm of `x` to base 10.
641 Throws [`DomainError`](@ref) for negative [`Real`](@ref) arguments.
642
643 # Examples
644 ```jldoctest; filter = r"Stacktrace:(\\n \\[[0-9]+\\].*)*"
645 julia> log10(100)
646 2.0
647
648 julia> log10(2)
649 0.3010299956639812
650
651 julia> log10(-2)
652 ERROR: DomainError with -2.0:
653 log10 was called with a negative real argument but will only return a complex result if called with a complex argument. Try log10(Complex(x)).
654 Stacktrace:
655 [1] throw_complex_domainerror(f::Symbol, x::Float64) at ./math.jl:31
656 [...]
657 ```
658 """
659 log10(x)
660
661 """
662 log1p(x)
663
664 Accurate natural logarithm of `1+x`. Throws [`DomainError`](@ref) for [`Real`](@ref)
665 arguments less than -1.
666
667 # Examples
668 ```jldoctest; filter = r"Stacktrace:(\\n \\[[0-9]+\\].*)*"
669 julia> log1p(-0.5)
670 -0.6931471805599453
671
672 julia> log1p(0)
673 0.0
674
675 julia> log1p(-2)
676 ERROR: DomainError with -2.0:
677 log1p was called with a real argument < -1 but will only return a complex result if called with a complex argument. Try log1p(Complex(x)).
678 Stacktrace:
679 [1] throw_complex_domainerror(::Symbol, ::Float64) at ./math.jl:31
680 [...]
681 ```
682 """
683 log1p(x)
684
685 @inline function sqrt(x::Union{Float32,Float64})
686 x < zero(x) && throw_complex_domainerror(:sqrt, x)
687 sqrt_llvm(x)
688 end
689
690 """
691 sqrt(x)
692
693 Return ``\\sqrt{x}``. Throws [`DomainError`](@ref) for negative [`Real`](@ref) arguments.
694 Use complex negative arguments instead. The prefix operator `√` is equivalent to `sqrt`.
695
696 See also: [`hypot`](@ref).
697
698 # Examples
699 ```jldoctest; filter = r"Stacktrace:(\\n \\[[0-9]+\\].*)*"
700 julia> sqrt(big(81))
701 9.0
702
703 julia> sqrt(big(-81))
704 ERROR: DomainError with -81.0:
705 NaN result for non-NaN input.
706 Stacktrace:
707 [1] sqrt(::BigFloat) at ./mpfr.jl:501
708 [...]
709
710 julia> sqrt(big(complex(-81)))
711 0.0 + 9.0im
712
713 julia> .√(1:4)
714 4-element Vector{Float64}:
715 1.0
716 1.4142135623730951
717 1.7320508075688772
718 2.0
719 ```
720 """
721 sqrt(x)
722
723 """
724 fourthroot(x)
725
726 Return the fourth root of `x` by applying `sqrt` twice successively.
727 """
728 fourthroot(x::Number) = sqrt(sqrt(x))
729
730 """
731 hypot(x, y)
732
733 Compute the hypotenuse ``\\sqrt{|x|^2+|y|^2}`` avoiding overflow and underflow.
734
735 This code is an implementation of the algorithm described in:
736 An Improved Algorithm for `hypot(a,b)`
737 by Carlos F. Borges
738 The article is available online at arXiv at the link
739 https://arxiv.org/abs/1904.09481
740
741 hypot(x...)
742
743 Compute the hypotenuse ``\\sqrt{\\sum |x_i|^2}`` avoiding overflow and underflow.
744
745 See also `norm` in the [`LinearAlgebra`](@ref man-linalg) standard library.
746
747 # Examples
748 ```jldoctest; filter = r"Stacktrace:(\\n \\[[0-9]+\\].*)*"
749 julia> a = Int64(10)^10;
750
751 julia> hypot(a, a)
752 1.4142135623730951e10
753
754 julia> √(a^2 + a^2) # a^2 overflows
755 ERROR: DomainError with -2.914184810805068e18:
756 sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
757 Stacktrace:
758 [...]
759
760 julia> hypot(3, 4im)
761 5.0
762
763 julia> hypot(-5.7)
764 5.7
765
766 julia> hypot(3, 4im, 12.0)
767 13.0
768
769 julia> using LinearAlgebra
770
771 julia> norm([a, a, a, a]) == hypot(a, a, a, a)
772 true
773 ```
774 """
775 hypot(x::Number) = abs(float(x))
776 hypot(x::Number, y::Number) = _hypot(promote(float(x), y)...)
777 hypot(x::Number, y::Number, xs::Number...) = _hypot(promote(float(x), y, xs...))
778 function _hypot(x, y)
779 # preserves unit
780 axu = abs(x)
781 ayu = abs(y)
782
783 # unitless
784 ax = axu / oneunit(axu)
785 ay = ayu / oneunit(ayu)
786
787 # Return Inf if either or both inputs is Inf (Compliance with IEEE754)
788 if isinf(ax) || isinf(ay)
789 return typeof(axu)(Inf)
790 end
791
792 # Order the operands
793 if ay > ax
794 axu, ayu = ayu, axu
795 ax, ay = ay, ax
796 end
797
798 # Widely varying operands
799 if ay <= ax*sqrt(eps(typeof(ax))/2) #Note: This also gets ay == 0
800 return axu
801 end
802
803 # Operands do not vary widely
804 scale = eps(typeof(ax))*sqrt(floatmin(ax)) #Rescaling constant
805 if ax > sqrt(floatmax(ax)/2)
806 ax = ax*scale
807 ay = ay*scale
808 scale = inv(scale)
809 elseif ay < sqrt(floatmin(ax))
810 ax = ax/scale
811 ay = ay/scale
812 else
813 scale = oneunit(scale)
814 end
815 h = sqrt(muladd(ax, ax, ay*ay))
816 # This branch is correctly rounded but requires a native hardware fma.
817 if Core.Intrinsics.have_fma(typeof(h))
818 hsquared = h*h
819 axsquared = ax*ax
820 h -= (fma(-ay, ay, hsquared-axsquared) + fma(h, h,-hsquared) - fma(ax, ax, -axsquared))/(2*h)
821 # This branch is within one ulp of correctly rounded.
822 else
823 if h <= 2*ay
824 delta = h-ay
825 h -= muladd(delta, delta-2*(ax-ay), ax*(2*delta - ax))/(2*h)
826 else
827 delta = h-ax
828 h -= muladd(delta, delta, muladd(ay, (4*delta - ay), 2*delta*(ax - 2*ay)))/(2*h)
829 end
830 end
831 return h*scale*oneunit(axu)
832 end
833 @inline function _hypot(x::Float32, y::Float32)
834 if isinf(x) || isinf(y)
835 return Inf32
836 end
837 _x, _y = Float64(x), Float64(y)
838 return Float32(sqrt(muladd(_x, _x, _y*_y)))
839 end
840 @inline function _hypot(x::Float16, y::Float16)
841 if isinf(x) || isinf(y)
842 return Inf16
843 end
844 _x, _y = Float32(x), Float32(y)
845 return Float16(sqrt(muladd(_x, _x, _y*_y)))
846 end
847 _hypot(x::ComplexF16, y::ComplexF16) = Float16(_hypot(ComplexF32(x), ComplexF32(y)))
848
849 function _hypot(x::NTuple{N,<:Number}) where {N}
850 maxabs = maximum(abs, x)
851 if isnan(maxabs) && any(isinf, x)
852 return typeof(maxabs)(Inf)
853 elseif (iszero(maxabs) || isinf(maxabs))
854 return maxabs
855 else
856 return maxabs * sqrt(sum(y -> abs2(y / maxabs), x))
857 end
858 end
859
860 function _hypot(x::NTuple{N,<:IEEEFloat}) where {N}
861 T = eltype(x)
862 infT = convert(T, Inf)
863 x = abs.(x) # doesn't change result but enables computational shortcuts
864 # note: any() was causing this to not inline for N=3 but mapreduce() was not
865 mapreduce(==(infT), |, x) && return infT # return Inf even if an argument is NaN
866 maxabs = reinterpret(T, maximum(z -> reinterpret(Signed, z), x)) # for abs(::IEEEFloat) values, a ::BitInteger cast does not change the result
867 maxabs > zero(T) || return maxabs # catch NaN before the @fastmath below, but also shortcut 0 since we can (remove if no more @fastmath below)
868 scale,invscale = scaleinv(maxabs)
869 # @fastmath(+) to allow reassociation (see #48129)
870 add_fast(x, y) = Core.Intrinsics.add_float_fast(x, y) # @fastmath is not available during bootstrap
871 return scale * sqrt(mapreduce(y -> abs2(y * invscale), add_fast, x))
872 end
873
874 atan(y::Real, x::Real) = atan(promote(float(y),float(x))...)
875 atan(y::T, x::T) where {T<:AbstractFloat} = Base.no_op_err("atan", T)
876
877 _isless(x::T, y::T) where {T<:AbstractFloat} = (x < y) || (signbit(x) > signbit(y))
878 min(x::T, y::T) where {T<:AbstractFloat} = isnan(x) || ~isnan(y) && _isless(x, y) ? x : y
879 max(x::T, y::T) where {T<:AbstractFloat} = isnan(x) || ~isnan(y) && _isless(y, x) ? x : y
880 minmax(x::T, y::T) where {T<:AbstractFloat} = min(x, y), max(x, y)
881
882 _isless(x::Float16, y::Float16) = signbit(widen(x) - widen(y))
883
884 const has_native_fminmax = Sys.ARCH === :aarch64
885 @static if has_native_fminmax
886 @eval begin
887 Base.@assume_effects :total @inline llvm_min(x::Float64, y::Float64) = ccall("llvm.minimum.f64", llvmcall, Float64, (Float64, Float64), x, y)
888 Base.@assume_effects :total @inline llvm_min(x::Float32, y::Float32) = ccall("llvm.minimum.f32", llvmcall, Float32, (Float32, Float32), x, y)
889 Base.@assume_effects :total @inline llvm_max(x::Float64, y::Float64) = ccall("llvm.maximum.f64", llvmcall, Float64, (Float64, Float64), x, y)
890 Base.@assume_effects :total @inline llvm_max(x::Float32, y::Float32) = ccall("llvm.maximum.f32", llvmcall, Float32, (Float32, Float32), x, y)
891 end
892 end
893
894 function min(x::T, y::T) where {T<:Union{Float32,Float64}}
895 @static if has_native_fminmax
896 return llvm_min(x,y)
897 end
898 diff = x - y
899 argmin = ifelse(signbit(diff), x, y)
900 anynan = isnan(x)|isnan(y)
901 return ifelse(anynan, diff, argmin)
902 end
903
904 function max(x::T, y::T) where {T<:Union{Float32,Float64}}
905 @static if has_native_fminmax
906 return llvm_max(x,y)
907 end
908 diff = x - y
909 argmax = ifelse(signbit(diff), y, x)
910 anynan = isnan(x)|isnan(y)
911 return ifelse(anynan, diff, argmax)
912 end
913
914 function minmax(x::T, y::T) where {T<:Union{Float32,Float64}}
915 @static if has_native_fminmax
916 return llvm_min(x, y), llvm_max(x, y)
917 end
918 diff = x - y
919 sdiff = signbit(diff)
920 min, max = ifelse(sdiff, x, y), ifelse(sdiff, y, x)
921 anynan = isnan(x)|isnan(y)
922 return ifelse(anynan, diff, min), ifelse(anynan, diff, max)
923 end
924
925 """
926 ldexp(x, n)
927
928 Compute ``x \\times 2^n``.
929
930 # Examples
931 ```jldoctest
932 julia> ldexp(5., 2)
933 20.0
934 ```
935 """
936 function ldexp(x::T, e::Integer) where T<:IEEEFloat
937 xu = reinterpret(Unsigned, x)
938 xs = xu & ~sign_mask(T)
939 xs >= exponent_mask(T) && return x # NaN or Inf
940 k = (xs >> significand_bits(T)) % Int
941 if k == 0 # x is subnormal
942 xs == 0 && return x # +-0
943 m = leading_zeros(xs) - exponent_bits(T)
944 ys = xs << unsigned(m)
945 xu = ys | (xu & sign_mask(T))
946 k = 1 - m
947 # underflow, otherwise may have integer underflow in the following n + k
948 e < -50000 && return flipsign(T(0.0), x)
949 end
950 # For cases where e of an Integer larger than Int make sure we properly
951 # overflow/underflow; this is optimized away otherwise.
952 if e > typemax(Int)
953 return flipsign(T(Inf), x)
954 elseif e < typemin(Int)
955 return flipsign(T(0.0), x)
956 end
957 n = e % Int
958 k += n
959 # overflow, if k is larger than maximum possible exponent
960 if k >= exponent_raw_max(T)
961 return flipsign(T(Inf), x)
962 end
963 if k > 0 # normal case
964 xu = (xu & ~exponent_mask(T)) | (rem(k, uinttype(T)) << significand_bits(T))
965 return reinterpret(T, xu)
966 else # subnormal case
967 if k <= -significand_bits(T) # underflow
968 # overflow, for the case of integer overflow in n + k
969 e > 50000 && return flipsign(T(Inf), x)
970 return flipsign(T(0.0), x)
971 end
972 k += significand_bits(T)
973 # z = T(2.0) ^ (-significand_bits(T))
974 z = reinterpret(T, rem(exponent_bias(T)-significand_bits(T), uinttype(T)) << significand_bits(T))
975 xu = (xu & ~exponent_mask(T)) | (rem(k, uinttype(T)) << significand_bits(T))
976 return z*reinterpret(T, xu)
977 end
978 end
979 ldexp(x::Float16, q::Integer) = Float16(ldexp(Float32(x), q))
980
981 """
982 exponent(x) -> Int
983
984 Returns the largest integer `y` such that `2^y ≤ abs(x)`.
985 For a normalized floating-point number `x`, this corresponds to the exponent of `x`.
986
987 # Examples
988 ```jldoctest
989 julia> exponent(8)
990 3
991
992 julia> exponent(64//1)
993 6
994
995 julia> exponent(6.5)
996 2
997
998 julia> exponent(16.0)
999 4
1000
1001 julia> exponent(3.142e-4)
1002 -12
1003 ```
1004 """
1005 function exponent(x::T) where T<:IEEEFloat
1006 @noinline throw1(x) = throw(DomainError(x, "Cannot be NaN or Inf."))
1007 @noinline throw2(x) = throw(DomainError(x, "Cannot be ±0.0."))
1008 xs = reinterpret(Unsigned, x) & ~sign_mask(T)
1009 xs >= exponent_mask(T) && throw1(x)
1010 k = Int(xs >> significand_bits(T))
1011 if k == 0 # x is subnormal
1012 xs == 0 && throw2(x)
1013 m = leading_zeros(xs) - exponent_bits(T)
1014 k = 1 - m
1015 end
1016 return k - exponent_bias(T)
1017 end
1018
1019 # Like exponent, but assumes the nothrow precondition. For
1020 # internal use only. Could be written as
1021 # @assume_effects :nothrow exponent()
1022 # but currently this form is easier on the compiler.
1023 function _exponent_finite_nonzero(x::T) where T<:IEEEFloat
1024 # @precond :nothrow !isnan(x) && !isinf(x) && !iszero(x)
1025 xs = reinterpret(Unsigned, x) & ~sign_mask(T)
1026 k = rem(xs >> significand_bits(T), Int)
1027 if k == 0 # x is subnormal
1028 m = leading_zeros(xs) - exponent_bits(T)
1029 k = 1 - m
1030 end
1031 return k - exponent_bias(T)
1032 end
1033
1034 """
1035 significand(x)
1036
1037 Extract the significand (a.k.a. mantissa) of a floating-point number. If `x` is
1038 a non-zero finite number, then the result will be a number of the same type and
1039 sign as `x`, and whose absolute value is on the interval ``[1,2)``. Otherwise
1040 `x` is returned.
1041
1042 # Examples
1043 ```jldoctest
1044 julia> significand(15.2)
1045 1.9
1046
1047 julia> significand(-15.2)
1048 -1.9
1049
1050 julia> significand(-15.2) * 2^3
1051 -15.2
1052
1053 julia> significand(-Inf), significand(Inf), significand(NaN)
1054 (-Inf, Inf, NaN)
1055 ```
1056 """
1057 function significand(x::T) where T<:IEEEFloat
1058 xu = reinterpret(Unsigned, x)
1059 xs = xu & ~sign_mask(T)
1060 xs >= exponent_mask(T) && return x # NaN or Inf
1061 if xs <= (~exponent_mask(T) & ~sign_mask(T)) # x is subnormal
1062 xs == 0 && return x # +-0
1063 m = unsigned(leading_zeros(xs) - exponent_bits(T))
1064 xs <<= m
1065 xu = xs | (xu & sign_mask(T))
1066 end
1067 xu = (xu & ~exponent_mask(T)) | exponent_one(T)
1068 return reinterpret(T, xu)
1069 end
1070
1071 """
1072 frexp(val)
1073
1074 Return `(x,exp)` such that `x` has a magnitude in the interval ``[1/2, 1)`` or 0,
1075 and `val` is equal to ``x \\times 2^{exp}``.
1076 # Examples
1077 ```jldoctest
1078 julia> frexp(12.8)
1079 (0.8, 4)
1080 ```
1081 """
1082 function frexp(x::T) where T<:IEEEFloat
1083 xu = reinterpret(Unsigned, x)
1084 xs = xu & ~sign_mask(T)
1085 xs >= exponent_mask(T) && return x, 0 # NaN or Inf
1086 k = Int(xs >> significand_bits(T))
1087 if k == 0 # x is subnormal
1088 xs == 0 && return x, 0 # +-0
1089 m = leading_zeros(xs) - exponent_bits(T)
1090 xs <<= unsigned(m)
1091 xu = xs | (xu & sign_mask(T))
1092 k = 1 - m
1093 end
1094 k -= (exponent_bias(T) - 1)
1095 xu = (xu & ~exponent_mask(T)) | exponent_half(T)
1096 return reinterpret(T, xu), k
1097 end
1098
1099 """
1100 $(@__MODULE__).scaleinv(x)
1101
1102 Compute `(scale, invscale)` where `scale` and `invscale` are non-subnormal
1103 (https://en.wikipedia.org/wiki/Subnormal_number) finite powers of two such that
1104 `scale * invscale == 1` and `scale` is roughly on the order of `abs(x)`.
1105 Inf, NaN, and zero inputs also result in finite nonzero outputs.
1106 These values are useful to scale computations to prevent overflow and underflow
1107 without round-off errors or division.
1108
1109 UNSTABLE DETAIL: For `x isa IEEEFLoat`, `scale` is chosen to be the
1110 `prevpow(2,abs(x))` when possible, but is never less than floatmin(x) or greater
1111 than inv(floatmin(x)). `Inf` and `NaN` resolve to `inv(floatmin(x))`. This
1112 behavior is subject to change.
1113
1114 # Examples
1115 ```jldoctest
1116 julia> $(@__MODULE__).scaleinv(7.5)
1117 (4.0, 0.25)
1118 ```
1119 """
1120 function scaleinv(x::T) where T<:IEEEFloat
1121 # by removing the sign and significand and restricting values to a limited range,
1122 # we can invert a number using bit-twiddling instead of division
1123 U = uinttype(T)
1124 umin = reinterpret(U, floatmin(T))
1125 umax = reinterpret(U, inv(floatmin(T)))
1126 emask = exponent_mask(T) # used to strip sign and significand
1127 u = clamp(reinterpret(U, x) & emask, umin, umax)
1128 scale = reinterpret(T, u)
1129 invscale = reinterpret(T, umin + umax - u) # inv(scale)
1130 return scale, invscale
1131 end
1132
1133 # NOTE: This `rem` method is adapted from the msun `remainder` and `remainderf`
1134 # functions, which are under the following license:
1135 #
1136 # Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
1137 #
1138 # Developed at SunSoft, a Sun Microsystems, Inc. business.
1139 # Permission to use, copy, modify, and distribute this
1140 # software is freely granted, provided that this notice
1141 # is preserved.
1142 function rem(x::T, p::T, ::RoundingMode{:Nearest}) where T<:IEEEFloat
1143 (iszero(p) || !isfinite(x) || isnan(p)) && return T(NaN)
1144 x == p && return copysign(zero(T), x)
1145 oldx = x
1146 x = abs(rem(x, 2p)) # 2p may overflow but that's okay
1147 p = abs(p)
1148 if p < 2 * floatmin(T) # Check whether dividing p by 2 will underflow
1149 if 2x > p
1150 x -= p
1151 if 2x >= p
1152 x -= p
1153 end
1154 end
1155 else
1156 p_half = p / 2
1157 if x > p_half
1158 x -= p
1159 if x >= p_half
1160 x -= p
1161 end
1162 end
1163 end
1164 return flipsign(x, oldx)
1165 end
1166
1167
1168 """
1169 modf(x)
1170
1171 Return a tuple `(fpart, ipart)` of the fractional and integral parts of a number. Both parts
1172 have the same sign as the argument.
1173
1174 # Examples
1175 ```jldoctest
1176 julia> modf(3.5)
1177 (0.5, 3.0)
1178
1179 julia> modf(-3.5)
1180 (-0.5, -3.0)
1181 ```
1182 """
1183 modf(x) = isinf(x) ? (flipsign(zero(x), x), x) : (rem(x, one(x)), trunc(x))
1184
1185 function modf(x::T) where T<:IEEEFloat
1186 isinf(x) && return (copysign(zero(T), x), x)
1187 ix = trunc(x)
1188 rx = copysign(x - ix, x)
1189 return (rx, ix)
1190 end
1191
1192 # @constprop aggressive to help the compiler see the switch between the integer and float
1193 # variants for callers with constant `y`
1194 @constprop :aggressive function ^(x::Float64, y::Float64)
1195 xu = reinterpret(UInt64, x)
1196 xu == reinterpret(UInt64, 1.0) && return 1.0
1197 # Exponents greater than this will always overflow or underflow.
1198 # Note that NaN can pass through this, but that will end up fine.
1199 if !(abs(y)<0x1.8p62)
1200 isnan(y) && return y
1201 y = sign(y)*0x1.8p62
1202 end
1203 yint = unsafe_trunc(Int64, y) # This is actually safe since julia freezes the result
1204 y == yint && return @noinline x^yint
1205 2*xu==0 && return abs(y)*Inf*(!(y>0)) # if x==0
1206 x<0 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer
1207 !isfinite(x) && return x*(y>0 || isnan(x)) # x is inf or NaN
1208 if xu < (UInt64(1)<<52) # x is subnormal
1209 xu = reinterpret(UInt64, x * 0x1p52) # normalize x
1210 xu &= ~sign_mask(Float64)
1211 xu -= UInt64(52) << 52 # mess with the exponent
1212 end
1213 return pow_body(xu, y)
1214 end
1215
1216 @inline function pow_body(xu::UInt64, y::Float64)
1217 logxhi,logxlo = Base.Math._log_ext(xu)
1218 xyhi, xylo = two_mul(logxhi,y)
1219 xylo = muladd(logxlo, y, xylo)
1220 hi = xyhi+xylo
1221 return Base.Math.exp_impl(hi, xylo-(hi-xyhi), Val(:ℯ))
1222 end
1223
1224 @constprop :aggressive function ^(x::T, y::T) where T <: Union{Float16, Float32}
1225 x == 1 && return one(T)
1226 # Exponents greater than this will always overflow or underflow.
1227 # Note that NaN can pass through this, but that will end up fine.
1228 max_exp = T == Float16 ? T(3<<14) : T(0x1.Ap30)
1229 if !(abs(y)<max_exp)
1230 isnan(y) && return y
1231 y = sign(y)*max_exp
1232 end
1233 yint = unsafe_trunc(Int32, y) # This is actually safe since julia freezes the result
1234 y == yint && return x^yint
1235 x < 0 && throw_exp_domainerror(x)
1236 !isfinite(x) && return x*(y>0 || isnan(x))
1237 x==0 && return abs(y)*T(Inf)*(!(y>0))
1238 return pow_body(x, y)
1239 end
1240
1241 @inline function pow_body(x::T, y::T) where T <: Union{Float16, Float32}
1242 return T(exp2(log2(abs(widen(x))) * y))
1243 end
1244
1245 # compensated power by squaring
1246 @constprop :aggressive @inline function ^(x::Float64, n::Integer)
1247 n == 0 && return one(x)
1248 return pow_body(x, n)
1249 end
1250
1251 @assume_effects :terminates_locally @noinline function pow_body(x::Float64, n::Integer)
1252 y = 1.0
1253 xnlo = ynlo = 0.0
1254 n == 3 && return x*x*x # keep compatibility with literal_pow
1255 if n < 0
1256 rx = inv(x)
1257 n==-2 && return rx*rx #keep compatibility with literal_pow
1258 isfinite(x) && (xnlo = -fma(x, rx, -1.) * rx)
1259 x = rx
1260 n = -n
1261 end
1262 while n > 1
1263 if n&1 > 0
1264 err = muladd(y, xnlo, x*ynlo)
1265 y, ynlo = two_mul(x,y)
1266 ynlo += err
1267 end
1268 err = x*2*xnlo
1269 x, xnlo = two_mul(x, x)
1270 xnlo += err
1271 n >>>= 1
1272 end
1273 err = muladd(y, xnlo, x*ynlo)
1274 return ifelse(isfinite(x) & isfinite(err), muladd(x, y, err), x*y)
1275 end
1276
1277 function ^(x::Float32, n::Integer)
1278 n == -2 && return (i=inv(x); i*i)
1279 n == 3 && return x*x*x #keep compatibility with literal_pow
1280 n < 0 && return Float32(Base.power_by_squaring(inv(Float64(x)),-n))
1281 Float32(Base.power_by_squaring(Float64(x),n))
1282 end
1283 @inline ^(x::Float16, y::Integer) = Float16(Float32(x) ^ y)
1284 @inline literal_pow(::typeof(^), x::Float16, ::Val{p}) where {p} = Float16(literal_pow(^,Float32(x),Val(p)))
1285
1286 ## rem2pi-related calculations ##
1287
1288 function add22condh(xh::Float64, xl::Float64, yh::Float64, yl::Float64)
1289 # This algorithm, due to Dekker, computes the sum of two
1290 # double-double numbers and returns the high double. References:
1291 # [1] http://www.digizeitschriften.de/en/dms/img/?PID=GDZPPN001170007
1292 # [2] https://doi.org/10.1007/BF01397083
1293 r = xh+yh
1294 s = (abs(xh) > abs(yh)) ? (xh-r+yh+yl+xl) : (yh-r+xh+xl+yl)
1295 zh = r+s
1296 return zh
1297 end
1298
1299 # multiples of pi/2, as double-double (ie with "tail")
1300 const pi1o2_h = 1.5707963267948966 # convert(Float64, pi * BigFloat(1/2))
1301 const pi1o2_l = 6.123233995736766e-17 # convert(Float64, pi * BigFloat(1/2) - pi1o2_h)
1302
1303 const pi2o2_h = 3.141592653589793 # convert(Float64, pi * BigFloat(1))
1304 const pi2o2_l = 1.2246467991473532e-16 # convert(Float64, pi * BigFloat(1) - pi2o2_h)
1305
1306 const pi3o2_h = 4.71238898038469 # convert(Float64, pi * BigFloat(3/2))
1307 const pi3o2_l = 1.8369701987210297e-16 # convert(Float64, pi * BigFloat(3/2) - pi3o2_h)
1308
1309 const pi4o2_h = 6.283185307179586 # convert(Float64, pi * BigFloat(2))
1310 const pi4o2_l = 2.4492935982947064e-16 # convert(Float64, pi * BigFloat(2) - pi4o2_h)
1311
1312 """
1313 rem2pi(x, r::RoundingMode)
1314
1315 Compute the remainder of `x` after integer division by `2π`, with the quotient rounded
1316 according to the rounding mode `r`. In other words, the quantity
1317
1318 x - 2π*round(x/(2π),r)
1319
1320 without any intermediate rounding. This internally uses a high precision approximation of
1321 2π, and so will give a more accurate result than `rem(x,2π,r)`
1322
1323 - if `r == RoundNearest`, then the result is in the interval ``[-π, π]``. This will generally
1324 be the most accurate result. See also [`RoundNearest`](@ref).
1325
1326 - if `r == RoundToZero`, then the result is in the interval ``[0, 2π]`` if `x` is positive,.
1327 or ``[-2π, 0]`` otherwise. See also [`RoundToZero`](@ref).
1328
1329 - if `r == RoundDown`, then the result is in the interval ``[0, 2π]``.
1330 See also [`RoundDown`](@ref).
1331 - if `r == RoundUp`, then the result is in the interval ``[-2π, 0]``.
1332 See also [`RoundUp`](@ref).
1333
1334 # Examples
1335 ```jldoctest
1336 julia> rem2pi(7pi/4, RoundNearest)
1337 -0.7853981633974485
1338
1339 julia> rem2pi(7pi/4, RoundDown)
1340 5.497787143782138
1341 ```
1342 """
1343 function rem2pi end
1344 function rem2pi(x::Float64, ::RoundingMode{:Nearest})
1345 isnan(x) && return x
1346 isinf(x) && return NaN
1347
1348 abs(x) < pi && return x
1349
1350 n,y = rem_pio2_kernel(x)
1351
1352 if iseven(n)
1353 if n & 2 == 2 # n % 4 == 2: add/subtract pi
1354 if y.hi <= 0
1355 return add22condh(y.hi,y.lo,pi2o2_h,pi2o2_l)
1356 else
1357 return add22condh(y.hi,y.lo,-pi2o2_h,-pi2o2_l)
1358 end
1359 else # n % 4 == 0: add 0
1360 return y.hi+y.lo
1361 end
1362 else
1363 if n & 2 == 2 # n % 4 == 3: subtract pi/2
1364 return add22condh(y.hi,y.lo,-pi1o2_h,-pi1o2_l)
1365 else # n % 4 == 1: add pi/2
1366 return add22condh(y.hi,y.lo,pi1o2_h,pi1o2_l)
1367 end
1368 end
1369 end
1370 function rem2pi(x::Float64, ::RoundingMode{:ToZero})
1371 isnan(x) && return x
1372 isinf(x) && return NaN
1373
1374 ax = abs(x)
1375 ax <= 2*Float64(pi,RoundDown) && return x
1376
1377 n,y = rem_pio2_kernel(ax)
1378
1379 if iseven(n)
1380 if n & 2 == 2 # n % 4 == 2: add pi
1381 z = add22condh(y.hi,y.lo,pi2o2_h,pi2o2_l)
1382 else # n % 4 == 0: add 0 or 2pi
1383 if y.hi > 0
1384 z = y.hi+y.lo
1385 else # negative: add 2pi
1386 z = add22condh(y.hi,y.lo,pi4o2_h,pi4o2_l)
1387 end
1388 end
1389 else
1390 if n & 2 == 2 # n % 4 == 3: add 3pi/2
1391 z = add22condh(y.hi,y.lo,pi3o2_h,pi3o2_l)
1392 else # n % 4 == 1: add pi/2
1393 z = add22condh(y.hi,y.lo,pi1o2_h,pi1o2_l)
1394 end
1395 end
1396 copysign(z,x)
1397 end
1398 function rem2pi(x::Float64, ::RoundingMode{:Down})
1399 isnan(x) && return x
1400 isinf(x) && return NaN
1401
1402 if x < pi4o2_h
1403 if x >= 0
1404 return x
1405 elseif x > -pi4o2_h
1406 return add22condh(x,0.0,pi4o2_h,pi4o2_l)
1407 end
1408 end
1409
1410 n,y = rem_pio2_kernel(x)
1411
1412 if iseven(n)
1413 if n & 2 == 2 # n % 4 == 2: add pi
1414 return add22condh(y.hi,y.lo,pi2o2_h,pi2o2_l)
1415 else # n % 4 == 0: add 0 or 2pi
1416 if y.hi > 0
1417 return y.hi+y.lo
1418 else # negative: add 2pi
1419 return add22condh(y.hi,y.lo,pi4o2_h,pi4o2_l)
1420 end
1421 end
1422 else
1423 if n & 2 == 2 # n % 4 == 3: add 3pi/2
1424 return add22condh(y.hi,y.lo,pi3o2_h,pi3o2_l)
1425 else # n % 4 == 1: add pi/2
1426 return add22condh(y.hi,y.lo,pi1o2_h,pi1o2_l)
1427 end
1428 end
1429 end
1430 function rem2pi(x::Float64, ::RoundingMode{:Up})
1431 isnan(x) && return x
1432 isinf(x) && return NaN
1433
1434 if x > -pi4o2_h
1435 if x <= 0
1436 return x
1437 elseif x < pi4o2_h
1438 return add22condh(x,0.0,-pi4o2_h,-pi4o2_l)
1439 end
1440 end
1441
1442 n,y = rem_pio2_kernel(x)
1443
1444 if iseven(n)
1445 if n & 2 == 2 # n % 4 == 2: sub pi
1446 return add22condh(y.hi,y.lo,-pi2o2_h,-pi2o2_l)
1447 else # n % 4 == 0: sub 0 or 2pi
1448 if y.hi < 0
1449 return y.hi+y.lo
1450 else # positive: sub 2pi
1451 return add22condh(y.hi,y.lo,-pi4o2_h,-pi4o2_l)
1452 end
1453 end
1454 else
1455 if n & 2 == 2 # n % 4 == 3: sub pi/2
1456 return add22condh(y.hi,y.lo,-pi1o2_h,-pi1o2_l)
1457 else # n % 4 == 1: sub 3pi/2
1458 return add22condh(y.hi,y.lo,-pi3o2_h,-pi3o2_l)
1459 end
1460 end
1461 end
1462
1463 rem2pi(x::Float32, r::RoundingMode) = Float32(rem2pi(Float64(x), r))
1464 rem2pi(x::Float16, r::RoundingMode) = Float16(rem2pi(Float64(x), r))
1465 rem2pi(x::Int32, r::RoundingMode) = rem2pi(Float64(x), r)
1466 function rem2pi(x::Int64, r::RoundingMode)
1467 fx = Float64(x)
1468 fx == x || throw(ArgumentError("Int64 argument to rem2pi is too large: $x"))
1469 rem2pi(fx, r)
1470 end
1471
1472 """
1473 mod2pi(x)
1474
1475 Modulus after division by `2π`, returning in the range ``[0,2π)``.
1476
1477 This function computes a floating point representation of the modulus after division by
1478 numerically exact `2π`, and is therefore not exactly the same as `mod(x,2π)`, which would
1479 compute the modulus of `x` relative to division by the floating-point number `2π`.
1480
1481 !!! note
1482 Depending on the format of the input value, the closest representable value to 2π may
1483 be less than 2π. For example, the expression `mod2pi(2π)` will not return `0`, because
1484 the intermediate value of `2*π` is a `Float64` and `2*Float64(π) < 2*big(π)`. See
1485 [`rem2pi`](@ref) for more refined control of this behavior.
1486
1487 # Examples
1488 ```jldoctest
1489 julia> mod2pi(9*pi/4)
1490 0.7853981633974481
1491 ```
1492 """
1493 mod2pi(x) = rem2pi(x,RoundDown)
1494
1495 # generic fallback; for number types, promotion.jl does promotion
1496
1497 """
1498 muladd(x, y, z)
1499
1500 Combined multiply-add: computes `x*y+z`, but allowing the add and multiply to be merged
1501 with each other or with surrounding operations for performance.
1502 For example, this may be implemented as an [`fma`](@ref) if the hardware supports it
1503 efficiently.
1504 The result can be different on different machines and can also be different on the same machine
1505 due to constant propagation or other optimizations.
1506 See [`fma`](@ref).
1507
1508 # Examples
1509 ```jldoctest
1510 julia> muladd(3, 2, 1)
1511 7
1512
1513 julia> 3 * 2 + 1
1514 7
1515 ```
1516 """
1517 muladd(x,y,z) = x*y+z
1518
1519
1520 # helper functions for Libm functionality
1521
1522 """
1523 highword(x)
1524
1525 Return the high word of `x` as a `UInt32`.
1526 """
1527 @inline highword(x::Float64) = highword(reinterpret(UInt64, x))
1528 @inline highword(x::UInt64) = (x >>> 32) % UInt32
1529 @inline highword(x::Float32) = reinterpret(UInt32, x)
1530
1531 @inline fromhighword(::Type{Float64}, u::UInt32) = reinterpret(Float64, UInt64(u) << 32)
1532 @inline fromhighword(::Type{Float32}, u::UInt32) = reinterpret(Float32, u)
1533
1534
1535 """
1536 poshighword(x)
1537
1538 Return positive part of the high word of `x` as a `UInt32`.
1539 """
1540 @inline poshighword(x::Float64) = poshighword(reinterpret(UInt64, x))
1541 @inline poshighword(x::UInt64) = highword(x) & 0x7fffffff
1542 @inline poshighword(x::Float32) = highword(x) & 0x7fffffff
1543
1544 # More special functions
1545 include("special/cbrt.jl")
1546 include("special/exp.jl")
1547 include("special/hyperbolic.jl")
1548 include("special/trig.jl")
1549 include("special/rem_pio2.jl")
1550 include("special/log.jl")
1551
1552
1553 # Float16 definitions
1554
1555 for func in (:sin,:cos,:tan,:asin,:acos,:atan,:cosh,:tanh,:asinh,:acosh,
1556 :atanh,:log,:log2,:log10,:sqrt,:fourthroot,:log1p)
1557 @eval begin
1558 $func(a::Float16) = Float16($func(Float32(a)))
1559 $func(a::ComplexF16) = ComplexF16($func(ComplexF32(a)))
1560 end
1561 end
1562
1563 for func in (:exp,:exp2,:exp10,:sinh)
1564 @eval $func(a::ComplexF16) = ComplexF16($func(ComplexF32(a)))
1565 end
1566
1567
1568 atan(a::Float16,b::Float16) = Float16(atan(Float32(a),Float32(b)))
1569 sincos(a::Float16) = Float16.(sincos(Float32(a)))
1570
1571 for f in (:sin, :cos, :tan, :asin, :atan, :acos,
1572 :sinh, :cosh, :tanh, :asinh, :acosh, :atanh,
1573 :exp, :exp2, :exp10, :expm1, :log, :log2, :log10, :log1p,
1574 :exponent, :sqrt, :cbrt, :sinpi, :cospi, :sincospi, :tanpi)
1575 @eval function ($f)(x::Real)
1576 xf = float(x)
1577 x === xf && throw(MethodError($f, (x,)))
1578 return ($f)(xf)
1579 end
1580 @eval $(f)(::Missing) = missing
1581 end
1582
1583 for f in (:atan, :hypot, :log)
1584 @eval $(f)(::Missing, ::Missing) = missing
1585 @eval $(f)(::Number, ::Missing) = missing
1586 @eval $(f)(::Missing, ::Number) = missing
1587 end
1588
1589 exp2(x::AbstractFloat) = 2^x
1590 exp10(x::AbstractFloat) = 10^x
1591 clamp(::Missing, lo, hi) = missing
1592 fourthroot(::Missing) = missing
1593
1594 end # module