```   1/*
2
3  fp_trig.c: floating-point math routines for the Linux-m68k
4  floating point emulator.
5
6  Copyright (c) 1998-1999 David Huggins-Daines / Roman Zippel.
7
8  I hereby give permission, free of charge, to copy, modify, and
9  redistribute this software, in source or binary form, provided that
10  the above copyright notice and the following disclaimer are included
11  in all such copies.
12
13  THIS SOFTWARE IS PROVIDED "AS IS", WITH ABSOLUTELY NO WARRANTY, REAL
14  OR IMPLIED.
15
16*/
17
18#include "fp_emu.h"
19
20static const struct fp_ext fp_one =
21{
22        .exp = 0x3fff,
23};
24
25extern struct fp_ext *fp_fadd(struct fp_ext *dest, const struct fp_ext *src);
26extern struct fp_ext *fp_fdiv(struct fp_ext *dest, const struct fp_ext *src);
27
28struct fp_ext *
29fp_fsqrt(struct fp_ext *dest, struct fp_ext *src)
30{
31        struct fp_ext tmp, src2;
32        int i, exp;
33
34        dprint(PINSTR, "fsqrt\n");
35
37
38        if (IS_ZERO(dest))
39                return dest;
40
41        if (dest->sign) {
42                fp_set_nan(dest);
43                return dest;
44        }
45        if (IS_INF(dest))
46                return dest;
47
48        /*
49         *               sqrt(m) * 2^(p)        , if e = 2*p
50         * sqrt(m*2^e) =
51         *               sqrt(2*m) * 2^(p)      , if e = 2*p + 1
52         *
53         * So we use the last bit of the exponent to decide whether to
54         * use the m or 2*m.
55         *
56         * Since only the fractional part of the mantissa is stored and
57         * the integer part is assumed to be one, we place a 1 or 2 into
58         * the fixed point representation.
59         */
60        exp = dest->exp;
61        dest->exp = 0x3FFF;
62        if (!(exp & 1))         /* lowest bit of exponent is set */
63                dest->exp++;
64        fp_copy_ext(&src2, dest);
65
66        /*
67         * The taylor row around a for sqrt(x) is:
68         *      sqrt(x) = sqrt(a) + 1/(2*sqrt(a))*(x-a) + R
69         * With a=1 this gives:
70         *      sqrt(x) = 1 + 1/2*(x-1)
71         *              = 1/2*(1+x)
72         */
74        dest->exp--;            /* * 1/2 */
75
76        /*
77         * We now apply the newton rule to the function
78         *      f(x) := x^2 - r
79         * which has a null point on x = sqrt(r).
80         *
81         * It gives:
82         *      x' := x - f(x)/f'(x)
83         *          = x - (x^2 -r)/(2*x)
84         *          = x - (x - r/x)/2
85         *          = (2*x - x + r/x)/2
86         *          = (x + r/x)/2
87         */
88        for (i = 0; i < 9; i++) {
89                fp_copy_ext(&tmp, &src2);
90
91                fp_fdiv(&tmp, dest);
93                dest->exp--;
94        }
95
96        dest->exp += (exp - 0x3FFF) / 2;
97
98        return dest;
99}
100
101struct fp_ext *
102fp_fetoxm1(struct fp_ext *dest, struct fp_ext *src)
103{
104        uprint("fetoxm1\n");
105
107
108        return dest;
109}
110
111struct fp_ext *
112fp_fetox(struct fp_ext *dest, struct fp_ext *src)
113{
114        uprint("fetox\n");
115
117
118        return dest;
119}
120
121struct fp_ext *
122fp_ftwotox(struct fp_ext *dest, struct fp_ext *src)
123{
124        uprint("ftwotox\n");
125
127
128        return dest;
129}
130
131struct fp_ext *
132fp_ftentox(struct fp_ext *dest, struct fp_ext *src)
133{
134        uprint("ftentox\n");
135
137
138        return dest;
139}
140
141struct fp_ext *
142fp_flogn(struct fp_ext *dest, struct fp_ext *src)
143{
144        uprint("flogn\n");
145
147
148        return dest;
149}
150
151struct fp_ext *
152fp_flognp1(struct fp_ext *dest, struct fp_ext *src)
153{
154        uprint("flognp1\n");
155
157
158        return dest;
159}
160
161struct fp_ext *
162fp_flog10(struct fp_ext *dest, struct fp_ext *src)
163{
164        uprint("flog10\n");
165
167
168        return dest;
169}
170
171struct fp_ext *
172fp_flog2(struct fp_ext *dest, struct fp_ext *src)
173{
174        uprint("flog2\n");
175
177
178        return dest;
179}
180
181struct fp_ext *
182fp_fgetexp(struct fp_ext *dest, struct fp_ext *src)
183{
184        dprint(PINSTR, "fgetexp\n");
185
187
188        if (IS_INF(dest)) {
189                fp_set_nan(dest);
190                return dest;
191        }
192        if (IS_ZERO(dest))
193                return dest;
194
195        fp_conv_long2ext(dest, (int)dest->exp - 0x3FFF);
196
197        fp_normalize_ext(dest);
198
199        return dest;
200}
201
202struct fp_ext *
203fp_fgetman(struct fp_ext *dest, struct fp_ext *src)
204{
205        dprint(PINSTR, "fgetman\n");
206
208
209        if (IS_ZERO(dest))
210                return dest;
211
212        if (IS_INF(dest))
213                return dest;
214
215        dest->exp = 0x3FFF;
216
217        return dest;
218}
219
220```
lxr.linux.no kindly hosted by Redpill Linpro AS, provider of Linux consulting and operations services since 1995.