r̲e̲a̲l̲ p̲r̲o̲c̲e̲d̲u̲r̲e̲ Expint(x);
   r̲e̲a̲l̲ x;
b̲e̲g̲i̲n̲ r̲e̲a̲l̲ y, w, z;
   i̲f̲ x < 1 t̲h̲e̲n̲
     z := ((((·00107857 × x - ·00976004) × x
     + ·05519968) × x - ·24991055) × x
     + ·99999193) × x - ·57721566 - ln(x)
   e̲l̲s̲e̲
   b̲e̲g̲i̲n̲
      y := (((x + 8·5733287401) × x
        + 18·059016973) × x + 8·6347608925) × x
        + ·2677737343 ;
      w := ((( x + 9·5733223454) × x
        + 25·6329561486) × x
        + 21·0996530827) × x + 3·9584969228;
      z := exp(-x) / x × (y/w)
   e̲n̲d̲;
   Expint := z
e̲n̲d̲;