1 /*-
2 * Copyright (c) 2015 Dag-Erling Smørgrav
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
7 * are met:
8 * 1. Redistributions of source code must retain the above copyright
9 * notice, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 * notice, this list of conditions and the following disclaimer in the
12 * documentation and/or other materials provided with the distribution.
13 *
14 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
15 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
16 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
17 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
18 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
19 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
20 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
21 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
22 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
23 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
24 * SUCH DAMAGE.
25 */
26
27 #ifdef _KERNEL
28 #include <sys/libkern.h>
29 #else
30 #include <stdio.h>
31 #include <strings.h>
32 #endif
33
34 #include "fp16.h"
35
36 /*
37 * Compute the quare root of x, using Newton's method with 2^(log2(x)/2)
38 * as the initial estimate.
39 */
40 fp16_t
fp16_sqrt(fp16_t x)41 fp16_sqrt(fp16_t x)
42 {
43 fp16_t y, delta;
44 signed int log2x;
45
46 /* special case */
47 if (x == 0)
48 return (0);
49
50 /* shift toward 0 by half the logarithm */
51 log2x = flsl(x) - 1;
52 if (log2x >= 16) {
53 y = x >> (log2x - 16) / 2;
54 } else {
55 #if 0
56 y = x << (16 - log2x) / 2;
57 #else
58 /* XXX for now, return 0 for anything < 1 */
59 return (0);
60 #endif
61 }
62 while (y > 0) {
63 /* delta = y^2 / 2y */
64 delta = fp16_div(fp16_sub(fp16_mul(y, y), x), y * 2);
65 if (delta == 0)
66 break;
67 y = fp16_sub(y, delta);
68 }
69 return (y);
70 }
71
72 static fp16_t fp16_sin_table[256] = {
73 0, 402, 804, 1206, 1608, 2010, 2412, 2814,
74 3215, 3617, 4018, 4420, 4821, 5222, 5622, 6023,
75 6423, 6823, 7223, 7623, 8022, 8421, 8819, 9218,
76 9616, 10013, 10410, 10807, 11204, 11600, 11995, 12390,
77 12785, 13179, 13573, 13966, 14359, 14751, 15142, 15533,
78 15923, 16313, 16702, 17091, 17479, 17866, 18253, 18638,
79 19024, 19408, 19792, 20175, 20557, 20938, 21319, 21699,
80 22078, 22456, 22833, 23210, 23586, 23960, 24334, 24707,
81 25079, 25450, 25820, 26189, 26557, 26925, 27291, 27656,
82 28020, 28383, 28745, 29105, 29465, 29824, 30181, 30538,
83 30893, 31247, 31600, 31952, 32302, 32651, 32999, 33346,
84 33692, 34036, 34379, 34721, 35061, 35400, 35738, 36074,
85 36409, 36743, 37075, 37406, 37736, 38064, 38390, 38716,
86 39039, 39362, 39682, 40002, 40319, 40636, 40950, 41263,
87 41575, 41885, 42194, 42501, 42806, 43110, 43412, 43712,
88 44011, 44308, 44603, 44897, 45189, 45480, 45768, 46055,
89 46340, 46624, 46906, 47186, 47464, 47740, 48015, 48288,
90 48558, 48828, 49095, 49360, 49624, 49886, 50146, 50403,
91 50660, 50914, 51166, 51416, 51665, 51911, 52155, 52398,
92 52639, 52877, 53114, 53348, 53581, 53811, 54040, 54266,
93 54491, 54713, 54933, 55152, 55368, 55582, 55794, 56004,
94 56212, 56417, 56621, 56822, 57022, 57219, 57414, 57606,
95 57797, 57986, 58172, 58356, 58538, 58718, 58895, 59070,
96 59243, 59414, 59583, 59749, 59913, 60075, 60235, 60392,
97 60547, 60700, 60850, 60998, 61144, 61288, 61429, 61568,
98 61705, 61839, 61971, 62100, 62228, 62353, 62475, 62596,
99 62714, 62829, 62942, 63053, 63162, 63268, 63371, 63473,
100 63571, 63668, 63762, 63854, 63943, 64030, 64115, 64197,
101 64276, 64353, 64428, 64501, 64571, 64638, 64703, 64766,
102 64826, 64884, 64939, 64992, 65043, 65091, 65136, 65179,
103 65220, 65258, 65294, 65327, 65358, 65386, 65412, 65436,
104 65457, 65475, 65491, 65505, 65516, 65524, 65531, 65534,
105 };
106
107 /*
108 * Compute the sine of theta.
109 */
110 fp16_t
fp16_sin(fp16_t theta)111 fp16_sin(fp16_t theta)
112 {
113 unsigned int i;
114
115 i = 1024 * (theta % FP16_2PI) / FP16_2PI;
116 switch (i / 256) {
117 case 0:
118 return (fp16_sin_table[i % 256]);
119 case 1:
120 return (fp16_sin_table[255 - i % 256]);
121 case 2:
122 return (-fp16_sin_table[i % 256]);
123 case 3:
124 return (-fp16_sin_table[255 - i % 256]);
125 default:
126 /* inconceivable! */
127 return (0);
128 }
129 }
130
131 /*
132 * Compute the cosine of theta.
133 */
134 fp16_t
fp16_cos(fp16_t theta)135 fp16_cos(fp16_t theta)
136 {
137 unsigned int i;
138
139 i = 1024 * (theta % FP16_2PI) / FP16_2PI;
140 switch (i / 256) {
141 case 0:
142 return (fp16_sin_table[255 - i % 256]);
143 case 1:
144 return (-fp16_sin_table[i % 256]);
145 case 2:
146 return (-fp16_sin_table[255 - i % 256]);
147 case 3:
148 return (fp16_sin_table[i % 256]);
149 default:
150 /* inconceivable! */
151 return (0);
152 }
153 }
154