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