
Math library functions that seem unnecessary - fogus
http://www.johndcook.com/blog/2010/06/07/math-library-functions-that-seem-unnecessary/
======
RiderOfGiraffes
When computing various values related to navigation you sometimes need
_1-cos(a)_ for very small values of _a._ If _a_ is very close to zero then
_cos(a)_ is very close to 1, and as _a_ varies (very sensibly) over modest
ranges, the value of _1-cos(a)_ as computed by the naive expression bounces
all over the place.

So one needs to detect when _a_ is small and use _1-cos(a) ~ a^2/2-a^4/24_

It makes a difference in some safety critical software I've written, and was
removed _despite the comment_ by a subsequent maintenance programmer. I was
then contracted to fix the major problem that ensued.

A similar problem is when using the naive solution to the quadractic equation.
If you have ax^2+bx+c=0 and you use x=(-b+/-sqrt(b^2-4ac))/(2a) then you run
into stability problems if ac is very small and "a" is small. Then the sqrt is
close to b, so for one choice on the numerator you get the subtraction of
nearly equal values, giving rubbish. Then you divide by a small number and get
big rubbish.

The solution in this case is to get the root where the numerator is close to
2b, and derive the second root from the first. More details if you want, but
after seeing why it's hard, it's no longer interesting.

~~~
marchdown
As far as I understand, mathematically von Neumann's computers simply use some
specific algebra that is decidedly different from naïve arithmetic. It is
actually very useful in certain applications, such as cryptography, that rely
on its algebraic properties. It may or may not be good approximation for some
numeric problems, but it is generally faster then working out elaborate
mathematical structures.

But I'm only an ignorant student, and I wonder if these solutions are really
significantly more effective than implementing the number tower or switching
to bignums? What are the best practices and the consensus among professionals
that have to deal with inadequacies of floating-point number in real-world
numeric applications?

~~~
lmkg
First of all, the problems with fixnums vs bignums are _very_ different from
the problems with floating point precision. Fixnums are a slightly different
mathematical structure from bignums, but they're both commutative groups, and
are identical for most inputs. Floating point arithmetic can basically be
treated as non-deterministic, which is an entirely different kettle of fish
(RoG can probably elaborate better than me how this characterization is and is
not valid).

If you're using a naive algorithm for numeric computing, you can get some
really atrocious numerical stability issues. If you happen to have exact
representations of your numbers, precision isn't an issue, but it's not common
that exact representations are even possible. A lot of common numerical
operations result in irrationals, including all root, logarithmic,
exponential, and trigonometric functions, and that would take an infinite
amount of space to represent faithfully. So really, you're always stuck with
deciding how much is "enough" for whatever you're doing.

The amount of space that is "enough" for a naive algorithm varies wildly
depending on the particular naive algorithm, and some of them are disgustingly
high. So, in order to have your naive algorithm output an acceptable answer,
you've have to do a lot of numerical stability analysis. But, a similar amount
of numerical stability analysis can usually give you a different algorithm
that can get "enough" precision with the standard sizes of floats. Using
standard floats gives all sorts of other nice benefits like being supported in
modern languages, better speed and space performance, and hardware
acceleration.

Additionally, most of the common pitfalls of floats are well-known by now, and
have well-known solutions, and there are a number of variants and libraries
for most 'workhorse' algorithms that trade between speed, space, and numerical
stability. It comes down to being another case of knowing your tools.

(PS--Numerical stability doesn't have to do with von Neumann architecture,
only with representation of numeric values)

~~~
eru
> (PS--Numerical stability doesn't have to do with von Neumann architecture,
> only with representation of numeric values)

Which you may view as part of the von Neumann architecture. But only if you
compare it to something radically different, like analogue computers.

------
pkaler
This is stuff you should have learned in school in a numerical methods or
numerical analysis class. If you didn't, take a read of _What Every Computer
Scientist Should Know About Floating-Point Arithmetic_.

<http://docs.sun.com/source/806-3568/ncg_goldberg.html>

------
lutorm
I saw the title and immediately thought: that's going to be about expm1 et
al... Anyone who's done numerical analysis connected to Blackbody functions or
mean-free-path calculations will have learned this the hard way, like I did.
;-)

~~~
copper
I think the authors point was that you shouldn't have to learn these things
the hard way :)

------
rmundo
Numerical Methods is one of the few classes I've taken in college that has
direct application to what I'm doing now. Chances are anyone majoring in a
science/engineering field would find it useful later on.

------
someone_here
Why is the "random()" function in libraries always not very random? Why do I
always have to use some other function for more random stuff? Shouldn't the
poor implimentation be called "fast_random()" or something similar?

~~~
sgift
Principle of least surprise.

If programmers want a "good enough" random value they will try random() first.
If you really have a need for a better random value (e.g. usable for strong
cryptography) you can be assumed to already know that random() is not good
enough and can be trusted to search for a better substitute.

~~~
zokier
I'd say that random() producing good pseudorandom numbers is less surprising
than not.

If programmers want a random number they will try random() first. If they find
it's performance inadequate, they can then search for faster alternatives.

~~~
sgift
In theory you are correct, in practice it is more complicated. There is no
such thing as a good pseudo-random number (good meaning good for all use
cases), only good enough. For random() to be useful it has to be generate
"good" number for the average programming task, but must be fast enough that
the average task can use it. Both properties are in conflict and so
programming language implementors have to balance them out.

~~~
mturmon
I disagree with this. For non-cryptographic purposes, the Mersenne Twister is
very pseudorandom, has a huge period, and is competitively fast. There are
also related shift-register techniques that have these properties.

The state of the art in PRNG's is still advancing, but it takes a fairly
esoteric example to find a blind spot of the MT (if initialized properly, as
wrappers should do). Way outside the "generate a random 512x512 bitmap"
variety.

~~~
jacobolus
Which is why Python’s random.random() (and friends) all just use the Mersenne
Twister, and programmers can happily use it for their Monte Carlo simulations
(or whatever), blissfully ignorant of the details of PRNGs.

------
qwzybug
Not quite what the article is talking about, but one of my greatest "aha!"
moments in math/code was when I realized how handy the seemingly redundant
atan2 function is.

~~~
orborde
That is precisely what the article is talking about: math library functions
that seem at first glance to be superfluous, until you happen to need exactly
the thing they are designed to do.

~~~
RiderOfGiraffes
It's much, much more than that - these are math library functions that do what
you want, and are such that the obvious implementations using the obvious math
functions are dangerously unstable or inaccurate.

------
jheriko
Maybe this is why Intel give us F2XM1 (2^x - 1) as an opcode instead of, what
would seem to be more useful, 2^x? Still - it requires a lot of thought from
the programmer to realise he might gain some precision by removing pows and
exps and splitting them up to avoid precision loss - and even more to make the
jump to assembler in order to do it.

~~~
Daniel_Newby
Another reason is that the operation may already be needed as a component of a
trigonometric function, and it doesn't cost much to promote it to a separate
opcode.

~~~
jheriko
There is no 2^x opcode btw...

------
dlsspy
I was assuming I'd find this:

[http://en.wikipedia.org/wiki/C%2B%2B_Technical_Report_1#Math...](http://en.wikipedia.org/wiki/C%2B%2B_Technical_Report_1#Mathematical_special_functions)

Parts of C++ that don't seem to be quite generally useful enough to warrant
inclusion in a standard library.

------
nkassis
why do they all define pi when it's just 22/7?

just kidding ;p

