Home | History | Annotate | Download | only in base
      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