libera/#sicl - IRC Chatlog
Search
13:41:45
kingcons
beach and/or prov: saw in logs that you're working on elementary functions (sin/cos/etc) and thought this might be of interest if you haven't seen it - https://blog.sigplan.org/2022/04/28/one-polynomial-approximation-to-produce-correctly-rounded-results-for-multiple-representations-and-rounding-modes/
16:19:04
beach
I came up with an idea for computing SIN and COS that I would like to try. It is very experimental, and the code has not been optimized. I might try to optimize it for SBCL at some point, just to see whether the technique is competitive with what SBCL does. For now, I just verify that the result is reasonably accurate.
16:19:10
beach
Here is the code with explanatory comments: https://github.com/robert-strandh/SICL/blob/master/Code/Arithmetic/sin-cos.lisp
16:19:11
beach
In half an hour or so, my (admittedly small) family will announce that dinner is ready, but before that I have a little time to discuss it. But feel free to examine it in greater detail, and give some thought to how well it might perform when optimized.
16:22:14
mfiano
Still analyzing it but my first thought is the runtime lookup of the special variables seems like it'd be easy to make a read or load-time table.
16:23:13
beach
Nothing has been optimized in fact. I am just counting arithmetic operations, and looking at accuracy.
16:23:42
mfiano
Also type stability may be a problem for efficiency, but there isn't much you can do about that
16:24:27
beach
I don't see that it would be more problematic here than in other methods. Am I wrong?
16:24:42
beach
The tables have 100 elements each, but it might be advantageous to have a lot more if that means decreasing the degree of the approximations. Caches are pretty big these days.
16:25:23
mfiano
Oh I just meant that if the return value's type has to be coerced from input types, that could be a severe price to pay.
16:30:34
beach
Well, I can't change the world, so I am mainly interested in comparisons between existing methods.
16:32:42
mfiano
I would have to study the math in detail to give any worthwhile advice, and I am much too sleepy for that.
16:33:43
Bike
bit of a typo in the comment on line 9. don't think i'm enough of a numerical analyst to comment on the actual algorithm, except i think the usual CORDIC method also works on some kind of angle sum trigonometric identity
16:34:47
Bike
sbcl i think calls out to libc, and glibc sin is... pretty involved, now that i'm looking at it
16:34:49
mfiano
If I was writing a CL compiler, I would want to look into standard functions being implemented in terms of generic functions with a JITer to generate efficient type combination functions at runtime as needed. That has been my dream for a long time, but I am too dumb and have too much on my plate to even see if it holds water.
16:34:58
beach
I briefly scanned the description of CORDIC and it looked strange to me. But I may be wrong.
16:35:21
mfiano
not necessariloy generic functions in terms of the standard, but a type-based dispatch variant
16:37:47
Bike
"This code uses some numerical hacks I've never seen before, though for all I know they might be well-known among floating-point experts. [...] The way this is done without division or branching is rather clever. But there's no comment at all!" good ol glibc.
16:38:37
Bike
no, it's on a stackoverflow answer to someone asking what glibc does, which is more informative than glibc itself
16:39:56
Bike
you could also just use the process instruction on x86, but i guess it's not very good
16:40:41
beach
Right, we had a long-ish discussion about that some time ago, and the conclusion is that one should not use such instructions, and they don't exist in most processors.
16:41:53
Bike
oh, see. i've been skimming through the math discussions here but i guess i missed that.
16:47:05
beach
My (admittedly small) family just announced that dinner is served. I'll be back tomorrow.
19:21:21
moon-child
mfiano: my understanding of 'type stability' as introduced by julia is that there is a strong correlation between a given function's parameter types and its result type. Which would seem to be the case here. No?
19:53:18
mfiano
In a more specific example than what was given, ASIN or ACOS can result in a different return type than either the inputs or the inferred types of intermediary calculations.
19:54:37
mfiano
The fact that there is only one function for the compiler to optimize, instead of several for different types, or using something like SBCL's defknown/deftransform introduces a problem
20:03:36
mfiano
If the compiler could infer the return type from the parametric values, it could dispatch at compile time to a more specialized function. Of course, this is not always possible with CL, and is subject to debate on whether it should be punted to the user to write such a compiler macro. The case I am talking about here is when they can return a compound object (of type COMPLEX). But this idea also
20:03:39
mfiano
extends to more specialized functions for the different floating point formats. But again, this is up to debate, for reasons of numerical stability for one, since it is often a good idea to perform intermediary calculations in a wider bit width like double-float for the added IEEE rounding bits.
20:05:25
moon-child
this is where dynamic analysis beats static analysis, imo. You may notice that asin _usually_ returns a non-complex number when given a non-complex number
20:06:17
moon-child
(doesn't have to be dynamic, even. Add a branch hint and you're off to the races)
20:09:18
mfiano
I overthink these things a lot. I have ran into issues where I had to introduce a bit array to get the rounding precision I needed for many chained double-float operations. Numerical stability compounds pretty fast in game development with lots of chained math.
20:09:54
mfiano
Implementing a subset IEEE whenever I need it is not my idea of fun, but sometimes I have to.
20:11:40
mfiano
Some older PC's atually did IEEE754 double calculations in 128bit FPU registers for the added rounding bits. I'm not sure if that is still the case today.
20:16:49
Bike
the glibc implementation of sine actually does like four different things depending on the magnitude of the argument
20:17:15
mfiano
Anyway, you can probably just ignore me. And yeah, I've been using Julia on and off for 7 years now, ever since Tomas Papp switched to it from CL.
20:20:33
mfiano
But maybe I did subconsciously? I did start a little Julia project about a week ago.
20:25:09
moon-child
Bike: it seems to me that it 'does' one thing, and just has a few minor variations on how it does range reduction
20:25:24
moon-child
(and also has a couple of paths where it doesn't do anything: very small, and very large)
20:27:06
moon-child
I am also a bit sceptical of code which relies heavily on lookup tables. They are difficult to benchmark; they may perform very well when you call them in a tight loop, but in a larger application, there are frequently knock-on effects and you can end up thrashing cache
20:27:30
Bike
all i'm thinking is that with a range analysis in the compiler that kind of thing could be avoided sometimes.
20:27:39
moon-child
(which is not to say they are universally bad, only a bit suspicious and to be treated with care)
20:28:13
moon-child
Bike: sounds like the sort of thing partial inlining would take care of. But I'm not sure how common it is to have bounds on floating-point numbers
20:29:06
Bike
sometimes they could be derived, like (sin (the real x)) is going to be between -1 and 1. though that's probably not going to be fed to sin again