1 /****************************************************************************
2 *
3 * fttrigon.c
4 *
5 * FreeType trigonometric functions (body).
6 *
7 * Copyright (C) 2001-2019 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 = { 1 << 24, 0 };
329
330
331 ft_trig_pseudo_rotate( &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 /* documentation is in fttrigon.h */
376
377 FT_EXPORT_DEF( void )
378 FT_Vector_Rotate( FT_Vector* vec,
379 FT_Angle angle )
380 {
381 FT_Int shift;
382 FT_Vector v;
383
384
385 if ( !vec || !angle )
386 return;
387
388 v = *vec;
389
390 if ( v.x == 0 && v.y == 0 )
391 return;
392
393 shift = ft_trig_prenorm( &v );
394 ft_trig_pseudo_rotate( &v, angle );
395 v.x = ft_trig_downscale( v.x );
396 v.y = ft_trig_downscale( v.y );
397
398 if ( shift > 0 )
399 {
400 FT_Int32 half = (FT_Int32)1L << ( shift - 1 );
401
402
403 vec->x = ( v.x + half - ( v.x < 0 ) ) >> shift;
404 vec->y = ( v.y + half - ( v.y < 0 ) ) >> shift;
405 }
406 else
407 {
408 shift = -shift;
409 vec->x = (FT_Pos)( (FT_ULong)v.x << shift );
410 vec->y = (FT_Pos)( (FT_ULong)v.y << shift );
411 }
412 }
413
414
415 /* documentation is in fttrigon.h */
416
417 FT_EXPORT_DEF( FT_Fixed )
418 FT_Vector_Length( FT_Vector* vec )
419 {
420 FT_Int shift;
421 FT_Vector v;
422
423
424 if ( !vec )
425 return 0;
426
427 v = *vec;
428
429 /* handle trivial cases */
430 if ( v.x == 0 )
431 {
432 return FT_ABS( v.y );
433 }
434 else if ( v.y == 0 )
435 {
436 return FT_ABS( v.x );
437 }
438
439 /* general case */
440 shift = ft_trig_prenorm( &v );
441 ft_trig_pseudo_polarize( &v );
442
443 v.x = ft_trig_downscale( v.x );
444
445 if ( shift > 0 )
446 return ( v.x + ( 1L << ( shift - 1 ) ) ) >> shift;
447
448 return (FT_Fixed)( (FT_UInt32)v.x << -shift );
449 }
450
451
452 /* documentation is in fttrigon.h */
453
454 FT_EXPORT_DEF( void )
455 FT_Vector_Polarize( FT_Vector* vec,
456 FT_Fixed *length,
457 FT_Angle *angle )
458 {
459 FT_Int shift;
460 FT_Vector v;
461
462
463 if ( !vec || !length || !angle )
464 return;
465
466 v = *vec;
467
468 if ( v.x == 0 && v.y == 0 )
469 return;
470
471 shift = ft_trig_prenorm( &v );
472 ft_trig_pseudo_polarize( &v );
473
474 v.x = ft_trig_downscale( v.x );
475
476 *length = shift >= 0 ? ( v.x >> shift )
477 : (FT_Fixed)( (FT_UInt32)v.x << -shift );
478 *angle = v.y;
479 }
480
481
482 /* documentation is in fttrigon.h */
483
484 FT_EXPORT_DEF( void )
485 FT_Vector_From_Polar( FT_Vector* vec,
486 FT_Fixed length,
487 FT_Angle angle )
488 {
489 if ( !vec )
490 return;
491
492 vec->x = length;
493 vec->y = 0;
494
495 FT_Vector_Rotate( vec, angle );
496 }
497
498
499 /* documentation is in fttrigon.h */
500
501 FT_EXPORT_DEF( FT_Angle )
502 FT_Angle_Diff( FT_Angle angle1,
503 FT_Angle angle2 )
504 {
505 FT_Angle delta = angle2 - angle1;
506
507
508 while ( delta <= -FT_ANGLE_PI )
509 delta += FT_ANGLE_2PI;
510
511 while ( delta > FT_ANGLE_PI )
512 delta -= FT_ANGLE_2PI;
513
514 return delta;
515 }
516
517
518 /* END */