1 | /**************************************************************************** |
2 | * |
3 | * fttrigon.c |
4 | * |
5 | * FreeType trigonometric functions (body). |
6 | * |
7 | * Copyright (C) 2001-2023 by |
8 | * David Turner, Robert Wilhelm, and Werner Lemberg. |
9 | * |
10 | * This file is part of the FreeType project, and may only be used, |
11 | * modified, and distributed under the terms of the FreeType project |
12 | * license, LICENSE.TXT. By continuing to use, modify, or distribute |
13 | * this file you indicate that you have read the license and |
14 | * understand and accept it fully. |
15 | * |
16 | */ |
17 | |
18 | /************************************************************************** |
19 | * |
20 | * This is a fixed-point CORDIC implementation of trigonometric |
21 | * functions as well as transformations between Cartesian and polar |
22 | * coordinates. The angles are represented as 16.16 fixed-point values |
23 | * in degrees, i.e., the angular resolution is 2^-16 degrees. Note that |
24 | * only vectors longer than 2^16*180/pi (or at least 22 bits) on a |
25 | * discrete Cartesian grid can have the same or better angular |
26 | * resolution. Therefore, to maintain this precision, some functions |
27 | * require an interim upscaling of the vectors, whereas others operate |
28 | * with 24-bit long vectors directly. |
29 | * |
30 | */ |
31 | |
32 | #include <freetype/internal/ftobjs.h> |
33 | #include <freetype/internal/ftcalc.h> |
34 | #include <freetype/fttrigon.h> |
35 | |
36 | |
37 | /* the Cordic shrink factor 0.858785336480436 * 2^32 */ |
38 | #define FT_TRIG_SCALE 0xDBD95B16UL |
39 | |
40 | /* the highest bit in overflow-safe vector components, */ |
41 | /* MSB of 0.858785336480436 * sqrt(0.5) * 2^30 */ |
42 | #define FT_TRIG_SAFE_MSB 29 |
43 | |
44 | /* this table was generated for FT_PI = 180L << 16, i.e. degrees */ |
45 | #define FT_TRIG_MAX_ITERS 23 |
46 | |
47 | static const FT_Angle |
48 | ft_trig_arctan_table[] = |
49 | { |
50 | 1740967L, 919879L, 466945L, 234379L, 117304L, 58666L, 29335L, |
51 | 14668L, 7334L, 3667L, 1833L, 917L, 458L, 229L, 115L, |
52 | 57L, 29L, 14L, 7L, 4L, 2L, 1L |
53 | }; |
54 | |
55 | |
56 | #ifdef FT_INT64 |
57 | |
58 | /* multiply a given value by the CORDIC shrink factor */ |
59 | static FT_Fixed |
60 | ft_trig_downscale( FT_Fixed val ) |
61 | { |
62 | FT_Int s = 1; |
63 | |
64 | |
65 | if ( val < 0 ) |
66 | { |
67 | val = -val; |
68 | s = -1; |
69 | } |
70 | |
71 | /* 0x40000000 comes from regression analysis between true */ |
72 | /* and CORDIC hypotenuse, so it minimizes the error */ |
73 | val = (FT_Fixed)( |
74 | ( (FT_UInt64)val * FT_TRIG_SCALE + 0x40000000UL ) >> 32 ); |
75 | |
76 | return s < 0 ? -val : val; |
77 | } |
78 | |
79 | #else /* !FT_INT64 */ |
80 | |
81 | /* multiply a given value by the CORDIC shrink factor */ |
82 | static FT_Fixed |
83 | ft_trig_downscale( FT_Fixed val ) |
84 | { |
85 | FT_Int s = 1; |
86 | FT_UInt32 lo1, hi1, lo2, hi2, lo, hi, i1, i2; |
87 | |
88 | |
89 | if ( val < 0 ) |
90 | { |
91 | val = -val; |
92 | s = -1; |
93 | } |
94 | |
95 | lo1 = (FT_UInt32)val & 0x0000FFFFU; |
96 | hi1 = (FT_UInt32)val >> 16; |
97 | lo2 = FT_TRIG_SCALE & 0x0000FFFFU; |
98 | hi2 = FT_TRIG_SCALE >> 16; |
99 | |
100 | lo = lo1 * lo2; |
101 | i1 = lo1 * hi2; |
102 | i2 = lo2 * hi1; |
103 | hi = hi1 * hi2; |
104 | |
105 | /* Check carry overflow of i1 + i2 */ |
106 | i1 += i2; |
107 | hi += (FT_UInt32)( i1 < i2 ) << 16; |
108 | |
109 | hi += i1 >> 16; |
110 | i1 = i1 << 16; |
111 | |
112 | /* Check carry overflow of i1 + lo */ |
113 | lo += i1; |
114 | hi += ( lo < i1 ); |
115 | |
116 | /* 0x40000000 comes from regression analysis between true */ |
117 | /* and CORDIC hypotenuse, so it minimizes the error */ |
118 | |
119 | /* Check carry overflow of lo + 0x40000000 */ |
120 | lo += 0x40000000UL; |
121 | hi += ( lo < 0x40000000UL ); |
122 | |
123 | val = (FT_Fixed)hi; |
124 | |
125 | return s < 0 ? -val : val; |
126 | } |
127 | |
128 | #endif /* !FT_INT64 */ |
129 | |
130 | |
131 | /* undefined and never called for zero vector */ |
132 | static FT_Int |
133 | ft_trig_prenorm( FT_Vector* vec ) |
134 | { |
135 | FT_Pos x, y; |
136 | FT_Int shift; |
137 | |
138 | |
139 | x = vec->x; |
140 | y = vec->y; |
141 | |
142 | shift = FT_MSB( (FT_UInt32)( FT_ABS( x ) | FT_ABS( y ) ) ); |
143 | |
144 | if ( shift <= FT_TRIG_SAFE_MSB ) |
145 | { |
146 | shift = FT_TRIG_SAFE_MSB - shift; |
147 | vec->x = (FT_Pos)( (FT_ULong)x << shift ); |
148 | vec->y = (FT_Pos)( (FT_ULong)y << shift ); |
149 | } |
150 | else |
151 | { |
152 | shift -= FT_TRIG_SAFE_MSB; |
153 | vec->x = x >> shift; |
154 | vec->y = y >> shift; |
155 | shift = -shift; |
156 | } |
157 | |
158 | return shift; |
159 | } |
160 | |
161 | |
162 | static void |
163 | ft_trig_pseudo_rotate( FT_Vector* vec, |
164 | FT_Angle theta ) |
165 | { |
166 | FT_Int i; |
167 | FT_Fixed x, y, xtemp, b; |
168 | const FT_Angle *arctanptr; |
169 | |
170 | |
171 | x = vec->x; |
172 | y = vec->y; |
173 | |
174 | /* Rotate inside [-PI/4,PI/4] sector */ |
175 | while ( theta < -FT_ANGLE_PI4 ) |
176 | { |
177 | xtemp = y; |
178 | y = -x; |
179 | x = xtemp; |
180 | theta += FT_ANGLE_PI2; |
181 | } |
182 | |
183 | while ( theta > FT_ANGLE_PI4 ) |
184 | { |
185 | xtemp = -y; |
186 | y = x; |
187 | x = xtemp; |
188 | theta -= FT_ANGLE_PI2; |
189 | } |
190 | |
191 | arctanptr = ft_trig_arctan_table; |
192 | |
193 | /* Pseudorotations, with right shifts */ |
194 | for ( i = 1, b = 1; i < FT_TRIG_MAX_ITERS; b <<= 1, i++ ) |
195 | { |
196 | if ( theta < 0 ) |
197 | { |
198 | xtemp = x + ( ( y + b ) >> i ); |
199 | y = y - ( ( x + b ) >> i ); |
200 | x = xtemp; |
201 | theta += *arctanptr++; |
202 | } |
203 | else |
204 | { |
205 | xtemp = x - ( ( y + b ) >> i ); |
206 | y = y + ( ( x + b ) >> i ); |
207 | x = xtemp; |
208 | theta -= *arctanptr++; |
209 | } |
210 | } |
211 | |
212 | vec->x = x; |
213 | vec->y = y; |
214 | } |
215 | |
216 | |
217 | static void |
218 | ft_trig_pseudo_polarize( FT_Vector* vec ) |
219 | { |
220 | FT_Angle theta; |
221 | FT_Int i; |
222 | FT_Fixed x, y, xtemp, b; |
223 | const FT_Angle *arctanptr; |
224 | |
225 | |
226 | x = vec->x; |
227 | y = vec->y; |
228 | |
229 | /* Get the vector into [-PI/4,PI/4] sector */ |
230 | if ( y > x ) |
231 | { |
232 | if ( y > -x ) |
233 | { |
234 | theta = FT_ANGLE_PI2; |
235 | xtemp = y; |
236 | y = -x; |
237 | x = xtemp; |
238 | } |
239 | else |
240 | { |
241 | theta = y > 0 ? FT_ANGLE_PI : -FT_ANGLE_PI; |
242 | x = -x; |
243 | y = -y; |
244 | } |
245 | } |
246 | else |
247 | { |
248 | if ( y < -x ) |
249 | { |
250 | theta = -FT_ANGLE_PI2; |
251 | xtemp = -y; |
252 | y = x; |
253 | x = xtemp; |
254 | } |
255 | else |
256 | { |
257 | theta = 0; |
258 | } |
259 | } |
260 | |
261 | arctanptr = ft_trig_arctan_table; |
262 | |
263 | /* Pseudorotations, with right shifts */ |
264 | for ( i = 1, b = 1; i < FT_TRIG_MAX_ITERS; b <<= 1, i++ ) |
265 | { |
266 | if ( y > 0 ) |
267 | { |
268 | xtemp = x + ( ( y + b ) >> i ); |
269 | y = y - ( ( x + b ) >> i ); |
270 | x = xtemp; |
271 | theta += *arctanptr++; |
272 | } |
273 | else |
274 | { |
275 | xtemp = x - ( ( y + b ) >> i ); |
276 | y = y + ( ( x + b ) >> i ); |
277 | x = xtemp; |
278 | theta -= *arctanptr++; |
279 | } |
280 | } |
281 | |
282 | /* round theta to acknowledge its error that mostly comes */ |
283 | /* from accumulated rounding errors in the arctan table */ |
284 | if ( theta >= 0 ) |
285 | theta = FT_PAD_ROUND( theta, 16 ); |
286 | else |
287 | theta = -FT_PAD_ROUND( -theta, 16 ); |
288 | |
289 | vec->x = x; |
290 | vec->y = theta; |
291 | } |
292 | |
293 | |
294 | /* documentation is in fttrigon.h */ |
295 | |
296 | FT_EXPORT_DEF( FT_Fixed ) |
297 | FT_Cos( FT_Angle angle ) |
298 | { |
299 | FT_Vector v; |
300 | |
301 | |
302 | FT_Vector_Unit( &v, angle ); |
303 | |
304 | return v.x; |
305 | } |
306 | |
307 | |
308 | /* documentation is in fttrigon.h */ |
309 | |
310 | FT_EXPORT_DEF( FT_Fixed ) |
311 | FT_Sin( FT_Angle angle ) |
312 | { |
313 | FT_Vector v; |
314 | |
315 | |
316 | FT_Vector_Unit( &v, angle ); |
317 | |
318 | return v.y; |
319 | } |
320 | |
321 | |
322 | /* documentation is in fttrigon.h */ |
323 | |
324 | FT_EXPORT_DEF( FT_Fixed ) |
325 | FT_Tan( FT_Angle angle ) |
326 | { |
327 | FT_Vector v = { 1 << 24, 0 }; |
328 | |
329 | |
330 | ft_trig_pseudo_rotate( &v, angle ); |
331 | |
332 | return FT_DivFix( v.y, v.x ); |
333 | } |
334 | |
335 | |
336 | /* documentation is in fttrigon.h */ |
337 | |
338 | FT_EXPORT_DEF( FT_Angle ) |
339 | FT_Atan2( FT_Fixed dx, |
340 | FT_Fixed dy ) |
341 | { |
342 | FT_Vector v; |
343 | |
344 | |
345 | if ( dx == 0 && dy == 0 ) |
346 | return 0; |
347 | |
348 | v.x = dx; |
349 | v.y = dy; |
350 | ft_trig_prenorm( &v ); |
351 | ft_trig_pseudo_polarize( &v ); |
352 | |
353 | return v.y; |
354 | } |
355 | |
356 | |
357 | /* documentation is in fttrigon.h */ |
358 | |
359 | FT_EXPORT_DEF( void ) |
360 | FT_Vector_Unit( FT_Vector* vec, |
361 | FT_Angle angle ) |
362 | { |
363 | if ( !vec ) |
364 | return; |
365 | |
366 | vec->x = FT_TRIG_SCALE >> 8; |
367 | vec->y = 0; |
368 | ft_trig_pseudo_rotate( vec, angle ); |
369 | vec->x = ( vec->x + 0x80L ) >> 8; |
370 | vec->y = ( vec->y + 0x80L ) >> 8; |
371 | } |
372 | |
373 | |
374 | /* documentation is in fttrigon.h */ |
375 | |
376 | FT_EXPORT_DEF( void ) |
377 | FT_Vector_Rotate( FT_Vector* vec, |
378 | FT_Angle angle ) |
379 | { |
380 | FT_Int shift; |
381 | FT_Vector v; |
382 | |
383 | |
384 | if ( !vec || !angle ) |
385 | return; |
386 | |
387 | v = *vec; |
388 | |
389 | if ( v.x == 0 && v.y == 0 ) |
390 | return; |
391 | |
392 | shift = ft_trig_prenorm( &v ); |
393 | ft_trig_pseudo_rotate( &v, angle ); |
394 | v.x = ft_trig_downscale( v.x ); |
395 | v.y = ft_trig_downscale( v.y ); |
396 | |
397 | if ( shift > 0 ) |
398 | { |
399 | FT_Int32 half = (FT_Int32)1L << ( shift - 1 ); |
400 | |
401 | |
402 | vec->x = ( v.x + half - ( v.x < 0 ) ) >> shift; |
403 | vec->y = ( v.y + half - ( v.y < 0 ) ) >> shift; |
404 | } |
405 | else |
406 | { |
407 | shift = -shift; |
408 | vec->x = (FT_Pos)( (FT_ULong)v.x << shift ); |
409 | vec->y = (FT_Pos)( (FT_ULong)v.y << shift ); |
410 | } |
411 | } |
412 | |
413 | |
414 | /* documentation is in fttrigon.h */ |
415 | |
416 | FT_EXPORT_DEF( FT_Fixed ) |
417 | FT_Vector_Length( FT_Vector* vec ) |
418 | { |
419 | FT_Int shift; |
420 | FT_Vector v; |
421 | |
422 | |
423 | if ( !vec ) |
424 | return 0; |
425 | |
426 | v = *vec; |
427 | |
428 | /* handle trivial cases */ |
429 | if ( v.x == 0 ) |
430 | { |
431 | return FT_ABS( v.y ); |
432 | } |
433 | else if ( v.y == 0 ) |
434 | { |
435 | return FT_ABS( v.x ); |
436 | } |
437 | |
438 | /* general case */ |
439 | shift = ft_trig_prenorm( &v ); |
440 | ft_trig_pseudo_polarize( &v ); |
441 | |
442 | v.x = ft_trig_downscale( v.x ); |
443 | |
444 | if ( shift > 0 ) |
445 | return ( v.x + ( 1L << ( shift - 1 ) ) ) >> shift; |
446 | |
447 | return (FT_Fixed)( (FT_UInt32)v.x << -shift ); |
448 | } |
449 | |
450 | |
451 | /* documentation is in fttrigon.h */ |
452 | |
453 | FT_EXPORT_DEF( void ) |
454 | FT_Vector_Polarize( FT_Vector* vec, |
455 | FT_Fixed *length, |
456 | FT_Angle *angle ) |
457 | { |
458 | FT_Int shift; |
459 | FT_Vector v; |
460 | |
461 | |
462 | if ( !vec || !length || !angle ) |
463 | return; |
464 | |
465 | v = *vec; |
466 | |
467 | if ( v.x == 0 && v.y == 0 ) |
468 | return; |
469 | |
470 | shift = ft_trig_prenorm( &v ); |
471 | ft_trig_pseudo_polarize( &v ); |
472 | |
473 | v.x = ft_trig_downscale( v.x ); |
474 | |
475 | *length = shift >= 0 ? ( v.x >> shift ) |
476 | : (FT_Fixed)( (FT_UInt32)v.x << -shift ); |
477 | *angle = v.y; |
478 | } |
479 | |
480 | |
481 | /* documentation is in fttrigon.h */ |
482 | |
483 | FT_EXPORT_DEF( void ) |
484 | FT_Vector_From_Polar( FT_Vector* vec, |
485 | FT_Fixed length, |
486 | FT_Angle angle ) |
487 | { |
488 | if ( !vec ) |
489 | return; |
490 | |
491 | vec->x = length; |
492 | vec->y = 0; |
493 | |
494 | FT_Vector_Rotate( vec, angle ); |
495 | } |
496 | |
497 | |
498 | /* documentation is in fttrigon.h */ |
499 | |
500 | FT_EXPORT_DEF( FT_Angle ) |
501 | FT_Angle_Diff( FT_Angle angle1, |
502 | FT_Angle angle2 ) |
503 | { |
504 | FT_Angle delta = angle2 - angle1; |
505 | |
506 | |
507 | while ( delta <= -FT_ANGLE_PI ) |
508 | delta += FT_ANGLE_2PI; |
509 | |
510 | while ( delta > FT_ANGLE_PI ) |
511 | delta -= FT_ANGLE_2PI; |
512 | |
513 | return delta; |
514 | } |
515 | |
516 | |
517 | /* END */ |
518 | |