sdk-world v2026.11.05
PS2 SDK - world tier (copyleft-allowed, non-essential packages)
Loading...
Searching...
No Matches
ftcalc.c
1/****************************************************************************
2 *
3 * ftcalc.c
4 *
5 * Arithmetic computations (body).
6 *
7 * Copyright (C) 1996-2025 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 * Support for 1-complement arithmetic has been totally dropped in this
21 * release. You can still write your own code if you need it.
22 *
23 */
24
25 /**************************************************************************
26 *
27 * Implementing basic computation routines.
28 *
29 * FT_MulDiv(), FT_MulFix(), FT_DivFix(), FT_RoundFix(), FT_CeilFix(),
30 * and FT_FloorFix() are declared in freetype.h.
31 *
32 */
33
34
35#include <freetype/ftglyph.h>
36#include <freetype/fttrigon.h>
37#include <freetype/internal/ftcalc.h>
38#include <freetype/internal/ftdebug.h>
39#include <freetype/internal/ftobjs.h>
40
41 /* cancel inlining macro from internal/ftcalc.h */
42#ifdef FT_MulFix
43# undef FT_MulFix
44#endif
45
46
47 /**************************************************************************
48 *
49 * The macro FT_COMPONENT is used in trace mode. It is an implicit
50 * parameter of the FT_TRACE() and FT_ERROR() macros, used to print/log
51 * messages during execution.
52 */
53#undef FT_COMPONENT
54#define FT_COMPONENT calc
55
56
57 /* transfer sign, leaving a positive number; */
58 /* we need an unsigned value to safely negate INT_MIN (or LONG_MIN) */
59#define FT_MOVE_SIGN( utype, x, x_unsigned, s ) \
60 FT_BEGIN_STMNT \
61 if ( x < 0 ) \
62 { \
63 x_unsigned = 0U - (utype)x; \
64 s = -s; \
65 } \
66 else \
67 x_unsigned = (utype)x; \
68 FT_END_STMNT
69
70 /* The following three functions are available regardless of whether */
71 /* FT_INT64 is defined. */
72
73 /* documentation is in freetype.h */
74
75 FT_EXPORT_DEF( FT_Fixed )
76 FT_RoundFix( FT_Fixed a )
77 {
78 return ADD_LONG( a, 0x8000L - ( a < 0 ) ) & ~0xFFFFL;
79 }
80
81
82 /* documentation is in freetype.h */
83
84 FT_EXPORT_DEF( FT_Fixed )
85 FT_CeilFix( FT_Fixed a )
86 {
87 return ADD_LONG( a, 0xFFFFL ) & ~0xFFFFL;
88 }
89
90
91 /* documentation is in freetype.h */
92
93 FT_EXPORT_DEF( FT_Fixed )
94 FT_FloorFix( FT_Fixed a )
95 {
96 return a & ~0xFFFFL;
97 }
98
99#ifndef FT_MSB
100
101 FT_BASE_DEF( FT_Int )
102 FT_MSB( FT_UInt32 z )
103 {
104 FT_Int shift = 0;
105
106
107 /* determine msb bit index in `shift' */
108 if ( z & 0xFFFF0000UL )
109 {
110 z >>= 16;
111 shift += 16;
112 }
113 if ( z & 0x0000FF00UL )
114 {
115 z >>= 8;
116 shift += 8;
117 }
118 if ( z & 0x000000F0UL )
119 {
120 z >>= 4;
121 shift += 4;
122 }
123 if ( z & 0x0000000CUL )
124 {
125 z >>= 2;
126 shift += 2;
127 }
128 if ( z & 0x00000002UL )
129 {
130 /* z >>= 1; */
131 shift += 1;
132 }
133
134 return shift;
135 }
136
137#endif /* !FT_MSB */
138
139
140 /* documentation is in ftcalc.h */
141
142 FT_BASE_DEF( FT_Fixed )
143 FT_Hypot( FT_Fixed x,
144 FT_Fixed y )
145 {
146 FT_Vector v;
147
148
149 v.x = x;
150 v.y = y;
151
152 return FT_Vector_Length( &v );
153 }
154
155
156#ifdef FT_INT64
157
158
159 /* documentation is in freetype.h */
160
161 FT_EXPORT_DEF( FT_Long )
162 FT_MulDiv( FT_Long a_,
163 FT_Long b_,
164 FT_Long c_ )
165 {
166 FT_Int s = 1;
167 FT_UInt64 a, b, c, d;
168 FT_Long d_;
169
170
171 FT_MOVE_SIGN( FT_UInt64, a_, a, s );
172 FT_MOVE_SIGN( FT_UInt64, b_, b, s );
173 FT_MOVE_SIGN( FT_UInt64, c_, c, s );
174
175 d = c > 0 ? ( a * b + ( c >> 1 ) ) / c
176 : 0x7FFFFFFFUL;
177
178 d_ = (FT_Long)d;
179
180 return s < 0 ? NEG_LONG( d_ ) : d_;
181 }
182
183
184 /* documentation is in ftcalc.h */
185
186 FT_BASE_DEF( FT_Long )
187 FT_MulDiv_No_Round( FT_Long a_,
188 FT_Long b_,
189 FT_Long c_ )
190 {
191 FT_Int s = 1;
192 FT_UInt64 a, b, c, d;
193 FT_Long d_;
194
195
196 FT_MOVE_SIGN( FT_UInt64, a_, a, s );
197 FT_MOVE_SIGN( FT_UInt64, b_, b, s );
198 FT_MOVE_SIGN( FT_UInt64, c_, c, s );
199
200 d = c > 0 ? a * b / c
201 : 0x7FFFFFFFUL;
202
203 d_ = (FT_Long)d;
204
205 return s < 0 ? NEG_LONG( d_ ) : d_;
206 }
207
208
209 /* documentation is in freetype.h */
210
211 FT_EXPORT_DEF( FT_Long )
212 FT_MulFix( FT_Long a_,
213 FT_Long b_ )
214 {
215#ifdef FT_CONFIG_OPTION_INLINE_MULFIX
216
217 return FT_MulFix_64( a_, b_ );
218
219#else
220
221 FT_Int64 ab = MUL_INT64( a_, b_ );
222
223 /* this requires arithmetic right shift of signed numbers */
224 return (FT_Long)( ( ab + 0x8000L + ( ab >> 63 ) ) >> 16 );
225
226#endif /* FT_CONFIG_OPTION_INLINE_MULFIX */
227 }
228
229
230 /* documentation is in freetype.h */
231
232 FT_EXPORT_DEF( FT_Long )
233 FT_DivFix( FT_Long a_,
234 FT_Long b_ )
235 {
236 FT_Int s = 1;
237 FT_UInt64 a, b, q;
238 FT_Long q_;
239
240
241 FT_MOVE_SIGN( FT_UInt64, a_, a, s );
242 FT_MOVE_SIGN( FT_UInt64, b_, b, s );
243
244 q = b > 0 ? ( ( a << 16 ) + ( b >> 1 ) ) / b
245 : 0x7FFFFFFFUL;
246
247 q_ = (FT_Long)q;
248
249 return s < 0 ? NEG_LONG( q_ ) : q_;
250 }
251
252
253#else /* !FT_INT64 */
254
255
256 static void
257 ft_multo64( FT_UInt32 x,
258 FT_UInt32 y,
259 FT_Int64 *z )
260 {
261 FT_UInt32 lo1, hi1, lo2, hi2, lo, hi, i1, i2;
262
263
264 lo1 = x & 0x0000FFFFU; hi1 = x >> 16;
265 lo2 = y & 0x0000FFFFU; hi2 = y >> 16;
266
267 lo = lo1 * lo2;
268 i1 = lo1 * hi2;
269 i2 = lo2 * hi1;
270 hi = hi1 * hi2;
271
272 /* Check carry overflow of i1 + i2 */
273 i1 += i2;
274 hi += (FT_UInt32)( i1 < i2 ) << 16;
275
276 hi += i1 >> 16;
277 i1 = i1 << 16;
278
279 /* Check carry overflow of i1 + lo */
280 lo += i1;
281 hi += ( lo < i1 );
282
283 z->lo = lo;
284 z->hi = hi;
285 }
286
287
288 static FT_UInt32
289 ft_div64by32( FT_UInt32 hi,
290 FT_UInt32 lo,
291 FT_UInt32 y )
292 {
293 FT_UInt32 r, q;
294 FT_Int i;
295
296
297 if ( hi >= y )
298 return (FT_UInt32)0x7FFFFFFFL;
299
300 /* We shift as many bits as we can into the high register, perform */
301 /* 32-bit division with modulo there, then work through the remaining */
302 /* bits with long division. This optimization is especially noticeable */
303 /* for smaller dividends that barely use the high register. */
304
305 i = 31 - FT_MSB( hi );
306 r = ( hi << i ) | ( lo >> ( 32 - i ) ); lo <<= i; /* left 64-bit shift */
307 q = r / y;
308 r -= q * y; /* remainder */
309
310 i = 32 - i; /* bits remaining in low register */
311 do
312 {
313 q <<= 1;
314 r = ( r << 1 ) | ( lo >> 31 ); lo <<= 1;
315
316 if ( r >= y )
317 {
318 r -= y;
319 q |= 1;
320 }
321 } while ( --i );
322
323 return q;
324 }
325
326
327 static void
328 FT_Add64( FT_Int64* x,
329 FT_Int64* y,
330 FT_Int64 *z )
331 {
332 FT_UInt32 lo, hi;
333
334
335 lo = x->lo + y->lo;
336 hi = x->hi + y->hi + ( lo < x->lo );
337
338 z->lo = lo;
339 z->hi = hi;
340 }
341
342
343 /* The FT_MulDiv function has been optimized thanks to ideas from */
344 /* Graham Asher and Alexei Podtelezhnikov. The trick is to optimize */
345 /* a rather common case when everything fits within 32-bits. */
346 /* */
347 /* We compute 'a*b+c/2', then divide it by 'c' (all positive values). */
348 /* */
349 /* The product of two positive numbers never exceeds the square of */
350 /* its mean values. Therefore, we always avoid the overflow by */
351 /* imposing */
352 /* */
353 /* (a + b) / 2 <= sqrt(X - c/2) , */
354 /* */
355 /* where X = 2^32 - 1, the maximum unsigned 32-bit value, and using */
356 /* unsigned arithmetic. Now we replace `sqrt' with a linear function */
357 /* that is smaller or equal for all values of c in the interval */
358 /* [0;X/2]; it should be equal to sqrt(X) and sqrt(3X/4) at the */
359 /* endpoints. Substituting the linear solution and explicit numbers */
360 /* we get */
361 /* */
362 /* a + b <= 131071.99 - c / 122291.84 . */
363 /* */
364 /* In practice, we should use a faster and even stronger inequality */
365 /* */
366 /* a + b <= 131071 - (c >> 16) */
367 /* */
368 /* or, alternatively, */
369 /* */
370 /* a + b <= 129894 - (c >> 17) . */
371 /* */
372 /* FT_MulFix, on the other hand, is optimized for a small value of */
373 /* the first argument, when the second argument can be much larger. */
374 /* This can be achieved by scaling the second argument and the limit */
375 /* in the above inequalities. For example, */
376 /* */
377 /* a + (b >> 8) <= (131071 >> 4) */
378 /* */
379 /* covers the practical range of use. The actual test below is a bit */
380 /* tighter to avoid the border case overflows. */
381 /* */
382 /* In the case of FT_DivFix, the exact overflow check */
383 /* */
384 /* a << 16 <= X - c/2 */
385 /* */
386 /* is scaled down by 2^16 and we use */
387 /* */
388 /* a <= 65535 - (c >> 17) . */
389
390 /* documentation is in freetype.h */
391
392 FT_EXPORT_DEF( FT_Long )
393 FT_MulDiv( FT_Long a_,
394 FT_Long b_,
395 FT_Long c_ )
396 {
397 FT_Int s = 1;
398 FT_UInt32 a, b, c;
399
400
401 /* XXX: this function does not allow 64-bit arguments */
402
403 FT_MOVE_SIGN( FT_UInt32, a_, a, s );
404 FT_MOVE_SIGN( FT_UInt32, b_, b, s );
405 FT_MOVE_SIGN( FT_UInt32, c_, c, s );
406
407 if ( c == 0 )
408 a = 0x7FFFFFFFUL;
409
410 else if ( a + b <= 129894UL - ( c >> 17 ) )
411 a = ( a * b + ( c >> 1 ) ) / c;
412
413 else
414 {
415 FT_Int64 temp, temp2;
416
417
418 ft_multo64( a, b, &temp );
419
420 temp2.hi = 0;
421 temp2.lo = c >> 1;
422
423 FT_Add64( &temp, &temp2, &temp );
424
425 /* last attempt to ditch long division */
426 a = ( temp.hi == 0 ) ? temp.lo / c
427 : ft_div64by32( temp.hi, temp.lo, c );
428 }
429
430 a_ = (FT_Long)a;
431
432 return s < 0 ? NEG_LONG( a_ ) : a_;
433 }
434
435
436 FT_BASE_DEF( FT_Long )
437 FT_MulDiv_No_Round( FT_Long a_,
438 FT_Long b_,
439 FT_Long c_ )
440 {
441 FT_Int s = 1;
442 FT_UInt32 a, b, c;
443
444
445 /* XXX: this function does not allow 64-bit arguments */
446
447 FT_MOVE_SIGN( FT_UInt32, a_, a, s );
448 FT_MOVE_SIGN( FT_UInt32, b_, b, s );
449 FT_MOVE_SIGN( FT_UInt32, c_, c, s );
450
451 if ( c == 0 )
452 a = 0x7FFFFFFFUL;
453
454 else if ( a + b <= 131071UL )
455 a = a * b / c;
456
457 else
458 {
459 FT_Int64 temp;
460
461
462 ft_multo64( a, b, &temp );
463
464 /* last attempt to ditch long division */
465 a = ( temp.hi == 0 ) ? temp.lo / c
466 : ft_div64by32( temp.hi, temp.lo, c );
467 }
468
469 a_ = (FT_Long)a;
470
471 return s < 0 ? NEG_LONG( a_ ) : a_;
472 }
473
474
475 /* documentation is in freetype.h */
476
477 FT_EXPORT_DEF( FT_Long )
478 FT_MulFix( FT_Long a_,
479 FT_Long b_ )
480 {
481#ifdef FT_MULFIX_ASSEMBLER
482
483 return FT_MULFIX_ASSEMBLER( a_, b_ );
484
485#elif 0
486
487 /*
488 * This code is nonportable. See comment below.
489 *
490 * However, on a platform where right-shift of a signed quantity fills
491 * the leftmost bits by copying the sign bit, it might be faster.
492 */
493
494 FT_Long sa, sb;
495 FT_UInt32 a, b;
496
497
498 /*
499 * This is a clever way of converting a signed number `a' into its
500 * absolute value (stored back into `a') and its sign. The sign is
501 * stored in `sa'; 0 means `a' was positive or zero, and -1 means `a'
502 * was negative. (Similarly for `b' and `sb').
503 *
504 * Unfortunately, it doesn't work (at least not portably).
505 *
506 * It makes the assumption that right-shift on a negative signed value
507 * fills the leftmost bits by copying the sign bit. This is wrong.
508 * According to K&R 2nd ed, section `A7.8 Shift Operators' on page 206,
509 * the result of right-shift of a negative signed value is
510 * implementation-defined. At least one implementation fills the
511 * leftmost bits with 0s (i.e., it is exactly the same as an unsigned
512 * right shift). This means that when `a' is negative, `sa' ends up
513 * with the value 1 rather than -1. After that, everything else goes
514 * wrong.
515 */
516 sa = ( a_ >> ( sizeof ( a_ ) * 8 - 1 ) );
517 a = ( a_ ^ sa ) - sa;
518 sb = ( b_ >> ( sizeof ( b_ ) * 8 - 1 ) );
519 b = ( b_ ^ sb ) - sb;
520
521 a = (FT_UInt32)a_;
522 b = (FT_UInt32)b_;
523
524 if ( a + ( b >> 8 ) <= 8190UL )
525 a = ( a * b + 0x8000U ) >> 16;
526 else
527 {
528 FT_UInt32 al = a & 0xFFFFUL;
529
530
531 a = ( a >> 16 ) * b + al * ( b >> 16 ) +
532 ( ( al * ( b & 0xFFFFUL ) + 0x8000UL ) >> 16 );
533 }
534
535 sa ^= sb;
536 a = ( a ^ sa ) - sa;
537
538 return (FT_Long)a;
539
540#else /* 0 */
541
542 FT_Int s = 1;
543 FT_UInt32 a, b;
544
545
546 /* XXX: this function does not allow 64-bit arguments */
547
548 FT_MOVE_SIGN( FT_UInt32, a_, a, s );
549 FT_MOVE_SIGN( FT_UInt32, b_, b, s );
550
551 if ( a + ( b >> 8 ) <= 8190UL )
552 a = ( a * b + 0x8000UL ) >> 16;
553 else
554 {
555 FT_UInt32 al = a & 0xFFFFUL;
556
557
558 a = ( a >> 16 ) * b + al * ( b >> 16 ) +
559 ( ( al * ( b & 0xFFFFUL ) + 0x8000UL ) >> 16 );
560 }
561
562 a_ = (FT_Long)a;
563
564 return s < 0 ? NEG_LONG( a_ ) : a_;
565
566#endif /* 0 */
567
568 }
569
570
571 /* documentation is in freetype.h */
572
573 FT_EXPORT_DEF( FT_Long )
574 FT_DivFix( FT_Long a_,
575 FT_Long b_ )
576 {
577 FT_Int s = 1;
578 FT_UInt32 a, b, q;
579 FT_Long q_;
580
581
582 /* XXX: this function does not allow 64-bit arguments */
583
584 FT_MOVE_SIGN( FT_UInt32, a_, a, s );
585 FT_MOVE_SIGN( FT_UInt32, b_, b, s );
586
587 if ( b == 0 )
588 {
589 /* check for division by 0 */
590 q = 0x7FFFFFFFUL;
591 }
592 else if ( a <= 65535UL - ( b >> 17 ) )
593 {
594 /* compute result directly */
595 q = ( ( a << 16 ) + ( b >> 1 ) ) / b;
596 }
597 else
598 {
599 /* we need more bits; we have to do it by hand */
600 FT_Int64 temp, temp2;
601
602
603 temp.hi = a >> 16;
604 temp.lo = a << 16;
605 temp2.hi = 0;
606 temp2.lo = b >> 1;
607
608 FT_Add64( &temp, &temp2, &temp );
609 q = ft_div64by32( temp.hi, temp.lo, b );
610 }
611
612 q_ = (FT_Long)q;
613
614 return s < 0 ? NEG_LONG( q_ ) : q_;
615 }
616
617
618#endif /* !FT_INT64 */
619
620
621 /* documentation is in ftglyph.h */
622
623 FT_EXPORT_DEF( void )
624 FT_Matrix_Multiply( const FT_Matrix* a,
625 FT_Matrix *b )
626 {
627 FT_Fixed xx, xy, yx, yy;
628
629
630 if ( !a || !b )
631 return;
632
633 xx = ADD_LONG( FT_MulFix( a->xx, b->xx ),
634 FT_MulFix( a->xy, b->yx ) );
635 xy = ADD_LONG( FT_MulFix( a->xx, b->xy ),
636 FT_MulFix( a->xy, b->yy ) );
637 yx = ADD_LONG( FT_MulFix( a->yx, b->xx ),
638 FT_MulFix( a->yy, b->yx ) );
639 yy = ADD_LONG( FT_MulFix( a->yx, b->xy ),
640 FT_MulFix( a->yy, b->yy ) );
641
642 b->xx = xx;
643 b->xy = xy;
644 b->yx = yx;
645 b->yy = yy;
646 }
647
648
649 /* documentation is in ftglyph.h */
650
651 FT_EXPORT_DEF( FT_Error )
652 FT_Matrix_Invert( FT_Matrix* matrix )
653 {
654 FT_Pos delta, xx, yy;
655
656
657 if ( !matrix )
658 return FT_THROW( Invalid_Argument );
659
660 /* compute discriminant */
661 delta = FT_MulFix( matrix->xx, matrix->yy ) -
662 FT_MulFix( matrix->xy, matrix->yx );
663
664 if ( !delta )
665 return FT_THROW( Invalid_Argument ); /* matrix can't be inverted */
666
667 matrix->xy = -FT_DivFix( matrix->xy, delta );
668 matrix->yx = -FT_DivFix( matrix->yx, delta );
669
670 xx = matrix->xx;
671 yy = matrix->yy;
672
673 matrix->xx = FT_DivFix( yy, delta );
674 matrix->yy = FT_DivFix( xx, delta );
675
676 return FT_Err_Ok;
677 }
678
679
680 /* documentation is in ftcalc.h */
681
682 FT_BASE_DEF( void )
683 FT_Matrix_Multiply_Scaled( const FT_Matrix* a,
684 FT_Matrix *b,
685 FT_Long scaling )
686 {
687 FT_Fixed xx, xy, yx, yy;
688
689 FT_Long val = 0x10000L * scaling;
690
691
692 if ( !a || !b )
693 return;
694
695 xx = ADD_LONG( FT_MulDiv( a->xx, b->xx, val ),
696 FT_MulDiv( a->xy, b->yx, val ) );
697 xy = ADD_LONG( FT_MulDiv( a->xx, b->xy, val ),
698 FT_MulDiv( a->xy, b->yy, val ) );
699 yx = ADD_LONG( FT_MulDiv( a->yx, b->xx, val ),
700 FT_MulDiv( a->yy, b->yx, val ) );
701 yy = ADD_LONG( FT_MulDiv( a->yx, b->xy, val ),
702 FT_MulDiv( a->yy, b->yy, val ) );
703
704 b->xx = xx;
705 b->xy = xy;
706 b->yx = yx;
707 b->yy = yy;
708 }
709
710
711 /* documentation is in ftcalc.h */
712
713 FT_BASE_DEF( FT_Bool )
714 FT_Matrix_Check( const FT_Matrix* matrix )
715 {
716 FT_Fixed xx, xy, yx, yy;
717 FT_Fixed val;
718 FT_Int shift;
719 FT_ULong temp1, temp2;
720
721
722 if ( !matrix )
723 return 0;
724
725 xx = matrix->xx;
726 xy = matrix->xy;
727 yx = matrix->yx;
728 yy = matrix->yy;
729 val = FT_ABS( xx ) | FT_ABS( xy ) | FT_ABS( yx ) | FT_ABS( yy );
730
731 /* we only handle non-zero 32-bit values */
732 if ( !val || val > 0x7FFFFFFFL )
733 return 0;
734
735 /* Scale matrix to avoid the temp1 overflow, which is */
736 /* more stringent than avoiding the temp2 overflow. */
737
738 shift = FT_MSB( val ) - 12;
739
740 if ( shift > 0 )
741 {
742 xx >>= shift;
743 xy >>= shift;
744 yx >>= shift;
745 yy >>= shift;
746 }
747
748 temp1 = 32U * (FT_ULong)FT_ABS( xx * yy - xy * yx );
749 temp2 = (FT_ULong)( xx * xx ) + (FT_ULong)( xy * xy ) +
750 (FT_ULong)( yx * yx ) + (FT_ULong)( yy * yy );
751
752 if ( temp1 <= temp2 )
753 return 0;
754
755 return 1;
756 }
757
758
759 /* documentation is in ftcalc.h */
760
761 FT_BASE_DEF( void )
762 FT_Vector_Transform_Scaled( FT_Vector* vector,
763 const FT_Matrix* matrix,
764 FT_Long scaling )
765 {
766 FT_Pos xz, yz;
767
768 FT_Long val = 0x10000L * scaling;
769
770
771 if ( !vector || !matrix )
772 return;
773
774 xz = ADD_LONG( FT_MulDiv( vector->x, matrix->xx, val ),
775 FT_MulDiv( vector->y, matrix->xy, val ) );
776 yz = ADD_LONG( FT_MulDiv( vector->x, matrix->yx, val ),
777 FT_MulDiv( vector->y, matrix->yy, val ) );
778
779 vector->x = xz;
780 vector->y = yz;
781 }
782
783
784 /* documentation is in ftcalc.h */
785
786 FT_BASE_DEF( FT_UInt32 )
787 FT_Vector_NormLen( FT_Vector* vector )
788 {
789 FT_Int32 x_ = vector->x;
790 FT_Int32 y_ = vector->y;
791 FT_Int32 b, z;
792 FT_UInt32 x, y, u, v, l;
793 FT_Int sx = 1, sy = 1, shift;
794
795
796 FT_MOVE_SIGN( FT_UInt32, x_, x, sx );
797 FT_MOVE_SIGN( FT_UInt32, y_, y, sy );
798
799 /* trivial cases */
800 if ( x == 0 )
801 {
802 if ( y > 0 )
803 vector->y = sy * 0x10000;
804 return y;
805 }
806 else if ( y == 0 )
807 {
808 if ( x > 0 )
809 vector->x = sx * 0x10000;
810 return x;
811 }
812
813 /* Estimate length and prenormalize by shifting so that */
814 /* the new approximate length is between 2/3 and 4/3. */
815 /* The magic constant 0xAAAAAAAAUL (2/3 of 2^32) helps */
816 /* achieve this in 16.16 fixed-point representation. */
817 l = x > y ? x + ( y >> 1 )
818 : y + ( x >> 1 );
819
820 shift = 31 - FT_MSB( l );
821 shift -= 15 + ( l >= ( 0xAAAAAAAAUL >> shift ) );
822
823 if ( shift > 0 )
824 {
825 x <<= shift;
826 y <<= shift;
827
828 /* re-estimate length for tiny vectors */
829 l = x > y ? x + ( y >> 1 )
830 : y + ( x >> 1 );
831 }
832 else
833 {
834 x >>= -shift;
835 y >>= -shift;
836 l >>= -shift;
837 }
838
839 /* lower linear approximation for reciprocal length minus one */
840 b = 0x10000 - (FT_Int32)l;
841
842 x_ = (FT_Int32)x;
843 y_ = (FT_Int32)y;
844
845 /* Newton's iterations */
846 do
847 {
848 u = (FT_UInt32)( x_ + ( x_ * b >> 16 ) );
849 v = (FT_UInt32)( y_ + ( y_ * b >> 16 ) );
850
851 /* Normalized squared length in the parentheses approaches 2^32. */
852 /* On two's complement systems, converting to signed gives the */
853 /* difference with 2^32 even if the expression wraps around. */
854 z = -(FT_Int32)( u * u + v * v ) / 0x200;
855 z = z * ( ( 0x10000 + b ) >> 8 ) / 0x10000;
856
857 b += z;
858
859 } while ( z > 0 );
860
861 vector->x = sx < 0 ? -(FT_Pos)u : (FT_Pos)u;
862 vector->y = sy < 0 ? -(FT_Pos)v : (FT_Pos)v;
863
864 /* Conversion to signed helps to recover from likely wrap around */
865 /* in calculating the prenormalized length, because it gives the */
866 /* correct difference with 2^32 on two's complement systems. */
867 l = (FT_UInt32)( 0x10000 + (FT_Int32)( u * x + v * y ) / 0x10000 );
868 if ( shift > 0 )
869 l = ( l + ( 1 << ( shift - 1 ) ) ) >> shift;
870 else
871 l <<= -shift;
872
873 return l;
874 }
875
876
877 /* documentation is in ftcalc.h */
878
879 FT_BASE_DEF( FT_UInt32 )
880 FT_SqrtFixed( FT_UInt32 v )
881 {
882 if ( v == 0 )
883 return 0;
884
885#ifndef FT_INT64
886
887 /* Algorithm by Christophe Meessen (1993) with overflow fixed and */
888 /* rounding added. Any unsigned fixed 16.16 argument is acceptable. */
889 /* However, this algorithm is slower than the Babylonian method with */
890 /* a good initial guess. We only use it for large 32-bit values when */
891 /* 64-bit computations are not desirable. */
892 else if ( v > 0x10000U )
893 {
894 FT_UInt32 r = v >> 1;
895 FT_UInt32 q = ( v & 1 ) << 15;
896 FT_UInt32 b = 0x20000000;
897 FT_UInt32 t;
898
899
900 do
901 {
902 t = q + b;
903 if ( r >= t )
904 {
905 r -= t;
906 q = t + b; /* equivalent to q += 2*b */
907 }
908 r <<= 1;
909 b >>= 1;
910
911 } while ( b > 0x10 ); /* exactly 25 cycles */
912
913 return ( q + 0x40 ) >> 7;
914 }
915 else
916 {
917 FT_UInt32 r = ( v << 16 ) - 1;
918
919#else /* FT_INT64 */
920
921 else
922 {
923 FT_UInt64 r = ( (FT_UInt64)v << 16 ) - 1;
924
925#endif /* FT_INT64 */
926
927 FT_UInt32 q = 1 << ( ( 17 + FT_MSB( v ) ) >> 1 );
928 FT_UInt32 t;
929
930
931 /* Babylonian method with rounded-up division */
932 do
933 {
934 t = q;
935 q = ( t + (FT_UInt32)( r / t ) + 1 ) >> 1;
936
937 } while ( q != t ); /* less than 6 cycles */
938
939 return q;
940 }
941 }
942
943
944 /* documentation is in ftcalc.h */
945
946 FT_BASE_DEF( FT_Int )
947 ft_corner_orientation( FT_Pos in_x,
948 FT_Pos in_y,
949 FT_Pos out_x,
950 FT_Pos out_y )
951 {
952 /* we silently ignore overflow errors since such large values */
953 /* lead to even more (harmless) rendering errors later on */
954
955#ifdef FT_INT64
956
957 FT_Int64 delta = SUB_INT64( MUL_INT64( in_x, out_y ),
958 MUL_INT64( in_y, out_x ) );
959
960
961 return ( delta > 0 ) - ( delta < 0 );
962
963#else
964
965 FT_Int64 z1, z2;
966 FT_Int result;
967
968
969 if ( (FT_ULong)FT_ABS( in_x ) + (FT_ULong)FT_ABS( out_y ) <= 92681UL )
970 {
971 z1.lo = (FT_UInt32)in_x * (FT_UInt32)out_y;
972 z1.hi = (FT_UInt32)( (FT_Int32)z1.lo >> 31 ); /* sign-expansion */
973 }
974 else
975 ft_multo64( (FT_UInt32)in_x, (FT_UInt32)out_y, &z1 );
976
977 if ( (FT_ULong)FT_ABS( in_y ) + (FT_ULong)FT_ABS( out_x ) <= 92681UL )
978 {
979 z2.lo = (FT_UInt32)in_y * (FT_UInt32)out_x;
980 z2.hi = (FT_UInt32)( (FT_Int32)z2.lo >> 31 ); /* sign-expansion */
981 }
982 else
983 ft_multo64( (FT_UInt32)in_y, (FT_UInt32)out_x, &z2 );
984
985 if ( (FT_Int32)z1.hi > (FT_Int32)z2.hi )
986 result = +1;
987 else if ( (FT_Int32)z1.hi < (FT_Int32)z2.hi )
988 result = -1;
989 else if ( z1.lo > z2.lo )
990 result = +1;
991 else if ( z1.lo < z2.lo )
992 result = -1;
993 else
994 result = 0;
995
996 /* XXX: only the sign of return value, +1/0/-1 must be used */
997 return result;
998
999#endif
1000 }
1001
1002
1003 /* documentation is in ftcalc.h */
1004
1005 FT_BASE_DEF( FT_Int )
1006 ft_corner_is_flat( FT_Pos in_x,
1007 FT_Pos in_y,
1008 FT_Pos out_x,
1009 FT_Pos out_y )
1010 {
1011 FT_Pos ax = in_x + out_x;
1012 FT_Pos ay = in_y + out_y;
1013
1014 FT_Pos d_in, d_out, d_hypot;
1015
1016
1017 /* The idea of this function is to compare the length of the */
1018 /* hypotenuse with the `in' and `out' length. The `corner' */
1019 /* represented by `in' and `out' is flat if the hypotenuse's */
1020 /* length isn't too large. */
1021 /* */
1022 /* This approach has the advantage that the angle between */
1023 /* `in' and `out' is not checked. In case one of the two */
1024 /* vectors is `dominant', that is, much larger than the */
1025 /* other vector, we thus always have a flat corner. */
1026 /* */
1027 /* hypotenuse */
1028 /* x---------------------------x */
1029 /* \ / */
1030 /* \ / */
1031 /* in \ / out */
1032 /* \ / */
1033 /* o */
1034 /* Point */
1035
1036 d_in = FT_HYPOT( in_x, in_y );
1037 d_out = FT_HYPOT( out_x, out_y );
1038 d_hypot = FT_HYPOT( ax, ay );
1039
1040 /* now do a simple length comparison: */
1041 /* */
1042 /* d_in + d_out < 17/16 d_hypot */
1043
1044 return ( d_in + d_out - d_hypot ) < ( d_hypot >> 4 );
1045 }
1046
1047
1048/* END */