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
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 }
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--) {
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 }
|
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 (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21 /*
22 * Copyright 2009 Sun Microsystems, Inc. All rights reserved.
23 * Use is subject to license terms.
24 */
25
26 /*
27 * If compiled without -DRF_INLINE_MACROS then needs -lm at link time
28 * If compiled with -DRF_INLINE_MACROS then needs conv.il at compile time
29 * (i.e. cc <compiler_flags> -DRF_INLINE_MACROS conv.il mont_mulf.c )
30 */
31
32 #include <sys/types.h>
33 #include <math.h>
34
35 static const double TwoTo16 = 65536.0;
36 static const double TwoToMinus16 = 1.0/65536.0;
37 static const double Zero = 0.0;
38 static const double TwoTo32 = 65536.0 * 65536.0;
39 static const double TwoToMinus32 = 1.0 / (65536.0 * 65536.0);
40
41 #ifdef RF_INLINE_MACROS
42
43 double upper32(double);
44 double lower32(double, double);
45 double mod(double, double, double);
46
47 #else
48
49 static double
74 int i;
75 double tmp, tmp1, x, x1;
76
77 tmp = tmp1 = Zero;
78
79 for (i = 2 * from; i < 2 * tlen; i += 2) {
80 x = dt[i];
81 x1 = dt[i + 1];
82 dt[i] = lower32(x, Zero) + tmp;
83 dt[i + 1] = lower32(x1, Zero) + tmp1;
84 tmp = upper32(x);
85 tmp1 = upper32(x1);
86 }
87 }
88
89
90 void
91 conv_d16_to_i32(uint32_t *i32, double *d16, int64_t *tmp, int ilen)
92 {
93 int i;
94 int64_t t, t1, /* Using int64_t and not uint64_t */
95 a, b, c, d; /* because more efficient code is */
96 /* generated this way, and there */
97 /* is no overflow. */
98 t1 = 0;
99 a = (int64_t)d16[0];
100 b = (int64_t)d16[1];
101 for (i = 0; i < ilen - 1; i++) {
102 c = (int64_t)d16[2 * i + 2];
103 t1 += a & 0xffffffff;
104 t = (a >> 32);
105 d = (int64_t)d16[2 * i + 3];
106 t1 += (b & 0xffff) << 16;
107 t += (b >> 16) + (t1 >> 32);
108 i32[i] = t1 & 0xffffffff;
109 t1 = t;
110 a = c;
111 b = d;
112 }
113 t1 += a & 0xffffffff;
114 t = (a >> 32);
115 t1 += (b & 0xffff) << 16;
116 i32[i] = t1 & 0xffffffff;
117 }
181 result16[4] = (double)(c & 0xffff);
182 result16[5] = (double)(c >> 16);
183 result32[2] = (double)c;
184 result16[6] = (double)(d & 0xffff);
185 result16[7] = (double)(d >> 16);
186 result32[3] = (double)d;
187 }
188
189 #endif
190
191
192 void
193 conv_i32_to_d32_and_d16(double *d32, double *d16, uint32_t *i32, int len)
194 {
195 int i;
196 uint32_t a;
197
198 #pragma pipeloop(0)
199 for (i = 0; i < len - 3; i += 4) {
200 i16_to_d16_and_d32x4(&TwoToMinus16, &TwoTo16, &Zero,
201 &(d16[2*i]), &(d32[i]), (float *)(&(i32[i])));
202 }
203 for (; i < len; i++) {
204 a = i32[i];
205 d32[i] = (double)(i32[i]);
206 d16[2 * i] = (double)(a & 0xffff);
207 d16[2 * i + 1] = (double)(a >> 16);
208 }
209 }
210
211
212 static void
213 adjust_montf_result(uint32_t *i32, uint32_t *nint, int len)
214 {
215 int64_t acc;
216 int i;
217
218 if (i32[len] > 0)
219 i = -1;
220 else {
221 for (i = len - 1; i >= 0; i--) {
305 m2j = pdm2[j];
306 a = pdtj[0] + pdn_0 * digit;
307 b = pdtj[1] + pdm1_0 * pdm2[j + 1] + a * TwoToMinus16;
308 pdtj[1] = b;
309
310 pdtj[2] += pdm1[1] *m2j + pdn[1] * digit;
311 pdtj[4] += pdm1[2] *m2j + pdn[2] * digit;
312 pdtj[6] += pdm1[3] *m2j + pdn[3] * digit;
313 pdtj[8] += pdm1[4] *m2j + pdn[4] * digit;
314 pdtj[10] += pdm1[5] *m2j + pdn[5] * digit;
315 pdtj[12] += pdm1[6] *m2j + pdn[6] * digit;
316 pdtj[14] += pdm1[7] *m2j + pdn[7] * digit;
317 pdtj[16] += pdm1[8] *m2j + pdn[8] * digit;
318 pdtj[18] += pdm1[9] *m2j + pdn[9] * digit;
319 pdtj[20] += pdm1[10] *m2j + pdn[10] * digit;
320 pdtj[22] += pdm1[11] *m2j + pdn[11] * digit;
321 pdtj[24] += pdm1[12] *m2j + pdn[12] * digit;
322 pdtj[26] += pdm1[13] *m2j + pdn[13] * digit;
323 pdtj[28] += pdm1[14] *m2j + pdn[14] * digit;
324 pdtj[30] += pdm1[15] *m2j + pdn[15] * digit;
325 /* no need for cleanup, cannot overflow */
326 digit = mod(lower32(b, Zero) * dn0,
327 TwoToMinus16, TwoTo16);
328 }
329 }
330
331 conv_d16_to_i32(result, dt + 2 * nlen, (int64_t *)dt, nlen + 1);
332 adjust_montf_result(result, nint, nlen);
333 }
|