In GNU libm, the implementation of `sin`

is system-dependent. Therefore you can find the implementation, for each platform, somewhere in the appropriate subdirectory of sysdeps.

One directory includes an implementation in C, contributed by IBM. Since October 2011, this is the code that actually runs when you call `sin()`

on a typical x86-64 Linux system. It is apparently faster than the `fsin`

assembly instruction. Source code: sysdeps/ieee754/dbl-64/s_sin.c, look for `__sin (double x)`

.

This code is very complex. No one software algorithm is as fast as possible and also accurate over the whole range of *x* values, so the library implements several different algorithms, and its first job is to look at *x* and decide which algorithm to use.

When *x* is very *very* close to 0, `sin(x) == x`

is the right answer.

A bit further out, `sin(x)`

uses the familiar Taylor series. However, this is only accurate near 0, so...

When the angle is more than about 7°, a different algorithm is used, computing Taylor-series approximations for both sin(x) and cos(x), then using values from a precomputed table to refine the approximation.

When |*x*| > 2, none of the above algorithms would work, so the code starts by computing some value closer to 0 that can be fed to `sin`

or `cos`

instead.

There's yet another branch to deal with *x* being a NaN or infinity.

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. Sometimes a few lines of code would take several paragraphs to explain. For example, these two lines

```
double t = (x * hpinv + toint);
double xn = t - toint;
```

are used (sometimes) in reducing *x* to a value close to 0 that differs from *x* by a multiple of π/2, specifically `xn`

× π/2. The way this is done without division or branching is rather clever. But there's no comment at all!

Older 32-bit versions of GCC/glibc used the `fsin`

instruction, which is surprisingly inaccurate for some inputs. There's a fascinating blog post illustrating this with just 2 lines of code.

fdlibm's implementation of `sin`

in pure C is much simpler than glibc's and is nicely commented. Source code: fdlibm/s_sin.c and fdlibm/k_sin.c

`CLOCKS_PER_SEC`

is a constant which is declared in `<time.h>`

. To get the CPU time used by a task within a C application, use:

```
clock_t begin = clock();
/* here, do your time-consuming job */
clock_t end = clock();
double time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
```

Note that this returns the time as a floating point type. This can be more precise than a second (e.g. you measure 4.52 seconds). Precision depends on the architecture; on modern systems you easily get 10ms or lower, but on older Windows machines (from the Win98 era) it was closer to 60ms.

`clock()`

is standard C; it works "everywhere". There are system-specific functions, such as `getrusage()`

on Unix-like systems.

Java's `System.currentTimeMillis()`

does not measure the same thing. It is a "wall clock": it can help you measure how much time it took for the program to execute, but it does not tell you how much CPU time was used. On a multitasking systems (i.e. all of them), these can be widely different.

## Best Solution

Minimizing a function is like trying to find lowest point on a surface. Think of yourself walking on a hilly surface and that you are trying to get to the lowest point. You would find the direction that goes downhill and walk until it doesn't go downhill anymore. Then you would chose a new direction that goes downhill and walk in that direction until it doesn't go downhill anymore, and so on. Eventually (hopefully) you would reach a point where no direction goes downhill anymore. You would then be at a (local) minimum.

The LM algorithm, and many other minimization algorithms, use this scheme.

Suppose that the function being minimized is F and we are at the point x(n) in our iteration. We wish to find the next iterate x(n+1) such that F(x(n+1)) < F(x(n)), i.e. the function value is smaller. In order to chose x(n+1) we need two things, a direction from x(n) and a step size (how far to go in that direction). The LM algorithm determines these values as follows -

First, compute a linear approximation to F at the point x(n). It is easy to find out the downhill direction of a linear function, so we use the linear approximating function to determine the downhill direction. Next, we need to know how far we can go in this chosen direction. If our approximating linear function is a good approximation for F for a large area around x(n), then we can take a fairly large step. If it's a good approximation only very close to x(n), then we can take only a very small step.

This is what LM does - calculates a linear approximation to F at x(n), thus giving the downhill direction, then it figures out how big a step to take based on how well the linear function approximates F at x(n). LM figures out how good the approximating function is by basically taking a step in the direction thus determined and comparing how much the linear approximation to F decreased to the how much the the actual function F decreased. If they are close, the approximating function is good and we can take a little larger step. If they are not close then the approximation function is not good and we should back off and take a smaller step.