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 * $FreeBSD$ 27 */ 28 29 #ifdef _KERNEL 30 #include <sys/libkern.h> 31 #else 32 #include <stdio.h> 33 #include <strings.h> 34 #endif 35 36 #include "fp16.h" 37 38 /* 39 * Compute the quare root of x, using Newton's method with 2^(log2(x)/2) 40 * as the initial estimate. 41 */ 42 fp16_t 43 fp16_sqrt(fp16_t x) 44 { 45 fp16_t y, delta; 46 signed int log2x; 47 48 /* special case */ 49 if (x == 0) 50 return (0); 51 52 /* shift toward 0 by half the logarithm */ 53 log2x = flsl(x) - 1; 54 if (log2x >= 16) { 55 y = x >> (log2x - 16) / 2; 56 } else { 57 #if 0 58 y = x << (16 - log2x) / 2; 59 #else 60 /* XXX for now, return 0 for anything < 1 */ 61 return (0); 62 #endif 63 } 64 while (y > 0) { 65 /* delta = y^2 / 2y */ 66 delta = fp16_div(fp16_sub(fp16_mul(y, y), x), y * 2); 67 if (delta == 0) 68 break; 69 y = fp16_sub(y, delta); 70 } 71 return (y); 72 } 73 74 static fp16_t fp16_sin_table[256] = { 75 0, 402, 804, 1206, 1608, 2010, 2412, 2814, 76 3215, 3617, 4018, 4420, 4821, 5222, 5622, 6023, 77 6423, 6823, 7223, 7623, 8022, 8421, 8819, 9218, 78 9616, 10013, 10410, 10807, 11204, 11600, 11995, 12390, 79 12785, 13179, 13573, 13966, 14359, 14751, 15142, 15533, 80 15923, 16313, 16702, 17091, 17479, 17866, 18253, 18638, 81 19024, 19408, 19792, 20175, 20557, 20938, 21319, 21699, 82 22078, 22456, 22833, 23210, 23586, 23960, 24334, 24707, 83 25079, 25450, 25820, 26189, 26557, 26925, 27291, 27656, 84 28020, 28383, 28745, 29105, 29465, 29824, 30181, 30538, 85 30893, 31247, 31600, 31952, 32302, 32651, 32999, 33346, 86 33692, 34036, 34379, 34721, 35061, 35400, 35738, 36074, 87 36409, 36743, 37075, 37406, 37736, 38064, 38390, 38716, 88 39039, 39362, 39682, 40002, 40319, 40636, 40950, 41263, 89 41575, 41885, 42194, 42501, 42806, 43110, 43412, 43712, 90 44011, 44308, 44603, 44897, 45189, 45480, 45768, 46055, 91 46340, 46624, 46906, 47186, 47464, 47740, 48015, 48288, 92 48558, 48828, 49095, 49360, 49624, 49886, 50146, 50403, 93 50660, 50914, 51166, 51416, 51665, 51911, 52155, 52398, 94 52639, 52877, 53114, 53348, 53581, 53811, 54040, 54266, 95 54491, 54713, 54933, 55152, 55368, 55582, 55794, 56004, 96 56212, 56417, 56621, 56822, 57022, 57219, 57414, 57606, 97 57797, 57986, 58172, 58356, 58538, 58718, 58895, 59070, 98 59243, 59414, 59583, 59749, 59913, 60075, 60235, 60392, 99 60547, 60700, 60850, 60998, 61144, 61288, 61429, 61568, 100 61705, 61839, 61971, 62100, 62228, 62353, 62475, 62596, 101 62714, 62829, 62942, 63053, 63162, 63268, 63371, 63473, 102 63571, 63668, 63762, 63854, 63943, 64030, 64115, 64197, 103 64276, 64353, 64428, 64501, 64571, 64638, 64703, 64766, 104 64826, 64884, 64939, 64992, 65043, 65091, 65136, 65179, 105 65220, 65258, 65294, 65327, 65358, 65386, 65412, 65436, 106 65457, 65475, 65491, 65505, 65516, 65524, 65531, 65534, 107 }; 108 109 /* 110 * Compute the sine of theta. 111 */ 112 fp16_t 113 fp16_sin(fp16_t theta) 114 { 115 unsigned int i; 116 117 i = 1024 * (theta % FP16_2PI) / FP16_2PI; 118 switch (i / 256) { 119 case 0: 120 return (fp16_sin_table[i % 256]); 121 case 1: 122 return (fp16_sin_table[255 - i % 256]); 123 case 2: 124 return (-fp16_sin_table[i % 256]); 125 case 3: 126 return (-fp16_sin_table[255 - i % 256]); 127 default: 128 /* inconceivable! */ 129 return (0); 130 } 131 } 132 133 /* 134 * Compute the cosine of theta. 135 */ 136 fp16_t 137 fp16_cos(fp16_t theta) 138 { 139 unsigned int i; 140 141 i = 1024 * (theta % FP16_2PI) / FP16_2PI; 142 switch (i / 256) { 143 case 0: 144 return (fp16_sin_table[255 - i % 256]); 145 case 1: 146 return (-fp16_sin_table[i % 256]); 147 case 2: 148 return (-fp16_sin_table[255 - i % 256]); 149 case 3: 150 return (fp16_sin_table[i % 256]); 151 default: 152 /* inconceivable! */ 153 return (0); 154 } 155 } 156