ia64/xen-unstable

view extras/mini-os/lib/math.c @ 16838:945820bfedb6

minios: POSIX fixes
Fixes some functions which are POSIX. Also make them ifndef HAVE_LIBC.

Signed-off-by: Samuel Thibault <samuel.thibault@eu.citrix.com>
author Keir Fraser <keir.fraser@citrix.com>
date Tue Jan 22 14:20:22 2008 +0000 (2008-01-22)
parents e9e327c3e81b
children e6c3006fd9be
line source
1 /* -*- Mode:C; c-basic-offset:4; tab-width:4 -*-
2 ****************************************************************************
3 * (C) 2003 - Rolf Neugebauer - Intel Research Cambridge
4 ****************************************************************************
5 *
6 * File: math.c
7 * Author: Rolf Neugebauer (neugebar@dcs.gla.ac.uk)
8 * Changes:
9 *
10 * Date: Aug 2003
11 *
12 * Environment: Xen Minimal OS
13 * Description: Library functions for 64bit arith and other
14 * from freebsd, files in sys/libkern/ (qdivrem.c, etc)
15 *
16 ****************************************************************************
17 * $Id: c-insert.c,v 1.7 2002/11/08 16:04:34 rn Exp $
18 ****************************************************************************
19 *-
20 * Copyright (c) 1992, 1993
21 * The Regents of the University of California. All rights reserved.
22 *
23 * This software was developed by the Computer Systems Engineering group
24 * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
25 * contributed to Berkeley.
26 *
27 * Redistribution and use in source and binary forms, with or without
28 * modification, are permitted provided that the following conditions
29 * are met:
30 * 1. Redistributions of source code must retain the above copyright
31 * notice, this list of conditions and the following disclaimer.
32 * 2. Redistributions in binary form must reproduce the above copyright
33 * notice, this list of conditions and the following disclaimer in the
34 * documentation and/or other materials provided with the distribution.
35 * 3. All advertising materials mentioning features or use of this software
36 * must display the following acknowledgement:
37 * This product includes software developed by the University of
38 * California, Berkeley and its contributors.
39 * 4. Neither the name of the University nor the names of its contributors
40 * may be used to endorse or promote products derived from this software
41 * without specific prior written permission.
42 *
43 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
44 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
45 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
46 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
47 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
48 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
49 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
50 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
51 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
52 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
53 * SUCH DAMAGE.
54 *
55 * $FreeBSD: src/sys/libkern/divdi3.c,v 1.6 1999/08/28 00:46:31 peter Exp $
56 */
58 #include <types.h>
59 #include <lib.h>
60 #include <time.h>
62 /* On ia64 these functions lead to crashes. These are replaced by
63 * assembler functions. */
64 #if !defined(__ia64__)
66 /*
67 * Depending on the desired operation, we view a `long long' (aka quad_t) in
68 * one or more of the following formats.
69 */
70 union uu {
71 s64 q; /* as a (signed) quad */
72 s64 uq; /* as an unsigned quad */
73 long sl[2]; /* as two signed longs */
74 unsigned long ul[2]; /* as two unsigned longs */
75 };
76 /* XXX RN: Yuck hardcoded endianess :) */
77 #define _QUAD_HIGHWORD 1
78 #define _QUAD_LOWWORD 0
79 /*
80 * Define high and low longwords.
81 */
82 #define H _QUAD_HIGHWORD
83 #define L _QUAD_LOWWORD
85 /*
86 * Total number of bits in a quad_t and in the pieces that make it up.
87 * These are used for shifting, and also below for halfword extraction
88 * and assembly.
89 */
90 #ifndef HAVE_LIBC
91 #define CHAR_BIT 8 /* number of bits in a char */
92 #endif
93 #define QUAD_BITS (sizeof(s64) * CHAR_BIT)
94 #define LONG_BITS (sizeof(long) * CHAR_BIT)
95 #define HALF_BITS (sizeof(long) * CHAR_BIT / 2)
97 /*
98 * Extract high and low shortwords from longword, and move low shortword of
99 * longword to upper half of long, i.e., produce the upper longword of
100 * ((quad_t)(x) << (number_of_bits_in_long/2)). (`x' must actually be u_long.)
101 *
102 * These are used in the multiply code, to split a longword into upper
103 * and lower halves, and to reassemble a product as a quad_t, shifted left
104 * (sizeof(long)*CHAR_BIT/2).
105 */
106 #define HHALF(x) ((x) >> HALF_BITS)
107 #define LHALF(x) ((x) & ((1UL << HALF_BITS) - 1))
108 #define LHUP(x) ((x) << HALF_BITS)
110 /*
111 * Multiprecision divide. This algorithm is from Knuth vol. 2 (2nd ed),
112 * section 4.3.1, pp. 257--259.
113 */
114 #define B (1UL << HALF_BITS) /* digit base */
116 /* Combine two `digits' to make a single two-digit number. */
117 #define COMBINE(a, b) (((u_long)(a) << HALF_BITS) | (b))
119 /* select a type for digits in base B: use unsigned short if they fit */
120 #if ULONG_MAX == 0xffffffff && USHRT_MAX >= 0xffff
121 typedef unsigned short digit;
122 #else
123 typedef u_long digit;
124 #endif
127 /*
128 * Shift p[0]..p[len] left `sh' bits, ignoring any bits that
129 * `fall out' the left (there never will be any such anyway).
130 * We may assume len >= 0. NOTE THAT THIS WRITES len+1 DIGITS.
131 */
132 static void
133 shl(register digit *p, register int len, register int sh)
134 {
135 register int i;
137 for (i = 0; i < len; i++)
138 p[i] = LHALF(p[i] << sh) | (p[i + 1] >> (HALF_BITS - sh));
139 p[i] = LHALF(p[i] << sh);
140 }
142 /*
143 * __qdivrem(u, v, rem) returns u/v and, optionally, sets *rem to u%v.
144 *
145 * We do this in base 2-sup-HALF_BITS, so that all intermediate products
146 * fit within u_long. As a consequence, the maximum length dividend and
147 * divisor are 4 `digits' in this base (they are shorter if they have
148 * leading zeros).
149 */
150 u64
151 __qdivrem(u64 uq, u64 vq, u64 *arq)
152 {
153 union uu tmp;
154 digit *u, *v, *q;
155 register digit v1, v2;
156 u_long qhat, rhat, t;
157 int m, n, d, j, i;
158 digit uspace[5], vspace[5], qspace[5];
160 /*
161 * Take care of special cases: divide by zero, and u < v.
162 */
163 if (vq == 0) {
164 /* divide by zero. */
165 static volatile const unsigned int zero = 0;
167 tmp.ul[H] = tmp.ul[L] = 1 / zero;
168 if (arq)
169 *arq = uq;
170 return (tmp.q);
171 }
172 if (uq < vq) {
173 if (arq)
174 *arq = uq;
175 return (0);
176 }
177 u = &uspace[0];
178 v = &vspace[0];
179 q = &qspace[0];
181 /*
182 * Break dividend and divisor into digits in base B, then
183 * count leading zeros to determine m and n. When done, we
184 * will have:
185 * u = (u[1]u[2]...u[m+n]) sub B
186 * v = (v[1]v[2]...v[n]) sub B
187 * v[1] != 0
188 * 1 < n <= 4 (if n = 1, we use a different division algorithm)
189 * m >= 0 (otherwise u < v, which we already checked)
190 * m + n = 4
191 * and thus
192 * m = 4 - n <= 2
193 */
194 tmp.uq = uq;
195 u[0] = 0;
196 u[1] = HHALF(tmp.ul[H]);
197 u[2] = LHALF(tmp.ul[H]);
198 u[3] = HHALF(tmp.ul[L]);
199 u[4] = LHALF(tmp.ul[L]);
200 tmp.uq = vq;
201 v[1] = HHALF(tmp.ul[H]);
202 v[2] = LHALF(tmp.ul[H]);
203 v[3] = HHALF(tmp.ul[L]);
204 v[4] = LHALF(tmp.ul[L]);
205 for (n = 4; v[1] == 0; v++) {
206 if (--n == 1) {
207 u_long rbj; /* r*B+u[j] (not root boy jim) */
208 digit q1, q2, q3, q4;
210 /*
211 * Change of plan, per exercise 16.
212 * r = 0;
213 * for j = 1..4:
214 * q[j] = floor((r*B + u[j]) / v),
215 * r = (r*B + u[j]) % v;
216 * We unroll this completely here.
217 */
218 t = v[2]; /* nonzero, by definition */
219 q1 = u[1] / t;
220 rbj = COMBINE(u[1] % t, u[2]);
221 q2 = rbj / t;
222 rbj = COMBINE(rbj % t, u[3]);
223 q3 = rbj / t;
224 rbj = COMBINE(rbj % t, u[4]);
225 q4 = rbj / t;
226 if (arq)
227 *arq = rbj % t;
228 tmp.ul[H] = COMBINE(q1, q2);
229 tmp.ul[L] = COMBINE(q3, q4);
230 return (tmp.q);
231 }
232 }
234 /*
235 * By adjusting q once we determine m, we can guarantee that
236 * there is a complete four-digit quotient at &qspace[1] when
237 * we finally stop.
238 */
239 for (m = 4 - n; u[1] == 0; u++)
240 m--;
241 for (i = 4 - m; --i >= 0;)
242 q[i] = 0;
243 q += 4 - m;
245 /*
246 * Here we run Program D, translated from MIX to C and acquiring
247 * a few minor changes.
248 *
249 * D1: choose multiplier 1 << d to ensure v[1] >= B/2.
250 */
251 d = 0;
252 for (t = v[1]; t < B / 2; t <<= 1)
253 d++;
254 if (d > 0) {
255 shl(&u[0], m + n, d); /* u <<= d */
256 shl(&v[1], n - 1, d); /* v <<= d */
257 }
258 /*
259 * D2: j = 0.
260 */
261 j = 0;
262 v1 = v[1]; /* for D3 -- note that v[1..n] are constant */
263 v2 = v[2]; /* for D3 */
264 do {
265 register digit uj0, uj1, uj2;
267 /*
268 * D3: Calculate qhat (\^q, in TeX notation).
269 * Let qhat = min((u[j]*B + u[j+1])/v[1], B-1), and
270 * let rhat = (u[j]*B + u[j+1]) mod v[1].
271 * While rhat < B and v[2]*qhat > rhat*B+u[j+2],
272 * decrement qhat and increase rhat correspondingly.
273 * Note that if rhat >= B, v[2]*qhat < rhat*B.
274 */
275 uj0 = u[j + 0]; /* for D3 only -- note that u[j+...] change */
276 uj1 = u[j + 1]; /* for D3 only */
277 uj2 = u[j + 2]; /* for D3 only */
278 if (uj0 == v1) {
279 qhat = B;
280 rhat = uj1;
281 goto qhat_too_big;
282 } else {
283 u_long nn = COMBINE(uj0, uj1);
284 qhat = nn / v1;
285 rhat = nn % v1;
286 }
287 while (v2 * qhat > COMBINE(rhat, uj2)) {
288 qhat_too_big:
289 qhat--;
290 if ((rhat += v1) >= B)
291 break;
292 }
293 /*
294 * D4: Multiply and subtract.
295 * The variable `t' holds any borrows across the loop.
296 * We split this up so that we do not require v[0] = 0,
297 * and to eliminate a final special case.
298 */
299 for (t = 0, i = n; i > 0; i--) {
300 t = u[i + j] - v[i] * qhat - t;
301 u[i + j] = LHALF(t);
302 t = (B - HHALF(t)) & (B - 1);
303 }
304 t = u[j] - t;
305 u[j] = LHALF(t);
306 /*
307 * D5: test remainder.
308 * There is a borrow if and only if HHALF(t) is nonzero;
309 * in that (rare) case, qhat was too large (by exactly 1).
310 * Fix it by adding v[1..n] to u[j..j+n].
311 */
312 if (HHALF(t)) {
313 qhat--;
314 for (t = 0, i = n; i > 0; i--) { /* D6: add back. */
315 t += u[i + j] + v[i];
316 u[i + j] = LHALF(t);
317 t = HHALF(t);
318 }
319 u[j] = LHALF(u[j] + t);
320 }
321 q[j] = qhat;
322 } while (++j <= m); /* D7: loop on j. */
324 /*
325 * If caller wants the remainder, we have to calculate it as
326 * u[m..m+n] >> d (this is at most n digits and thus fits in
327 * u[m+1..m+n], but we may need more source digits).
328 */
329 if (arq) {
330 if (d) {
331 for (i = m + n; i > m; --i)
332 u[i] = (u[i] >> d) |
333 LHALF(u[i - 1] << (HALF_BITS - d));
334 u[i] = 0;
335 }
336 tmp.ul[H] = COMBINE(uspace[1], uspace[2]);
337 tmp.ul[L] = COMBINE(uspace[3], uspace[4]);
338 *arq = tmp.q;
339 }
341 tmp.ul[H] = COMBINE(qspace[1], qspace[2]);
342 tmp.ul[L] = COMBINE(qspace[3], qspace[4]);
343 return (tmp.q);
344 }
347 /*
348 * Divide two signed quads.
349 * ??? if -1/2 should produce -1 on this machine, this code is wrong
350 */
351 s64
352 __divdi3(s64 a, s64 b)
353 {
354 u64 ua, ub, uq;
355 int neg;
357 if (a < 0)
358 ua = -(u64)a, neg = 1;
359 else
360 ua = a, neg = 0;
361 if (b < 0)
362 ub = -(u64)b, neg ^= 1;
363 else
364 ub = b;
365 uq = __qdivrem(ua, ub, (u64 *)0);
366 return (neg ? -uq : uq);
367 }
369 /*
370 * Divide two unsigned quads.
371 */
372 u64
373 __udivdi3(u64 a, u64 b)
374 {
375 return (__qdivrem(a, b, (u64 *)0));
376 }
379 /*
380 * Return remainder after dividing two unsigned quads.
381 */
382 u_quad_t
383 __umoddi3(u_quad_t a, u_quad_t b)
384 {
385 u_quad_t r;
387 (void)__qdivrem(a, b, &r);
388 return (r);
389 }
391 #endif /* !defined(__ia64__) */
393 #ifndef HAVE_LIBC
394 /* Should be random enough for our uses */
395 int rand(void)
396 {
397 static unsigned int previous;
398 struct timeval tv;
399 gettimeofday(&tv, NULL);
400 previous += tv.tv_sec + tv.tv_usec;
401 previous *= RAND_MIX;
402 return previous;
403 }
404 #endif