src/base/ftcalc.c

   1 /***************************************************************************/
   2 /*                                                                         */
   3 /*  ftcalc.c                                                               */
   4 /*                                                                         */
   5 /*    Arithmetic computations (body).                                      */
   6 /*                                                                         */
   7 /*  Copyright 1996-2006, 2008, 2012-2014 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 #define FT_MOVE_SIGN( x, s ) \
  73   FT_BEGIN_STMNT             \
  74     if ( x < 0 )             \
  75     {                        \
  76       x = -x;                \
  77       s = -s;                \
  78     }                        \
  79   FT_END_STMNT
  80
  81   /* The following three functions are available regardless of whether */
  82   /* FT_LONG64 is defined.                                             */
  83
  84   /* documentation is in freetype.h */
  85
  86   FT_EXPORT_DEF( FT_Fixed )
  87   FT_RoundFix( FT_Fixed  a )
  88   {
  89     return a >= 0 ?   ( a + 0x8000L ) & ~0xFFFFL
  90                   : -((-a + 0x8000L ) & ~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 a >= 0 ?   ( a + 0xFFFFL ) & ~0xFFFFL
 100                   : -((-a + 0xFFFFL ) & ~0xFFFFL );
 101   }
 102
 103
 104   /* documentation is in freetype.h */
 105
 106   FT_EXPORT_DEF( FT_Fixed )
 107   FT_FloorFix( FT_Fixed  a )
 108   {
 109     return a >= 0 ?   a & ~0xFFFFL
 110                   : -((-a) & ~0xFFFFL );
 111   }
 112
 113 #ifndef FT_MSB
 114
 115   FT_BASE_DEF ( FT_Int )
 116   FT_MSB( FT_UInt32 z )
 117   {
 118     FT_Int  shift = 0;
 119
 120
 121     /* determine msb bit index in `shift' */
 122     if ( z & 0xFFFF0000UL )
 123     {
 124       z     >>= 16;
 125       shift  += 16;
 126     }
 127     if ( z & 0x0000FF00UL )
 128     {
 129       z     >>= 8;
 130       shift  += 8;
 131     }
 132     if ( z & 0x000000F0UL )
 133     {
 134       z     >>= 4;
 135       shift  += 4;
 136     }
 137     if ( z & 0x0000000CUL )
 138     {
 139       z     >>= 2;
 140       shift  += 2;
 141     }
 142     if ( z & 0x00000002UL )
 143     {
 144    /* z     >>= 1; */
 145       shift  += 1;
 146     }
 147
 148     return shift;
 149   }
 150
 151 #endif /* !FT_MSB */
 152
 153
 154   /* documentation is in ftcalc.h */
 155
 156   FT_BASE_DEF( FT_Fixed )
 157   FT_Hypot( FT_Fixed  x,
 158             FT_Fixed  y )
 159   {
 160     FT_Vector  v;
 161
 162
 163     v.x = x;
 164     v.y = y;
 165
 166     return FT_Vector_Length( &v );
 167   }
 168
 169
 170 #ifdef FT_LONG64
 171
 172
 173   /* documentation is in freetype.h */
 174
 175   FT_EXPORT_DEF( FT_Long )
 176   FT_MulDiv( FT_Long  a,
 177              FT_Long  b,
 178              FT_Long  c )
 179   {
 180     FT_Int   s = 1;
 181     FT_Long  d;
 182
 183
 184     FT_MOVE_SIGN( a, s );
 185     FT_MOVE_SIGN( b, s );
 186     FT_MOVE_SIGN( c, s );
 187
 188     d = (FT_Long)( c > 0 ? ( (FT_Int64)a * b + ( c >> 1 ) ) / c
 189                          : 0x7FFFFFFFL );
 190
 191     return s < 0 ? -d : d;
 192   }
 193
 194
 195   /* documentation is in ftcalc.h */
 196
 197   FT_BASE_DEF( FT_Long )
 198   FT_MulDiv_No_Round( FT_Long  a,
 199                       FT_Long  b,
 200                       FT_Long  c )
 201   {
 202     FT_Int   s = 1;
 203     FT_Long  d;
 204
 205
 206     FT_MOVE_SIGN( a, s );
 207     FT_MOVE_SIGN( b, s );
 208     FT_MOVE_SIGN( c, s );
 209
 210     d = (FT_Long)( c > 0 ? (FT_Int64)a * b / c
 211                          : 0x7FFFFFFFL );
 212
 213     return s < 0 ? -d : d;
 214   }
 215
 216
 217   /* documentation is in freetype.h */
 218
 219   FT_EXPORT_DEF( FT_Long )
 220   FT_MulFix( FT_Long  a,
 221              FT_Long  b )
 222   {
 223 #ifdef FT_MULFIX_ASSEMBLER
 224
 225     return FT_MULFIX_ASSEMBLER( a, b );
 226
 227 #else
 228
 229     FT_Int   s = 1;
 230     FT_Long  c;
 231
 232
 233     FT_MOVE_SIGN( a, s );
 234     FT_MOVE_SIGN( b, s );
 235
 236     c = (FT_Long)( ( (FT_Int64)a * b + 0x8000L ) >> 16 );
 237
 238     return s < 0 ? -c : c;
 239
 240 #endif /* FT_MULFIX_ASSEMBLER */
 241   }
 242
 243
 244   /* documentation is in freetype.h */
 245
 246   FT_EXPORT_DEF( FT_Long )
 247   FT_DivFix( FT_Long  a,
 248              FT_Long  b )
 249   {
 250     FT_Int   s = 1;
 251     FT_Long  q;
 252
 253
 254     FT_MOVE_SIGN( a, s );
 255     FT_MOVE_SIGN( b, s );
 256
 257     q = (FT_Long)( b > 0 ? ( ( (FT_UInt64)a << 16 ) + ( b >> 1 ) ) / b
 258                          : 0x7FFFFFFFL );
 259
 260     return s < 0 ? -q : q;
 261   }
 262
 263
 264 #else /* !FT_LONG64 */
 265
 266
 267   static void
 268   ft_multo64( FT_UInt32  x,
 269               FT_UInt32  y,
 270               FT_Int64  *z )
 271   {
 272     FT_UInt32  lo1, hi1, lo2, hi2, lo, hi, i1, i2;
 273
 274
 275     lo1 = x & 0x0000FFFFU;  hi1 = x >> 16;
 276     lo2 = y & 0x0000FFFFU;  hi2 = y >> 16;
 277
 278     lo = lo1 * lo2;
 279     i1 = lo1 * hi2;
 280     i2 = lo2 * hi1;
 281     hi = hi1 * hi2;
 282
 283     /* Check carry overflow of i1 + i2 */
 284     i1 += i2;
 285     hi += (FT_UInt32)( i1 < i2 ) << 16;
 286
 287     hi += i1 >> 16;
 288     i1  = i1 << 16;
 289
 290     /* Check carry overflow of i1 + lo */
 291     lo += i1;
 292     hi += ( lo < i1 );
 293
 294     z->lo = lo;
 295     z->hi = hi;
 296   }
 297
 298
 299   static FT_UInt32
 300   ft_div64by32( FT_UInt32  hi,
 301                 FT_UInt32  lo,
 302                 FT_UInt32  y )
 303   {
 304     FT_UInt32  r, q;
 305     FT_Int     i;
 306
 307
 308     if ( hi >= y )
 309       return (FT_UInt32)0x7FFFFFFFL;
 310
 311     /* We shift as many bits as we can into the high register, perform     */
 312     /* 32-bit division with modulo there, then work through the remaining  */
 313     /* bits with long division. This optimization is especially noticeable */
 314     /* for smaller dividends that barely use the high register.            */
 315
 316     i = 31 - FT_MSB( hi );
 317     r = ( hi << i ) | ( lo >> ( 32 - i ) ); lo <<= i; /* left 64-bit shift */
 318     q = r / y;
 319     r -= q * y;   /* remainder */
 320
 321     i = 32 - i;   /* bits remaining in low register */
 322     do
 323     {
 324       q <<= 1;
 325       r   = ( r << 1 ) | ( lo >> 31 ); lo <<= 1;
 326
 327       if ( r >= y )
 328       {
 329         r -= y;
 330         q |= 1;
 331       }
 332     } while ( --i );
 333
 334     return q;
 335   }
 336
 337
 338   static void
 339   FT_Add64( FT_Int64*  x,
 340             FT_Int64*  y,
 341             FT_Int64  *z )
 342   {
 343     FT_UInt32  lo, hi;
 344
 345
 346     lo = x->lo + y->lo;
 347     hi = x->hi + y->hi + ( lo < x->lo );
 348
 349     z->lo = lo;
 350     z->hi = hi;
 351   }
 352
 353
 354   /*  The FT_MulDiv function has been optimized thanks to ideas from     */
 355   /*  Graham Asher and Alexei Podtelezhnikov.  The trick is to optimize  */
 356   /*  a rather common case when everything fits within 32-bits.          */
 357   /*                                                                     */
 358   /*  We compute 'a*b+c/2', then divide it by 'c' (all positive values). */
 359   /*                                                                     */
 360   /*  The product of two positive numbers never exceeds the square of    */
 361   /*  its mean values.  Therefore, we always avoid the overflow by       */
 362   /*  imposing                                                           */
 363   /*                                                                     */
 364   /*    (a + b) / 2 <= sqrt(X - c/2)    ,                                */
 365   /*                                                                     */
 366   /*  where X = 2^32 - 1, the maximum unsigned 32-bit value, and using   */
 367   /*  unsigned arithmetic.  Now we replace `sqrt' with a linear function */
 368   /*  that is smaller or equal for all values of c in the interval       */
 369   /*  [0;X/2]; it should be equal to sqrt(X) and sqrt(3X/4) at the       */
 370   /*  endpoints.  Substituting the linear solution and explicit numbers  */
 371   /*  we get                                                             */
 372   /*                                                                     */
 373   /*    a + b <= 131071.99 - c / 122291.84    .                          */
 374   /*                                                                     */
 375   /*  In practice, we should use a faster and even stronger inequality   */
 376   /*                                                                     */
 377   /*    a + b <= 131071 - (c >> 16)                                      */
 378   /*                                                                     */
 379   /*  or, alternatively,                                                 */
 380   /*                                                                     */
 381   /*    a + b <= 129894 - (c >> 17)    .                                 */
 382   /*                                                                     */
 383   /*  FT_MulFix, on the other hand, is optimized for a small value of    */
 384   /*  the first argument, when the second argument can be much larger.   */
 385   /*  This can be achieved by scaling the second argument and the limit  */
 386   /*  in the above inequalities.  For example,                           */
 387   /*                                                                     */
 388   /*    a + (b >> 8) <= (131071 >> 4)                                    */
 389   /*                                                                     */
 390   /*  covers the practical range of use. The actual test below is a bit  */
 391   /*  tighter to avoid the border case overflows.                        */
 392   /*                                                                     */
 393   /*  In the case of FT_DivFix, the exact overflow check                 */
 394   /*                                                                     */
 395   /*    a << 16 <= X - c/2                                               */
 396   /*                                                                     */
 397   /*  is scaled down by 2^16 and we use                                  */
 398   /*                                                                     */
 399   /*    a <= 65535 - (c >> 17)    .                                      */
 400
 401   /* documentation is in freetype.h */
 402
 403   FT_EXPORT_DEF( FT_Long )
 404   FT_MulDiv( FT_Long  a,
 405              FT_Long  b,
 406              FT_Long  c )
 407   {
 408     FT_Int  s = 1;
 409
 410
 411     /* XXX: this function does not allow 64-bit arguments */
 412     if ( a == 0 || b == c )
 413       return a;
 414
 415     FT_MOVE_SIGN( a, s );
 416     FT_MOVE_SIGN( b, s );
 417     FT_MOVE_SIGN( c, s );
 418
 419     if ( c == 0 )
 420       a = 0x7FFFFFFFL;
 421
 422     else if ( (FT_ULong)a + b <= 129894UL - ( c >> 17 ) )
 423       a = ( (FT_ULong)a * b + ( c >> 1 ) ) / c;
 424
 425     else
 426     {
 427       FT_Int64  temp, temp2;
 428
 429
 430       ft_multo64( a, b, &temp );
 431
 432       temp2.hi = 0;
 433       temp2.lo = c >> 1;
 434
 435       FT_Add64( &temp, &temp2, &temp );
 436
 437       /* last attempt to ditch long division */
 438       a = temp.hi == 0 ? temp.lo / c
 439                        : ft_div64by32( temp.hi, temp.lo, c );
 440     }
 441
 442     return s < 0 ? -a : a;
 443   }
 444
 445
 446   FT_BASE_DEF( FT_Long )
 447   FT_MulDiv_No_Round( FT_Long  a,
 448                       FT_Long  b,
 449                       FT_Long  c )
 450   {
 451     FT_Int  s = 1;
 452
 453
 454     if ( a == 0 || b == c )
 455       return a;
 456
 457     FT_MOVE_SIGN( a, s );
 458     FT_MOVE_SIGN( b, s );
 459     FT_MOVE_SIGN( c, s );
 460
 461     if ( c == 0 )
 462       a = 0x7FFFFFFFL;
 463
 464     else if ( (FT_ULong)a + b <= 131071UL )
 465       a = (FT_ULong)a * b / c;
 466
 467     else
 468     {
 469       FT_Int64  temp;
 470
 471
 472       ft_multo64( a, b, &temp );
 473
 474       /* last attempt to ditch long division */
 475       a = temp.hi == 0 ? temp.lo / c
 476                        : ft_div64by32( temp.hi, temp.lo, c );
 477     }
 478
 479     return s < 0 ? -a : a;
 480   }
 481
 482
 483   /* documentation is in freetype.h */
 484
 485   FT_EXPORT_DEF( FT_Long )
 486   FT_MulFix( FT_Long  a,
 487              FT_Long  b )
 488   {
 489 #ifdef FT_MULFIX_ASSEMBLER
 490
 491     return FT_MULFIX_ASSEMBLER( a, b );
 492
 493 #elif 0
 494
 495     /*
 496      *  This code is nonportable.  See comment below.
 497      *
 498      *  However, on a platform where right-shift of a signed quantity fills
 499      *  the leftmost bits by copying the sign bit, it might be faster.
 500      */
 501
 502     FT_Long   sa, sb;
 503     FT_ULong  ua, ub;
 504
 505
 506     if ( a == 0 || b == 0x10000L )
 507       return a;
 508
 509     /*
 510      *  This is a clever way of converting a signed number `a' into its
 511      *  absolute value (stored back into `a') and its sign.  The sign is
 512      *  stored in `sa'; 0 means `a' was positive or zero, and -1 means `a'
 513      *  was negative.  (Similarly for `b' and `sb').
 514      *
 515      *  Unfortunately, it doesn't work (at least not portably).
 516      *
 517      *  It makes the assumption that right-shift on a negative signed value
 518      *  fills the leftmost bits by copying the sign bit.  This is wrong.
 519      *  According to K&R 2nd ed, section `A7.8 Shift Operators' on page 206,
 520      *  the result of right-shift of a negative signed value is
 521      *  implementation-defined.  At least one implementation fills the
 522      *  leftmost bits with 0s (i.e., it is exactly the same as an unsigned
 523      *  right shift).  This means that when `a' is negative, `sa' ends up
 524      *  with the value 1 rather than -1.  After that, everything else goes
 525      *  wrong.
 526      */
 527     sa = ( a >> ( sizeof ( a ) * 8 - 1 ) );
 528     a  = ( a ^ sa ) - sa;
 529     sb = ( b >> ( sizeof ( b ) * 8 - 1 ) );
 530     b  = ( b ^ sb ) - sb;
 531
 532     ua = (FT_ULong)a;
 533     ub = (FT_ULong)b;
 534
 535     if ( ua + ( ub >> 8 ) <= 8190UL )
 536       ua = ( ua * ub + 0x8000U ) >> 16;
 537     else
 538     {
 539       FT_ULong  al = ua & 0xFFFFU;
 540
 541
 542       ua = ( ua >> 16 ) * ub +  al * ( ub >> 16 ) +
 543            ( ( al * ( ub & 0xFFFFU ) + 0x8000U ) >> 16 );
 544     }
 545
 546     sa ^= sb,
 547     ua  = (FT_ULong)(( ua ^ sa ) - sa);
 548
 549     return (FT_Long)ua;
 550
 551 #else /* 0 */
 552
 553     FT_Int    s = 1;
 554     FT_ULong  ua, ub;
 555
 556
 557     if ( a == 0 || b == 0x10000L )
 558       return a;
 559
 560     FT_MOVE_SIGN( a, s );
 561     FT_MOVE_SIGN( b, s );
 562
 563     ua = (FT_ULong)a;
 564     ub = (FT_ULong)b;
 565
 566     if ( ua + ( ub >> 8 ) <= 8190UL )
 567       ua = ( ua * ub + 0x8000UL ) >> 16;
 568     else
 569     {
 570       FT_ULong  al = ua & 0xFFFFUL;
 571
 572
 573       ua = ( ua >> 16 ) * ub +  al * ( ub >> 16 ) +
 574            ( ( al * ( ub & 0xFFFFUL ) + 0x8000UL ) >> 16 );
 575     }
 576
 577     return s < 0 ? -(FT_Long)ua : (FT_Long)ua;
 578
 579 #endif /* 0 */
 580
 581   }
 582
 583
 584   /* documentation is in freetype.h */
 585
 586   FT_EXPORT_DEF( FT_Long )
 587   FT_DivFix( FT_Long  a,
 588              FT_Long  b )
 589   {
 590     FT_Int   s = 1;
 591     FT_Long  q;
 592
 593
 594     /* XXX: this function does not allow 64-bit arguments */
 595
 596     FT_MOVE_SIGN( a, s );
 597     FT_MOVE_SIGN( b, s );
 598
 599     if ( b == 0 )
 600     {
 601       /* check for division by 0 */
 602       q = 0x7FFFFFFFL;
 603     }
 604     else if ( a <= 65535L - ( b >> 17 ) )
 605     {
 606       /* compute result directly */
 607       q = (FT_Long)( ( ( (FT_ULong)a << 16 ) + ( b >> 1 ) ) / b );
 608     }
 609     else
 610     {
 611       /* we need more bits; we have to do it by hand */
 612       FT_Int64  temp, temp2;
 613
 614
 615       temp.hi  = a >> 16;
 616       temp.lo  = a << 16;
 617       temp2.hi = 0;
 618       temp2.lo = b >> 1;
 619
 620       FT_Add64( &temp, &temp2, &temp );
 621       q = (FT_Long)ft_div64by32( temp.hi, temp.lo, b );
 622     }
 623
 624     return s < 0 ? -q : q;
 625   }
 626
 627
 628 #endif /* FT_LONG64 */
 629
 630
 631   /* documentation is in ftglyph.h */
 632
 633   FT_EXPORT_DEF( void )
 634   FT_Matrix_Multiply( const FT_Matrix*  a,
 635                       FT_Matrix        *b )
 636   {
 637     FT_Fixed  xx, xy, yx, yy;
 638
 639
 640     if ( !a || !b )
 641       return;
 642
 643     xx = FT_MulFix( a->xx, b->xx ) + FT_MulFix( a->xy, b->yx );
 644     xy = FT_MulFix( a->xx, b->xy ) + FT_MulFix( a->xy, b->yy );
 645     yx = FT_MulFix( a->yx, b->xx ) + FT_MulFix( a->yy, b->yx );
 646     yy = FT_MulFix( a->yx, b->xy ) + FT_MulFix( a->yy, b->yy );
 647
 648     b->xx = xx;  b->xy = xy;
 649     b->yx = yx;  b->yy = yy;
 650   }
 651
 652
 653   /* documentation is in ftglyph.h */
 654
 655   FT_EXPORT_DEF( FT_Error )
 656   FT_Matrix_Invert( FT_Matrix*  matrix )
 657   {
 658     FT_Pos  delta, xx, yy;
 659
 660
 661     if ( !matrix )
 662       return FT_THROW( Invalid_Argument );
 663
 664     /* compute discriminant */
 665     delta = FT_MulFix( matrix->xx, matrix->yy ) -
 666             FT_MulFix( matrix->xy, matrix->yx );
 667
 668     if ( !delta )
 669       return FT_THROW( Invalid_Argument );  /* matrix can't be inverted */
 670
 671     matrix->xy = - FT_DivFix( matrix->xy, delta );
 672     matrix->yx = - FT_DivFix( matrix->yx, delta );
 673
 674     xx = matrix->xx;
 675     yy = matrix->yy;
 676
 677     matrix->xx = FT_DivFix( yy, delta );
 678     matrix->yy = FT_DivFix( xx, delta );
 679
 680     return FT_Err_Ok;
 681   }
 682
 683
 684   /* documentation is in ftcalc.h */
 685
 686   FT_BASE_DEF( void )
 687   FT_Matrix_Multiply_Scaled( const FT_Matrix*  a,
 688                              FT_Matrix        *b,
 689                              FT_Long           scaling )
 690   {
 691     FT_Fixed  xx, xy, yx, yy;
 692
 693     FT_Long   val = 0x10000L * scaling;
 694
 695
 696     if ( !a || !b )
 697       return;
 698
 699     xx = FT_MulDiv( a->xx, b->xx, val ) + FT_MulDiv( a->xy, b->yx, val );
 700     xy = FT_MulDiv( a->xx, b->xy, val ) + FT_MulDiv( a->xy, b->yy, val );
 701     yx = FT_MulDiv( a->yx, b->xx, val ) + FT_MulDiv( a->yy, b->yx, val );
 702     yy = FT_MulDiv( a->yx, b->xy, val ) + FT_MulDiv( a->yy, b->yy, val );
 703
 704     b->xx = xx;  b->xy = xy;
 705     b->yx = yx;  b->yy = yy;
 706   }
 707
 708
 709   /* documentation is in ftcalc.h */
 710
 711   FT_BASE_DEF( void )
 712   FT_Vector_Transform_Scaled( FT_Vector*        vector,
 713                               const FT_Matrix*  matrix,
 714                               FT_Long           scaling )
 715   {
 716     FT_Pos   xz, yz;
 717
 718     FT_Long  val = 0x10000L * scaling;
 719
 720
 721     if ( !vector || !matrix )
 722       return;
 723
 724     xz = FT_MulDiv( vector->x, matrix->xx, val ) +
 725          FT_MulDiv( vector->y, matrix->xy, val );
 726
 727     yz = FT_MulDiv( vector->x, matrix->yx, val ) +
 728          FT_MulDiv( vector->y, matrix->yy, val );
 729
 730     vector->x = xz;
 731     vector->y = yz;
 732   }
 733
 734
 735 #if 0
 736
 737   /* documentation is in ftcalc.h */
 738
 739   FT_BASE_DEF( FT_Int32 )
 740   FT_SqrtFixed( FT_Int32  x )
 741   {
 742     FT_UInt32  root, rem_hi, rem_lo, test_div;
 743     FT_Int     count;
 744
 745
 746     root = 0;
 747
 748     if ( x > 0 )
 749     {
 750       rem_hi = 0;
 751       rem_lo = x;
 752       count  = 24;
 753       do
 754       {
 755         rem_hi   = ( rem_hi << 2 ) | ( rem_lo >> 30 );
 756         rem_lo <<= 2;
 757         root   <<= 1;
 758         test_div = ( root << 1 ) + 1;
 759
 760         if ( rem_hi >= test_div )
 761         {
 762           rem_hi -= test_div;
 763           root   += 1;
 764         }
 765       } while ( --count );
 766     }
 767
 768     return (FT_Int32)root;
 769   }
 770
 771 #endif /* 0 */
 772
 773
 774   /* documentation is in ftcalc.h */
 775
 776   FT_BASE_DEF( FT_Int )
 777   ft_corner_orientation( FT_Pos  in_x,
 778                          FT_Pos  in_y,
 779                          FT_Pos  out_x,
 780                          FT_Pos  out_y )
 781   {
 782     FT_Long  result; /* avoid overflow on 16-bit system */
 783
 784
 785     /* deal with the trivial cases quickly */
 786     if ( in_y == 0 )
 787     {
 788       if ( in_x >= 0 )
 789         result = out_y;
 790       else
 791         result = -out_y;
 792     }
 793     else if ( in_x == 0 )
 794     {
 795       if ( in_y >= 0 )
 796         result = -out_x;
 797       else
 798         result = out_x;
 799     }
 800     else if ( out_y == 0 )
 801     {
 802       if ( out_x >= 0 )
 803         result = in_y;
 804       else
 805         result = -in_y;
 806     }
 807     else if ( out_x == 0 )
 808     {
 809       if ( out_y >= 0 )
 810         result = -in_x;
 811       else
 812         result =  in_x;
 813     }
 814     else /* general case */
 815     {
 816 #ifdef FT_LONG64
 817
 818       FT_Int64  delta = (FT_Int64)in_x * out_y - (FT_Int64)in_y * out_x;
 819
 820
 821       if ( delta == 0 )
 822         result = 0;
 823       else
 824         result = 1 - 2 * ( delta < 0 );
 825
 826 #else
 827
 828       FT_Int64  z1, z2;
 829
 830
 831       /* XXX: this function does not allow 64-bit arguments */
 832       ft_multo64( (FT_Int32)in_x, (FT_Int32)out_y, &z1 );
 833       ft_multo64( (FT_Int32)in_y, (FT_Int32)out_x, &z2 );
 834
 835       if ( z1.hi > z2.hi )
 836         result = +1;
 837       else if ( z1.hi < z2.hi )
 838         result = -1;
 839       else if ( z1.lo > z2.lo )
 840         result = +1;
 841       else if ( z1.lo < z2.lo )
 842         result = -1;
 843       else
 844         result = 0;
 845
 846 #endif
 847     }
 848
 849     /* XXX: only the sign of return value, +1/0/-1 must be used */
 850     return (FT_Int)result;
 851   }
 852
 853
 854   /* documentation is in ftcalc.h */
 855
 856   FT_BASE_DEF( FT_Int )
 857   ft_corner_is_flat( FT_Pos  in_x,
 858                      FT_Pos  in_y,
 859                      FT_Pos  out_x,
 860                      FT_Pos  out_y )
 861   {
 862     FT_Pos  ax = in_x + out_x;
 863     FT_Pos  ay = in_y + out_y;
 864
 865     FT_Pos  d_in, d_out, d_hypot;
 866
 867
 868     /* The idea of this function is to compare the length of the */
 869     /* hypotenuse with the `in' and `out' length.  The `corner'  */
 870     /* represented by `in' and `out' is flat if the hypotenuse's */
 871     /* length isn't too large.                                   */
 872     /*                                                           */
 873     /* This approach has the advantage that the angle between    */
 874     /* `in' and `out' is not checked.  In case one of the two    */
 875     /* vectors is `dominant', this is, much larger than the      */
 876     /* other vector, we thus always have a flat corner.          */
 877     /*                                                           */
 878     /*                hypotenuse                                 */
 879     /*       x---------------------------x                       */
 880     /*        \                      /                           */
 881     /*         \                /                                */
 882     /*      in  \          /  out                                */
 883     /*           \    /                                          */
 884     /*            o                                              */
 885     /*              Point                                        */
 886
 887     d_in    = FT_HYPOT(  in_x,  in_y );
 888     d_out   = FT_HYPOT( out_x, out_y );
 889     d_hypot = FT_HYPOT(    ax,    ay );
 890
 891     /* now do a simple length comparison: */
 892     /*                                    */
 893     /*   d_in + d_out < 17/16 d_hypot     */
 894
 895     return ( d_in + d_out - d_hypot ) < ( d_hypot >> 4 );
 896   }
 897
 898
 899 /* END */