00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #include "scs.h"
00020 #include "scs_private.h"
00021
00022
00023
00024
00025
00026
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
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
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
00064 if(val != 0){
00065
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
00083
00084
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
00093 if(val != 0){
00094
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 }