title: Range Reduction Poor Precision For Single Precision Floating Point
body: <p>I am trying to implement range reduction as the first step of implementing the sine function.</p>
  
  <p>I am following the method described in the paper <a href="http://www.derekroconnor.net/Software/Ng--ArgReduction.pdf" rel="noreferrer">"ARGUMENT REDUCTION FOR HUGE ARGUMENTS" by K.C. NG</a></p>
  
  <p>I am getting error as large as 0.002339146 when using the input range of x from 0 to 20000.  My error obviously shouldn't be that large, and I'm not sure how I can reduce it.  I noticed that the error magnitude is associated with the input theta magnitude to cosine/sine.</p>
  
  <p>I was able to obtain the nearpi.c code that the paper mentions, but I'm not sure how to utilize the code for single precision floating point.  If anyone is interested, the nearpi.c file can be found at this link: <a href="http://www.cs.berkeley.edu/~wkahan/testpi/nearpi.c" rel="noreferrer">nearpi.c</a></p>
  
  <p>Here is my MATLAB code:</p>
  
  <pre><code>x = 0:0.1:20000;
  
  % Perform range reduction
  % Store constant 2/pi
  twooverpi = single(2/pi);
  
  % Compute y
  y = (x.*twooverpi);
  
  % Compute k (round to nearest integer
  k = round(y);
  
  % Solve for f
  f = single(y-k);
  
  % Solve for r
  r = single(f*single(pi/2));
  
  % Find last two bits of k
  n = bitand(fi(k,1,32,0),fi(3,1,32,0));
  n = single(n);
  
  % Preallocate for speed
  z(length(x)) = 0;
  for i = 1:length(x)
  
      switch(n(i))
          case 0
              z(i)=sin(r(i));
          case 1
              z(i) = single(cos(r(i)));
          case 2
              z(i) = -sin(r(i));
          case 3
              z(i) = single(-cos(r(i)));
          otherwise
      end
  
  end
  
  maxerror = max(abs(single(z - single(sin(single(x))))))
  minerror = min(abs(single(z - single(sin(single(x))))))
  </code></pre>
  
  <p>I have edited the program nearpi.c so that it compiles.  However I am not sure how to interpret the output.  Also the file expects an input, which I had to input by hand, also I am not sure of the significance of the input.</p>
  
  <p>Here is the working nearpi.c: </p>
  
  <pre><code>/*
   ============================================================================
   Name        : nearpi.c
   Author      : 
   Version     :
   Copyright   : Your copyright notice
   Description : Hello World in C, Ansi-style
   ============================================================================
   */
  
  #include &lt;stdio.h&gt;
  #include &lt;stdlib.h&gt;
  #include &lt;math.h&gt;
  
  
  /*
   * Global macro definitions.
   */
  
  # define hex( double )  *(1 + ((long *) &amp;double)), *((long *) &amp;double)
  # define sgn(a)         (a &gt;= 0 ? 1 : -1)
  # define MAX_k          2500
  # define D              56
  # define MAX_EXP        127
  # define THRESHOLD      2.22e-16
  
  /*
   *  Global Variables
   */
  
  int     CFlength,               /* length of CF including terminator */
          binade;
  double  e,
          f;                      /* [e,f] range of D-bit unsigned int of f;
                                     form 1X...X */
  
  // Function Prototypes
  int dbleCF (double i[], double j[]);
  void input (double i[]);
  void nearPiOver2 (double i[]);
  
  
  /*
   *  This is the start of the main program.
   */
  
  int main (void)
  {
      int     k;                  /* subscript variable */
      double  i[MAX_k],
              j[MAX_k];           /* i and j are continued fractions
                                     (coeffs) */
  
  
     // fp = fopen("/src/cfpi.txt", "r");
  
  
  /*
   *  Compute global variables e and f, where
   *
   *      e = 2 ^ (D-1), i.e. the D bit number 10...0
   *  and
   *      f = 2 ^ D - 1, i.e. the D bit number 11...1  .
   */
  
      e = 1;
      for (k = 2; k &lt;= D; k = k + 1)
          e = 2 * e;
      f = 2 * e - 1;
  
   /*
    *  Compute the continued fraction for  (2/e)/(pi/2)  , i.e.
    *  q's starting value for the first binade, given the continued
    *  fraction for  pi  as input; set the global variable CFlength
    *  to the length of the resulting continued fraction (including
    *  its negative valued terminator).  One should use as many
    *  partial coefficients of  pi  as necessary to resolve numbers
    *  of the width of the underflow plus the overflow threshold.
    *  A rule of thumb is 0.97 partial coefficients are generated
    *  for every decimal digit of  pi .
    *
    *  Note: for radix B machines, subroutine  input  should compute
    *  the continued fraction for  (B/e)/(pi/2)  where  e = B ^ (D - 1).
    */
  
      input (i);
  
  /*
   *  Begin main loop over all binades:
   *  For each binade, find the nearest multiples of pi/2 in that binade.
   *
   *  [ Note: for hexadecimal machines ( B = 16 ), the rest of the main
   *  program simplifies(!) to
   *
   *                      B_ade = 1;
   *                      while (B_ade &lt; MAX_EXP)
   *                      {
   *                          dbleCF (i, j);
   *                          dbleCF (j, i);
   *                          dbleCF (i, j);
   *                          CFlength = dbleCF (j, i);
   *                          B_ade = B_ade + 1;
   *                      }
   *                  }
   *
   *  because the alternation of source &amp; destination are no longer necessary. ]
   */
  
      binade = 1;
      while (binade &lt; MAX_EXP)
      {
  
  /*
   *  For the current (odd) binade, find the nearest multiples of pi/2.
   */
  
          nearPiOver2 (i);
  
  /*
   *  Double the continued fraction to get to the next (even) binade.
   *  To save copying arrays, i and j will alternate as the source
   *  and destination for the continued fractions.
   */
  
          CFlength = dbleCF (i, j);
          binade = binade + 1;
  
  /*
   *  Check for main loop termination again because of the
   *  alternation.
   */
  
          if (binade &gt;= MAX_EXP)
              break;
  
  /*
   *  For the current (even) binade, find the nearest multiples of pi/2.
   */
  
          nearPiOver2 (j);
  
  /*
   *  Double the continued fraction to get to the next (odd) binade.
   */
  
          CFlength = dbleCF (j, i);
          binade = binade + 1;
      }
  
      return 0;
  }                               /* end of Main Program */
  
  /*
   *  Subroutine  DbleCF  doubles a continued fraction whose partial
   *  coefficients are i[] into a continued fraction j[], where both
   *  arrays are of a type sufficient to do D-bit integer arithmetic.
   *
   *  In my case ( D = 56 ) , I am forced to treat integers as double
   *  precision reals because my machine does not have integers of
   *  sufficient width to handle D-bit integer arithmetic.
   *
   *  Adapted from a Basic program written by W. Kahan.
   *
   *  Algorithm based on Hurwitz's method of doubling continued
   *  fractions (see Knuth Vol. 3, p.360).
   *
   *  A negative value terminates the last partial quotient.
   *
   *  Note:  for the non-C programmers, the statement  break
   *  exits a loop and the statement  continue  skips to the next
   *  case in the same loop.
   *
   *  The call  modf ( l / 2, &amp;l0 )  assigns the integer portion of
   *  half of L to L0.
   */
  
  int dbleCF (double i[], double j[])
  {
        double k,
                      l,
                      l0,
                      j0;
        int    n,
                      m;
      n = 1;
      m = 0;
      j0 = i[0] + i[0];
      l = i[n];
      while (1)
      {
          if (l &lt; 0)
          {
              j[m] = j0;
              break;
          };
          modf (l / 2, &amp;l0);
          l = l - l0 - l0;
          k = i[n + 1];
          if (l0 &gt; 0)
          {
              j[m] = j0;
              j[m + 1] = l0;
              j0 = 0;
              m = m + 2;
          };
          if (l == 0) {
  /*
   *  Even case.
   */
              if (k &lt; 0)
              {
                  m = m - 1;
                  break;
              }
              else
              {
                  j0 = j0 + k + k;
                  n = n + 2;
                  l = i[n];
                  continue;
              };
          }
  /*
   *  Odd case.
   */
          if (k &lt; 0)
          {
              j[m] = j0 + 2;
              break;
          };
          if (k == 0)
          {
              n = n + 2;
              l = l + i[n];
              continue;
          };
          j[m] = j0 + 1;
          m = m + 1;
          j0 = 1;
          l = k - 1;
          n = n + 1;
          continue;
      };
      m = m + 1;
      j[m] = -99999;
      return (m);
  }
  
  /*
   *  Subroutine  input  computes the continued fraction for
   *  (2/e) / (pi/2) , where  e = 2 ^ (D-1) , given  pi 's
   *  continued fraction as input.  That is, double the continued
   *  fraction of  pi   D-3  times and place a zero at the front.
   *
   *  One should use as many partial coefficients of  pi  as
   *  necessary to resolve numbers of the width of the underflow
   *  plus the overflow threshold.  A rule of thumb is  0.97
   *  partial coefficients are generated for every decimal digit
   *  of  pi .  The last coefficient of  pi  is terminated by a
   *  negative number.
   *
   *  I'll be happy to supply anyone with the partial coefficients
   *  of  pi .  My ARPA address is  mcdonald@ucbdali.BERKELEY.ARPA .
   *
   *  I computed the partial coefficients of  pi  using a method of
   *  Bill Gosper's.  I need only compute with integers, albeit
   *  large ones.  After writing the program in  bc  and  Vaxima  ,
   *  Prof. Fateman suggested  FranzLisp .  To my surprise, FranzLisp
   *  ran the fastest!  the reason?   FranzLisp's  Bignum  package is
   *  hand coded in assembler.  Also,  FranzLisp  can be compiled.
   *
   *
   *  Note: for radix B machines, subroutine  input  should compute
   *  the continued fraction for  (B/e)/(pi/2)  where  e = B ^ (D - 1).
   *  In the case of hexadecimal ( B = 16 ), this is done by repeated
   *  doubling the appropriate number of times.
   */
  
  void input (double i[])
  {
      int     k;
      double  j[MAX_k];
  
  /*
   *  Read in the partial coefficients of  pi  from a precalculated file
   *  until a negative value is encountered.
   */
  
      k = -1;
      do
      {
          k = k + 1;
          scanf ("%lE", &amp;i[k]);
          printf("hello\n");
          printf("%d", k);
      } while (i[k] &gt;= 0);
  
  /*
   *  Double the continued fraction for  pi  D-3  times using
   *  i  and  j  alternately as source and destination.  On my
   *  machine  D = 56  so  D-3  is odd; hence the following code:
   *
   *  Double twice  (D-3)/2  times,
   */
      for (k = 1; k &lt;= (D - 3) / 2; k = k + 1)
      {
          dbleCF (i, j);
          dbleCF (j, i);
      };
  /*
   *  then double once more.
   */
      dbleCF (i, j);
  
  /*
   *  Now append a zero on the front (reciprocate the continued
   *  fraction) and the return the coefficients in  i .
   */
  
      i[0] = 0;
      k = -1;
      do
      {
          k = k + 1;
          i[k + 1] = j[k];
      } while (j[k] &gt;= 0);
  
  /*
   *  Return the length of the continued fraction, including its
   *  terminator and initial zero, in the global variable CFlength.
   */
  
      CFlength = k;
  }
  
  /*
   *  Given a continued fraction's coefficients in an array  i ,
   *  subroutine  nearPiOver2  finds all machine representable
   *  values near a integer multiple of  pi/2  in the current binade.
   */
  
  void nearPiOver2 (double i[])
  {
      int     k,                  /* subscript for recurrences    (see
                                     handout) */
              K;                  /* like  k , but used during cancel. elim.
                                     */
      double  p[MAX_k],           /* product of the q's           (see
                                     handout) */
              q[MAX_k],           /* successive tail evals of CF  (see
                                     handout) */
              j[MAX_k],           /* like convergent numerators   (see
                                     handout) */
              tmp,                /* temporary used during cancellation
                                     elim. */
              mk0,                /* m[k - 1]                     (see
                                     handout) */
              mk,                 /* m[k] is one of the few ints  (see
                                     handout) */
              mkAbs,              /* absolute value of m sub k
                                  */
              mK0,                /* like  mk0 , but used during cancel.
                                     elim. */
              mK,                 /* like  mk  , but used during cancel.
                                     elim. */
              z,                  /* the object of our quest (the argument)
                                  */
              m0,                 /* the mantissa of z as a D-bit integer
                                  */
              x,                  /* the reduced argument         (see
                                     handout) */
              ldexp (),           /* sys routine to multiply by a power of
                                     two  */
              fabs (),            /* sys routine to compute FP absolute
                                     value   */
              floor (),           /* sys routine to compute greatest int &lt;=
                                     value   */
              ceil ();            /* sys routine to compute least int &gt;=
                                     value   */
  
   /*
    *  Compute the q's by evaluating the continued fraction from
    *  bottom up.
    *
    *  Start evaluation with a big number in the terminator position.
    */
  
      q[CFlength] = 1.0 + 30;
  
      for (k = CFlength - 1; k &gt;= 0; k = k - 1)
          q[k] = i[k] + 1 / q[k + 1];
  
  /*
   *  Let  THRESHOLD  be the biggest  | x |  that we are interesed in
   *  seeing.
   *
   *  Compute the p's and j's by the recurrences from the top down.
   *
   *  Stop when
   *
   *        1                          1
   *      -----   &gt;=  THRESHOLD  &gt;   ------    .
   *      2 |j |                     2 |j  |
   *          k                          k+1
   */
  
      p[0] = 1;
      j[0] = 0;
      j[1] = 1;
      k = 0;
      do
      {
          p[k + 1] = -q[k + 1] * p[k];
          if (k &gt; 0)
              j[1 + k] = j[k - 1] - i[k] * j[k];
          k = k + 1;
      } while (1 / (2 * fabs (j[k])) &gt;= THRESHOLD);
  
  /*
   *  Then  mk  runs through the integers between
   *
   *                  k        +                   k        +
   *              (-1)  e / p  -  1/2     &amp;    (-1)  f / p  -  1/2  .
   *                         k                            k
   */
  
      for (mkAbs = floor (e / fabs (p[k]));
              mkAbs &lt;= ceil (f / fabs (p[k])); mkAbs = mkAbs + 1)
      {
  
          mk = mkAbs * sgn (p[k]);
  
  /*
   *  For each  mk ,  mk0  runs through integers between
   *
   *                    +
   *              m  q  -  p  THRESHOLD  .
   *               k  k     k
   */
  
          for (mk0 = floor (mk * q[k] - fabs (p[k]) * THRESHOLD);
                  mk0 &lt;= ceil (mk * q[k] + fabs (p[k]) * THRESHOLD);
                  mk0 = mk0 + 1)
          {
  
  /*
   *  For each pair  { mk ,  mk0 } , check that
   *
   *                             k
   *              m       =  (-1)  ( j   m  - j  m   )
   *               0                  k-1 k    k  k-1
   */
              m0 = (k &amp; 1 ? -1 : 1) * (j[k - 1] * mk - j[k] * mk0);
  
  /*
   *  lies between  e  and  f .
   */
              if (e &lt;= fabs (m0) &amp;&amp; fabs (m0) &lt;= f)
              {
  
  /*
   *  If so, then we have found an
   *
   *                              k
   *              x       =  ((-1)  m  / p  - m ) / j
   *                                 0    k    k     k
   *
   *                      =  ( m  q  - m   ) / p  .
   *                            k  k    k-1     k
   *
   *  But this later formula can suffer cancellation.  Therefore,
   *  run the recurrence for the  mk 's  to get  mK  with minimal
   *   | mK | + | mK0 |  in the hope  mK  is  0  .
   */
                  K = k;
                  mK = mk;
                  mK0 = mk0;
                  while (fabs (mK) &gt; 0)
                  {
                      p[K + 1] = -q[K + 1] * p[K];
                      tmp = mK0 - i[K] * mK;
                      if (fabs (tmp) &gt; fabs (mK0))
                          break;
                      mK0 = mK;
                      mK = tmp;
                      K = K + 1;
                  };
  
  /*
   *  Then
   *              x       =  ( m  q  - m   ) / p
   *                            K  K    K-1     K
   *
   *  as accurately as one could hope.
   */
                  x = (mK * q[K] - mK0) / p[K];
  
  /*
   *  To return  z  and  m0  as positive numbers,
   *   x  must take the sign of  m0  .
   */
                  x = x * sgn (m0);
                  m0 = fabs (m0);
  
  /*d
   *  Set  z = m0 * 2 ^ (binade+1-D) .
   */
                  z = ldexp (m0, binade + 1 - D);
  
  /*
   *  Print  z (hex),  z (dec),  m0 (dec),  binade+1-D,  x (hex), x (dec).
   */
  
                  printf ("%08lx %08lx    Z=%22.16E    M=%17.17G    L+1-%d=%3d    %08lx %08lx    x=%23.16E\n", hex (z), z, m0, D, binade + 1 - D, hex (x), x);
  
              }
          }
      }
  }
  </code></pre>