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