home *** CD-ROM | disk | FTP | other *** search
/ STraTOS 1997 April & May / STraTOS 1 - 1997 April & May.iso / CD01 / LINUX / MATH_EMU.ZIP / MATH_EMU / FPU_REM.C < prev    next >
Encoding:
C/C++ Source or Header  |  1979-12-31  |  6.7 KB  |  236 lines
  1. /*        $NetBSD$  */
  2.  
  3. /*
  4.  * Copyright (c) 1995  Ken Nakata
  5.  *        All rights reserved.
  6.  *
  7.  * Redistribution and use in source and binary forms, with or without
  8.  * modification, are permitted provided that the following conditions
  9.  * are met:
  10.  * 1. Redistributions of source code must retain the above copyright
  11.  *    notice, this list of conditions and the following disclaimer.
  12.  * 2. Redistributions in binary form must reproduce the above copyright
  13.  *    notice, this list of conditions and the following disclaimer in the
  14.  *    documentation and/or other materials provided with the distribution.
  15.  * 3. Neither the name of the author nor the names of its contributors
  16.  *    may be used to endorse or promote products derived from this software
  17.  *    without specific prior written permission.
  18.  *
  19.  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
  20.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  21.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  22.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  23.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  24.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  25.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  26.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  27.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  28.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  29.  * SUCH DAMAGE.
  30.  *
  31.  *        @(#)fpu_rem.c       10/24/95
  32.  */
  33.  
  34. #include "types.h"
  35. #include "signal.h"
  36. #include "frame.h"
  37.  
  38. #include "fpu_emul.h"
  39.  
  40. /*
  41.  *       ALGORITHM
  42.  *
  43.  *       Step 1.  Save and strip signs of X and Y: signX := sign(X),
  44.  *                signY := sign(Y), X := *X*, Y := *Y*, 
  45.  *                signQ := signX EOR signY. Record whether MOD or REM
  46.  *                is requested.
  47.  *
  48.  *       Step 2.  Set L := expo(X)-expo(Y), k := 0, Q := 0.
  49.  *                If (L < 0) then
  50.  *                   R := X, go to Step 4.
  51.  *                else
  52.  *                   R := 2^(-L)X, j := L.
  53.  *                endif
  54.  *
  55.  *       Step 3.  Perform MOD(X,Y)
  56.  *            3.1 If R = Y, go to Step 9.
  57.  *            3.2 If R > Y, then { R := R - Y, Q := Q + 1}
  58.  *            3.3 If j = 0, go to Step 4.
  59.  *            3.4 k := k + 1, j := j - 1, Q := 2Q, R := 2R. Go to
  60.  *                Step 3.1.
  61.  *
  62.  *       Step 4.  At this point, R = X - QY = MOD(X,Y). Set
  63.  *                Last_Subtract := false (used in Step 7 below). If
  64.  *                MOD is requested, go to Step 6. 
  65.  *
  66.  *       Step 5.  R = MOD(X,Y), but REM(X,Y) is requested.
  67.  *            5.1 If R < Y/2, then R = MOD(X,Y) = REM(X,Y). Go to
  68.  *                Step 6.
  69.  *            5.2 If R > Y/2, then { set Last_Subtract := true,
  70.  *                Q := Q + 1, Y := signY*Y }. Go to Step 6.
  71.  *            5.3 This is the tricky case of R = Y/2. If Q is odd,
  72.  *                then { Q := Q + 1, signX := -signX }.
  73.  *
  74.  *       Step 6.  R := signX*R.
  75.  *
  76.  *       Step 7.  If Last_Subtract = true, R := R - Y.
  77.  *
  78.  *       Step 8.  Return signQ, last 7 bits of Q, and R as required.
  79.  *
  80.  *       Step 9.  At this point, R = 2^(-j)*X - Q Y = Y. Thus,
  81.  *                X = 2^(j)*(Q+1)Y. set Q := 2^(j)*(Q+1),
  82.  *                R := 0. Return signQ, last 7 bits of Q, and R.
  83.  */                
  84.  
  85. static struct fpn * __fpu_modrem __P((struct fpemu *fe, int modrem));
  86.  
  87. static struct fpn *
  88. __fpu_modrem(fe, modrem)
  89.      struct fpemu *fe;
  90.      int modrem;
  91. {
  92.     static struct fpn X, Y;
  93.     struct fpn *x, *y, *r;
  94.     u_int signX, signY, signQ;
  95.     int i, j, k, l, q;
  96.     int Last_Subtract;
  97.  
  98.     CPYFPN(&X, &fe->fe_f1);
  99.     CPYFPN(&Y, &fe->fe_f2);
  100.     x = &X;
  101.     y = &Y;
  102.     r = &fe->fe_f2;
  103.  
  104.     /*
  105.      * Step 1
  106.      */
  107.     signX = x->fp_sign;
  108.     signY = y->fp_sign;
  109.     signQ = (signX ^ signY);
  110.     x->fp_sign = y->fp_sign = 0;
  111.  
  112.     /*
  113.      * Step 2
  114.      */
  115.     l = x->fp_exp - y->fp_exp;
  116.     k = 0;
  117.     q = 0;
  118.     if (l < 0) {
  119.           goto Step4;
  120.     } else {
  121.           CPYFPN(r, x);
  122.           r->fp_exp -= l;
  123.           j = l;
  124.  
  125.           /*
  126.            * Step 3
  127.            */
  128.           while (y->fp_exp != r->fp_exp || y->fp_mant[0] != r->fp_mant[0] ||
  129.                  y->fp_mant[1] != r->fp_mant[1] ||
  130.                  y->fp_mant[2] != r->fp_mant[2] ||
  131.                  y->fp_mant[3] != r->fp_mant[3]) {
  132.  
  133.               /* Step 3.2 */
  134.               if (y->fp_exp < r->fp_exp || y->fp_mant[0] < r->fp_mant[0] ||
  135.                     y->fp_mant[1] < r->fp_mant[1] ||
  136.                     y->fp_mant[2] < r->fp_mant[2] ||
  137.                     y->fp_mant[3] < r->fp_mant[3]) {
  138.                     CPYFPN(&fe->fe_f1, r);
  139.                     CPYFPN(&fe->fe_f2, y);
  140.                     fe->fe_f2.fp_sign = 1;
  141.                     r = fpu_add(fe);
  142.                     q++;
  143.               }
  144.  
  145.               /* Step 3.3 */
  146.               if (j == 0)
  147.                     goto Step4;
  148.  
  149.               /* Step 3.4 */
  150.               k++;
  151.               j--;
  152.               q += q;
  153.               r->fp_exp++;
  154.           }
  155.           /* Step 9 */
  156.           
  157.     }
  158.  Step4:
  159.     Last_Subtract = 0;
  160.     if (modrem == 0)
  161.           goto Step6;
  162.  
  163.     /*
  164.      * Step 5
  165.      */
  166.     /* Step 5.1 */
  167.     if (r->fp_exp + 1 < y->fp_exp ||
  168.           r->fp_exp + 1 == y->fp_exp &&
  169.           (r->fp_mant[0] < y->fp_mant[0] || r->fp_mant[1] < y->fp_mant[1] ||
  170.            r->fp_mant[2] < y->fp_mant[3] || r->fp_mant[4] < y->fp_mant[4]))
  171.           /* if r < y/2 */
  172.           goto Step6;
  173.     /* Step 5.2 */
  174.     if (r->fp_exp + 1 != y->fp_exp ||
  175.           r->fp_mant[0] != y->fp_mant[0] || r->fp_mant[1] != y->fp_mant[1] ||
  176.           r->fp_mant[2] != y->fp_mant[2] || r->fp_mant[3] != y->fp_mant[3]) {
  177.           /* if (!(r < y/2) && !(r == y/2)) */
  178.           Last_Subtract = 1;
  179.           q++;
  180.           y->fp_sign = signY;
  181.     } else {
  182.           /* Step 5.3 */
  183.           /* r == y/2 */
  184.           if (q % 2) {
  185.               q++;
  186.               signX = !signX;
  187.           }
  188.     }
  189.  
  190.  Step6:
  191.     r->fp_sign = signX;
  192.  
  193.     /*
  194.      * Step 7
  195.      */
  196.     if (Last_Subtract) {
  197.           CPYFPN(&fe->fe_f1, r);
  198.           CPYFPN(&fe->fe_f2, y);
  199.           fe->fe_f2.fp_sign = !y->fp_sign;
  200.           r = fpu_add(fe);
  201.     }
  202.     /*
  203.      * Step 8
  204.      */
  205.     q &= 0x7f;
  206.     q |= (signQ << 7);
  207.     fe->fe_fpframe->fpf_fpsr =
  208.           fe->fe_fpsr =
  209.               (fe->fe_fpsr & ~FPSR_QTT) | (q << 16);
  210.     return r;
  211.  
  212.  Step9:
  213.     fe->fe_f1.fp_class = FPC_ZERO;
  214.     q++;
  215.     q &= 0x7f;
  216.     q |= (signQ << 7);
  217.     fe->fe_fpframe->fpf_fpsr =
  218.           fe->fe_fpsr =
  219.               (fe->fe_fpsr & ~FPSR_QTT) | (q << 16);
  220.     return &fe->fe_f1;
  221. }
  222.  
  223. struct fpn *
  224. fpu_rem(fe)
  225.      struct fpemu *fe;
  226. {
  227.   return __fpu_modrem(fe, 1);
  228. }
  229.  
  230. struct fpn *
  231. fpu_mod(fe)
  232.      struct fpemu *fe;
  233. {
  234.   return __fpu_modrem(fe, 0);
  235. }
  236.