1 /***************************************************************************/
5 /* Arithmetic computations (body). */
7 /* Copyright 1996-2016 by */
8 /* David Turner, Robert Wilhelm, and Werner Lemberg. */
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. */
16 /***************************************************************************/
18 /*************************************************************************/
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. */
23 /*************************************************************************/
25 /*************************************************************************/
27 /* Implementing basic computation routines. */
29 /* FT_MulDiv(), FT_MulFix(), FT_DivFix(), FT_RoundFix(), FT_CeilFix(), */
30 /* and FT_FloorFix() are declared in freetype.h. */
32 /*************************************************************************/
37 #include FT_TRIGONOMETRY_H
38 #include FT_INTERNAL_CALC_H
39 #include FT_INTERNAL_DEBUG_H
40 #include FT_INTERNAL_OBJECTS_H
43 #ifdef FT_MULFIX_ASSEMBLER
47 /* we need to emulate a 64-bit data type if a real one isn't available */
51 typedef struct FT_Int64_
58 #endif /* !FT_LONG64 */
61 /*************************************************************************/
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. */
68 #define FT_COMPONENT trace_calc
71 /* transfer sign leaving a positive number */
72 #define FT_MOVE_SIGN( x, s ) \
81 /* The following three functions are available regardless of whether */
82 /* FT_LONG64 is defined. */
84 /* documentation is in freetype.h */
86 FT_EXPORT_DEF( FT_Fixed
)
87 FT_RoundFix( FT_Fixed a
)
89 return ( a
+ 0x8000L
- ( a
< 0 ) ) & ~0xFFFFL
;
93 /* documentation is in freetype.h */
95 FT_EXPORT_DEF( FT_Fixed
)
96 FT_CeilFix( FT_Fixed a
)
98 return ( a
+ 0xFFFFL
) & ~0xFFFFL
;
102 /* documentation is in freetype.h */
104 FT_EXPORT_DEF( FT_Fixed
)
105 FT_FloorFix( FT_Fixed a
)
112 FT_BASE_DEF ( FT_Int
)
113 FT_MSB( FT_UInt32 z
)
118 /* determine msb bit index in `shift' */
119 if ( z
& 0xFFFF0000UL
)
124 if ( z
& 0x0000FF00UL
)
129 if ( z
& 0x000000F0UL
)
134 if ( z
& 0x0000000CUL
)
139 if ( z
& 0x00000002UL
)
151 /* documentation is in ftcalc.h */
153 FT_BASE_DEF( FT_Fixed
)
154 FT_Hypot( FT_Fixed x
,
163 return FT_Vector_Length( &v
);
170 /* documentation is in freetype.h */
172 FT_EXPORT_DEF( FT_Long
)
173 FT_MulDiv( FT_Long a_
,
178 FT_UInt64 a
, b
, c
, d
;
182 FT_MOVE_SIGN( a_
, s
);
183 FT_MOVE_SIGN( b_
, s
);
184 FT_MOVE_SIGN( c_
, s
);
190 d
= c
> 0 ? ( a
* b
+ ( c
>> 1 ) ) / c
195 return s
< 0 ? -d_
: d_
;
199 /* documentation is in ftcalc.h */
201 FT_BASE_DEF( FT_Long
)
202 FT_MulDiv_No_Round( FT_Long a_
,
207 FT_UInt64 a
, b
, c
, d
;
211 FT_MOVE_SIGN( a_
, s
);
212 FT_MOVE_SIGN( b_
, s
);
213 FT_MOVE_SIGN( c_
, s
);
219 d
= c
> 0 ? a
* b
/ c
224 return s
< 0 ? -d_
: d_
;
228 /* documentation is in freetype.h */
230 FT_EXPORT_DEF( FT_Long
)
231 FT_MulFix( FT_Long a_
,
234 #ifdef FT_MULFIX_ASSEMBLER
236 return FT_MULFIX_ASSEMBLER( (FT_Int32
)a_
, (FT_Int32
)b_
);
240 FT_Int64 ab
= (FT_Int64
)a_
* (FT_Int64
)b_
;
242 /* this requires arithmetic right shift of signed numbers */
243 return (FT_Long
)( ( ab
+ 0x8000L
- ( ab
< 0 ) ) >> 16 );
245 #endif /* FT_MULFIX_ASSEMBLER */
249 /* documentation is in freetype.h */
251 FT_EXPORT_DEF( FT_Long
)
252 FT_DivFix( FT_Long a_
,
260 FT_MOVE_SIGN( a_
, s
);
261 FT_MOVE_SIGN( b_
, s
);
266 q
= b
> 0 ? ( ( a
<< 16 ) + ( b
>> 1 ) ) / b
271 return s
< 0 ? -q_
: q_
;
275 #else /* !FT_LONG64 */
279 ft_multo64( FT_UInt32 x
,
283 FT_UInt32 lo1
, hi1
, lo2
, hi2
, lo
, hi
, i1
, i2
;
286 lo1
= x
& 0x0000FFFFU
; hi1
= x
>> 16;
287 lo2
= y
& 0x0000FFFFU
; hi2
= y
>> 16;
294 /* Check carry overflow of i1 + i2 */
296 hi
+= (FT_UInt32
)( i1
< i2
) << 16;
301 /* Check carry overflow of i1 + lo */
311 ft_div64by32( FT_UInt32 hi
,
320 return (FT_UInt32
)0x7FFFFFFFL
;
322 /* We shift as many bits as we can into the high register, perform */
323 /* 32-bit division with modulo there, then work through the remaining */
324 /* bits with long division. This optimization is especially noticeable */
325 /* for smaller dividends that barely use the high register. */
327 i
= 31 - FT_MSB( hi
);
328 r
= ( hi
<< i
) | ( lo
>> ( 32 - i
) ); lo
<<= i
; /* left 64-bit shift */
330 r
-= q
* y
; /* remainder */
332 i
= 32 - i
; /* bits remaining in low register */
336 r
= ( r
<< 1 ) | ( lo
>> 31 ); lo
<<= 1;
350 FT_Add64( FT_Int64
* x
,
358 hi
= x
->hi
+ y
->hi
+ ( lo
< x
->lo
);
365 /* The FT_MulDiv function has been optimized thanks to ideas from */
366 /* Graham Asher and Alexei Podtelezhnikov. The trick is to optimize */
367 /* a rather common case when everything fits within 32-bits. */
369 /* We compute 'a*b+c/2', then divide it by 'c' (all positive values). */
371 /* The product of two positive numbers never exceeds the square of */
372 /* its mean values. Therefore, we always avoid the overflow by */
375 /* (a + b) / 2 <= sqrt(X - c/2) , */
377 /* where X = 2^32 - 1, the maximum unsigned 32-bit value, and using */
378 /* unsigned arithmetic. Now we replace `sqrt' with a linear function */
379 /* that is smaller or equal for all values of c in the interval */
380 /* [0;X/2]; it should be equal to sqrt(X) and sqrt(3X/4) at the */
381 /* endpoints. Substituting the linear solution and explicit numbers */
384 /* a + b <= 131071.99 - c / 122291.84 . */
386 /* In practice, we should use a faster and even stronger inequality */
388 /* a + b <= 131071 - (c >> 16) */
390 /* or, alternatively, */
392 /* a + b <= 129894 - (c >> 17) . */
394 /* FT_MulFix, on the other hand, is optimized for a small value of */
395 /* the first argument, when the second argument can be much larger. */
396 /* This can be achieved by scaling the second argument and the limit */
397 /* in the above inequalities. For example, */
399 /* a + (b >> 8) <= (131071 >> 4) */
401 /* covers the practical range of use. The actual test below is a bit */
402 /* tighter to avoid the border case overflows. */
404 /* In the case of FT_DivFix, the exact overflow check */
406 /* a << 16 <= X - c/2 */
408 /* is scaled down by 2^16 and we use */
410 /* a <= 65535 - (c >> 17) . */
412 /* documentation is in freetype.h */
414 FT_EXPORT_DEF( FT_Long
)
415 FT_MulDiv( FT_Long a_
,
423 /* XXX: this function does not allow 64-bit arguments */
425 FT_MOVE_SIGN( a_
, s
);
426 FT_MOVE_SIGN( b_
, s
);
427 FT_MOVE_SIGN( c_
, s
);
436 else if ( a
+ b
<= 129894UL - ( c
>> 17 ) )
437 a
= ( a
* b
+ ( c
>> 1 ) ) / c
;
441 FT_Int64 temp
, temp2
;
444 ft_multo64( a
, b
, &temp
);
449 FT_Add64( &temp
, &temp2
, &temp
);
451 /* last attempt to ditch long division */
452 a
= ( temp
.hi
== 0 ) ? temp
.lo
/ c
453 : ft_div64by32( temp
.hi
, temp
.lo
, c
);
458 return s
< 0 ? -a_
: a_
;
462 FT_BASE_DEF( FT_Long
)
463 FT_MulDiv_No_Round( FT_Long a_
,
471 /* XXX: this function does not allow 64-bit arguments */
473 FT_MOVE_SIGN( a_
, s
);
474 FT_MOVE_SIGN( b_
, s
);
475 FT_MOVE_SIGN( c_
, s
);
484 else if ( a
+ b
<= 131071UL )
492 ft_multo64( a
, b
, &temp
);
494 /* last attempt to ditch long division */
495 a
= ( temp
.hi
== 0 ) ? temp
.lo
/ c
496 : ft_div64by32( temp
.hi
, temp
.lo
, c
);
501 return s
< 0 ? -a_
: a_
;
505 /* documentation is in freetype.h */
507 FT_EXPORT_DEF( FT_Long
)
508 FT_MulFix( FT_Long a_
,
511 #ifdef FT_MULFIX_ASSEMBLER
513 return FT_MULFIX_ASSEMBLER( a_
, b_
);
518 * This code is nonportable. See comment below.
520 * However, on a platform where right-shift of a signed quantity fills
521 * the leftmost bits by copying the sign bit, it might be faster.
529 * This is a clever way of converting a signed number `a' into its
530 * absolute value (stored back into `a') and its sign. The sign is
531 * stored in `sa'; 0 means `a' was positive or zero, and -1 means `a'
532 * was negative. (Similarly for `b' and `sb').
534 * Unfortunately, it doesn't work (at least not portably).
536 * It makes the assumption that right-shift on a negative signed value
537 * fills the leftmost bits by copying the sign bit. This is wrong.
538 * According to K&R 2nd ed, section `A7.8 Shift Operators' on page 206,
539 * the result of right-shift of a negative signed value is
540 * implementation-defined. At least one implementation fills the
541 * leftmost bits with 0s (i.e., it is exactly the same as an unsigned
542 * right shift). This means that when `a' is negative, `sa' ends up
543 * with the value 1 rather than -1. After that, everything else goes
546 sa
= ( a_
>> ( sizeof ( a_
) * 8 - 1 ) );
547 a
= ( a_
^ sa
) - sa
;
548 sb
= ( b_
>> ( sizeof ( b_
) * 8 - 1 ) );
549 b
= ( b_
^ sb
) - sb
;
554 if ( a
+ ( b
>> 8 ) <= 8190UL )
555 a
= ( a
* b
+ 0x8000U
) >> 16;
558 FT_UInt32 al
= a
& 0xFFFFUL
;
561 a
= ( a
>> 16 ) * b
+ al
* ( b
>> 16 ) +
562 ( ( al
* ( b
& 0xFFFFUL
) + 0x8000UL
) >> 16 );
576 /* XXX: this function does not allow 64-bit arguments */
578 FT_MOVE_SIGN( a_
, s
);
579 FT_MOVE_SIGN( b_
, s
);
584 if ( a
+ ( b
>> 8 ) <= 8190UL )
585 a
= ( a
* b
+ 0x8000UL
) >> 16;
588 FT_UInt32 al
= a
& 0xFFFFUL
;
591 a
= ( a
>> 16 ) * b
+ al
* ( b
>> 16 ) +
592 ( ( al
* ( b
& 0xFFFFUL
) + 0x8000UL
) >> 16 );
597 return s
< 0 ? -a_
: a_
;
604 /* documentation is in freetype.h */
606 FT_EXPORT_DEF( FT_Long
)
607 FT_DivFix( FT_Long a_
,
615 /* XXX: this function does not allow 64-bit arguments */
617 FT_MOVE_SIGN( a_
, s
);
618 FT_MOVE_SIGN( b_
, s
);
625 /* check for division by 0 */
628 else if ( a
<= 65535UL - ( b
>> 17 ) )
630 /* compute result directly */
631 q
= ( ( a
<< 16 ) + ( b
>> 1 ) ) / b
;
635 /* we need more bits; we have to do it by hand */
636 FT_Int64 temp
, temp2
;
644 FT_Add64( &temp
, &temp2
, &temp
);
645 q
= ft_div64by32( temp
.hi
, temp
.lo
, b
);
650 return s
< 0 ? -q_
: q_
;
654 #endif /* !FT_LONG64 */
657 /* documentation is in ftglyph.h */
659 FT_EXPORT_DEF( void )
660 FT_Matrix_Multiply( const FT_Matrix
* a
,
663 FT_Fixed xx
, xy
, yx
, yy
;
669 xx
= FT_MulFix( a
->xx
, b
->xx
) + FT_MulFix( a
->xy
, b
->yx
);
670 xy
= FT_MulFix( a
->xx
, b
->xy
) + FT_MulFix( a
->xy
, b
->yy
);
671 yx
= FT_MulFix( a
->yx
, b
->xx
) + FT_MulFix( a
->yy
, b
->yx
);
672 yy
= FT_MulFix( a
->yx
, b
->xy
) + FT_MulFix( a
->yy
, b
->yy
);
674 b
->xx
= xx
; b
->xy
= xy
;
675 b
->yx
= yx
; b
->yy
= yy
;
679 /* documentation is in ftglyph.h */
681 FT_EXPORT_DEF( FT_Error
)
682 FT_Matrix_Invert( FT_Matrix
* matrix
)
684 FT_Pos delta
, xx
, yy
;
688 return FT_THROW( Invalid_Argument
);
690 /* compute discriminant */
691 delta
= FT_MulFix( matrix
->xx
, matrix
->yy
) -
692 FT_MulFix( matrix
->xy
, matrix
->yx
);
695 return FT_THROW( Invalid_Argument
); /* matrix can't be inverted */
697 matrix
->xy
= - FT_DivFix( matrix
->xy
, delta
);
698 matrix
->yx
= - FT_DivFix( matrix
->yx
, delta
);
703 matrix
->xx
= FT_DivFix( yy
, delta
);
704 matrix
->yy
= FT_DivFix( xx
, delta
);
710 /* documentation is in ftcalc.h */
713 FT_Matrix_Multiply_Scaled( const FT_Matrix
* a
,
717 FT_Fixed xx
, xy
, yx
, yy
;
719 FT_Long val
= 0x10000L
* scaling
;
725 xx
= FT_MulDiv( a
->xx
, b
->xx
, val
) + FT_MulDiv( a
->xy
, b
->yx
, val
);
726 xy
= FT_MulDiv( a
->xx
, b
->xy
, val
) + FT_MulDiv( a
->xy
, b
->yy
, val
);
727 yx
= FT_MulDiv( a
->yx
, b
->xx
, val
) + FT_MulDiv( a
->yy
, b
->yx
, val
);
728 yy
= FT_MulDiv( a
->yx
, b
->xy
, val
) + FT_MulDiv( a
->yy
, b
->yy
, val
);
730 b
->xx
= xx
; b
->xy
= xy
;
731 b
->yx
= yx
; b
->yy
= yy
;
735 /* documentation is in ftcalc.h */
738 FT_Vector_Transform_Scaled( FT_Vector
* vector
,
739 const FT_Matrix
* matrix
,
744 FT_Long val
= 0x10000L
* scaling
;
747 if ( !vector
|| !matrix
)
750 xz
= FT_MulDiv( vector
->x
, matrix
->xx
, val
) +
751 FT_MulDiv( vector
->y
, matrix
->xy
, val
);
753 yz
= FT_MulDiv( vector
->x
, matrix
->yx
, val
) +
754 FT_MulDiv( vector
->y
, matrix
->yy
, val
);
761 /* documentation is in ftcalc.h */
763 FT_BASE_DEF( FT_UInt32
)
764 FT_Vector_NormLen( FT_Vector
* vector
)
766 FT_Int32 x_
= vector
->x
;
767 FT_Int32 y_
= vector
->y
;
769 FT_UInt32 x
, y
, u
, v
, l
;
770 FT_Int sx
= 1, sy
= 1, shift
;
773 FT_MOVE_SIGN( x_
, sx
);
774 FT_MOVE_SIGN( y_
, sy
);
783 vector
->y
= sy
* 0x10000;
789 vector
->x
= sx
* 0x10000;
793 /* Estimate length and prenormalize by shifting so that */
794 /* the new approximate length is between 2/3 and 4/3. */
795 /* The magic constant 0xAAAAAAAAUL (2/3 of 2^32) helps */
796 /* achieve this in 16.16 fixed-point representation. */
797 l
= x
> y
? x
+ ( y
>> 1 )
800 shift
= 31 - FT_MSB( l
);
801 shift
-= 15 + ( l
>= ( 0xAAAAAAAAUL
>> shift
) );
808 /* re-estimate length for tiny vectors */
809 l
= x
> y
? x
+ ( y
>> 1 )
819 /* lower linear approximation for reciprocal length minus one */
820 b
= 0x10000 - (FT_Int32
)l
;
825 /* Newton's iterations */
828 u
= (FT_UInt32
)( x_
+ ( x_
* b
>> 16 ) );
829 v
= (FT_UInt32
)( y_
+ ( y_
* b
>> 16 ) );
831 /* Normalized squared length in the parentheses approaches 2^32. */
832 /* On two's complement systems, converting to signed gives the */
833 /* difference with 2^32 even if the expression wraps around. */
834 z
= -(FT_Int32
)( u
* u
+ v
* v
) / 0x200;
835 z
= z
* ( ( 0x10000 + b
) >> 8 ) / 0x10000;
841 vector
->x
= sx
< 0 ? -(FT_Pos
)u
: (FT_Pos
)u
;
842 vector
->y
= sy
< 0 ? -(FT_Pos
)v
: (FT_Pos
)v
;
844 /* Conversion to signed helps to recover from likely wrap around */
845 /* in calculating the prenormalized length, because it gives the */
846 /* correct difference with 2^32 on two's complement systems. */
847 l
= (FT_UInt32
)( 0x10000 + (FT_Int32
)( u
* x
+ v
* y
) / 0x10000 );
849 l
= ( l
+ ( 1 << ( shift
- 1 ) ) ) >> shift
;
859 /* documentation is in ftcalc.h */
861 FT_BASE_DEF( FT_Int32
)
862 FT_SqrtFixed( FT_Int32 x
)
864 FT_UInt32 root
, rem_hi
, rem_lo
, test_div
;
873 rem_lo
= (FT_UInt32
)x
;
877 rem_hi
= ( rem_hi
<< 2 ) | ( rem_lo
>> 30 );
880 test_div
= ( root
<< 1 ) + 1;
882 if ( rem_hi
>= test_div
)
890 return (FT_Int32
)root
;
896 /* documentation is in ftcalc.h */
898 FT_BASE_DEF( FT_Int
)
899 ft_corner_orientation( FT_Pos in_x
,
906 FT_Int64 delta
= (FT_Int64
)in_x
* out_y
- (FT_Int64
)in_y
* out_x
;
909 return ( delta
> 0 ) - ( delta
< 0 );
916 if ( (FT_ULong
)FT_ABS( in_x
) + (FT_ULong
)FT_ABS( out_y
) <= 131071UL &&
917 (FT_ULong
)FT_ABS( in_y
) + (FT_ULong
)FT_ABS( out_x
) <= 131071UL )
919 FT_Long z1
= in_x
* out_y
;
920 FT_Long z2
= in_y
* out_x
;
930 else /* products might overflow 32 bits */
935 /* XXX: this function does not allow 64-bit arguments */
936 ft_multo64( (FT_UInt32
)in_x
, (FT_UInt32
)out_y
, &z1
);
937 ft_multo64( (FT_UInt32
)in_y
, (FT_UInt32
)out_x
, &z2
);
941 else if ( z1
.hi
< z2
.hi
)
943 else if ( z1
.lo
> z2
.lo
)
945 else if ( z1
.lo
< z2
.lo
)
951 /* XXX: only the sign of return value, +1/0/-1 must be used */
958 /* documentation is in ftcalc.h */
960 FT_BASE_DEF( FT_Int
)
961 ft_corner_is_flat( FT_Pos in_x
,
966 FT_Pos ax
= in_x
+ out_x
;
967 FT_Pos ay
= in_y
+ out_y
;
969 FT_Pos d_in
, d_out
, d_hypot
;
972 /* The idea of this function is to compare the length of the */
973 /* hypotenuse with the `in' and `out' length. The `corner' */
974 /* represented by `in' and `out' is flat if the hypotenuse's */
975 /* length isn't too large. */
977 /* This approach has the advantage that the angle between */
978 /* `in' and `out' is not checked. In case one of the two */
979 /* vectors is `dominant', this is, much larger than the */
980 /* other vector, we thus always have a flat corner. */
983 /* x---------------------------x */
991 d_in
= FT_HYPOT( in_x
, in_y
);
992 d_out
= FT_HYPOT( out_x
, out_y
);
993 d_hypot
= FT_HYPOT( ax
, ay
);
995 /* now do a simple length comparison: */
997 /* d_in + d_out < 17/16 d_hypot */
999 return ( d_in
+ d_out
- d_hypot
) < ( d_hypot
>> 4 );