1 /* $OpenBSD: vecast.c,v 1.1 2022/10/22 17:50:28 gkoehler Exp $ */
2
3 /*
4 * Copyright (c) 2022 George Koehler <gkoehler@openbsd.org>
5 *
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18
19 #include <altivec.h>
20 #include <err.h>
21 #include <math.h>
22 #include <stdio.h>
23
24 struct double4 {
25 double d[4];
26 };
27
28 union vu {
29 vector float vf;
30 vector int vi;
31 vector unsigned vu;
32 float f[4];
33 int i[4];
34 unsigned u[4];
35 };
36
37 #define AD(a, b, c, d) (struct double4){a, b, c, d}
38 #define VF(a, b, c, d) (vector float)(a, b, c, d)
39 #define VI(a, b, c, d) (vector int)(a, b, c, d)
40 #define VU(a, b, c, d) (vector unsigned)(a, b, c, d)
41 #define rsqrt(f) (1.0 / sqrt(f))
42
43 int fail;
44
45 void
ck_equal(const char * what,vector float out,vector float answer)46 ck_equal(const char *what, vector float out, vector float answer)
47 {
48 if (vec_any_ne(out, answer)) {
49 union vu a, b;
50
51 a.vf = out;
52 b.vf = answer;
53 warnx("%s: {%a, %a, %a, %a} should be {%a, %a, %a, %a}",
54 what, a.f[0], a.f[1], a.f[2], a.f[3],
55 b.f[0], b.f[1], b.f[2], b.f[3]);
56 fail = 1;
57 }
58 }
59
60 void
ck_equal_i(const char * what,vector int out,vector int answer)61 ck_equal_i(const char *what, vector int out, vector int answer)
62 {
63 if (vec_any_ne(out, answer)) {
64 union vu a, b;
65
66 a.vi = out;
67 b.vi = answer;
68 warnx("%s: {%d, %d, %d, %d} should be {%d, %d, %d, %d}",
69 what, a.i[0], a.i[1], a.i[2], a.i[3],
70 b.i[0], b.i[1], b.i[2], b.i[3]);
71 fail = 1;
72 }
73 }
74
75 void
ck_equal_u(const char * what,vector unsigned out,vector unsigned answer)76 ck_equal_u(const char *what, vector unsigned out, vector unsigned answer)
77 {
78 if (vec_any_ne(out, answer)) {
79 union vu a, b;
80
81 a.vi = out;
82 b.vi = answer;
83 warnx("%s: {%u, %u, %u, %u} should be {%u, %u, %u, %u}",
84 what, a.u[0], a.u[1], a.u[2], a.u[3],
85 b.u[0], b.u[1], b.u[2], b.u[3]);
86 fail = 1;
87 }
88 }
89
90 enum error_check {REL_1_IN, ABS_1_IN};
91
92 /* Checks that error is at most 1 in err_den. */
93 void
ck_estimate(const char * what,vector float out,struct double4 answer,enum error_check which,double err_den)94 ck_estimate(const char *what, vector float out, struct double4 answer,
95 enum error_check which, double err_den)
96 {
97 union vu u;
98 int i, warned = 0;
99
100 u.vf = out;
101 for (i = 0; i < 4; i++) {
102 double estimate = u.f[i];
103 double target = answer.d[i];
104 double error;
105
106 switch (which) {
107 case REL_1_IN: /* relative error */
108 error = fabs(target / (estimate - target));
109 break;
110 case ABS_1_IN: /* absolute error */
111 error = fabs(1 / (estimate - target));
112 break;
113 default:
114 errx(1, "invalid check");
115 }
116
117 if (error < err_den) {
118 if (!warned) {
119 warnx("%s: {%a, %a, %a, %a} should be "
120 "near {%a, %a, %a, %a} (1/%g)",
121 what, u.f[0], u.f[1], u.f[2], u.f[3],
122 answer.d[0], answer.d[1], answer.d[2],
123 answer.d[3], err_den);
124 warned = 1;
125 fail = 1;
126 }
127 warnx("%a is off %a by 1/%g", estimate, target,
128 error);
129 }
130 }
131 }
132
133 /*
134 * Tries altivec with denormal or subnormal floats.
135 * These are single-precision floats f, where
136 * 0 < |f| < 2**-126 = 0x1p-126 = 0x10p-130 = 1.17549435E-38F
137 */
138 int
main(void)139 main(void)
140 {
141 struct double4 dan;
142 volatile vector float in1, in2, in3;
143 vector float ans;
144 vector int ian;
145 vector unsigned uan;
146
147 /* in1 + in2 */
148 in1 = VF(10, 0x10p-140, 0x20p-130, -0x2000p-134);
149 in2 = VF( 4, 0x5p-140, -0x1p-130, 0x1fffp-134);
150 ans = VF(14, 0x15p-140, 0x1fp-130, -0x1p-134);
151 ck_equal("vec_add", vec_add(in1, in2), ans);
152
153 /* in1 - in2 */
154 in1 = VF(0x4000p-134, 10, 0x10p-140, 0x3p-130);
155 in2 = VF(0x3ffep-134, 4, 0x5p-140, 0x40p-130);
156 ans = VF( 0x2p-134, 6, 0xbp-140, -0x3dp-130);
157 ck_equal("vec_sub", vec_sub(in1, in2), ans);
158
159 /* in1 * in2 + in3 */
160 in1 = VF( 0x6p-70, 0x6p-140, 6, 0x6p-100);
161 in2 = VF( 0x7p-70, 0x7p50, 7, 0x7p-30);
162 in3 = VF( 0, 0, 1, -0x20p-130);
163 ans = VF(0x2ap-140, 0x2ap-90, 43, 0xap-130);
164 ck_equal("vec_madd", vec_madd(in1, in2, in3), ans);
165
166 /* in3 - in1 * in2 */
167 in1 = VF( 0xbp-30, 0xbp-70, 0xbp44, 11);
168 in2 = VF( 0x3p-100, 0x3p-70, -0x3p-138, 3);
169 in3 = VF(0x25p-130, 0, 0, 35);
170 ans = VF( 0x4p-130, -0x21p-140, 0x21p-94, 2);
171 ck_equal("vec_nmsub", vec_nmsub(in1, in2, in3), ans);
172
173 /* 1 / in1 */
174 in1 = VF( 3, 0x3p126, 0x3p-126, 0x1p127);
175 dan = AD(1.0 / 3, 1.0 / 0x3p126, 1.0 / 0x3p-126, 0x1p-127);
176 ck_estimate("vec_re", vec_re(in1), dan, REL_1_IN, 4096);
177
178 /* 1 / sqrt(in1) */
179 in1 = VF(1, 2, 0x1p-128, 0x5p-135);
180 dan = AD(1, rsqrt(2), rsqrt(0x1p-128), rsqrt(0x5p-135));
181 ck_estimate("vec_rsqrt", vec_rsqrte(in1), dan, REL_1_IN, 4096);
182
183 /* log2(in1) */
184 in1 = VF(0x1p-130, 0x1p-149, 32, 0x1p-10);
185 dan = AD( -130, -149, 5, -10);
186 ck_estimate("vec_loge", vec_loge(in1), dan, ABS_1_IN, 32);
187 in1 = VF( 0x123p-139, 0xabcp-145, 1, 1);
188 dan = AD(log2(0x123p-139), log2(0xabcp-145), 0, 0);
189 ck_estimate("vec_loge", vec_loge(in1), dan, ABS_1_IN, 32);
190
191 /* 2**in1 */
192 in1 = VF( -149, -138, -127, 10);
193 ans = VF(0x1p-149, 0x1p-138, 0x1p-127, 1024);
194 ck_equal("vec_expte", vec_expte(in1), ans);
195 in1 = VF( -10, -145.3, -136.9, -127.1);
196 dan = AD(0x1p-10, exp2(-145.3), exp2(-136.9), exp2(-127.1));
197 ck_estimate("vec_expte", vec_expte(in1), dan, REL_1_IN, 16);
198
199 /* (int)(in1 * 2**exponent) */
200 in1 = VF(0x1p-127, 2.34, -0xfedp-140, -19.8);
201 ian = VI( 0, 2, 0, -19);
202 ck_equal_i("vec_cts", vec_cts(in1, 0), ian);
203 in1 = VF(0x1p-113, -1, -0xabcp-143, 0x1fp-10);
204 ian = VI( 0, -1024, 0, 0x1f);
205 ck_equal_i("vec_cts", vec_cts(in1, 10), ian);
206
207 /* (unsigned)(in1 * 2**exponent) */
208 in1 = VF(0x1.ap-130, 0x1.ep-140, 24000012, 0);
209 uan = VU( 0, 0, 3072001536, 0);
210 ck_equal_u("vec_ctu", vec_ctu(in1, 7), uan);
211
212 return fail;
213 }
214