mirror of
https://github.com/classilla/tenfourfox.git
synced 2024-08-10 13:29:05 +00:00
296 lines
6.6 KiB
C
296 lines
6.6 KiB
C
/* This Source Code Form is subject to the terms of the Mozilla Public
|
|
* License, v. 2.0. If a copy of the MPL was not distributed with this
|
|
* file, You can obtain one at http://mozilla.org/MPL/2.0/. */
|
|
|
|
#ifdef SOLARIS
|
|
#define RF_INLINE_MACROS 1
|
|
#endif
|
|
|
|
static const double TwoTo16=65536.0;
|
|
static const double TwoToMinus16=1.0/65536.0;
|
|
static const double Zero=0.0;
|
|
static const double TwoTo32=65536.0*65536.0;
|
|
static const double TwoToMinus32=1.0/(65536.0*65536.0);
|
|
|
|
#ifdef RF_INLINE_MACROS
|
|
|
|
double upper32(double);
|
|
double lower32(double, double);
|
|
double mod(double, double, double);
|
|
|
|
void i16_to_d16_and_d32x4(const double * /*1/(2^16)*/,
|
|
const double * /* 2^16*/,
|
|
const double * /* 0 */,
|
|
double * /*result16*/,
|
|
double * /* result32 */,
|
|
float * /*source - should be unsigned int*
|
|
converted to float* */);
|
|
|
|
#else
|
|
#ifdef MP_USE_FLOOR
|
|
#include <math.h>
|
|
#else
|
|
#define floor(d) ((double)((unsigned long long)(d)))
|
|
#endif
|
|
|
|
static double upper32(double x)
|
|
{
|
|
return floor(x*TwoToMinus32);
|
|
}
|
|
|
|
static double lower32(double x, double y)
|
|
{
|
|
return x-TwoTo32*floor(x*TwoToMinus32);
|
|
}
|
|
|
|
static double mod(double x, double oneoverm, double m)
|
|
{
|
|
return x-m*floor(x*oneoverm);
|
|
}
|
|
|
|
#endif
|
|
|
|
|
|
static void cleanup(double *dt, int from, int tlen)
|
|
{
|
|
int i;
|
|
double tmp,tmp1,x,x1;
|
|
|
|
tmp=tmp1=Zero;
|
|
/* original code **
|
|
for(i=2*from;i<2*tlen-2;i++)
|
|
{
|
|
x=dt[i];
|
|
dt[i]=lower32(x,Zero)+tmp1;
|
|
tmp1=tmp;
|
|
tmp=upper32(x);
|
|
}
|
|
dt[tlen-2]+=tmp1;
|
|
dt[tlen-1]+=tmp;
|
|
**end original code ***/
|
|
/* new code ***/
|
|
for(i=2*from;i<2*tlen;i+=2)
|
|
{
|
|
x=dt[i];
|
|
x1=dt[i+1];
|
|
dt[i]=lower32(x,Zero)+tmp;
|
|
dt[i+1]=lower32(x1,Zero)+tmp1;
|
|
tmp=upper32(x);
|
|
tmp1=upper32(x1);
|
|
}
|
|
/** end new code **/
|
|
}
|
|
|
|
|
|
void conv_d16_to_i32(unsigned int *i32, double *d16, long long *tmp, int ilen)
|
|
{
|
|
int i;
|
|
long long t, t1, a, b, c, d;
|
|
|
|
t1=0;
|
|
a=(long long)d16[0];
|
|
b=(long long)d16[1];
|
|
for(i=0; i<ilen-1; i++)
|
|
{
|
|
c=(long long)d16[2*i+2];
|
|
t1+=(unsigned int)a;
|
|
t=(a>>32);
|
|
d=(long long)d16[2*i+3];
|
|
t1+=(b&0xffff)<<16;
|
|
t+=(b>>16)+(t1>>32);
|
|
i32[i]=(unsigned int)t1;
|
|
t1=t;
|
|
a=c;
|
|
b=d;
|
|
}
|
|
t1+=(unsigned int)a;
|
|
t=(a>>32);
|
|
t1+=(b&0xffff)<<16;
|
|
i32[i]=(unsigned int)t1;
|
|
}
|
|
|
|
void conv_i32_to_d32(double *d32, unsigned int *i32, int len)
|
|
{
|
|
int i;
|
|
|
|
#pragma pipeloop(0)
|
|
for(i=0;i<len;i++) d32[i]=(double)(i32[i]);
|
|
}
|
|
|
|
|
|
void conv_i32_to_d16(double *d16, unsigned int *i32, int len)
|
|
{
|
|
int i;
|
|
unsigned int a;
|
|
|
|
#pragma pipeloop(0)
|
|
for(i=0;i<len;i++)
|
|
{
|
|
a=i32[i];
|
|
d16[2*i]=(double)(a&0xffff);
|
|
d16[2*i+1]=(double)(a>>16);
|
|
}
|
|
}
|
|
|
|
|
|
void conv_i32_to_d32_and_d16(double *d32, double *d16,
|
|
unsigned int *i32, int len)
|
|
{
|
|
int i = 0;
|
|
unsigned int a;
|
|
|
|
#pragma pipeloop(0)
|
|
#ifdef RF_INLINE_MACROS
|
|
for(;i<len-3;i+=4)
|
|
{
|
|
i16_to_d16_and_d32x4(&TwoToMinus16, &TwoTo16, &Zero,
|
|
&(d16[2*i]), &(d32[i]), (float *)(&(i32[i])));
|
|
}
|
|
#endif
|
|
for(;i<len;i++)
|
|
{
|
|
a=i32[i];
|
|
d32[i]=(double)(i32[i]);
|
|
d16[2*i]=(double)(a&0xffff);
|
|
d16[2*i+1]=(double)(a>>16);
|
|
}
|
|
}
|
|
|
|
|
|
void adjust_montf_result(unsigned int *i32, unsigned int *nint, int len)
|
|
{
|
|
long long acc;
|
|
int i;
|
|
|
|
if(i32[len]>0) i=-1;
|
|
else
|
|
{
|
|
for(i=len-1; i>=0; i--)
|
|
{
|
|
if(i32[i]!=nint[i]) break;
|
|
}
|
|
}
|
|
if((i<0)||(i32[i]>nint[i]))
|
|
{
|
|
acc=0;
|
|
for(i=0;i<len;i++)
|
|
{
|
|
acc=acc+(unsigned long long)(i32[i])-(unsigned long long)(nint[i]);
|
|
i32[i]=(unsigned int)acc;
|
|
acc=acc>>32;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
|
|
/*
|
|
** the lengths of the input arrays should be at least the following:
|
|
** result[nlen+1], dm1[nlen], dm2[2*nlen+1], dt[4*nlen+2], dn[nlen], nint[nlen]
|
|
** all of them should be different from one another
|
|
**
|
|
*/
|
|
void mont_mulf_noconv(unsigned int *result,
|
|
double *dm1, double *dm2, double *dt,
|
|
double *dn, unsigned int *nint,
|
|
int nlen, double dn0)
|
|
{
|
|
int i, j, jj;
|
|
int tmp;
|
|
double digit, m2j, nextm2j, a, b;
|
|
double *dptmp, *pdm1, *pdm2, *pdn, *pdtj, pdn_0, pdm1_0;
|
|
|
|
pdm1=&(dm1[0]);
|
|
pdm2=&(dm2[0]);
|
|
pdn=&(dn[0]);
|
|
pdm2[2*nlen]=Zero;
|
|
|
|
if (nlen!=16)
|
|
{
|
|
for(i=0;i<4*nlen+2;i++) dt[i]=Zero;
|
|
|
|
a=dt[0]=pdm1[0]*pdm2[0];
|
|
digit=mod(lower32(a,Zero)*dn0,TwoToMinus16,TwoTo16);
|
|
|
|
pdtj=&(dt[0]);
|
|
for(j=jj=0;j<2*nlen;j++,jj++,pdtj++)
|
|
{
|
|
m2j=pdm2[j];
|
|
a=pdtj[0]+pdn[0]*digit;
|
|
b=pdtj[1]+pdm1[0]*pdm2[j+1]+a*TwoToMinus16;
|
|
pdtj[1]=b;
|
|
|
|
#pragma pipeloop(0)
|
|
for(i=1;i<nlen;i++)
|
|
{
|
|
pdtj[2*i]+=pdm1[i]*m2j+pdn[i]*digit;
|
|
}
|
|
if((jj==30)) {cleanup(dt,j/2+1,2*nlen+1); jj=0;}
|
|
|
|
digit=mod(lower32(b,Zero)*dn0,TwoToMinus16,TwoTo16);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
a=dt[0]=pdm1[0]*pdm2[0];
|
|
|
|
dt[65]= dt[64]= dt[63]= dt[62]= dt[61]= dt[60]=
|
|
dt[59]= dt[58]= dt[57]= dt[56]= dt[55]= dt[54]=
|
|
dt[53]= dt[52]= dt[51]= dt[50]= dt[49]= dt[48]=
|
|
dt[47]= dt[46]= dt[45]= dt[44]= dt[43]= dt[42]=
|
|
dt[41]= dt[40]= dt[39]= dt[38]= dt[37]= dt[36]=
|
|
dt[35]= dt[34]= dt[33]= dt[32]= dt[31]= dt[30]=
|
|
dt[29]= dt[28]= dt[27]= dt[26]= dt[25]= dt[24]=
|
|
dt[23]= dt[22]= dt[21]= dt[20]= dt[19]= dt[18]=
|
|
dt[17]= dt[16]= dt[15]= dt[14]= dt[13]= dt[12]=
|
|
dt[11]= dt[10]= dt[ 9]= dt[ 8]= dt[ 7]= dt[ 6]=
|
|
dt[ 5]= dt[ 4]= dt[ 3]= dt[ 2]= dt[ 1]=Zero;
|
|
|
|
pdn_0=pdn[0];
|
|
pdm1_0=pdm1[0];
|
|
|
|
digit=mod(lower32(a,Zero)*dn0,TwoToMinus16,TwoTo16);
|
|
pdtj=&(dt[0]);
|
|
|
|
for(j=0;j<32;j++,pdtj++)
|
|
{
|
|
|
|
m2j=pdm2[j];
|
|
a=pdtj[0]+pdn_0*digit;
|
|
b=pdtj[1]+pdm1_0*pdm2[j+1]+a*TwoToMinus16;
|
|
pdtj[1]=b;
|
|
|
|
/**** this loop will be fully unrolled:
|
|
for(i=1;i<16;i++)
|
|
{
|
|
pdtj[2*i]+=pdm1[i]*m2j+pdn[i]*digit;
|
|
}
|
|
*************************************/
|
|
pdtj[2]+=pdm1[1]*m2j+pdn[1]*digit;
|
|
pdtj[4]+=pdm1[2]*m2j+pdn[2]*digit;
|
|
pdtj[6]+=pdm1[3]*m2j+pdn[3]*digit;
|
|
pdtj[8]+=pdm1[4]*m2j+pdn[4]*digit;
|
|
pdtj[10]+=pdm1[5]*m2j+pdn[5]*digit;
|
|
pdtj[12]+=pdm1[6]*m2j+pdn[6]*digit;
|
|
pdtj[14]+=pdm1[7]*m2j+pdn[7]*digit;
|
|
pdtj[16]+=pdm1[8]*m2j+pdn[8]*digit;
|
|
pdtj[18]+=pdm1[9]*m2j+pdn[9]*digit;
|
|
pdtj[20]+=pdm1[10]*m2j+pdn[10]*digit;
|
|
pdtj[22]+=pdm1[11]*m2j+pdn[11]*digit;
|
|
pdtj[24]+=pdm1[12]*m2j+pdn[12]*digit;
|
|
pdtj[26]+=pdm1[13]*m2j+pdn[13]*digit;
|
|
pdtj[28]+=pdm1[14]*m2j+pdn[14]*digit;
|
|
pdtj[30]+=pdm1[15]*m2j+pdn[15]*digit;
|
|
/* no need for cleenup, cannot overflow */
|
|
digit=mod(lower32(b,Zero)*dn0,TwoToMinus16,TwoTo16);
|
|
}
|
|
}
|
|
|
|
conv_d16_to_i32(result,dt+2*nlen,(long long *)dt,nlen+1);
|
|
|
|
adjust_montf_result(result,nint,nlen);
|
|
|
|
}
|
|
|