text: #include <math.h>
  #include "float_lib.h"
  #include "exp2.h"
  
  float rlibm_fast_exp10(float x) {
    float_x fx;
    fx.f = x;
    
    // Take care of special cases
    if (0x421a209b <= fx.x && fx.x <= 0xb25e5bd8) {
      if (fx.x <= 0x7F800000) return 1.0/0.0;
      if (fx.x < 0x80000000) return 0.0/0.0;
      return 1.0;
    }
    
    if (fx.x <= 0x32de5bd8) {
      return 1.0;
    }
    
    if (fx.x >= 0xc2349e36) {
      if (fx.x <= 0xFF800000) return 0.0;
      return 0.0/0.0;
    }
    
    // Perform range reduction
    double xp = x * 2.12603398072791179629348334856331348419189453125e+02;
    int N = (int)xp;
    int N2 = N % 64;
    if (N2 < 0) N2 += 64;
    int N1 = N - N2;
    
    int M = N1 / 64;
    int J = N2;
    double R = x - N *
    4.703593682249706219022922226713490090332925319671630859375e-03;
    
    double y;
    if (R == -4.673187077026297320259118350804783403873443603515625e-03)
      y = 0.98929727522680910301033918585744686424732208251953125;
    else if (R == -4.655308259524655600802844901409116573631763458251953125e-03)
      y = 0.98933800295939189339833319536410272121429443359375;
    else if (R == -1.31447054445743560791015625e-03)
      y = 0.9969778954982757568359375;
    else {
      y = 0.575063325143780001980076121981255710124969482421875;
      y *= R;
      y += 1.1712319945670899290490751809556968510150909423828125;
      y *= R;
      y += 2.034677894860589475456436048261821269989013671875;
      y *= R;
      y += 2.6509490556773211977770188241265714168548583984375;
      y *= R;
      y += 2.302585092996466631376506484230048954486846923828125;
      y *= R;
      y += 0.9999999999999993338661852249060757458209991455078125;
    }
      
    // Perform output compensation
    return y * ldexp(exp2JBy64[J], M);
  }
  