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