Home | History | Annotate | Download | only in bignum
      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 /*
     23  * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
     24  * Use is subject to license terms.
     25  */
     26 
     27 #pragma ident	"%Z%%M%	%I%	%E% SMI"
     28 
     29 /*
     30  * If compiled without -DRF_INLINE_MACROS then needs -lm at link time
     31  * If compiled with -DRF_INLINE_MACROS then needs conv.il at compile time
     32  * (i.e. cc <compileer_flags> -DRF_INLINE_MACROS conv.il mont_mulf.c )
     33  */
     34 
     35 #include <sys/types.h>
     36 #include <math.h>
     37 
     38 static const double TwoTo16 = 65536.0;
     39 static const double TwoToMinus16 = 1.0/65536.0;
     40 static const double Zero = 0.0;
     41 static const double TwoTo32 = 65536.0 * 65536.0;
     42 static const double TwoToMinus32 = 1.0 / (65536.0 * 65536.0);
     43 
     44 #ifdef RF_INLINE_MACROS
     45 
     46 double upper32(double);
     47 double lower32(double, double);
     48 double mod(double, double, double);
     49 
     50 #else
     51 
     52 static double
     53 upper32(double x)
     54 {
     55 	return (floor(x * TwoToMinus32));
     56 }
     57 
     58 
     59 static double
     60 lower32(double x, double y)
     61 {
     62 	return (x - TwoTo32 * floor(x * TwoToMinus32));
     63 }
     64 
     65 static double
     66 mod(double x, double oneoverm, double m)
     67 {
     68 	return (x - m * floor(x * oneoverm));
     69 }
     70 
     71 #endif
     72 
     73 
     74 static void
     75 cleanup(double *dt, int from, int tlen)
     76 {
     77 	int i;
     78 	double tmp, tmp1, x, x1;
     79 
     80 	tmp = tmp1 = Zero;
     81 
     82 	for (i = 2 * from; i < 2 * tlen; i += 2) {
     83 		x = dt[i];
     84 		x1 = dt[i + 1];
     85 		dt[i] = lower32(x, Zero) + tmp;
     86 		dt[i + 1] = lower32(x1, Zero) + tmp1;
     87 		tmp = upper32(x);
     88 		tmp1 = upper32(x1);
     89 	}
     90 }
     91 
     92 
     93 void
     94 conv_d16_to_i32(uint32_t *i32, double *d16, int64_t *tmp, int ilen)
     95 {
     96 	int i;
     97 	int64_t t, t1,		/* using int64_t and not uint64_t */
     98 		a, b, c, d;	/* because more efficient code is */
     99 				/* generated this way, and there  */
    100 				/* is no overflow  */
    101 	t1 = 0;
    102 	a = (int64_t)d16[0];
    103 	b = (int64_t)d16[1];
    104 	for (i = 0; i < ilen - 1; i++) {
    105 		c = (int64_t)d16[2 * i + 2];
    106 		t1 += a & 0xffffffff;
    107 		t = (a >> 32);
    108 		d = (int64_t)d16[2 * i + 3];
    109 		t1 += (b & 0xffff) << 16;
    110 		t += (b >> 16) + (t1 >> 32);
    111 		i32[i] = t1 & 0xffffffff;
    112 		t1 = t;
    113 		a = c;
    114 		b = d;
    115 	}
    116 	t1 += a & 0xffffffff;
    117 	t = (a >> 32);
    118 	t1 += (b & 0xffff) << 16;
    119 	i32[i] = t1 & 0xffffffff;
    120 }
    121 
    122 void
    123 conv_i32_to_d32(double *d32, uint32_t *i32, int len)
    124 {
    125 	int i;
    126 
    127 #pragma pipeloop(0)
    128 	for (i = 0; i < len; i++)
    129 		d32[i] = (double)(i32[i]);
    130 }
    131 
    132 
    133 void
    134 conv_i32_to_d16(double *d16, uint32_t *i32, int len)
    135 {
    136 	int i;
    137 	uint32_t a;
    138 
    139 #pragma pipeloop(0)
    140 	for (i = 0; i < len; i++) {
    141 		a = i32[i];
    142 		d16[2 * i] = (double)(a & 0xffff);
    143 		d16[2 * i + 1] = (double)(a >> 16);
    144 	}
    145 }
    146 
    147 #ifdef RF_INLINE_MACROS
    148 
    149 void
    150 i16_to_d16_and_d32x4(const double *,	/* 1/(2^16) */
    151 			const double *,	/* 2^16 */
    152 			const double *,	/* 0 */
    153 			double *,	/* result16 */
    154 			double *,	/* result32 */
    155 			float *);	/* source - should be unsigned int* */
    156 					/* converted to float* */
    157 
    158 #else
    159 
    160 
    161 static void
    162 i16_to_d16_and_d32x4(const double *dummy1,	/* 1/(2^16) */
    163 			const double *dummy2,	/* 2^16 */
    164 			const double *dummy3,	/* 0 */
    165 			double *result16,
    166 			double *result32,
    167 			float *src)	/* source - should be unsigned int* */
    168 					/* converted to float* */
    169 {
    170 	uint32_t *i32;
    171 	uint32_t a, b, c, d;
    172 
    173 	i32 = (uint32_t *)src;
    174 	a = i32[0];
    175 	b = i32[1];
    176 	c = i32[2];
    177 	d = i32[3];
    178 	result16[0] = (double)(a & 0xffff);
    179 	result16[1] = (double)(a >> 16);
    180 	result32[0] = (double)a;
    181 	result16[2] = (double)(b & 0xffff);
    182 	result16[3] = (double)(b >> 16);
    183 	result32[1] = (double)b;
    184 	result16[4] = (double)(c & 0xffff);
    185 	result16[5] = (double)(c >> 16);
    186 	result32[2] = (double)c;
    187 	result16[6] = (double)(d & 0xffff);
    188 	result16[7] = (double)(d >> 16);
    189 	result32[3] = (double)d;
    190 }
    191 
    192 #endif
    193 
    194 
    195 void
    196 conv_i32_to_d32_and_d16(double *d32, double *d16, uint32_t *i32, int len)
    197 {
    198 	int i;
    199 	uint32_t a;
    200 
    201 #pragma pipeloop(0)
    202 	for (i = 0; i < len - 3; i += 4) {
    203 		i16_to_d16_and_d32x4(&TwoToMinus16, &TwoTo16, &Zero,
    204 					&(d16[2*i]), &(d32[i]),
    205 					(float *)(&(i32[i])));
    206 	}
    207 	for (; i < len; i++) {
    208 		a = i32[i];
    209 		d32[i] = (double)(i32[i]);
    210 		d16[2 * i] = (double)(a & 0xffff);
    211 		d16[2 * i + 1] = (double)(a >> 16);
    212 	}
    213 }
    214 
    215 
    216 static void
    217 adjust_montf_result(uint32_t *i32, uint32_t *nint, int len)
    218 {
    219 	int64_t acc;
    220 	int i;
    221 
    222 	if (i32[len] > 0)
    223 		i = -1;
    224 	else {
    225 		for (i = len - 1; i >= 0; i--) {
    226 			if (i32[i] != nint[i]) break;
    227 		}
    228 	}
    229 	if ((i < 0) || (i32[i] > nint[i])) {
    230 		acc = 0;
    231 		for (i = 0; i < len; i++) {
    232 			acc = acc + (uint64_t)(i32[i]) - (uint64_t)(nint[i]);
    233 			i32[i] = acc & 0xffffffff;
    234 			acc = acc >> 32;
    235 		}
    236 	}
    237 }
    238 
    239 
    240 /*
    241  * the lengths of the input arrays should be at least the following:
    242  * result[nlen+1], dm1[nlen], dm2[2*nlen+1], dt[4*nlen+2], dn[nlen], nint[nlen]
    243  * all of them should be different from one another
    244  */
    245 void mont_mulf_noconv(uint32_t *result,
    246 			double *dm1, double *dm2, double *dt,
    247 			double *dn, uint32_t *nint,
    248 			int nlen, double dn0)
    249 {
    250 	int i, j, jj;
    251 	double digit, m2j, a, b;
    252 	double *pdm1, *pdm2, *pdn, *pdtj, pdn_0, pdm1_0;
    253 
    254 	pdm1 = &(dm1[0]);
    255 	pdm2 = &(dm2[0]);
    256 	pdn = &(dn[0]);
    257 	pdm2[2 * nlen] = Zero;
    258 
    259 	if (nlen != 16) {
    260 		for (i = 0; i < 4 * nlen + 2; i++)
    261 			dt[i] = Zero;
    262 		a = dt[0] = pdm1[0] * pdm2[0];
    263 		digit = mod(lower32(a, Zero) * dn0, TwoToMinus16, TwoTo16);
    264 
    265 		pdtj = &(dt[0]);
    266 		for (j = jj = 0; j < 2 * nlen; j++, jj++, pdtj++) {
    267 			m2j = pdm2[j];
    268 			a = pdtj[0] + pdn[0] * digit;
    269 			b = pdtj[1] + pdm1[0] * pdm2[j + 1] + a * TwoToMinus16;
    270 			pdtj[1] = b;
    271 
    272 #pragma pipeloop(0)
    273 			for (i = 1; i < nlen; i++) {
    274 				pdtj[2 * i] += pdm1[i] * m2j + pdn[i] * digit;
    275 			}
    276 			if (jj == 30) {
    277 				cleanup(dt, j / 2 + 1, 2 * nlen + 1);
    278 				jj = 0;
    279 			}
    280 
    281 			digit = mod(lower32(b, Zero) * dn0,
    282 				    TwoToMinus16, TwoTo16);
    283 		}
    284 	} else {
    285 		a = dt[0] = pdm1[0] * pdm2[0];
    286 
    287 		dt[65] = dt[64] = dt[63] = dt[62] = dt[61] = dt[60] =
    288 			dt[59] = dt[58] = dt[57] = dt[56] = dt[55] =
    289 			dt[54] = dt[53] = dt[52] = dt[51] = dt[50] =
    290 			dt[49] = dt[48] = dt[47] = dt[46] = dt[45] =
    291 			dt[44] = dt[43] = dt[42] = dt[41] = dt[40] =
    292 			dt[39] = dt[38] = dt[37] = dt[36] = dt[35] =
    293 			dt[34] = dt[33] = dt[32] = dt[31] = dt[30] =
    294 			dt[29] = dt[28] = dt[27] = dt[26] = dt[25] =
    295 			dt[24] = dt[23] = dt[22] = dt[21] = dt[20] =
    296 			dt[19] = dt[18] = dt[17] = dt[16] = dt[15] =
    297 			dt[14] = dt[13] = dt[12] = dt[11] = dt[10] =
    298 			dt[9] = dt[8] = dt[7] = dt[6] = dt[5] = dt[4] =
    299 			dt[3] = dt[2] = dt[1] = Zero;
    300 
    301 		pdn_0 = pdn[0];
    302 		pdm1_0 = pdm1[0];
    303 
    304 		digit = mod(lower32(a, Zero) * dn0, TwoToMinus16, TwoTo16);
    305 		pdtj = &(dt[0]);
    306 
    307 		for (j = 0; j < 32; j++, pdtj++) {
    308 
    309 			m2j = pdm2[j];
    310 			a = pdtj[0] + pdn_0 * digit;
    311 			b = pdtj[1] + pdm1_0 * pdm2[j + 1] + a * TwoToMinus16;
    312 			pdtj[1] = b;
    313 
    314 			pdtj[2] += pdm1[1] *m2j + pdn[1] * digit;
    315 			pdtj[4] += pdm1[2] *m2j + pdn[2] * digit;
    316 			pdtj[6] += pdm1[3] *m2j + pdn[3] * digit;
    317 			pdtj[8] += pdm1[4] *m2j + pdn[4] * digit;
    318 			pdtj[10] += pdm1[5] *m2j + pdn[5] * digit;
    319 			pdtj[12] += pdm1[6] *m2j + pdn[6] * digit;
    320 			pdtj[14] += pdm1[7] *m2j + pdn[7] * digit;
    321 			pdtj[16] += pdm1[8] *m2j + pdn[8] * digit;
    322 			pdtj[18] += pdm1[9] *m2j + pdn[9] * digit;
    323 			pdtj[20] += pdm1[10] *m2j + pdn[10] * digit;
    324 			pdtj[22] += pdm1[11] *m2j + pdn[11] * digit;
    325 			pdtj[24] += pdm1[12] *m2j + pdn[12] * digit;
    326 			pdtj[26] += pdm1[13] *m2j + pdn[13] * digit;
    327 			pdtj[28] += pdm1[14] *m2j + pdn[14] * digit;
    328 			pdtj[30] += pdm1[15] *m2j + pdn[15] * digit;
    329 			/* no need for cleenup, cannot overflow */
    330 			digit = mod(lower32(b, Zero) * dn0,
    331 				    TwoToMinus16, TwoTo16);
    332 		}
    333 	}
    334 
    335 	conv_d16_to_i32(result, dt + 2 * nlen, (int64_t *)dt, nlen + 1);
    336 	adjust_montf_result(result, nint, nlen);
    337 }
    338