Faster `voigt`, `exp`
- Langage dominant
- Julia
- Étoiles
- 67
- Forks
- 15
- Merge moyen
- 2 j 6 h
- PR mergées (30 j)
- 6
Description
As promised, I'm posting my faster exponential function. I'm curious to hear if using this inside of `voigt` buys you any speed.
Tonight, I also put together a branchless version of `voigt`. Unfortunately, all of the extra work that is being done by the branchless voigt makes it slower than the unvectorized version. I suspect that if someone spent a lot of time refactoring, they might be able to get it to be slightly faster. But I don't think that's worth the time. As you had suggested a while back, I think a better use of time would be to replace the current approximation with a better approach.
I did some googling and the consensus seems to be that we should leverage the relationship of the Voigt function with the real function. According to [eqn 9 of Thompson (93)](https://aip.scitation.org/doi/abs/10.1063/1.4823236), `V(u,a) = Re[exp(z²) erfc(z)]`, where `z = a + i*u`. I suspect that nobody will have time to pursue this for a while, so I'll briefly document what I read. I've seen 2 ways to do this (there are probably more):
1. I think that [Turbospectrum's implementation](https://github.com/bertrandplez/Turbospectrum2019/blob/master/source-v19.1/newvoigt.f) is based on the same underlying approach as [Kuntz (1997)](https://www.sciencedirect.com/science/article/pii/S0022407396001628). (Interestingly, I think that turbospectrum and the article might be suggesting that for `V(a,v)`, `z = a - iv`...). Under this approach, I think that they directly approximate the whole expression `Re[exp(z²) erfc(z)]`, but it introduces at least some level of branching. I think this article talks a little about [variants of this approach](https://elib.dlr.de/120657/1/FSchreier_Voigt_1806.11560.pdf).
2. At a glance, I suspect this might potentially be slightly faster to adopt a strategy like the one used in [this article](https://opensky.ucar.edu/islandora/object/conference%3A828/datastream/PDF/view) (if we target a particular accuracy). They seemed to use an approximation for `exp(z²) erfc(z) ≈ (Σⱼᵖ aⱼ zʲ)/(z⁽ᵖ⁺¹⁾ + Σⱼᵖ bⱼ zʲ)` with coefficients from [here](https://ui.adsabs.harvard.edu/abs/1978JQSRT..19..509H/abstract), and directly computed `exp(z²)`. The benefit of this approach is that there is basically no-branching. However, plugging `Re(exp(-1*(x + i*y)^2) * (a + i b))` into wolfram alpha suggests that you also need to evaluate sin and cos (which I think boils down to 2 more polynomial evaluations)...
If either of these approaches are adopted, I believe some care is required to avoid using Julia's complex datatype (since I think ForwardDiff.jl is incompatible with complex numbers)
Guide de contribution
Aucun guide de contribution indexé pour ce dépôt
Évaluation
Cette issue n'a pas encore été évaluée.