Main Page   Data Structures   File List   Data Fields   Globals  

poly_fct.c

00001 /*
00002  * Author  : Defour David
00003  * Contact : David.Defour@ens-lyon.fr
00004  *
00005  * This program is free software; you can redistribute it and/or modify
00006  * it under the terms of the GNU Lesser General Public License as published by
00007  * the Free Software Foundation; either version 2 of the License, or 
00008  * (at your option) any later version.
00009  * 
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  *
00015  * You should have received a copy of the GNU Lesser General Public License
00016  * along with this program; if not, write to the Free Software
00017  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  
00018  */
00019 #include "scs.h"
00020 #include "scs_private.h"
00021 
00022 
00023 /*
00024  * result = z + (x * y)
00025  */
00026 /* z->sign = X_SGN . Y_SGN */
00027 void scs_fma(scs_ptr result,  scs_ptr x,  scs_ptr y,  scs_ptr z){
00028   unsigned long long int RES[2*SCS_NB_WORDS];
00029   unsigned long long int val, tmp;
00030   int i, j, ind, Diff;    
00031 
00032   ind = X_IND + Y_IND; 
00033 
00034   for(i=0; i<=SCS_NB_WORDS+1; i++)
00035     RES[i]=0;
00036   
00037   for(i=0 ; i<SCS_NB_WORDS; i++){
00038     for(j=0; j<(SCS_NB_WORDS-i); j++){
00039       RES[i+j] += (unsigned long long int)X_HW[i] * Y_HW[j];
00040     }}
00041 
00042   /* if we can perform an add */
00043   if (z->sign == (X_SGN * Y_SGN)){
00044     Diff = z->index - ind;
00045     if (Diff >= 0){
00046       for(i=(SCS_NB_WORDS-1), j=(SCS_NB_WORDS-Diff); j>=0; i--, j--)
00047         RES[i] = z->h_word[i] + RES[j];     
00048       for(  ; i>=0; i--)
00049         RES[i] = z->h_word[i]; 
00050     }else {    
00051       for(i=(SCS_NB_WORDS+Diff), j=(SCS_NB_WORDS-1); i>=0; i--, j--)
00052         RES[j] = z->h_word[i] + RES[j];     
00053     }
00054 
00055     /* Carry propagate */
00056     RES[SCS_NB_WORDS-1] += (RES[SCS_NB_WORDS]>>SCS_NB_BITS);
00057     for(i=(SCS_NB_WORDS-1); i>0; i--)
00058       {tmp = RES[i]>>SCS_NB_BITS;  RES[i-1] += tmp;  RES[i] -= (tmp<<SCS_NB_BITS);}
00059     
00060     val = RES[0] >> SCS_NB_BITS;
00061     R_IND = X_IND + Y_IND;
00062     
00063     /* Store the result */
00064     if(val != 0){
00065       /* shift all the digits ! */     
00066       R_HW[0] = val;
00067       R_HW[1] = RES[0] - (val<<SCS_NB_BITS);
00068       for(i=2; i<SCS_NB_WORDS; i++)
00069         R_HW[i] = RES[i-1];
00070       
00071       R_IND += 1;
00072     }
00073     else {
00074       for(i=0; i<SCS_NB_WORDS; i++)
00075         R_HW[i] = RES[i];
00076     }
00077     
00078     R_EXP = (z->exception.d + (X_EXP * Y_EXP)) - 1;
00079     R_SGN = X_SGN * Y_SGN;
00080   
00081   }else {
00082     /* we have to do a sub */
00083 
00084     /* Carry propagate */
00085     RES[SCS_NB_WORDS-1] += (RES[SCS_NB_WORDS]>>SCS_NB_BITS);
00086     for(i=(SCS_NB_WORDS-1); i>0; i--)
00087       {tmp = RES[i]>>SCS_NB_BITS;  RES[i-1] += tmp;  RES[i] -= (tmp<<SCS_NB_BITS);}
00088     
00089     val = RES[0] >> SCS_NB_BITS;
00090     R_IND = X_IND + Y_IND;
00091     
00092     /* Store the result */
00093     if(val != 0){
00094       /* shift all the digits ! */     
00095       R_HW[0] = val;
00096       R_HW[1] = RES[0] - (val<<SCS_NB_BITS);
00097       for(i=2; i<SCS_NB_WORDS; i++)
00098         R_HW[i] = RES[i-1];
00099       
00100       R_IND += 1;
00101     }
00102     else {
00103       for(i=0; i<SCS_NB_WORDS; i++)
00104         R_HW[i] = RES[i];
00105     }
00106     
00107     R_EXP = (X_EXP * Y_EXP);
00108     R_SGN = X_SGN * Y_SGN;
00109 
00110     scs_add(result, result, z);
00111   }
00112 }

Generated on Tue Jun 17 10:15:51 2003 for SCSLib by doxygen1.2.15