1/***************************************************************************/
2/* */
3/* ftcalc.c */
4/* */
5/* Arithmetic computations (body). */
6/* */
7/* Copyright 1996-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 /* 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 trace_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( void )
751 FT_Vector_Transform_Scaled( FT_Vector* vector,
752 const FT_Matrix* matrix,
753 FT_Long scaling )
754 {
755 FT_Pos xz, yz;
756
757 FT_Long val = 0x10000L * scaling;
758
759
760 if ( !vector || !matrix )
761 return;
762
763 xz = ADD_LONG( FT_MulDiv( vector->x, matrix->xx, val ),
764 FT_MulDiv( vector->y, matrix->xy, val ) );
765 yz = ADD_LONG( FT_MulDiv( vector->x, matrix->yx, val ),
766 FT_MulDiv( vector->y, matrix->yy, val ) );
767
768 vector->x = xz;
769 vector->y = yz;
770 }
771
772
773 /* documentation is in ftcalc.h */
774
775 FT_BASE_DEF( FT_UInt32 )
776 FT_Vector_NormLen( FT_Vector* vector )
777 {
778 FT_Int32 x_ = vector->x;
779 FT_Int32 y_ = vector->y;
780 FT_Int32 b, z;
781 FT_UInt32 x, y, u, v, l;
782 FT_Int sx = 1, sy = 1, shift;
783
784
785 x = (FT_UInt32)x_;
786 y = (FT_UInt32)y_;
787
788 FT_MOVE_SIGN( x_, x, sx );
789 FT_MOVE_SIGN( y_, y, sy );
790
791 /* trivial cases */
792 if ( x == 0 )
793 {
794 if ( y > 0 )
795 vector->y = sy * 0x10000;
796 return y;
797 }
798 else if ( y == 0 )
799 {
800 if ( x > 0 )
801 vector->x = sx * 0x10000;
802 return x;
803 }
804
805 /* Estimate length and prenormalize by shifting so that */
806 /* the new approximate length is between 2/3 and 4/3. */
807 /* The magic constant 0xAAAAAAAAUL (2/3 of 2^32) helps */
808 /* achieve this in 16.16 fixed-point representation. */
809 l = x > y ? x + ( y >> 1 )
810 : y + ( x >> 1 );
811
812 shift = 31 - FT_MSB( l );
813 shift -= 15 + ( l >= ( 0xAAAAAAAAUL >> shift ) );
814
815 if ( shift > 0 )
816 {
817 x <<= shift;
818 y <<= shift;
819
820 /* re-estimate length for tiny vectors */
821 l = x > y ? x + ( y >> 1 )
822 : y + ( x >> 1 );
823 }
824 else
825 {
826 x >>= -shift;
827 y >>= -shift;
828 l >>= -shift;
829 }
830
831 /* lower linear approximation for reciprocal length minus one */
832 b = 0x10000 - (FT_Int32)l;
833
834 x_ = (FT_Int32)x;
835 y_ = (FT_Int32)y;
836
837 /* Newton's iterations */
838 do
839 {
840 u = (FT_UInt32)( x_ + ( x_ * b >> 16 ) );
841 v = (FT_UInt32)( y_ + ( y_ * b >> 16 ) );
842
843 /* Normalized squared length in the parentheses approaches 2^32. */
844 /* On two's complement systems, converting to signed gives the */
845 /* difference with 2^32 even if the expression wraps around. */
846 z = -(FT_Int32)( u * u + v * v ) / 0x200;
847 z = z * ( ( 0x10000 + b ) >> 8 ) / 0x10000;
848
849 b += z;
850
851 } while ( z > 0 );
852
853 vector->x = sx < 0 ? -(FT_Pos)u : (FT_Pos)u;
854 vector->y = sy < 0 ? -(FT_Pos)v : (FT_Pos)v;
855
856 /* Conversion to signed helps to recover from likely wrap around */
857 /* in calculating the prenormalized length, because it gives the */
858 /* correct difference with 2^32 on two's complement systems. */
859 l = (FT_UInt32)( 0x10000 + (FT_Int32)( u * x + v * y ) / 0x10000 );
860 if ( shift > 0 )
861 l = ( l + ( 1 << ( shift - 1 ) ) ) >> shift;
862 else
863 l <<= -shift;
864
865 return l;
866 }
867
868
869#if 0
870
871 /* documentation is in ftcalc.h */
872
873 FT_BASE_DEF( FT_Int32 )
874 FT_SqrtFixed( FT_Int32 x )
875 {
876 FT_UInt32 root, rem_hi, rem_lo, test_div;
877 FT_Int count;
878
879
880 root = 0;
881
882 if ( x > 0 )
883 {
884 rem_hi = 0;
885 rem_lo = (FT_UInt32)x;
886 count = 24;
887 do
888 {
889 rem_hi = ( rem_hi << 2 ) | ( rem_lo >> 30 );
890 rem_lo <<= 2;
891 root <<= 1;
892 test_div = ( root << 1 ) + 1;
893
894 if ( rem_hi >= test_div )
895 {
896 rem_hi -= test_div;
897 root += 1;
898 }
899 } while ( --count );
900 }
901
902 return (FT_Int32)root;
903 }
904
905#endif /* 0 */
906
907
908 /* documentation is in ftcalc.h */
909
910 FT_BASE_DEF( FT_Int )
911 ft_corner_orientation( FT_Pos in_x,
912 FT_Pos in_y,
913 FT_Pos out_x,
914 FT_Pos out_y )
915 {
916#ifdef FT_LONG64
917
918 FT_Int64 delta = (FT_Int64)in_x * out_y - (FT_Int64)in_y * out_x;
919
920
921 return ( delta > 0 ) - ( delta < 0 );
922
923#else
924
925 FT_Int result;
926
927
928 /* we silently ignore overflow errors, since such large values */
929 /* lead to even more (harmless) rendering errors later on */
930 if ( ADD_LONG( FT_ABS( in_x ), FT_ABS( out_y ) ) <= 131071L &&
931 ADD_LONG( FT_ABS( in_y ), FT_ABS( out_x ) ) <= 131071L )
932 {
933 FT_Long z1 = MUL_LONG( in_x, out_y );
934 FT_Long z2 = MUL_LONG( in_y, out_x );
935
936
937 if ( z1 > z2 )
938 result = +1;
939 else if ( z1 < z2 )
940 result = -1;
941 else
942 result = 0;
943 }
944 else /* products might overflow 32 bits */
945 {
946 FT_Int64 z1, z2;
947
948
949 /* XXX: this function does not allow 64-bit arguments */
950 ft_multo64( (FT_UInt32)in_x, (FT_UInt32)out_y, &z1 );
951 ft_multo64( (FT_UInt32)in_y, (FT_UInt32)out_x, &z2 );
952
953 if ( z1.hi > z2.hi )
954 result = +1;
955 else if ( z1.hi < z2.hi )
956 result = -1;
957 else if ( z1.lo > z2.lo )
958 result = +1;
959 else if ( z1.lo < z2.lo )
960 result = -1;
961 else
962 result = 0;
963 }
964
965 /* XXX: only the sign of return value, +1/0/-1 must be used */
966 return result;
967
968#endif
969 }
970
971
972 /* documentation is in ftcalc.h */
973
974 FT_BASE_DEF( FT_Int )
975 ft_corner_is_flat( FT_Pos in_x,
976 FT_Pos in_y,
977 FT_Pos out_x,
978 FT_Pos out_y )
979 {
980 FT_Pos ax = in_x + out_x;
981 FT_Pos ay = in_y + out_y;
982
983 FT_Pos d_in, d_out, d_hypot;
984
985
986 /* The idea of this function is to compare the length of the */
987 /* hypotenuse with the `in' and `out' length. The `corner' */
988 /* represented by `in' and `out' is flat if the hypotenuse's */
989 /* length isn't too large. */
990 /* */
991 /* This approach has the advantage that the angle between */
992 /* `in' and `out' is not checked. In case one of the two */
993 /* vectors is `dominant', this is, much larger than the */
994 /* other vector, we thus always have a flat corner. */
995 /* */
996 /* hypotenuse */
997 /* x---------------------------x */
998 /* \ / */
999 /* \ / */
1000 /* in \ / out */
1001 /* \ / */
1002 /* o */
1003 /* Point */
1004
1005 d_in = FT_HYPOT( in_x, in_y );
1006 d_out = FT_HYPOT( out_x, out_y );
1007 d_hypot = FT_HYPOT( ax, ay );
1008
1009 /* now do a simple length comparison: */
1010 /* */
1011 /* d_in + d_out < 17/16 d_hypot */
1012
1013 return ( d_in + d_out - d_hypot ) < ( d_hypot >> 4 );
1014 }
1015
1016
1017/* END */
1018