![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to fpu_trig.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Isaki's NoNo m68k/m88k emulator] / nono / fpe|
Return to fpu_trig.c CVS log | Up to [Isaki's NoNo m68k/m88k emulator] / nono / fpe |
1.1.1.2 ! root 1: /* $NetBSD: fpu_trig.c,v 1.18 2017/01/16 12:05:40 isaki Exp $ */
1.1 root 2:
3: /*
4: * Copyright (c) 1995 Ken Nakata
5: * All rights reserved.
6: *
7: * Redistribution and use in source and binary forms, with or without
8: * modification, are permitted provided that the following conditions
9: * are met:
10: * 1. Redistributions of source code must retain the above copyright
11: * notice, this list of conditions and the following disclaimer.
12: * 2. Redistributions in binary form must reproduce the above copyright
13: * notice, this list of conditions and the following disclaimer in the
14: * documentation and/or other materials provided with the distribution.
15: * 3. Neither the name of the author nor the names of its contributors
16: * may be used to endorse or promote products derived from this software
17: * without specific prior written permission.
18: *
19: * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
20: * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21: * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22: * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
23: * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
24: * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
25: * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
26: * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
27: * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
28: * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
29: * SUCH DAMAGE.
30: *
31: * @(#)fpu_trig.c 10/24/95
32: */
33:
34: /*
35: * Copyright (c) 2011 Tetsuya Isaki. All rights reserved.
36: *
37: * Redistribution and use in source and binary forms, with or without
38: * modification, are permitted provided that the following conditions
39: * are met:
40: * 1. Redistributions of source code must retain the above copyright
41: * notice, this list of conditions and the following disclaimer.
42: * 2. Redistributions in binary form must reproduce the above copyright
43: * notice, this list of conditions and the following disclaimer in the
44: * documentation and/or other materials provided with the distribution.
45: *
46: * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
47: * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
48: * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
49: * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
50: * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
51: * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
52: * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
53: * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
54: * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
55: * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
56: * SUCH DAMAGE.
57: */
58:
59: #include "fpu_emulate.h"
60:
61: /*
62: * arccos(x) = pi/2 - arcsin(x)
63: */
64: struct fpn *
65: fpu_acos(struct fpemu *fe)
66: {
67: struct fpn *r;
68:
69: if (ISNAN(&fe->fe_f2))
70: return &fe->fe_f2;
71: if (ISINF(&fe->fe_f2))
72: return fpu_newnan(fe);
73:
74: r = fpu_asin(fe);
75: if (ISNAN(r)) {
76: return r;
77: }
78: CPYFPN(&fe->fe_f2, r);
79:
80: /* pi/2 - asin(x) */
81: fpu_const(&fe->fe_f1, FPU_CONST_PI);
82: fe->fe_f1.fp_exp--;
83: fe->fe_f2.fp_sign = !fe->fe_f2.fp_sign;
84: r = fpu_add(fe);
85:
86: /*
87: * small hack
88: */
89: r->fp_sign = 0;
90:
91: return r;
92: }
93:
94: /*
95: * x
96: * arcsin(x) = arctan(---------------)
97: * sqrt(1 - x^2)
98: */
99: struct fpn *
100: fpu_asin(struct fpemu *fe)
101: {
102: struct fpn x;
103: struct fpn *r;
104:
105: if (ISNAN(&fe->fe_f2))
106: return &fe->fe_f2;
107: if (ISZERO(&fe->fe_f2))
108: return &fe->fe_f2;
109:
110: if (ISINF(&fe->fe_f2))
111: return fpu_newnan(fe);
112:
113: /*
114: * +pi/2 if x == 1
115: * -pi/2 if x == -1
116: */
117: r = &fe->fe_f2;
118: if (r->fp_exp == 0 && r->fp_mant[0] == FP_1 &&
119: r->fp_mant[1] == 0 && r->fp_mant[2] == 0) {
120:
121: /* f1 <- pi/2 */
122: fpu_const(&fe->fe_f1, FPU_CONST_PI);
123: fe->fe_f1.fp_exp--;
124:
125: fe->fe_f1.fp_sign = r->fp_sign;
126: return &fe->fe_f1;
127: }
128:
129: /* NAN if |x| > 1 */
130: if (fe->fe_f2.fp_exp >= 0)
131: return fpu_newnan(fe);
132:
133: CPYFPN(&x, &fe->fe_f2);
134:
135: /* x^2 */
136: CPYFPN(&fe->fe_f1, &fe->fe_f2);
137: r = fpu_mul(fe);
138:
139: /* 1 - x^2 */
140: CPYFPN(&fe->fe_f2, r);
141: fe->fe_f2.fp_sign = 1;
142: fpu_const(&fe->fe_f1, FPU_CONST_1);
143: r = fpu_add(fe);
144:
145: /* sqrt(1-x^2) */
146: CPYFPN(&fe->fe_f2, r);
147: r = fpu_sqrt(fe);
148:
149: /* x/sqrt */
150: CPYFPN(&fe->fe_f2, r);
151: CPYFPN(&fe->fe_f1, &x);
152: r = fpu_div(fe);
153:
154: /* arctan */
155: CPYFPN(&fe->fe_f2, r);
156: return fpu_atan(fe);
157: }
158:
159: /*
160: * arctan(x):
161: *
162: * if (x < 0) {
163: * x = abs(x);
164: * sign = 1;
165: * }
166: * y = arctan(x);
167: * if (sign) {
168: * y = -y;
169: * }
170: */
171: struct fpn *
172: fpu_atan(struct fpemu *fe)
173: {
174: struct fpn a;
175: struct fpn x;
176: struct fpn v;
177:
178: if (ISNAN(&fe->fe_f2))
179: return &fe->fe_f2;
180: if (ISZERO(&fe->fe_f2))
181: return &fe->fe_f2;
182:
183: CPYFPN(&a, &fe->fe_f2);
184:
185: if (ISINF(&fe->fe_f2)) {
186: /* f2 <- pi/2 */
187: fpu_const(&fe->fe_f2, FPU_CONST_PI);
188: fe->fe_f2.fp_exp--;
189:
190: fe->fe_f2.fp_sign = a.fp_sign;
191: return &fe->fe_f2;
192: }
193:
194: fpu_const(&x, FPU_CONST_1);
195: fpu_const(&fe->fe_f2, FPU_CONST_0);
196: CPYFPN(&v, &fe->fe_f2);
197: fpu_cordit1(fe, &x, &a, &fe->fe_f2, &v);
198:
199: return &fe->fe_f2;
200: }
201:
202:
203: /*
204: * fe_f1 := sin(in)
205: * fe_f2 := cos(in)
206: */
207: static void
208: __fpu_sincos_cordic(struct fpemu *fe, const struct fpn *in)
209: {
210: struct fpn a;
211: struct fpn v;
212:
213: CPYFPN(&a, in);
214: fpu_const(&fe->fe_f1, FPU_CONST_0);
215: CPYFPN(&fe->fe_f2, &fpu_cordic_inv_gain1);
216: fpu_const(&v, FPU_CONST_1);
217: v.fp_sign = 1;
218: fpu_cordit1(fe, &fe->fe_f2, &fe->fe_f1, &a, &v);
219: }
220:
221: /*
222: * cos(x):
223: *
224: * if (x < 0) {
225: * x = abs(x);
226: * }
227: * if (x >= 2*pi) {
228: * x %= 2*pi;
229: * }
230: * if (x >= pi) {
231: * x -= pi;
232: * sign inverse;
233: * }
234: * if (x >= pi/2) {
235: * y = sin(x - pi/2);
236: * sign inverse;
237: * } else {
238: * y = cos(x);
239: * }
240: * if (sign) {
241: * y = -y;
242: * }
243: */
244: struct fpn *
245: fpu_cos(struct fpemu *fe)
246: {
247: struct fpn x;
248: struct fpn p;
249: struct fpn *r;
250: int sign;
251:
252: if (ISNAN(&fe->fe_f2))
253: return &fe->fe_f2;
254: if (ISINF(&fe->fe_f2))
255: return fpu_newnan(fe);
256:
257: /* if x is +0/-0, return 1 */
258: if (ISZERO(&fe->fe_f2)) {
259: fpu_const(&fe->fe_f2, FPU_CONST_1);
260: return &fe->fe_f2;
261: }
262:
263: /* x = abs(input) */
264: sign = 0;
265: CPYFPN(&x, &fe->fe_f2);
266: x.fp_sign = 0;
267:
268: /* p <- 2*pi */
269: fpu_const(&p, FPU_CONST_PI);
270: p.fp_exp++;
271:
272: /*
273: * if (x >= 2*pi)
274: * cos(x) is cos(x % 2*pi)
275: */
276: CPYFPN(&fe->fe_f1, &x);
277: CPYFPN(&fe->fe_f2, &p);
278: r = fpu_cmp(fe);
279: if (r->fp_sign == 0) {
280: CPYFPN(&fe->fe_f1, &x);
281: CPYFPN(&fe->fe_f2, &p);
282: r = fpu_modrem(fe, 1);
283: CPYFPN(&x, r);
284: }
285:
286: /* p <- pi */
287: p.fp_exp--;
288:
289: /*
290: * if (x >= pi)
291: * cos(x) is -cos(x - pi)
292: */
293: CPYFPN(&fe->fe_f1, &x);
294: CPYFPN(&fe->fe_f2, &p);
295: fe->fe_f2.fp_sign = 1;
296: r = fpu_add(fe);
297: if (r->fp_sign == 0) {
298: CPYFPN(&x, r);
299: sign ^= 1;
300: }
301:
302: /* p <- pi/2 */
303: p.fp_exp--;
304:
305: /*
306: * if (x >= pi/2)
307: * cos(x) is -sin(x - pi/2)
308: * else
309: * cos(x)
310: */
311: CPYFPN(&fe->fe_f1, &x);
312: CPYFPN(&fe->fe_f2, &p);
313: fe->fe_f2.fp_sign = 1;
314: r = fpu_add(fe);
315: if (r->fp_sign == 0) {
316: __fpu_sincos_cordic(fe, r);
317: r = &fe->fe_f1;
318: sign ^= 1;
319: } else {
320: __fpu_sincos_cordic(fe, &x);
321: r = &fe->fe_f2;
322: }
323: r->fp_sign = sign;
324: return r;
325: }
326:
327: /*
328: * sin(x):
329: *
330: * if (x < 0) {
331: * x = abs(x);
332: * sign = 1;
333: * }
334: * if (x >= 2*pi) {
335: * x %= 2*pi;
336: * }
337: * if (x >= pi) {
338: * x -= pi;
339: * sign inverse;
340: * }
341: * if (x >= pi/2) {
342: * y = cos(x - pi/2);
343: * } else {
344: * y = sin(x);
345: * }
346: * if (sign) {
347: * y = -y;
348: * }
349: */
350: struct fpn *
351: fpu_sin(struct fpemu *fe)
352: {
353: struct fpn x;
354: struct fpn p;
355: struct fpn *r;
356: int sign;
357:
358: if (ISNAN(&fe->fe_f2))
359: return &fe->fe_f2;
360: if (ISINF(&fe->fe_f2))
361: return fpu_newnan(fe);
362:
363: /* if x is +0/-0, return +0/-0 */
364: if (ISZERO(&fe->fe_f2))
365: return &fe->fe_f2;
366:
367: /* x = abs(input) */
368: sign = fe->fe_f2.fp_sign;
369: CPYFPN(&x, &fe->fe_f2);
370: x.fp_sign = 0;
371:
372: /* p <- 2*pi */
373: fpu_const(&p, FPU_CONST_PI);
374: p.fp_exp++;
375:
376: /*
377: * if (x >= 2*pi)
378: * sin(x) is sin(x % 2*pi)
379: */
380: CPYFPN(&fe->fe_f1, &x);
381: CPYFPN(&fe->fe_f2, &p);
382: r = fpu_cmp(fe);
383: if (r->fp_sign == 0) {
384: CPYFPN(&fe->fe_f1, &x);
385: CPYFPN(&fe->fe_f2, &p);
386: r = fpu_modrem(fe, 1);
387: CPYFPN(&x, r);
388: }
389:
390: /* p <- pi */
391: p.fp_exp--;
392:
393: /*
394: * if (x >= pi)
395: * sin(x) is -sin(x - pi)
396: */
397: CPYFPN(&fe->fe_f1, &x);
398: CPYFPN(&fe->fe_f2, &p);
399: fe->fe_f2.fp_sign = 1;
400: r = fpu_add(fe);
401: if (r->fp_sign == 0) {
402: CPYFPN(&x, r);
403: sign ^= 1;
404: }
405:
406: /* if x is small enough, sin(x) is x */
407: /* XXX the border value is not well considered... */
408: if (r->fp_exp < -25) {
409: r->fp_sign = sign;
410: r->fp_sticky = 1;
411: return r;
412: }
413:
414: /* p <- pi/2 */
415: p.fp_exp--;
416:
417: /*
418: * if (x >= pi/2)
419: * sin(x) is cos(x - pi/2)
420: * else
421: * sin(x)
422: */
423: CPYFPN(&fe->fe_f1, &x);
424: CPYFPN(&fe->fe_f2, &p);
425: fe->fe_f2.fp_sign = 1;
426: r = fpu_add(fe);
427: if (r->fp_sign == 0) {
428: __fpu_sincos_cordic(fe, r);
429: r = &fe->fe_f2;
430: } else {
431: __fpu_sincos_cordic(fe, &x);
432: r = &fe->fe_f1;
433: }
434: r->fp_sign = sign;
435: return r;
436: }
437:
438: /*
439: * tan(x):
440: *
441: * if (x < 0) {
442: * x = abs(x);
443: * sign = 1;
444: * }
445: * if (x >= pi) {
446: * x %= pi;
447: * }
448: * if (x >= pi/2) {
449: * y = tan(pi - x);
450: * sign inverse;
451: * } else {
452: * y = tan(x);
453: * }
454: */
455: struct fpn *
456: fpu_tan(struct fpemu *fe)
457: {
458: struct fpn x;
459: struct fpn p;
460: struct fpn *r;
461: int sign;
462:
463: if (ISNAN(&fe->fe_f2))
464: return &fe->fe_f2;
465: if (ISINF(&fe->fe_f2))
466: return fpu_newnan(fe);
467:
468: /* if x is +0/-0, return +0/-0 */
469: if (ISZERO(&fe->fe_f2))
470: return &fe->fe_f2;
471:
472: /* x = abs(input) */
473: sign = fe->fe_f2.fp_sign;
474: CPYFPN(&x, &fe->fe_f2);
475: x.fp_sign = 0;
476:
477: /* p <- pi */
478: fpu_const(&p, FPU_CONST_PI);
479:
480: /*
481: * if (x >= pi)
482: * tan(x) is tan(x % pi)
483: */
484: CPYFPN(&fe->fe_f1, &x);
485: CPYFPN(&fe->fe_f2, &p);
486: r = fpu_cmp(fe);
487: if (r->fp_sign == 0) {
488: CPYFPN(&fe->fe_f1, &x);
489: CPYFPN(&fe->fe_f2, &p);
490: r = fpu_modrem(fe, 1);
491: CPYFPN(&x, r);
492: }
493:
494: /*
495: * if (x >= pi/2)
496: * tan(x) is -tan(pi - x)
497: * else
498: * tan(x)
499: */
500: CPYFPN(&fe->fe_f1, &x);
501: CPYFPN(&fe->fe_f2, &p);
502: fe->fe_f2.fp_exp--;
503: r = fpu_cmp(fe);
504: if (r->fp_sign == 0) {
505: CPYFPN(&fe->fe_f1, &p);
506: CPYFPN(&fe->fe_f2, &x);
507: fe->fe_f2.fp_sign = 1;
508: r = fpu_add(fe);
509: sign ^= 1;
510: } else {
511: r = &x;
512: }
513:
514: /* tan(x) = sin(x)/cos(x) */
515: __fpu_sincos_cordic(fe, r);
516: r = fpu_div(fe);
517:
518: if (sign) {
519: r->fp_sign = !r->fp_sign;
520: }
521: return r;
522: }
523:
524: struct fpn *
525: fpu_sincos(struct fpemu *fe, int regc)
526: {
527: struct fpn x;
528: struct fpn p;
529: struct fpn *r;
530:
531: if (ISNAN(&fe->fe_f2)) {
532: fpu_implode(fe, &fe->fe_f2, FTYPE_EXT,
533: &fe->fe_fpframe->fpf_regs[regc * 3]);
534: return &fe->fe_f2;
535: }
536: if (ISINF(&fe->fe_f2)) {
537: r = fpu_newnan(fe);
538: fpu_implode(fe, r, FTYPE_EXT,
539: &fe->fe_fpframe->fpf_regs[regc * 3]);
540: return r;
541: }
542:
543: /* if x is +0/-0, cos(x) is 1 and sin(x) is +0/-0 */
544: if (ISZERO(&fe->fe_f2)) {
545: fpu_const(&fe->fe_f1, FPU_CONST_1);
546: fpu_implode(fe, &fe->fe_f1, FTYPE_EXT,
547: &fe->fe_fpframe->fpf_regs[regc * 3]);
548: return &fe->fe_f2;
549: }
550:
551: /* p <- 2*pi */
552: fpu_const(&p, FPU_CONST_PI);
553: p.fp_exp++;
554:
555: /*
556: * if (x >= 2*pi)
557: * x %= 2*pi
558: */
559: CPYFPN(&x, &fe->fe_f2);
560: CPYFPN(&fe->fe_f1, &x);
561: CPYFPN(&fe->fe_f2, &p);
562: r = fpu_cmp(fe);
563: if (r->fp_sign == 0) {
564: CPYFPN(&fe->fe_f1, &x);
565: CPYFPN(&fe->fe_f2, &p);
566: r = fpu_modrem(fe, 1);
567: CPYFPN(&x, r);
568: }
569:
570: __fpu_sincos_cordic(fe, &x);
571:
572: /* cos(x) */
573: fpu_round_prec(fe, &fe->fe_f2);
574: fpu_implode(fe, &fe->fe_f2, FTYPE_EXT, &fe->fe_fpframe->fpf_regs[regc * 3]);
575:
576: /* sin(x) */
577: return &fe->fe_f1;
578: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.