real procedure Expint(x);
   real x;
begin real y, w, z;
   if x < 1 then
     z ≔ ((((·00107857 × x - ·00976004) × x
     + ·05519968) × x - ·24991055) × x
     + ·99999193) × x - ·57721566 - ln(x)
   else
   begin
      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)
   end;
   Expint ≔ z
end;