I hope I didn’t make it too complex. It can be as simple or complex as you want. I still maintain that Julia is the simplest language to get code right (for science at least; Clojure has good properties that most Julia programmers to not emulate, they could kill performance).
It’s hard in all languages to get the most speed, including in Julia, and type instability is NOT a bug per se, while can be considered a performance bug, even lowering it to under Python’s so I find it very important that Julia’s defaults are good, e.g. in this case, or easy to “fix” type instability and debug/track down. How did you find it, were tools good enough for that?
My emphasis was wrong with your code, sorry about that, I was agreeing with you! That the type stability was there, but also considering a possible “usability bug” in Julia.
Now to make this even more complex. You opened quite a can of worms, and I thank you for opening this interesting thread. The division trick I knew of is maybe the most obscure Julia trick I know of.
Even I had to reread stuff, including mine, I was thinking zero(r) that’s what this is for, but I forgot the integer (and Rational) type-instability (I think stable for all other (float) types). I thought actually I would get away with 0 (but not 0.0):
func(r) = r == 0 ? zero(r) : r/3
@code_native func(1.0) # it's actually slightly larger (still type-stable I think), because of cast I thought would be optimized away.
func(r) = r == Bool(0) ? zero(r) : r/3 # this seems as good but back into obscurity.
Now I realize this is only an MWE, and the check not needed, the more interesting case is 1/r or even a/b:
func(r) = r/3 # simplest version, but you can absolutely not do r/3.0, if you care about rationals (or Float32 or lower).
How do you extend that type-stability /1 trick to a/b? Note, a/b does NOT guarantee Float64 result, it.s. e.g. Float32 if both are, but you will get the larger type of the two (a question when float and e.g. Dec64 or Posit), or Rationals when applies, then the larger type of the two, and I don’t want to go there if num and dem are not the same type (not possible in Julia’s type, but possible in theory, I’ve been thinking of making my own (vectorized) rational type, and that’s one thing I’m considering…).
You want and trust that the /1 goes away in optimization, and it seems to, i.e. this sufficient:
func(a, b) = iszero(a) ? zero(a)/one(b) : a/b
@code_native func(1.0, 1) # only a screen-full, still much better than 3 for div(1.0, 1)
func2(a, b) = iszero(a) ? zero(a)/one(b) : div(a, b)
@code_native func2(1.0, 1) # is 5 screen-fulls of assembly...
At the assembly level you get a floating-point exception on most CPUs if you divide by zero (I don’t know of any that do not for floating-point).
First for integer divide:
@code_native div(1, 0) # I do not just get the slow idiv, but 20+-instruction overhead, for ijl_throw
The div-by-0 avoidance above was very costly. Some CPUs define idiv equivalent to return 0 for it. x86 does not, it has an interrupt, and I suppose it crashes your code, note that interrupt (“exception”) isn’t the same as a Julia exception, but can it be caught and handled by Julia? I.e. turn it into a Julia exception? Doing what Julia does natively, adding it’s own check, to prevent the interrupt ever happening, might be better than the user code doing it. And then catching that Julia exception and checking after the fact.
For a float div, you get a +/- Inf (or Inf32 etc.), but be careful or NaN for 0/0, so do not check for only one case after the fact for floting point. You CAN check for all cases with isfinite(0/0) and substitute your 0.0 (of right type…) then. But note, say this was exp (or part of your expression), and you can expect Inf, then do not substitute blindly. Division has high latency, so you may want to do it first, not a check first, (unpredictable) branches are also costly on top of that. If division-by-zero is likely rare, it might be worth it to check after.
You see I managed to made it even more complex, than just about type-stability…
If that somefunction is sin, i.e. you’re trying to implement sinc, then, intentionally showing the wrong version:
my_naive_sinc(r) = (r == 0) ? one(r) : sin(r)/r # sinc in base is the other definition
In such an example you can’t simply return r, and be careful not to return r + 1, either. For type stability you might think r works if r + one(r), and then your back to just one(r), so rather useless, but neither that is type-stable!
That’s even though I see one used in Julia Base:
_sinc(x::Number) = iszero(x) ? one(x) : isinf_real(x) ? zero(x) : sinpi(x)/(pi*x)
Note there the integer case is handled separately. For my one-liner it needs to be:
my_naive_sinc(r) = (r == zero(r)) ? one(r)/1 : sin(r)/r # sinc in base is the other definition
So the lesson is avoid literal constants for type stability, but sadly, seeming one function (also zero I guess) is NOT a guarantee for success. It most often is enough though… And when it does seem to apply to have constants, e.g. for E=mc^2 that’s a separate (overflow-causing) issue.
The good thing is that perfect is the enemy of the good. Try to be type-stable, even if you only do it for floats, then it might be enough for you. In you non-package code you don’t even need to be fully generic, or even in part, just work with e.g. Float64. But in packages it’s best to be fully generic, if you fail you’re just one bug report away from that. You still don’t get incorrect results. The power ^ issue is an exception, one I’ve complained about and want to fix, get Float64 always… make it harder to opt into integer result.
https://www.felixcloutier.com/x86/fdiv:fdivp:fidiv
https://groups.google.com/g/llvm-dev/c/1FOyTcEJ3FI
Per C and C++, integer division by 0 is undefined. That means, if it happens, the compiler is free to do whatever it wants. It is perfectly legal for LLVM to define r to be, say, 42 in this code; it is not required to preserve the fact that the idiv instruction on x86 andx86-64 will trap.
Note, while those languages are standardized, and don’t handle in the standard (I think most compiler implementations do something sensible), and Julia is NOT standardized, it’s till guarantees throwing its Julia exception, for integers only, otherwise following float-standard (maybe it should say why?):
julia> div(0,0)
ERROR: DivideError: integer division error
- […] The Virtualization Extensions introduce the requirement for an ARMv7-A implementation to include
SDIVandUDIV.- The ARMv7-M profile also includes the
SDIVandUDIVinstructions.In the ARMv7-R profile, the SCTLR.DZ bit enables divide by zero fault detection:
SCTLR.DZ == 0
Divide-by-zero returns a zero result.
SCTLR.DZ == 1
SDIVandUDIVgenerate an Undefined Instruction exception on a divide-by-zero.
Julia supports ARMv7-A not ARMv7-R profile I think, but I’m just showing an example where you can get 0 out. I don’t know what the default flag is (for ARMv7-A), likely “divide by zero fault detection” enabled. I would like to know the best way to get single-instruction divide in machine code, and no Julia boilerplate, also for x86:
https://pdos.csail.mit.edu/6.828/2008/readings/i386/IDIV.htm
Interrupt 0 if the quotient is too large to fit in the designated register (AL or AX), or if the divisor is 0