# Best way to use/implement Debye function in Julia?
**URL:** https://discourse.julialang.org/t/best-way-to-use-implement-debye-function-in-julia/111090
**Category:** General Usage
**Tags:** numerics
**Created:** [March 3, 2024, 12:36pm UTC](https://discourse.julialang.org/t/best-way-to-use-implement-debye-function-in-julia/111090 "2024-03-03T12:36:06Z")
**Posts on this page:** 10
**Page:** 2
### Author:  [@trg818](https://discourse.julialang.org/u/trg818)
#### Post date: [March 5, 2024, 9:52am UTC](https://discourse.julialang.org/t/best-way-to-use-implement-debye-function-in-julia/111090/21 "2024-03-05T09:52:59Z")
Right, improving this one cut computation time by a factor 2 for the `x < pi` branch of the if-clause. The other branch was also sped up a bit by substituting the powers, but that improvement was rather minor. I did have to introduce a check whether to use `big` or not, because for certain x it was needed to achieve the desired accuracy.
Still, performance is a far cry from TOMS or quadrature. One thing I don’t understand is why the benchmark tool reports several “allocations” (29 in the first, 237 in the second branch of the if block), but none in the TOMS and quadrature functions. What is it that is being allocated here, and may this be the cause for the poor performance?
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [March 5, 2024, 1:35pm UTC](https://discourse.julialang.org/t/best-way-to-use-implement-debye-function-in-julia/111090/22 "2024-03-05T13:35:14Z")
> [@trg818](#):
>
> Still, performance is a far cry from TOMS or quadrature. One thing I don’t understand is why the benchmark tool reports several “allocations”
I haven’t looked at your code, but are you aware that most `BigFloat` operations allocate a new `BigFloat` object? If you want to avoid this, use MutableArithmetics.jl.
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [March 15, 2024, 9:57am UTC](https://discourse.julialang.org/t/best-way-to-use-implement-debye-function-in-julia/111090/23 "2024-03-15T09:57:49Z")
> [@trg818](#):
>
> What are the pros and cons of the different implementations?
The quadrature-based solution can be arbitrarily accurate, but gets slower the further you go from zero. `quadgk` from QuadGK.jl takes tolerance parameters which are used to set the desired accuracy.
An important thing that was missed in the discussion above was that one should use `expm1(x)` instead of `exp(x) - 1` in the integrand, as that’s more accurate for tiny `x` and allows `quadgk` to be faster.
Polynomial approximations may be useful as an even faster alternative to QuadGK.jl, depending on what interval you’re targeting and what’s your desired accuracy.
Could you give an estimate for an upper bound for your interval? The interval doesn’t necessarily need to include _all_ your inputs. If, say, 99% of your inputs were solved using a polynomial approximation, and the other 1% using quadrature, that could still provide a speed-up compared to using quadrature for _all_ inputs.
Also, what is your desired accuracy? Usually that’s expressed in terms of maximum allowed relative error, or, alternatively, as the required number of accurate bits in the significand of the floating-point number.
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [March 15, 2024, 10:13am UTC](https://discourse.julialang.org/t/best-way-to-use-implement-debye-function-in-julia/111090/24 "2024-03-15T10:13:04Z")
Here’s a quadrature-based solution for D\_1, D\_2, D\_3 and D\_4:
```julia
module DebyeQuadrature
import QuadGK, MutableArithmetics
pow(x, ::Val{1}) = x
pow(x, ::Val{2}) = x*x
pow(x, ::Val{3}) = x^3
pow(x, ::Val{4}) = x^4
function pow(x::BigFloat, ::Val{3})
o! = MutableArithmetics.operate!
y = x*x
o!(*, y, x)::BigFloat
end
function pow(x::BigFloat, ::Val{4})
o! = MutableArithmetics.operate!
y = x*x
o!(*, y, y)::BigFloat
end
function evaluate(
x, dord::Val{DOrd}, relative_tolerance, quadrature_order::Int,
) where {DOrd}
DOrd::Int
f = t -> pow(t, dord)/expm1(t)
if iszero(x)
one(x)
else
let z = false
i = QuadGK.quadgk(
f, z, x, maxevals = Inf, atol = z, rtol = relative_tolerance,
order = quadrature_order,
)
DOrd/pow(x, dord)*first(i)
end
end
end
end
```
The quadrature order affects performance, a small value like 7 should be enough for `Float64`, while `BigFloat` performance benefits from larger values like 35, or even in the hundreds. Larger values increase compilation time.
---
### Author:  [@trg818](https://discourse.julialang.org/u/trg818)
#### Post date: [March 15, 2024, 11:20am UTC](https://discourse.julialang.org/t/best-way-to-use-implement-debye-function-in-julia/111090/25 "2024-03-15T11:20:01Z")
For my purposes, I rarely expect values for x above 7; in most cases, they should barely exceed 2 and indeed often be smaller than 1. This is because my applications do not generally target cryogenic temperatures but rather the intermediate-to-high-T realm.
---
### Author:  [@Eben60](https://discourse.julialang.org/u/Eben60)
#### Post date: [March 15, 2024, 2:41pm UTC](https://discourse.julialang.org/t/best-way-to-use-implement-debye-function-in-julia/111090/26 "2024-03-15T14:41:51Z")
> [@nsajko](#):
>
> The quadrature-based solution can be arbitrarily accurate, but gets slower the further you go from zero.
What about pre-computing sufficiently much intermeadiate points to maximal desired accuracy, and integrate from the nearest one?
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [March 15, 2024, 2:43pm UTC](https://discourse.julialang.org/t/best-way-to-use-implement-debye-function-in-julia/111090/27 "2024-03-15T14:43:49Z")
generally you can get better results just by finding the series expansion around infinity and using that for large values.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [March 15, 2024, 5:06pm UTC](https://discourse.julialang.org/t/best-way-to-use-implement-debye-function-in-julia/111090/28 "2024-03-15T17:06:25Z")
Yes, for a given fixed accuracy (e.g. `Float64`), typically high-performance implementations of special functions stitch together a set of different approximations — Taylor series near zeros, minimax polynomials or rational approximants in finite real intervals, continued-fractions and/or asymptotic series for large arguments. It takes a lot of tuning to find the right set of approximations and to find good crossover points to switch between approximations.
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [March 15, 2024, 11:51pm UTC](https://discourse.julialang.org/t/best-way-to-use-implement-debye-function-in-julia/111090/29 "2024-03-15T23:51:25Z")
A Julia package for evaluating Debye functions is being registered, it will live at DebyeFunctions.jl once it’s registered:
>
>
> \- Registering package: DebyeFunctions
> \- Repository: https://gitlab.com/nsajko/De…byeFunctions.jl
> \- Created by: https://gitlab.com/nsajko
> \- Version: v0.0.1
> \- Commit: 215b3bc05c752486fc0b1a17567692a50ba3230e
> \- Git reference: HEAD
> \- Description: Evaluate Debye functions
Git repo at Gitlab: [Neven Sajko / DebyeFunctions.jl · GitLab](https://gitlab.com/nsajko/DebyeFunctions.jl)
Please test it out and maybe suggest some improvements.
EDIT: well, an obvious and easy improvement would be to add a few more polynomial approximations (for each Debye function order). The trouble is choosing when to stop adding more polynomials, I guess.
---
### Author:  [@longemen3000](https://discourse.julialang.org/u/longemen3000)
#### Post date: [November 18, 2024, 8:56pm UTC](https://discourse.julialang.org/t/best-way-to-use-implement-debye-function-in-julia/111090/30 "2024-11-18T20:56:46Z")
for orders 1,2,3, there are specific formulations that are defined in terms of the polylogarithm and the riemann’s zeta function:
```julia
import PolyLog,SpecialFunctions,LogExpFunctions
function debye3(x)
N = 3
ζ = SpecialFunctions.zeta(N+1)
c1 = LogExpFunctions.log1mexp(x) #log(1-exp(x))
c2 = 3*PolyLog.reli2(ex)
c3 = -6*PolyLog.reli3(ex)
c4 = 6*PolyLog.reli4(ex)
ex = exp(x)
return x*x*x*(-0.25*x + c1) + c2 + c3 + c4 - 6*ζ
end
```
(from [https://mathworld.wolfram.com/DebyeFunctions.html](https://mathworld.wolfram.com/DebyeFunctions.html))
[Previous page](https://discourse.julialang.org/t/best-way-to-use-implement-debye-function-in-julia/111090.md?page=1)