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 */