|
|
1.1 root 1: /*
2: * linux/kernel/math/sqrt.c
3: *
1.1.1.3 ! root 4: * Copyright (C) 1991, 1992 Linus Torvalds
1.1 root 5: */
6:
7: /*
8: * simple and stupid temporary real fsqrt() routine
9: *
10: * There are probably better ways to do this, but this should work ok.
11: */
12:
13: #include <linux/math_emu.h>
14: #include <linux/sched.h>
15:
16: static void shift_right(int * c)
17: {
18: __asm__("shrl $1,12(%0) ; rcrl $1,8(%0) ; rcrl $1,4(%0) ; rcrl $1,(%0)"
19: ::"r" ((long) c));
20: }
21:
22: static int sqr64(unsigned long * a, unsigned long * b)
23: {
24: unsigned long tmp[4];
25:
26: __asm__("movl (%0),%%eax ; mull %%eax\n\t"
27: "movl %%eax,(%1) ; movl %%edx,4(%1)\n\t"
28: "movl 4(%0),%%eax ; mull %%eax\n\t"
29: "movl %%eax,8(%1) ; movl %%edx,12(%1)\n\t"
30: "movl (%0),%%eax ; mull 4(%0)\n\t"
31: "addl %%eax,%%eax ; adcl %%edx,%%edx\n\t"
32: "adcl $0,12(%1) ; addl %%eax,4(%1)\n\t"
33: "adcl %%edx,8(%1) ; adcl $0,12(%1)"
34: ::"b" ((long) a),"c" ((long) tmp)
35: :"ax","bx","cx","dx");
36: if (tmp[3] > b[3] ||
37: (tmp[3] == b[3] && (tmp[2] > b[2] ||
38: (tmp[2] == b[2] && (tmp[1] > b[1] ||
39: (tmp[1] == b[1] && tmp[0] > b[0]))))))
40: return 0;
41: return 1;
42: }
43:
44: void fsqrt(const temp_real * s, temp_real * d)
45: {
46: unsigned long src[4];
47: unsigned long res[2];
48: int exponent;
49: unsigned long mask, *c;
50: int i;
51:
52: exponent = s->exponent;
53: src[0] = src[1] = 0;
54: src[2] = s->a;
55: src[3] = s->b;
56: d->exponent = 0;
57: d->a = d->b = 0;
1.1.1.2 root 58: if (!exponent) /* fsqrt(0.0) = 0.0 */
1.1 root 59: return;
60: if (!src[2] && !src[3])
61: return;
62: if (exponent & 0x8000) {
63: send_sig(SIGFPE,current,0);
64: return;
65: }
66: if (exponent & 1) {
67: shift_right(src);
68: exponent++;
69: }
70: exponent >>= 1;
71: exponent += 0x1fff;
72: c = res + 2;
73: mask = 0;
74: for (i = 64 ; i > 0 ; i--) {
75: if (!(mask >>= 1)) {
76: c--;
77: mask = 0x80000000;
78: }
79: res[0] = d->a; res[1] = d->b;
80: *c |= mask;
81: if (sqr64(res,src)) {
82: d->a = res[0];
83: d->b = res[1];
84: }
85: }
86: if (!d->a && !d->b)
87: return;
88: while (!(d->b & 0x80000000)) {
89: __asm__("addl %%eax,%%eax ; adcl %%edx,%%edx"
90: :"=a" (d->a),"=d" (d->b)
91: :"0" (d->a),"1" (d->b));
92: exponent--;
93: }
94: d->exponent = exponent;
95: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.