1 /***************************************************************************/
5 /* Arithmetic computations (body). */
7 /* Copyright 1996-2006, 2008, 2012-2014 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 >= 0 ) ? ( a + 0x8000L ) & ~0xFFFFL
90 : -((-a + 0x8000L ) & ~0xFFFFL );
94 /* documentation is in freetype.h */
96 FT_EXPORT_DEF( FT_Fixed )
97 FT_CeilFix( FT_Fixed a )
99 return ( a >= 0 ) ? ( a + 0xFFFFL ) & ~0xFFFFL
100 : -((-a + 0xFFFFL ) & ~0xFFFFL );
104 /* documentation is in freetype.h */
106 FT_EXPORT_DEF( FT_Fixed )
107 FT_FloorFix( FT_Fixed a )
109 return ( a >= 0 ) ? a & ~0xFFFFL
110 : -((-a) & ~0xFFFFL );
115 FT_BASE_DEF ( FT_Int )
116 FT_MSB( FT_UInt32 z )
120 /* determine msb bit index in `shift' */
121 if ( z & 0xFFFF0000U )
126 if ( z & 0x0000FF00U )
131 if ( z & 0x000000F0U )
136 if ( z & 0x0000000CU )
141 if ( z & 0x00000002U )
153 /* documentation is in ftcalc.h */
155 FT_BASE_DEF( FT_Fixed )
156 FT_Hypot( FT_Fixed x,
165 return FT_Vector_Length( &v );
172 /* documentation is in freetype.h */
174 FT_EXPORT_DEF( FT_Long )
175 FT_MulDiv( FT_Long a,
183 FT_MOVE_SIGN( a, s );
184 FT_MOVE_SIGN( b, s );
185 FT_MOVE_SIGN( c, s );
187 d = (FT_Long)( c > 0 ? ( (FT_Int64)a * b + ( c >> 1 ) ) / c
190 return ( s > 0 ) ? d : -d;
194 /* documentation is in ftcalc.h */
196 FT_BASE_DEF( FT_Long )
197 FT_MulDiv_No_Round( FT_Long a,
205 FT_MOVE_SIGN( a, s );
206 FT_MOVE_SIGN( b, s );
207 FT_MOVE_SIGN( c, s );
209 d = (FT_Long)( c > 0 ? (FT_Int64)a * b / c
212 return ( s > 0 ) ? d : -d;
216 /* documentation is in freetype.h */
218 FT_EXPORT_DEF( FT_Long )
219 FT_MulFix( FT_Long a,
222 #ifdef FT_MULFIX_ASSEMBLER
224 return FT_MULFIX_ASSEMBLER( a, b );
232 FT_MOVE_SIGN( a, s );
233 FT_MOVE_SIGN( b, s );
235 c = (FT_Long)( ( (FT_Int64)a * b + 0x8000L ) >> 16 );
237 return ( s > 0 ) ? c : -c;
239 #endif /* FT_MULFIX_ASSEMBLER */
243 /* documentation is in freetype.h */
245 FT_EXPORT_DEF( FT_Long )
246 FT_DivFix( FT_Long a,
253 FT_MOVE_SIGN( a, s );
254 FT_MOVE_SIGN( b, s );
256 q = (FT_Long)( b > 0 ? ( ( (FT_UInt64)a << 16 ) + ( b >> 1 ) ) / b
259 return ( s < 0 ? -q : q );
263 #else /* !FT_LONG64 */
267 ft_multo64( FT_UInt32 x,
271 FT_UInt32 lo1, hi1, lo2, hi2, lo, hi, i1, i2;
274 lo1 = x & 0x0000FFFFU; hi1 = x >> 16;
275 lo2 = y & 0x0000FFFFU; hi2 = y >> 16;
282 /* Check carry overflow of i1 + i2 */
284 hi += (FT_UInt32)( i1 < i2 ) << 16;
289 /* Check carry overflow of i1 + lo */
299 ft_div64by32( FT_UInt32 hi,
308 return (FT_UInt32)0x7FFFFFFFL;
310 /* We shift as many bits as we can into the high register, perform */
311 /* 32-bit division with modulo there, then work through the remaining */
312 /* bits with long division. This optimization is especially noticeable */
313 /* for smaller dividends that barely use the high register. */
315 i = 31 - FT_MSB( hi );
316 r = ( hi << i ) | ( lo >> ( 32 - i ) ); lo <<= i; /* left 64-bit shift */
318 r -= q * y; /* remainder */
320 i = 32 - i; /* bits remaining in low register */
324 r = ( r << 1 ) | ( lo >> 31 ); lo <<= 1;
338 FT_Add64( FT_Int64* x,
346 hi = x->hi + y->hi + ( lo < x->lo );
353 /* The FT_MulDiv function has been optimized thanks to ideas from */
354 /* Graham Asher and Alexei Podtelezhnikov. The trick is to optimize */
355 /* a rather common case when everything fits within 32-bits. */
357 /* We compute 'a*b+c/2', then divide it by 'c' (all positive values). */
359 /* The product of two positive numbers never exceeds the square of */
360 /* its mean values. Therefore, we always avoid the overflow by */
363 /* (a + b) / 2 <= sqrt(X - c/2) , */
365 /* where X = 2^32 - 1, the maximum unsigned 32-bit value, and using */
366 /* unsigned arithmetic. Now we replace `sqrt' with a linear function */
367 /* that is smaller or equal for all values of c in the interval */
368 /* [0;X/2]; it should be equal to sqrt(X) and sqrt(3X/4) at the */
369 /* endpoints. Substituting the linear solution and explicit numbers */
372 /* a + b <= 131071.99 - c / 122291.84 . */
374 /* In practice, we should use a faster and even stronger inequality */
376 /* a + b <= 131071 - (c >> 16) */
378 /* or, alternatively, */
380 /* a + b <= 129894 - (c >> 17) . */
382 /* FT_MulFix, on the other hand, is optimized for a small value of */
383 /* the first argument, when the second argument can be much larger. */
384 /* This can be achieved by scaling the second argument and the limit */
385 /* in the above inequalities. For example, */
387 /* a + (b >> 8) <= (131071 >> 4) */
389 /* covers the practical range of use. The actual test below is a bit */
390 /* tighter to avoid the border case overflows. */
392 /* In the case of FT_DivFix, the exact overflow check */
394 /* a << 16 <= X - c/2 */
396 /* is scaled down by 2^16 and we use */
398 /* a <= 65535 - (c >> 17) . */
400 /* documentation is in freetype.h */
402 FT_EXPORT_DEF( FT_Long )
403 FT_MulDiv( FT_Long a,
410 /* XXX: this function does not allow 64-bit arguments */
411 if ( a == 0 || b == c )
414 FT_MOVE_SIGN( a, s );
415 FT_MOVE_SIGN( b, s );
416 FT_MOVE_SIGN( c, s );
421 else if ( (FT_ULong)a + b <= 129894UL - ( c >> 17 ) )
422 a = ( (FT_ULong)a * b + ( c >> 1 ) ) / c;
426 FT_Int64 temp, temp2;
429 ft_multo64( a, b, &temp );
434 FT_Add64( &temp, &temp2, &temp );
436 /* last attempt to ditch long division */
437 a = temp.hi == 0 ? temp.lo / c
438 : ft_div64by32( temp.hi, temp.lo, c );
441 return ( s < 0 ? -a : a );
445 FT_BASE_DEF( FT_Long )
446 FT_MulDiv_No_Round( FT_Long a,
453 if ( a == 0 || b == c )
456 FT_MOVE_SIGN( a, s );
457 FT_MOVE_SIGN( b, s );
458 FT_MOVE_SIGN( c, s );
463 else if ( (FT_ULong)a + b <= 131071UL )
464 a = (FT_ULong)a * b / c;
471 ft_multo64( a, b, &temp );
473 /* last attempt to ditch long division */
474 a = temp.hi == 0 ? temp.lo / c
475 : ft_div64by32( temp.hi, temp.lo, c );
478 return ( s < 0 ? -a : a );
482 /* documentation is in freetype.h */
484 FT_EXPORT_DEF( FT_Long )
485 FT_MulFix( FT_Long a,
488 #ifdef FT_MULFIX_ASSEMBLER
490 return FT_MULFIX_ASSEMBLER( a, b );
495 * This code is nonportable. See comment below.
497 * However, on a platform where right-shift of a signed quantity fills
498 * the leftmost bits by copying the sign bit, it might be faster.
505 if ( a == 0 || b == 0x10000L )
509 * This is a clever way of converting a signed number `a' into its
510 * absolute value (stored back into `a') and its sign. The sign is
511 * stored in `sa'; 0 means `a' was positive or zero, and -1 means `a'
512 * was negative. (Similarly for `b' and `sb').
514 * Unfortunately, it doesn't work (at least not portably).
516 * It makes the assumption that right-shift on a negative signed value
517 * fills the leftmost bits by copying the sign bit. This is wrong.
518 * According to K&R 2nd ed, section `A7.8 Shift Operators' on page 206,
519 * the result of right-shift of a negative signed value is
520 * implementation-defined. At least one implementation fills the
521 * leftmost bits with 0s (i.e., it is exactly the same as an unsigned
522 * right shift). This means that when `a' is negative, `sa' ends up
523 * with the value 1 rather than -1. After that, everything else goes
526 sa = ( a >> ( sizeof ( a ) * 8 - 1 ) );
528 sb = ( b >> ( sizeof ( b ) * 8 - 1 ) );
534 if ( ua + ( ub >> 8 ) <= 8190UL )
535 ua = ( ua * ub + 0x8000U ) >> 16;
538 FT_ULong al = ua & 0xFFFFU;
541 ua = ( ua >> 16 ) * ub + al * ( ub >> 16 ) +
542 ( ( al * ( ub & 0xFFFFU ) + 0x8000U ) >> 16 );
546 ua = (FT_ULong)(( ua ^ sa ) - sa);
556 if ( a == 0 || b == 0x10000L )
559 FT_MOVE_SIGN( a, s );
560 FT_MOVE_SIGN( b, s );
565 if ( ua + ( ub >> 8 ) <= 8190UL )
566 ua = ( ua * ub + 0x8000UL ) >> 16;
569 FT_ULong al = ua & 0xFFFFUL;
572 ua = ( ua >> 16 ) * ub + al * ( ub >> 16 ) +
573 ( ( al * ( ub & 0xFFFFUL ) + 0x8000UL ) >> 16 );
576 return ( s < 0 ? -(FT_Long)ua : (FT_Long)ua );
583 /* documentation is in freetype.h */
585 FT_EXPORT_DEF( FT_Long )
586 FT_DivFix( FT_Long a,
593 /* XXX: this function does not allow 64-bit arguments */
595 FT_MOVE_SIGN( a, s );
596 FT_MOVE_SIGN( b, s );
600 /* check for division by 0 */
603 else if ( a <= 65535L - ( b >> 17 ) )
605 /* compute result directly */
606 q = (FT_Long)( ( ( (FT_ULong)a << 16 ) + ( b >> 1 ) ) / b );
610 /* we need more bits; we have to do it by hand */
611 FT_Int64 temp, temp2;
619 FT_Add64( &temp, &temp2, &temp );
620 q = (FT_Long)ft_div64by32( temp.hi, temp.lo, b );
623 return ( s < 0 ? -q : q );
627 #endif /* FT_LONG64 */
630 /* documentation is in ftglyph.h */
632 FT_EXPORT_DEF( void )
633 FT_Matrix_Multiply( const FT_Matrix* a,
636 FT_Fixed xx, xy, yx, yy;
642 xx = FT_MulFix( a->xx, b->xx ) + FT_MulFix( a->xy, b->yx );
643 xy = FT_MulFix( a->xx, b->xy ) + FT_MulFix( a->xy, b->yy );
644 yx = FT_MulFix( a->yx, b->xx ) + FT_MulFix( a->yy, b->yx );
645 yy = FT_MulFix( a->yx, b->xy ) + FT_MulFix( a->yy, b->yy );
647 b->xx = xx; b->xy = xy;
648 b->yx = yx; b->yy = yy;
652 /* documentation is in ftglyph.h */
654 FT_EXPORT_DEF( FT_Error )
655 FT_Matrix_Invert( FT_Matrix* matrix )
657 FT_Pos delta, xx, yy;
661 return FT_THROW( Invalid_Argument );
663 /* compute discriminant */
664 delta = FT_MulFix( matrix->xx, matrix->yy ) -
665 FT_MulFix( matrix->xy, matrix->yx );
668 return FT_THROW( Invalid_Argument ); /* matrix can't be inverted */
670 matrix->xy = - FT_DivFix( matrix->xy, delta );
671 matrix->yx = - FT_DivFix( matrix->yx, delta );
676 matrix->xx = FT_DivFix( yy, delta );
677 matrix->yy = FT_DivFix( xx, delta );
683 /* documentation is in ftcalc.h */
686 FT_Matrix_Multiply_Scaled( const FT_Matrix* a,
690 FT_Fixed xx, xy, yx, yy;
692 FT_Long val = 0x10000L * scaling;
698 xx = FT_MulDiv( a->xx, b->xx, val ) + FT_MulDiv( a->xy, b->yx, val );
699 xy = FT_MulDiv( a->xx, b->xy, val ) + FT_MulDiv( a->xy, b->yy, val );
700 yx = FT_MulDiv( a->yx, b->xx, val ) + FT_MulDiv( a->yy, b->yx, val );
701 yy = FT_MulDiv( a->yx, b->xy, val ) + FT_MulDiv( a->yy, b->yy, val );
703 b->xx = xx; b->xy = xy;
704 b->yx = yx; b->yy = yy;
708 /* documentation is in ftcalc.h */
711 FT_Vector_Transform_Scaled( FT_Vector* vector,
712 const FT_Matrix* matrix,
717 FT_Long val = 0x10000L * scaling;
720 if ( !vector || !matrix )
723 xz = FT_MulDiv( vector->x, matrix->xx, val ) +
724 FT_MulDiv( vector->y, matrix->xy, val );
726 yz = FT_MulDiv( vector->x, matrix->yx, val ) +
727 FT_MulDiv( vector->y, matrix->yy, val );
736 /* documentation is in ftcalc.h */
738 FT_BASE_DEF( FT_Int32 )
739 FT_SqrtFixed( FT_Int32 x )
741 FT_UInt32 root, rem_hi, rem_lo, test_div;
754 rem_hi = ( rem_hi << 2 ) | ( rem_lo >> 30 );
757 test_div = ( root << 1 ) + 1;
759 if ( rem_hi >= test_div )
767 return (FT_Int32)root;
773 /* documentation is in ftcalc.h */
775 FT_BASE_DEF( FT_Int )
776 ft_corner_orientation( FT_Pos in_x,
781 FT_Long result; /* avoid overflow on 16-bit system */
784 /* deal with the trivial cases quickly */
792 else if ( in_x == 0 )
799 else if ( out_y == 0 )
806 else if ( out_x == 0 )
813 else /* general case */
817 FT_Int64 delta = (FT_Int64)in_x * out_y - (FT_Int64)in_y * out_x;
823 result = 1 - 2 * ( delta < 0 );
830 /* XXX: this function does not allow 64-bit arguments */
831 ft_multo64( (FT_Int32)in_x, (FT_Int32)out_y, &z1 );
832 ft_multo64( (FT_Int32)in_y, (FT_Int32)out_x, &z2 );
836 else if ( z1.hi < z2.hi )
838 else if ( z1.lo > z2.lo )
840 else if ( z1.lo < z2.lo )
848 /* XXX: only the sign of return value, +1/0/-1 must be used */
849 return (FT_Int)result;
853 /* documentation is in ftcalc.h */
855 FT_BASE_DEF( FT_Int )
856 ft_corner_is_flat( FT_Pos in_x,
861 FT_Pos ax = in_x + out_x;
862 FT_Pos ay = in_y + out_y;
864 FT_Pos d_in, d_out, d_hypot;
867 /* The idea of this function is to compare the length of the */
868 /* hypotenuse with the `in' and `out' length. The `corner' */
869 /* represented by `in' and `out' is flat if the hypotenuse's */
870 /* length isn't too large. */
872 /* This approach has the advantage that the angle between */
873 /* `in' and `out' is not checked. In case one of the two */
874 /* vectors is `dominant', this is, much larger than the */
875 /* other vector, we thus always have a flat corner. */
878 /* x---------------------------x */
886 d_in = FT_HYPOT( in_x, in_y );
887 d_out = FT_HYPOT( out_x, out_y );
888 d_hypot = FT_HYPOT( ax, ay );
890 /* now do a simple length comparison: */
892 /* d_in + d_out < 17/16 d_hypot */
894 return ( d_in + d_out - d_hypot ) < ( d_hypot >> 4 );