The Standard C Library

Author: P. J. Plauger
4.4
All Stack Overflow 31
This Year Stack Overflow 1
This Month Stack Overflow 1

Comments

by AviewAnew   2017-08-20

Reference Style - All Levels

Beginner

Intermediate

Above Intermediate

Uncategorized Additional C Programming Books

by anonymous   2017-08-20

Read Plauger's "The Standard C Library" for some insights into why various features of the (C89) standard library are as they are - and in particular why parts of the standard I/O library are as they are. One reason is that C runs on very diverse systems and with diverse media; devices such as tapes may well need to be handled somewhat differently from the disk drive you're accustomed to thinking of. Also, on Unix, consider your 'tty' device - it connects a keyboard and a mouse to a screen - three quite different bits of hardware. Coordinating between those is tricky enough; the rules in the standard make it easier.

by anonymous   2017-08-20

The functions from the standard C library that help most are frexp() and ldexp().

  • double frexp(double x, int *p)

    Returns either 0 or a fraction, f in the range [0.5 .. 1.0). Sets *p to the power of two such that the f * 2*p equals x.

  • double ldexp(double x, int p)

    Returns x * 2p

Hence, a decent first approximation to the square root of x can be found by:

int p;
double f = frexp(x, &p);
double r = ldexp(f, p / 2);

As this code shows, you need at most 5 Newton-Raphson 'divide and average' cycles to produce a fully accurate answer:

#include <stdio.h>
#include <math.h>

static double test_sqrt(double x)
{
    if (x <= 0.0)
        return 0.0;
    int p;
    double frac = frexp(x, &p);
    printf(" x = %22.15e;     f = %22.15e; p = %d\n", x, frac, p);
    double x0 = ldexp(frac, p / 2);
    printf("x0 = %22.15e; x0*x0 = %22.15e\n", x0, x0 * x0);
    for (int i = 0; i < 5; i++)
    {
        double x1 = (x / x0 + x0) / 2.0;
        printf("x1 = %22.15e; x1*x1 = %22.15e\n", x1, x1 * x1);
        x0 = x1;
    }
    printf(" r = %22.15e;   r*r = %22.15e;  x = %22.15e\n", x0, x0 * x0, x);
    return x0;
}

int main(void)
{
    double x = 1.23456E-125;
    while (x < 1E+125)
    {
        test_sqrt(x);
        x *= 987.654;
    }

    return 0;
}

Extracts from the output:

 x = 1.234560000000000e-125;     f =  5.223124843710012e-01; p = -414
x0 =  2.539342632858085e-63; x0*x0 = 6.448261007050633e-126
x1 =  3.700536659343551e-63; x1*x1 = 1.369397156714553e-125
x1 =  3.518350710213579e-63; x1*x1 = 1.237879172006039e-125
x1 =  3.513633767135102e-63; x1*x1 = 1.234562224955201e-125
x1 =  3.513630600961067e-63; x1*x1 = 1.234560000001003e-125
x1 =  3.513630600959640e-63; x1*x1 = 1.234560000000000e-125
 r =  3.513630600959640e-63;   r*r = 1.234560000000000e-125;  x = 1.234560000000000e-125
 x = 1.219318122240000e-122;     f =  5.037734516005438e-01; p = -404
x0 =  7.837474714727561e-62; x0*x0 = 6.142600990399386e-123
x1 =  1.169750645469021e-61; x1*x1 = 1.368316572575191e-122
x1 =  1.106062520605688e-61; x1*x1 = 1.223374299488608e-122
x1 =  1.104228909406935e-61; x1*x1 = 1.219321484370028e-122
x1 =  1.104227387018778e-61; x1*x1 = 1.219318122242317e-122
x1 =  1.104227387017728e-61; x1*x1 = 1.219318122240000e-122
 r =  1.104227387017728e-61;   r*r = 1.219318122240000e-122;  x = 1.219318122240000e-122

 x =  8.193538364768438e-27;     f =  6.339443253229819e-01; p = -86
x0 =  7.207112563753321e-14; x0*x0 =  5.194247150661096e-27
x1 =  9.287898167957462e-14; x1*x1 =  8.626505237834758e-27
x1 =  9.054816977123554e-14; x1*x1 =  8.198971048920493e-27
x1 =  9.051817090894058e-14; x1*x1 =  8.193539264700178e-27
x1 =  9.051816593794011e-14; x1*x1 =  8.193538364768461e-27
x1 =  9.051816593793999e-14; x1*x1 =  8.193538364768438e-27
 r =  9.051816593793999e-14;   r*r =  8.193538364768438e-27;  x =  8.193538364768438e-27
 x =  8.092380940117008e-24;     f =  6.114430162915473e-01; p = -76
x0 =  2.224416735012882e-12; x0*x0 =  4.948029811005370e-24
x1 =  2.931197771052297e-12; x1*x1 =  8.591920373021955e-24
x1 =  2.845986967837607e-12; x1*x1 =  8.099641821101499e-24
x1 =  2.844711332870451e-12; x1*x1 =  8.092382567361578e-24
x1 =  2.844711046858202e-12; x1*x1 =  8.092380940117088e-24
x1 =  2.844711046858188e-12; x1*x1 =  8.092380940117008e-24
 r =  2.844711046858188e-12;   r*r =  8.092380940117008e-24;  x =  8.092380940117008e-24

 x =  7.511130028860641e-06;     f =  9.844988351428219e-01; p = -17
x0 =  3.845698574776648e-03; x0*x0 =  1.478939752803914e-05
x1 =  2.899411787388324e-03; x1*x1 =  8.406588712846356e-06
x1 =  2.744991037655443e-03; x1*x1 =  7.534975796808706e-06
x1 =  2.740647532044505e-03; x1*x1 =  7.511148894901633e-06
x1 =  2.740644090149702e-03; x1*x1 =  7.511130028872487e-06
x1 =  2.740644090147541e-03; x1*x1 =  7.511130028860642e-06
 r =  2.740644090147541e-03;   r*r =  7.511130028860642e-06;  x =  7.511130028860641e-06
 x =  7.418397617524327e-03;     f =  9.495548950431139e-01; p = -7
x0 =  1.186943618803892e-01; x0*x0 =  1.408835154219280e-02
x1 =  9.059718094019462e-02; x1*x1 =  8.207849194310363e-03
x1 =  8.624024859090236e-02; x1*x1 =  7.437380477020637e-03
x1 =  8.613019058546711e-02; x1*x1 =  7.418409730288888e-03
x1 =  8.613012026886571e-02; x1*x1 =  7.418397617529272e-03
x1 =  8.613012026883701e-02; x1*x1 =  7.418397617524329e-03
 r =  8.613012026883701e-02;   r*r =  7.418397617524329e-03;  x =  7.418397617524327e-03
 x =  7.326810080538372e+00;     f =  9.158512600672964e-01; p = 3
x0 =  1.831702520134593e+00; x0*x0 =  3.355134122267418e+00
x1 =  2.915851260067297e+00; x1*x1 =  8.502188570836042e+00
x1 =  2.714301457717203e+00; x1*x1 =  7.367432403365734e+00
x1 =  2.706818441652081e+00; x1*x1 =  7.326866076067802e+00
x1 =  2.706808098230342e+00; x1*x1 =  7.326810080645360e+00
x1 =  2.706808098210579e+00; x1*x1 =  7.326810080538371e+00
 r =  2.706808098210579e+00;   r*r =  7.326810080538371e+00;  x =  7.326810080538372e+00
 x =  7.236353283284045e+03;     f =  8.833439066508844e-01; p = 13
x0 =  5.653401002565660e+01; x0*x0 =  3.196094289581041e+03
x1 =  9.226700501282829e+01; x1*x1 =  8.513200214037281e+03
x1 =  8.534770092045144e+01; x1*x1 =  7.284230052406829e+03
x1 =  8.506722020095668e+01; x1*x1 =  7.236431952718052e+03
x1 =  8.506675780525467e+01; x1*x1 =  7.236353283497856e+03
x1 =  8.506675780399794e+01; x1*x1 =  7.236353283284046e+03
 r =  8.506675780399794e+01;   r*r =  7.236353283284046e+03;  x =  7.236353283284045e+03
 x =  7.147013265648620e+06;     f =  8.519903738079810e-01; p = 23
x0 =  1.744876285558745e+03; x0*x0 =  3.044593251905283e+06
x1 =  2.920438142779372e+03; x1*x1 =  8.528958945800630e+06
x1 =  2.683839109930689e+03; x1*x1 =  7.202992367993553e+06
x1 =  2.673410187023612e+03; x1*x1 =  7.147122028081623e+06
x1 =  2.673389845507460e+03; x1*x1 =  7.147013266062398e+06
x1 =  2.673389845430072e+03; x1*x1 =  7.147013265648622e+06
 r =  2.673389845430072e+03;   r*r =  7.147013265648622e+06;  x =  7.147013265648620e+06

 x =  5.714937347533526e+60;     f =  8.891035606430090e-01; p = 202
x0 =  2.254145324628333e+30; x0*x0 =  5.081171144543771e+60
x1 =  2.394723262542396e+30; x1*x1 =  5.734699504161695e+60
x1 =  2.390597074573772e+30; x1*x1 =  5.714954372960678e+60
x1 =  2.390593513658524e+30; x1*x1 =  5.714937347546205e+60
x1 =  2.390593513655872e+30; x1*x1 =  5.714937347533526e+60
x1 =  2.390593513655872e+30; x1*x1 =  5.714937347533526e+60
 r =  2.390593513655872e+30;   r*r =  5.714937347533526e+60;  x =  5.714937347533526e+60
 x =  5.644380731040877e+63;     f =  8.575455938313579e-01; p = 212
x0 =  6.957236395157723e+31; x0*x0 =  4.840313825810723e+63
x1 =  7.535100118309196e+31; x1*x1 =  5.677773379294325e+63
x1 =  7.512942052902534e+31; x1*x1 =  5.644429829027134e+63
x1 =  7.512909377296952e+31; x1*x1 =  5.644380731147647e+63
x1 =  7.512909377225895e+31; x1*x1 =  5.644380731040878e+63
x1 =  7.512909377225895e+31; x1*x1 =  5.644380731040878e+63
 r =  7.512909377225895e+31;   r*r =  5.644380731040878e+63;  x =  5.644380731040877e+63

 x = 4.457671009241398e+120;     f =  8.631370354723219e-01; p = 401
x0 =  1.387007739709396e+60; x0*x0 = 1.923790470013767e+120
x1 =  2.300441914113688e+60; x1*x1 = 5.292033000211049e+120
x1 =  2.119093716219477e+60; x1*x1 = 4.490558178120875e+120
x1 =  2.111333991241823e+60; x1*x1 = 4.457731222573127e+120
x1 =  2.111319731695020e+60; x1*x1 = 4.457671009444733e+120
x1 =  2.111319731646867e+60; x1*x1 = 4.457671009241398e+120
 r =  2.111319731646867e+60;   r*r = 4.457671009241398e+120;  x = 4.457671009241398e+120
 x = 4.402636602961303e+123;     f =  8.325007281566217e-01; p = 411
x0 =  4.280886694234198e+61; x0*x0 = 1.832599088887140e+123
x1 =  7.282645088745868e+61; x1*x1 = 5.303691948863432e+123
x1 =  6.664013166606369e+61; x1*x1 = 4.440907148470304e+123
x1 =  6.635298828449913e+61; x1*x1 = 4.402719054282879e+123
x1 =  6.635236697622269e+61; x1*x1 = 4.402636603347328e+123
x1 =  6.635236697331380e+61; x1*x1 = 4.402636602961305e+123
 r =  6.635236697331380e+61;   r*r = 4.402636602961305e+123;  x = 4.402636602961303e+123

You can apply more sophisticated techniques if you wish (for example, improving the estimate of the fraction), but the frexp() and ldexp() functions are likely to be useful.

You can find a discussion of this in Plauger's The Standard C Library, and also in Bentley's More Programming Pearls.

by anonymous   2017-08-20

Most compilers provide source code for the library - however that source code is usually rather more complex and convoluted than you might expect.

A good resource for this is P.J. Plauger's book, "The Standard C Library", which can be had pretty cheaply if you go for a used one. The code presented is not always straight-forward, but Plauger explains it quite well (and gives the reasons why it can't always be straight-forward and still follow the standard).