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 */
42fp16_t
43fp16_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
74static 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 */
112fp16_t
113fp16_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 */
136fp16_t
137fp16_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