1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/security/nss/lib/freebl/mpi/montmulf.c Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,295 @@
1.4 +/* This Source Code Form is subject to the terms of the Mozilla Public
1.5 + * License, v. 2.0. If a copy of the MPL was not distributed with this
1.6 + * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
1.7 +
1.8 +#ifdef SOLARIS
1.9 +#define RF_INLINE_MACROS 1
1.10 +#endif
1.11 +
1.12 +static const double TwoTo16=65536.0;
1.13 +static const double TwoToMinus16=1.0/65536.0;
1.14 +static const double Zero=0.0;
1.15 +static const double TwoTo32=65536.0*65536.0;
1.16 +static const double TwoToMinus32=1.0/(65536.0*65536.0);
1.17 +
1.18 +#ifdef RF_INLINE_MACROS
1.19 +
1.20 +double upper32(double);
1.21 +double lower32(double, double);
1.22 +double mod(double, double, double);
1.23 +
1.24 +void i16_to_d16_and_d32x4(const double * /*1/(2^16)*/,
1.25 + const double * /* 2^16*/,
1.26 + const double * /* 0 */,
1.27 + double * /*result16*/,
1.28 + double * /* result32 */,
1.29 + float * /*source - should be unsigned int*
1.30 + converted to float* */);
1.31 +
1.32 +#else
1.33 +#ifdef MP_USE_FLOOR
1.34 +#include <math.h>
1.35 +#else
1.36 +#define floor(d) ((double)((unsigned long long)(d)))
1.37 +#endif
1.38 +
1.39 +static double upper32(double x)
1.40 +{
1.41 + return floor(x*TwoToMinus32);
1.42 +}
1.43 +
1.44 +static double lower32(double x, double y)
1.45 +{
1.46 + return x-TwoTo32*floor(x*TwoToMinus32);
1.47 +}
1.48 +
1.49 +static double mod(double x, double oneoverm, double m)
1.50 +{
1.51 + return x-m*floor(x*oneoverm);
1.52 +}
1.53 +
1.54 +#endif
1.55 +
1.56 +
1.57 +static void cleanup(double *dt, int from, int tlen)
1.58 +{
1.59 + int i;
1.60 + double tmp,tmp1,x,x1;
1.61 +
1.62 + tmp=tmp1=Zero;
1.63 + /* original code **
1.64 + for(i=2*from;i<2*tlen-2;i++)
1.65 + {
1.66 + x=dt[i];
1.67 + dt[i]=lower32(x,Zero)+tmp1;
1.68 + tmp1=tmp;
1.69 + tmp=upper32(x);
1.70 + }
1.71 + dt[tlen-2]+=tmp1;
1.72 + dt[tlen-1]+=tmp;
1.73 + **end original code ***/
1.74 + /* new code ***/
1.75 + for(i=2*from;i<2*tlen;i+=2)
1.76 + {
1.77 + x=dt[i];
1.78 + x1=dt[i+1];
1.79 + dt[i]=lower32(x,Zero)+tmp;
1.80 + dt[i+1]=lower32(x1,Zero)+tmp1;
1.81 + tmp=upper32(x);
1.82 + tmp1=upper32(x1);
1.83 + }
1.84 + /** end new code **/
1.85 +}
1.86 +
1.87 +
1.88 +void conv_d16_to_i32(unsigned int *i32, double *d16, long long *tmp, int ilen)
1.89 +{
1.90 +int i;
1.91 +long long t, t1, a, b, c, d;
1.92 +
1.93 + t1=0;
1.94 + a=(long long)d16[0];
1.95 + b=(long long)d16[1];
1.96 + for(i=0; i<ilen-1; i++)
1.97 + {
1.98 + c=(long long)d16[2*i+2];
1.99 + t1+=(unsigned int)a;
1.100 + t=(a>>32);
1.101 + d=(long long)d16[2*i+3];
1.102 + t1+=(b&0xffff)<<16;
1.103 + t+=(b>>16)+(t1>>32);
1.104 + i32[i]=(unsigned int)t1;
1.105 + t1=t;
1.106 + a=c;
1.107 + b=d;
1.108 + }
1.109 + t1+=(unsigned int)a;
1.110 + t=(a>>32);
1.111 + t1+=(b&0xffff)<<16;
1.112 + i32[i]=(unsigned int)t1;
1.113 +}
1.114 +
1.115 +void conv_i32_to_d32(double *d32, unsigned int *i32, int len)
1.116 +{
1.117 +int i;
1.118 +
1.119 +#pragma pipeloop(0)
1.120 + for(i=0;i<len;i++) d32[i]=(double)(i32[i]);
1.121 +}
1.122 +
1.123 +
1.124 +void conv_i32_to_d16(double *d16, unsigned int *i32, int len)
1.125 +{
1.126 +int i;
1.127 +unsigned int a;
1.128 +
1.129 +#pragma pipeloop(0)
1.130 + for(i=0;i<len;i++)
1.131 + {
1.132 + a=i32[i];
1.133 + d16[2*i]=(double)(a&0xffff);
1.134 + d16[2*i+1]=(double)(a>>16);
1.135 + }
1.136 +}
1.137 +
1.138 +
1.139 +void conv_i32_to_d32_and_d16(double *d32, double *d16,
1.140 + unsigned int *i32, int len)
1.141 +{
1.142 +int i = 0;
1.143 +unsigned int a;
1.144 +
1.145 +#pragma pipeloop(0)
1.146 +#ifdef RF_INLINE_MACROS
1.147 + for(;i<len-3;i+=4)
1.148 + {
1.149 + i16_to_d16_and_d32x4(&TwoToMinus16, &TwoTo16, &Zero,
1.150 + &(d16[2*i]), &(d32[i]), (float *)(&(i32[i])));
1.151 + }
1.152 +#endif
1.153 + for(;i<len;i++)
1.154 + {
1.155 + a=i32[i];
1.156 + d32[i]=(double)(i32[i]);
1.157 + d16[2*i]=(double)(a&0xffff);
1.158 + d16[2*i+1]=(double)(a>>16);
1.159 + }
1.160 +}
1.161 +
1.162 +
1.163 +void adjust_montf_result(unsigned int *i32, unsigned int *nint, int len)
1.164 +{
1.165 +long long acc;
1.166 +int i;
1.167 +
1.168 + if(i32[len]>0) i=-1;
1.169 + else
1.170 + {
1.171 + for(i=len-1; i>=0; i--)
1.172 + {
1.173 + if(i32[i]!=nint[i]) break;
1.174 + }
1.175 + }
1.176 + if((i<0)||(i32[i]>nint[i]))
1.177 + {
1.178 + acc=0;
1.179 + for(i=0;i<len;i++)
1.180 + {
1.181 + acc=acc+(unsigned long long)(i32[i])-(unsigned long long)(nint[i]);
1.182 + i32[i]=(unsigned int)acc;
1.183 + acc=acc>>32;
1.184 + }
1.185 + }
1.186 +}
1.187 +
1.188 +
1.189 +
1.190 +
1.191 +/*
1.192 +** the lengths of the input arrays should be at least the following:
1.193 +** result[nlen+1], dm1[nlen], dm2[2*nlen+1], dt[4*nlen+2], dn[nlen], nint[nlen]
1.194 +** all of them should be different from one another
1.195 +**
1.196 +*/
1.197 +void mont_mulf_noconv(unsigned int *result,
1.198 + double *dm1, double *dm2, double *dt,
1.199 + double *dn, unsigned int *nint,
1.200 + int nlen, double dn0)
1.201 +{
1.202 + int i, j, jj;
1.203 + int tmp;
1.204 + double digit, m2j, nextm2j, a, b;
1.205 + double *dptmp, *pdm1, *pdm2, *pdn, *pdtj, pdn_0, pdm1_0;
1.206 +
1.207 + pdm1=&(dm1[0]);
1.208 + pdm2=&(dm2[0]);
1.209 + pdn=&(dn[0]);
1.210 + pdm2[2*nlen]=Zero;
1.211 +
1.212 + if (nlen!=16)
1.213 + {
1.214 + for(i=0;i<4*nlen+2;i++) dt[i]=Zero;
1.215 +
1.216 + a=dt[0]=pdm1[0]*pdm2[0];
1.217 + digit=mod(lower32(a,Zero)*dn0,TwoToMinus16,TwoTo16);
1.218 +
1.219 + pdtj=&(dt[0]);
1.220 + for(j=jj=0;j<2*nlen;j++,jj++,pdtj++)
1.221 + {
1.222 + m2j=pdm2[j];
1.223 + a=pdtj[0]+pdn[0]*digit;
1.224 + b=pdtj[1]+pdm1[0]*pdm2[j+1]+a*TwoToMinus16;
1.225 + pdtj[1]=b;
1.226 +
1.227 +#pragma pipeloop(0)
1.228 + for(i=1;i<nlen;i++)
1.229 + {
1.230 + pdtj[2*i]+=pdm1[i]*m2j+pdn[i]*digit;
1.231 + }
1.232 + if((jj==30)) {cleanup(dt,j/2+1,2*nlen+1); jj=0;}
1.233 +
1.234 + digit=mod(lower32(b,Zero)*dn0,TwoToMinus16,TwoTo16);
1.235 + }
1.236 + }
1.237 + else
1.238 + {
1.239 + a=dt[0]=pdm1[0]*pdm2[0];
1.240 +
1.241 + dt[65]= dt[64]= dt[63]= dt[62]= dt[61]= dt[60]=
1.242 + dt[59]= dt[58]= dt[57]= dt[56]= dt[55]= dt[54]=
1.243 + dt[53]= dt[52]= dt[51]= dt[50]= dt[49]= dt[48]=
1.244 + dt[47]= dt[46]= dt[45]= dt[44]= dt[43]= dt[42]=
1.245 + dt[41]= dt[40]= dt[39]= dt[38]= dt[37]= dt[36]=
1.246 + dt[35]= dt[34]= dt[33]= dt[32]= dt[31]= dt[30]=
1.247 + dt[29]= dt[28]= dt[27]= dt[26]= dt[25]= dt[24]=
1.248 + dt[23]= dt[22]= dt[21]= dt[20]= dt[19]= dt[18]=
1.249 + dt[17]= dt[16]= dt[15]= dt[14]= dt[13]= dt[12]=
1.250 + dt[11]= dt[10]= dt[ 9]= dt[ 8]= dt[ 7]= dt[ 6]=
1.251 + dt[ 5]= dt[ 4]= dt[ 3]= dt[ 2]= dt[ 1]=Zero;
1.252 +
1.253 + pdn_0=pdn[0];
1.254 + pdm1_0=pdm1[0];
1.255 +
1.256 + digit=mod(lower32(a,Zero)*dn0,TwoToMinus16,TwoTo16);
1.257 + pdtj=&(dt[0]);
1.258 +
1.259 + for(j=0;j<32;j++,pdtj++)
1.260 + {
1.261 +
1.262 + m2j=pdm2[j];
1.263 + a=pdtj[0]+pdn_0*digit;
1.264 + b=pdtj[1]+pdm1_0*pdm2[j+1]+a*TwoToMinus16;
1.265 + pdtj[1]=b;
1.266 +
1.267 + /**** this loop will be fully unrolled:
1.268 + for(i=1;i<16;i++)
1.269 + {
1.270 + pdtj[2*i]+=pdm1[i]*m2j+pdn[i]*digit;
1.271 + }
1.272 + *************************************/
1.273 + pdtj[2]+=pdm1[1]*m2j+pdn[1]*digit;
1.274 + pdtj[4]+=pdm1[2]*m2j+pdn[2]*digit;
1.275 + pdtj[6]+=pdm1[3]*m2j+pdn[3]*digit;
1.276 + pdtj[8]+=pdm1[4]*m2j+pdn[4]*digit;
1.277 + pdtj[10]+=pdm1[5]*m2j+pdn[5]*digit;
1.278 + pdtj[12]+=pdm1[6]*m2j+pdn[6]*digit;
1.279 + pdtj[14]+=pdm1[7]*m2j+pdn[7]*digit;
1.280 + pdtj[16]+=pdm1[8]*m2j+pdn[8]*digit;
1.281 + pdtj[18]+=pdm1[9]*m2j+pdn[9]*digit;
1.282 + pdtj[20]+=pdm1[10]*m2j+pdn[10]*digit;
1.283 + pdtj[22]+=pdm1[11]*m2j+pdn[11]*digit;
1.284 + pdtj[24]+=pdm1[12]*m2j+pdn[12]*digit;
1.285 + pdtj[26]+=pdm1[13]*m2j+pdn[13]*digit;
1.286 + pdtj[28]+=pdm1[14]*m2j+pdn[14]*digit;
1.287 + pdtj[30]+=pdm1[15]*m2j+pdn[15]*digit;
1.288 + /* no need for cleenup, cannot overflow */
1.289 + digit=mod(lower32(b,Zero)*dn0,TwoToMinus16,TwoTo16);
1.290 + }
1.291 + }
1.292 +
1.293 + conv_d16_to_i32(result,dt+2*nlen,(long long *)dt,nlen+1);
1.294 +
1.295 + adjust_montf_result(result,nint,nlen);
1.296 +
1.297 +}
1.298 +
