OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libgcc/] [config/] [libbid/] [bid_sqrt_macros.h] - Blame information for rev 734

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 734 jeremybenn
/* Copyright (C) 2007, 2009  Free Software Foundation, Inc.
2
 
3
This file is part of GCC.
4
 
5
GCC is free software; you can redistribute it and/or modify it under
6
the terms of the GNU General Public License as published by the Free
7
Software Foundation; either version 3, or (at your option) any later
8
version.
9
 
10
GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11
WARRANTY; without even the implied warranty of MERCHANTABILITY or
12
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
13
for more details.
14
 
15
Under Section 7 of GPL version 3, you are granted additional
16
permissions described in the GCC Runtime Library Exception, version
17
3.1, as published by the Free Software Foundation.
18
 
19
You should have received a copy of the GNU General Public License and
20
a copy of the GCC Runtime Library Exception along with this program;
21
see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
22
<http://www.gnu.org/licenses/>.  */
23
 
24
#ifndef _SQRT_MACROS_H_
25
#define _SQRT_MACROS_H_
26
 
27
#define FENCE __fence
28
 
29
#if DOUBLE_EXTENDED_ON
30
 
31
extern BINARY80 SQRT80 (BINARY80);
32
 
33
 
34
__BID_INLINE__ UINT64
35
short_sqrt128 (UINT128 A10) {
36
  BINARY80 lx, ly, l64;
37
  int_float f64;
38
 
39
  // 2^64
40
  f64.i = 0x5f800000;
41
  l64 = (BINARY80) f64.d;
42
  lx = (BINARY80) A10.w[1] * l64 + (BINARY80) A10.w[0];
43
  ly = SQRT80 (lx);
44
  return (UINT64) ly;
45
}
46
 
47
 
48
__BID_INLINE__ void
49
long_sqrt128 (UINT128 * pCS, UINT256 C256) {
50
  UINT256 C4;
51
  UINT128 CS;
52
  UINT64 X;
53
  SINT64 SE;
54
  BINARY80 l64, lm64, l128, lxL, lx, ly, lS, lSH, lSL, lE, l3, l2,
55
    l1, l0, lp, lCl;
56
  int_float fx, f64, fm64;
57
  int *ple = (int *) &lx;
58
 
59
  // 2^64
60
  f64.i = 0x5f800000;
61
  l64 = (BINARY80) f64.d;
62
 
63
  l128 = l64 * l64;
64
  lx = l3 = (BINARY80) C256.w[3] * l64 * l128;
65
  l2 = (BINARY80) C256.w[2] * l128;
66
  lx = FENCE (lx + l2);
67
  l1 = (BINARY80) C256.w[1] * l64;
68
  lx = FENCE (lx + l1);
69
  l0 = (BINARY80) C256.w[0];
70
  lx = FENCE (lx + l0);
71
  // sqrt(C256)
72
  lS = SQRT80 (lx);
73
 
74
  // get coefficient
75
  // 2^(-64)
76
  fm64.i = 0x1f800000;
77
  lm64 = (BINARY80) fm64.d;
78
  CS.w[1] = (UINT64) (lS * lm64);
79
  CS.w[0] = (UINT64) (lS - (BINARY80) CS.w[1] * l64);
80
 
81
  ///////////////////////////////////////
82
  //  CAUTION!
83
  //  little endian code only
84
  //  add solution for big endian
85
  //////////////////////////////////////
86
  lSH = lS;
87
  *((UINT64 *) & lSH) &= 0xffffffff00000000ull;
88
 
89
  // correction for C256 rounding
90
  lCl = FENCE (l3 - lx);
91
  lCl = FENCE (lCl + l2);
92
  lCl = FENCE (lCl + l1);
93
  lCl = FENCE (lCl + l0);
94
 
95
  lSL = lS - lSH;
96
 
97
  //////////////////////////////////////////
98
  //   Watch for compiler re-ordering
99
  //
100
  /////////////////////////////////////////
101
  // C256-S^2
102
  lxL = FENCE (lx - lSH * lSH);
103
  lp = lSH * lSL;
104
  lp += lp;
105
  lxL = FENCE (lxL - lp);
106
  lSL *= lSL;
107
  lxL = FENCE (lxL - lSL);
108
  lCl += lxL;
109
 
110
  // correction term
111
  lE = lCl / (lS + lS);
112
 
113
  // get low part of coefficient
114
  X = CS.w[0];
115
  if (lCl >= 0) {
116
    SE = (SINT64) (lE);
117
    CS.w[0] += SE;
118
    if (CS.w[0] < X)
119
      CS.w[1]++;
120
  } else {
121
    SE = (SINT64) (-lE);
122
    CS.w[0] -= SE;
123
    if (CS.w[0] > X)
124
      CS.w[1]--;
125
  }
126
 
127
  pCS->w[0] = CS.w[0];
128
  pCS->w[1] = CS.w[1];
129
}
130
 
131
#else
132
 
133
extern double sqrt (double);
134
 
135
__BID_INLINE__ UINT64
136
short_sqrt128 (UINT128 A10) {
137
  UINT256 ARS, ARS0, AE0, AE, S;
138
 
139
  UINT64 MY, ES, CY;
140
  double lx, l64;
141
  int_double f64, ly;
142
  int ey, k;
143
 
144
  // 2^64
145
  f64.i = 0x43f0000000000000ull;
146
  l64 = f64.d;
147
  lx = (double) A10.w[1] * l64 + (double) A10.w[0];
148
  ly.d = 1.0 / sqrt (lx);
149
 
150
  MY = (ly.i & 0x000fffffffffffffull) | 0x0010000000000000ull;
151
  ey = 0x3ff - (ly.i >> 52);
152
 
153
  // A10*RS^2
154
  __mul_64x128_to_192 (ARS0, MY, A10);
155
  __mul_64x192_to_256 (ARS, MY, ARS0);
156
 
157
  // shr by 2*ey+40, to get a 64-bit value
158
  k = (ey << 1) + 104 - 64;
159
  if (k >= 128) {
160
    if (k > 128)
161
      ES = (ARS.w[2] >> (k - 128)) | (ARS.w[3] << (192 - k));
162
    else
163
      ES = ARS.w[2];
164
  } else {
165
    if (k >= 64) {
166
      ARS.w[0] = ARS.w[1];
167
      ARS.w[1] = ARS.w[2];
168
      k -= 64;
169
    }
170
    if (k) {
171
      __shr_128 (ARS, ARS, k);
172
    }
173
    ES = ARS.w[0];
174
  }
175
 
176
  ES = ((SINT64) ES) >> 1;
177
 
178
  if (((SINT64) ES) < 0) {
179
    ES = -ES;
180
 
181
    // A*RS*eps (scaled by 2^64)
182
    __mul_64x192_to_256 (AE0, ES, ARS0);
183
 
184
    AE.w[0] = AE0.w[1];
185
    AE.w[1] = AE0.w[2];
186
    AE.w[2] = AE0.w[3];
187
 
188
    __add_carry_out (S.w[0], CY, ARS0.w[0], AE.w[0]);
189
    __add_carry_in_out (S.w[1], CY, ARS0.w[1], AE.w[1], CY);
190
    S.w[2] = ARS0.w[2] + AE.w[2] + CY;
191
  } else {
192
    // A*RS*eps (scaled by 2^64)
193
    __mul_64x192_to_256 (AE0, ES, ARS0);
194
 
195
    AE.w[0] = AE0.w[1];
196
    AE.w[1] = AE0.w[2];
197
    AE.w[2] = AE0.w[3];
198
 
199
    __sub_borrow_out (S.w[0], CY, ARS0.w[0], AE.w[0]);
200
    __sub_borrow_in_out (S.w[1], CY, ARS0.w[1], AE.w[1], CY);
201
    S.w[2] = ARS0.w[2] - AE.w[2] - CY;
202
  }
203
 
204
  k = ey + 51;
205
 
206
  if (k >= 64) {
207
    if (k >= 128) {
208
      S.w[0] = S.w[2];
209
      S.w[1] = 0;
210
      k -= 128;
211
    } else {
212
      S.w[0] = S.w[1];
213
      S.w[1] = S.w[2];
214
    }
215
    k -= 64;
216
  }
217
  if (k) {
218
    __shr_128 (S, S, k);
219
  }
220
 
221
 
222
  return (UINT64) ((S.w[0] + 1) >> 1);
223
 
224
}
225
 
226
 
227
 
228
__BID_INLINE__ void
229
long_sqrt128 (UINT128 * pCS, UINT256 C256) {
230
  UINT512 ARS0, ARS;
231
  UINT256 ARS00, AE, AE2, S;
232
  UINT128 ES, ES2, ARS1;
233
  UINT64 ES32, CY, MY;
234
  double l64, l128, lx, l2, l1, l0;
235
  int_double f64, ly;
236
  int ey, k, k2;
237
 
238
  // 2^64
239
  f64.i = 0x43f0000000000000ull;
240
  l64 = f64.d;
241
 
242
  l128 = l64 * l64;
243
  lx = (double) C256.w[3] * l64 * l128;
244
  l2 = (double) C256.w[2] * l128;
245
  lx = FENCE (lx + l2);
246
  l1 = (double) C256.w[1] * l64;
247
  lx = FENCE (lx + l1);
248
  l0 = (double) C256.w[0];
249
  lx = FENCE (lx + l0);
250
  // sqrt(C256)
251
  ly.d = 1.0 / sqrt (lx);
252
 
253
  MY = (ly.i & 0x000fffffffffffffull) | 0x0010000000000000ull;
254
  ey = 0x3ff - (ly.i >> 52);
255
 
256
  // A10*RS^2, scaled by 2^(2*ey+104)
257
  __mul_64x256_to_320 (ARS0, MY, C256);
258
  __mul_64x320_to_384 (ARS, MY, ARS0);
259
 
260
  // shr by k=(2*ey+104)-128
261
  // expect k is in the range (192, 256) if result in [10^33, 10^34)
262
  // apply an additional signed shift by 1 at the same time (to get eps=eps0/2)
263
  k = (ey << 1) + 104 - 128 - 192;
264
  k2 = 64 - k;
265
  ES.w[0] = (ARS.w[3] >> (k + 1)) | (ARS.w[4] << (k2 - 1));
266
  ES.w[1] = (ARS.w[4] >> k) | (ARS.w[5] << k2);
267
  ES.w[1] = ((SINT64) ES.w[1]) >> 1;
268
 
269
  // A*RS >> 192 (for error term computation)
270
  ARS1.w[0] = ARS0.w[3];
271
  ARS1.w[1] = ARS0.w[4];
272
 
273
  // A*RS>>64
274
  ARS00.w[0] = ARS0.w[1];
275
  ARS00.w[1] = ARS0.w[2];
276
  ARS00.w[2] = ARS0.w[3];
277
  ARS00.w[3] = ARS0.w[4];
278
 
279
  if (((SINT64) ES.w[1]) < 0) {
280
    ES.w[0] = -ES.w[0];
281
    ES.w[1] = -ES.w[1];
282
    if (ES.w[0])
283
      ES.w[1]--;
284
 
285
    // A*RS*eps 
286
    __mul_128x128_to_256 (AE, ES, ARS1);
287
 
288
    __add_carry_out (S.w[0], CY, ARS00.w[0], AE.w[0]);
289
    __add_carry_in_out (S.w[1], CY, ARS00.w[1], AE.w[1], CY);
290
    __add_carry_in_out (S.w[2], CY, ARS00.w[2], AE.w[2], CY);
291
    S.w[3] = ARS00.w[3] + AE.w[3] + CY;
292
  } else {
293
    // A*RS*eps 
294
    __mul_128x128_to_256 (AE, ES, ARS1);
295
 
296
    __sub_borrow_out (S.w[0], CY, ARS00.w[0], AE.w[0]);
297
    __sub_borrow_in_out (S.w[1], CY, ARS00.w[1], AE.w[1], CY);
298
    __sub_borrow_in_out (S.w[2], CY, ARS00.w[2], AE.w[2], CY);
299
    S.w[3] = ARS00.w[3] - AE.w[3] - CY;
300
  }
301
 
302
  // 3/2*eps^2, scaled by 2^128
303
  ES32 = ES.w[1] + (ES.w[1] >> 1);
304
  __mul_64x64_to_128 (ES2, ES32, ES.w[1]);
305
  // A*RS*3/2*eps^2
306
  __mul_128x128_to_256 (AE2, ES2, ARS1);
307
 
308
  // result, scaled by 2^(ey+52-64)
309
  __add_carry_out (S.w[0], CY, S.w[0], AE2.w[0]);
310
  __add_carry_in_out (S.w[1], CY, S.w[1], AE2.w[1], CY);
311
  __add_carry_in_out (S.w[2], CY, S.w[2], AE2.w[2], CY);
312
  S.w[3] = S.w[3] + AE2.w[3] + CY;
313
 
314
  // k in (0, 64)
315
  k = ey + 51 - 128;
316
  k2 = 64 - k;
317
  S.w[0] = (S.w[1] >> k) | (S.w[2] << k2);
318
  S.w[1] = (S.w[2] >> k) | (S.w[3] << k2);
319
 
320
  // round to nearest
321
  S.w[0]++;
322
  if (!S.w[0])
323
    S.w[1]++;
324
 
325
  pCS->w[0] = (S.w[1] << 63) | (S.w[0] >> 1);
326
  pCS->w[1] = S.w[1] >> 1;
327
 
328
}
329
 
330
#endif
331
#endif

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.