usr/src/lib/libbc/libc/gen/common/_Qfmul.c

   1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License, Version 1.0 only
   6  * (the "License").  You may not use this file except in compliance
   7  * with the License.
   8  *
   9  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
  10  * or http://www.opensolaris.org/os/licensing.
  11  * See the License for the specific language governing permissions
  12  * and limitations under the License.
  13  *
  14  * When distributing Covered Code, include this CDDL HEADER in each
  15  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  16  * If applicable, add the following below this CDDL HEADER, with the
  17  * fields enclosed by brackets "[]" replaced with your own identifying
  18  * information: Portions Copyright [yyyy] [name of copyright owner]
  19  *
  20  * CDDL HEADER END
  21  */
  22 #pragma ident   "%Z%%M% %I%     %E% SMI"
  23
  24 /*
  25  * Copyright (c) 1988 by Sun Microsystems, Inc.
  26  */
  27
  28 #include "_Qquad.h"
  29 #include "_Qglobals.h"
  30
  31
  32 void
  33 _fp_mul(px, py, pz)
  34         unpacked       *px, *py, *pz;
  35
  36 {
  37         unpacked       *pt;
  38         unsigned       acc[4];          /* Product accumulator. */
  39         unsigned       i,j,y,*x,s,r,c;
  40
  41         if ((int) px->fpclass < (int) py->fpclass) {
  42                 pt = px;
  43                 px = py;
  44                 py = pt;
  45         }
  46         /* Now class(x) >= class(y).  */
  47
  48         *pz = *px;
  49         pz->sign = px->sign ^ py->sign;
  50
  51         switch (px->fpclass) {
  52         case fp_quiet:
  53         case fp_signaling:
  54         case fp_zero:
  55                 return;
  56         case fp_infinity:
  57                 if (py->fpclass == fp_zero) {
  58                         fpu_error_nan(pz);
  59                         pz->fpclass = fp_quiet;
  60                 }
  61                 return;
  62         case fp_normal:
  63                 if (py->fpclass == fp_zero) {
  64                         pz->fpclass = fp_zero;
  65                         return;
  66                 }
  67         }
  68
  69         /* Now x and y are both normal or subnormal. */
  70
  71         x = px->significand;    /* save typing */
  72
  73         s=r=acc[0]=acc[1]=acc[2]=acc[3]=0;      /* intialize acc to zero */
  74
  75         y = py->significand[3];         /* py->significand[3] * x */
  76         if(y!=0) {
  77             j=1;
  78             do {
  79                 s |= r;         /* shift acc right one bit */
  80                 r  = acc[3]&1;
  81                 acc[3] = ((acc[2]&1)<<31)|(acc[3]>>1);
  82                 acc[2] = ((acc[1]&1)<<31)|(acc[2]>>1);
  83                 acc[1] = ((acc[0]&1)<<31)|(acc[1]>>1);
  84                 acc[0] = (acc[0]>>1);
  85                 if(j&y) {               /* bit i of y != 0, add x to acc */
  86                         c = 0;
  87                         c = fpu_add3wc(&acc[3],acc[3],x[3],c);
  88                         c = fpu_add3wc(&acc[2],acc[2],x[2],c);
  89                         c = fpu_add3wc(&acc[1],acc[1],x[1],c);
  90                         c = fpu_add3wc(&acc[0],acc[0],x[0],c);
  91                 }
  92                 j += j;
  93             } while (j!=0);
  94         }
  95
  96         y = py->significand[2];         /* py->significand[2] * x */
  97         if(y!=0) {
  98             j=1;
  99             do {
 100                 s |= r;         /* shift acc right one bit */
 101                 r  = acc[3]&1;
 102                 acc[3] = ((acc[2]&1)<<31)|(acc[3]>>1);
 103                 acc[2] = ((acc[1]&1)<<31)|(acc[2]>>1);
 104                 acc[1] = ((acc[0]&1)<<31)|(acc[1]>>1);
 105                 acc[0] = (acc[0]>>1);
 106                 if(j&y) {               /* bit i of y != 0, add x to acc */
 107                         c = 0;
 108                         c = fpu_add3wc(&acc[3],acc[3],x[3],c);
 109                         c = fpu_add3wc(&acc[2],acc[2],x[2],c);
 110                         c = fpu_add3wc(&acc[1],acc[1],x[1],c);
 111                         c = fpu_add3wc(&acc[0],acc[0],x[0],c);
 112                 }
 113                 j += j;
 114             } while (j!=0);
 115         } else {
 116                 s |= r|(acc[3]&0x7fffffff);
 117                 r  = (acc[3]&0x80000000)>>31;
 118                 acc[3]=acc[2];acc[2]=acc[1];acc[1]=acc[0];acc[0]=0;
 119         }
 120
 121         y = py->significand[1];         /* py->significand[1] * x */
 122         if(y!=0) {
 123             j=1;
 124             do {
 125                 s |= r;         /* shift acc right one bit */
 126                 r  = acc[3]&1;
 127                 acc[3] = ((acc[2]&1)<<31)|(acc[3]>>1);
 128                 acc[2] = ((acc[1]&1)<<31)|(acc[2]>>1);
 129                 acc[1] = ((acc[0]&1)<<31)|(acc[1]>>1);
 130                 acc[0] = (acc[0]>>1);
 131                 if(j&y) {               /* bit i of y != 0, add x to acc */
 132                         c = 0;
 133                         c = fpu_add3wc(&acc[3],acc[3],x[3],c);
 134                         c = fpu_add3wc(&acc[2],acc[2],x[2],c);
 135                         c = fpu_add3wc(&acc[1],acc[1],x[1],c);
 136                         c = fpu_add3wc(&acc[0],acc[0],x[0],c);
 137                 }
 138                 j += j;
 139             } while (j!=0);
 140         } else {
 141                 s |= r|(acc[3]&0x7fffffff);
 142                 r  = (acc[3]&0x80000000)>>31;
 143                 acc[3]=acc[2];acc[2]=acc[1];acc[1]=acc[0];acc[0]=0;
 144         }
 145
 146                                         /* py->significand[0] * x */
 147         y = py->significand[0];         /* y is of form 0x0001???? */
 148         j=1;
 149         do {
 150                 s |= r;         /* shift acc right one bit */
 151                 r  = acc[3]&1;
 152                 acc[3] = ((acc[2]&1)<<31)|(acc[3]>>1);
 153                 acc[2] = ((acc[1]&1)<<31)|(acc[2]>>1);
 154                 acc[1] = ((acc[0]&1)<<31)|(acc[1]>>1);
 155                 acc[0] = (acc[0]>>1);
 156                 if(j&y) {               /* bit i of y != 0, add x to acc */
 157                         c = 0;
 158                         c = fpu_add3wc(&acc[3],acc[3],x[3],c);
 159                         c = fpu_add3wc(&acc[2],acc[2],x[2],c);
 160                         c = fpu_add3wc(&acc[1],acc[1],x[1],c);
 161                         c = fpu_add3wc(&acc[0],acc[0],x[0],c);
 162                 }
 163                 j += j;
 164         } while (j<=y);
 165
 166         if(acc[0]>=0x20000) {   /* right shift one bit to normalize */
 167                 pz->exponent = px->exponent + py->exponent + 1;
 168                 pz->sticky = s|r;
 169                 pz->rounded = acc[3]&1;
 170                 pz->significand[3]=((acc[2]&1)<<31)|(acc[3]>>1);
 171                 pz->significand[2]=((acc[1]&1)<<31)|(acc[2]>>1);
 172                 pz->significand[1]=((acc[0]&1)<<31)|(acc[1]>>1);
 173                 pz->significand[0]=(acc[0]>>1);
 174         } else {
 175                 pz->exponent = px->exponent + py->exponent;
 176                 pz->sticky = s;
 177                 pz->rounded = r;
 178                 pz->significand[3]=acc[3];
 179                 pz->significand[2]=acc[2];
 180                 pz->significand[1]=acc[1];
 181                 pz->significand[0]=acc[0];
 182         }
 183 }