Chebyshev Approximation
embeddedrelated.com
embeddedrelated.com
He has nice sets of coefficients derived by taking into account the discrete values of floats/doubles, and the performance gain by setting some coefficients to 0, 1 or 2.
The only thing he did not account for is precision loss when multiplying and adding values when actually calculating the polynomial, which unfortunately results in a large relative error around log(1). Eigen's vectorized implementation does some weird tricks with changing the range to sqrt(0.5) to sqrt(2), and does not exhibit this large relatively error. The unvectorized MSVC implementation does not exhibit this either.
Does anyone how Microsoft and Intel and the like derive the coefficients for their approximations of the elementary functions in math.h?
The standard Lagrange interpolation formula implemented in that blog post works too, but is O(n^2) vs. O(n) for each evaluation, and not numerically stable.
Also see http://www.maths.manchester.ac.uk/~higham/narep/narep440.pdf https://people.maths.ox.ac.uk/trefethen/publication/PDF/2011...
Inre the second blog post link, it’s generally a bad idea to evaluate polynomials on an interval in terms of the monomial basis.
But anyhow, on the subject of the Remez algorithm, see http://www.chebfun.org/publications/remez.pdf http://www.chebfun.org/examples/approx/BestApprox.html
Note that while a standard Chebyshev polynomial approximant isn’t quite “optimal” for any particular function, in an L∞ norm minimax sense, it’s much easier to find, and often better for most values of x. See Trefethen’s book for details, and also the Chebfun guide, http://www.chebfun.org/docs/guide/guide04.html or this paper has a nice example picture under “myth 3” https://people.maths.ox.ac.uk/trefethen/mythspaper.pdf
When I am trying to teach someone about algorithms, I turn to Python because it is more readable, quicker for me to write, it's platform-independent, and I can incorporate it into an IPython (now Jupyter) Notebook and postprocess that into my blog articles with a minimum of hiccups.
http://nim.community/answers/create_a_shared_library/
Although I haven't figured out the garbage collection bit for the shared object bit. It is a lovely language and it's reaching 1.0 sometime this year.
Unfortunately just because a language exists, doesn't mean it is usable on a particular embedded system architecture. The processor core and system libraries need support.
I avoided putting anything in linker scripts that didn't have to be there (IIRC, memory blocks and VTOR (for assembly startup routine that can handle flash or RAM)). Peripheral addresses stayed in Rust so that the optimizer could see them (eg writing to RA and then RB, the high u16 stays the same).
BTW F103C8 "dev boards" are available from China for $5.
EDIT: nevermind, I realize I misunderstood; you're referring to the optimizer being able to see that both reg addresses have the same high 16 bits. I thought you were referring to the actual interactions with the register contents. Yea, that makes sense.
EDIT 2: for what it's worth, i did a quick test just now and found that moving the addresses into the source actually made the generated code significantly worse. I'm just tinkering around on a lunch break and I need to get back to real work, so I can't investigate too deeply, but the 2 biggest factors I see are that (1) it simply doesn't perform the optimization you describe, and (2) it makes things even worse because rather than addressing relative to the peripheral base address it loads a complete new constant address for every operation.
(And if one isn't atomically writing to control registers when interrupts are enabled, you are most likely building on a foundation of nondeterministic bugs. Rust's linear types kind of inform you of this, but to do anything useful the temptation is to unsafely share across threads)
re edit 2: Really? I'd obviously looked at my generated code and saw the expected improvement, so we must be doing something different.
It's also possible something has changed in rustc since you were playing with it.
My basic register write is volatile_store(($base as * mut u8).offset(idx as isize) as * mut u32, v); $base is just a hex literal address. There is an extra cast because I defined 'idx' in terms of bytes.
I might not have been concerned about the absolute smallest code when I accessed $base+4 then $base+8 (use of a small offset vs loading a second constant), but I definitely observed eliminating redundant loading of the high half between different registers.
I was compiling with "-O -Z no-landing-pads -C lto".
It's sad to hear that they've dug their heels in even further and are now one of the few examples of bad-faith attacking the GPL through contracts. Especially because recently reading their ENC28J60 documentation was a nice break after wading through STM/ARM's stuff.
If I'm ever again in a position where I'm choosing chips for a production design, I'll be verifying what you said for myself (diligence) and then steering clear to less risky alternatives.
This is a quite serious concrete accusation. It would put Microchip in the company of Tivo and Sveasoft - not simply ignoring the GPL, but actively attacking it through bad-faith contractual restrictions.
But yes, I would definitely have to see the specific license for myself before forming any conclusions. HN is essentially just gossip, and GP could have very well interpreted a clause too broadly. But I'm not interested in sorting through what's new in the Microchip world just to preemptively answer this, especially because he could be referring to a single specific library that would be hard to find.
Microchip licenses to you the right to use,
modify, copy and distribute Software only
when embedded on a Microchip microcontroller
or digital signal controller that is integrated
into your product or third party product
(pursuant to the sublicense terms in the
accompanying license agreement).
Microchip has a number of open-source initiatives. You can use the internal guts of MPLAB X from Java via MDBcore, for example, as a scriptable debugger/simulator.
Um. Been there, done that. Well, you see, [: discussion of internal politics deleted :]
Details about a particular variable's type, or a languages semantics don't tell me anything about (in this case) chebyshev or taylor functions, but to follow the code I'd have to internalize the interactions of those things, along with the actual math for it all to make sense.
Consider it an abstraction layer for learning. Just like when dealing with files on disk I use open(), read(), etc rather than work in my fs + hdd controller.
Doing the same in C is going to be a lot more effort.
Also at Lloyd Trefethen’s book Approximation Theory and Approximation Practice, of which the first 6 chapters are available online, http://www.chebfun.org/ATAP/
The author of this blog post might as well have used Clenshaw’s algorithm as he mentioned was suggested by Numerical Recipes. It’s like 4–5 lines of code, fast, and has been proven numerically stable.
If bk1, and bk2 are the previous two terms of the recurrence, and c the coefficient array, code to find y = f(x) for x in [-1, 1] looks like:
x2 = 2*x
bk2 = bk1 = 0
for ck in c[:0:-1]:
bk2, bk1 = bk1, ck + x2*bk1 - bk2
y = c[0] + x * bk1 - bk2
Or if he wanted an easy Python tool, just used this (partial) python port of Chebfun (though I guess it dates from around the same time as his post), https://github.com/olivierverdier/pychebfun
I’m not convinced about his measurements on an arbitrary grid -> chebyshev coefficients via weighted least squares method. If his measurement points are reasonably clustered near the endpoints of his interval, he could just use an exactly interpolating polynomial and use the barycentric interpolation formula to get the values at the chebyshev points (see Trefethen’s book). His approach might also be fine, I haven’t thought deeply about it. Also see http://www.chebfun.org/examples/approx/RationalInterp.html http://www.chebfun.org/examples/approx/EquispacedData.html http://www.chebfun.org/examples/approx/BestL2Approximation.h...
The neat thing about Chebyshev polynomials is that it’s possible to convert from values at the Chebyshev points (either extrema or roots of the Chebyshev polynomial) to Chebyshev coefficients by using an FFT, which takes O(n log n) time. Along with various other efficient algorithms for differentiation, integration, root finding, solving differential equations, etc., Chebfun makes it practical to deal with polynomials of very high degree, and do lots of neat stuff with them.
This blog post’s statement that “higher degree polynomials can have problems with numerical accuracy” turns out to be a myth: https://people.maths.ox.ac.uk/trefethen/mythspaper.pdf
(author here) I should have qualified or corrected that statement in my article. There are a few reasons I rarely use n>5:
- when polynomial approximation is appropriate, high-degree polynomials are usually unnecessary on embedded systems (0.1% accuracy usually sufficient)
- it's usually a sign that polynomial approximation is being misapplied and there are better ways (via range reduction, transformation like 1/f(x), etc.)
- with least-squares, computation of the coefficients can have problems with numerical accuracy. When I've used least-squares fitting polynomials with n>10, I often get a complaint from my numerical tool (MATLAB or Python) that the matrix in question is ill-conditioned. That may be because I used the Vandermonde matrix rather than Chebyshev polynomials, or because I used a large number of points. Chebyshev approximation doesn't have the same problem, I don't think, and I should not have mapped my experience with least-squares fitting onto Chebyshev approximation.
Evaluation of the polynomials themselves should be fine (after all, we can easily compute a 20th-degree Chebyshev polynomial).
So is that essentially sort of like Horner's Rule but for Chebyshev polynomials? e.g. instead of a0 + (a1 * x) + (a2 * x * x) + (a3 * x * x * x) + ... , use a0 + x * (a1 + x * (a2 + x * (a3 + ... so that the coefficients are used in the reverse order, and instead of having to maintain a sum and a value of x^n, we just have a single accumulator; in Clenshaw's algorithm we have the two accumulators bk1 and bk2, but they look similar.
Such a workload is much more common than a case where you need to evaluate a polynomial of degree one million at one single point.
Note that you can use SIMD, multi-core, GPUs, etc. for your multiple point evaluations, and these algorithms can be turned into mostly FMA instructions.
EDIT, moving this comment inline: To see several suggested improvements to Horner’s rule, take a look at Knuth’s TAOCP volume 2, §4.6.4, which is (I assume) a decent summary of the state of the art in the late 90s, though there’s surely more literature since.
Anyway though, we’re mostly talking about evaluating Chebyshev polynomials here. If you know some papers about speeding up Clenshaw’s algorithm, I’d love to look them up. (I’m currently trying to figure out how fast I can make such code run on a GPU/CPU, for my own use projecting high-resolution raster satellite images onto different map projections, so speedups would be welcome.)
EDIT @jacobolus does it cover Estrin's ? Its certainly older than the 90s. The main problem is Horners does not have any parallelism. With Estrin's one can oarallelize along two axes: Use threads for one and simd for another. Although I believe one can play with the parse tree to get better speed than Estrins on your specific hardware.
Knuth, page 488:
Another attractive method for parallel computation has been suggested by G. Estrin [Proc. Western Joint Computing Conf. 17 (1960), 33–44]; for n = 7, Estrin’s method is:
Processor 1: Processor 2: Processor 3:
a1 = u7*x + u6 b1 = u5*x + u4 c1 = u3*x + u2
a2 = a1*x^2 + b1 c2 = c1*x^2 + d1
a3 = a2*x^4 + c2
Processor 4: Processor 5:
d1 = u1*x + u0 x^2
x^4
* * *
Personally this seems like a lot of implementation trouble and communication overhead, and I suspect it will lose out badly to a GPU cranking out evaluations at separate points.
I have this learned fear of a lot of operations on polynomials, like root-finding, that are based on using the wrong set of basis functions (i.e., the monomial powers x^n), and the linked paper helped me to put this superstition in context.
[1] https://github.com/daniel-levin/numerical-methods-3/blob/mas...
[1] http://siber.cankaya.edu.tr/NumericalComputations/Fall2004/c...
Somehow my children cannot grow up as mathematically stunted as I am.
- It's not as bad as that. Chebyshev approximation + other techniques for evaluating mathematical functions are pretty easy to learn.
- You're absolutely right. Welcome to mathematics! The mountain goes infinitely high. There are tons of technical papers I'll never be able to digest simply because they've compressed lots of thoughts into cryptic-looking but standard notation.
Suggest group theory or linear algebra instead. ;-)