ftcalc.c (26846B)
1 /**************************************************************************** 2 * 3 * ftcalc.c 4 * 5 * Arithmetic computations (body). 6 * 7 * Copyright (C) 1996-2025 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 <freetype/ftglyph.h> 36 #include <freetype/fttrigon.h> 37 #include <freetype/internal/ftcalc.h> 38 #include <freetype/internal/ftdebug.h> 39 #include <freetype/internal/ftobjs.h> 40 41 /* cancel inlining macro from internal/ftcalc.h */ 42 #ifdef FT_MulFix 43 # undef FT_MulFix 44 #endif 45 46 47 /************************************************************************** 48 * 49 * The macro FT_COMPONENT is used in trace mode. It is an implicit 50 * parameter of the FT_TRACE() and FT_ERROR() macros, used to print/log 51 * messages during execution. 52 */ 53 #undef FT_COMPONENT 54 #define FT_COMPONENT calc 55 56 57 /* transfer sign, leaving a positive number; */ 58 /* we need an unsigned value to safely negate INT_MIN (or LONG_MIN) */ 59 #define FT_MOVE_SIGN( utype, x, x_unsigned, s ) \ 60 FT_BEGIN_STMNT \ 61 if ( x < 0 ) \ 62 { \ 63 x_unsigned = 0U - (utype)x; \ 64 s = -s; \ 65 } \ 66 else \ 67 x_unsigned = (utype)x; \ 68 FT_END_STMNT 69 70 /* The following three functions are available regardless of whether */ 71 /* FT_INT64 is defined. */ 72 73 /* documentation is in freetype.h */ 74 75 FT_EXPORT_DEF( FT_Fixed ) 76 FT_RoundFix( FT_Fixed a ) 77 { 78 return ADD_LONG( a, 0x8000L - ( a < 0 ) ) & ~0xFFFFL; 79 } 80 81 82 /* documentation is in freetype.h */ 83 84 FT_EXPORT_DEF( FT_Fixed ) 85 FT_CeilFix( FT_Fixed a ) 86 { 87 return ADD_LONG( a, 0xFFFFL ) & ~0xFFFFL; 88 } 89 90 91 /* documentation is in freetype.h */ 92 93 FT_EXPORT_DEF( FT_Fixed ) 94 FT_FloorFix( FT_Fixed a ) 95 { 96 return a & ~0xFFFFL; 97 } 98 99 #ifndef FT_MSB 100 101 FT_BASE_DEF( FT_Int ) 102 FT_MSB( FT_UInt32 z ) 103 { 104 FT_Int shift = 0; 105 106 107 /* determine msb bit index in `shift' */ 108 if ( z & 0xFFFF0000UL ) 109 { 110 z >>= 16; 111 shift += 16; 112 } 113 if ( z & 0x0000FF00UL ) 114 { 115 z >>= 8; 116 shift += 8; 117 } 118 if ( z & 0x000000F0UL ) 119 { 120 z >>= 4; 121 shift += 4; 122 } 123 if ( z & 0x0000000CUL ) 124 { 125 z >>= 2; 126 shift += 2; 127 } 128 if ( z & 0x00000002UL ) 129 { 130 /* z >>= 1; */ 131 shift += 1; 132 } 133 134 return shift; 135 } 136 137 #endif /* !FT_MSB */ 138 139 140 /* documentation is in ftcalc.h */ 141 142 FT_BASE_DEF( FT_Fixed ) 143 FT_Hypot( FT_Fixed x, 144 FT_Fixed y ) 145 { 146 FT_Vector v; 147 148 149 v.x = x; 150 v.y = y; 151 152 return FT_Vector_Length( &v ); 153 } 154 155 156 #ifdef FT_INT64 157 158 159 /* documentation is in freetype.h */ 160 161 FT_EXPORT_DEF( FT_Long ) 162 FT_MulDiv( FT_Long a_, 163 FT_Long b_, 164 FT_Long c_ ) 165 { 166 FT_Int s = 1; 167 FT_UInt64 a, b, c, d; 168 FT_Long d_; 169 170 171 FT_MOVE_SIGN( FT_UInt64, a_, a, s ); 172 FT_MOVE_SIGN( FT_UInt64, b_, b, s ); 173 FT_MOVE_SIGN( FT_UInt64, c_, c, s ); 174 175 d = c > 0 ? ( a * b + ( c >> 1 ) ) / c 176 : 0x7FFFFFFFUL; 177 178 d_ = (FT_Long)d; 179 180 return s < 0 ? NEG_LONG( d_ ) : d_; 181 } 182 183 184 /* documentation is in ftcalc.h */ 185 186 FT_BASE_DEF( FT_Long ) 187 FT_MulDiv_No_Round( FT_Long a_, 188 FT_Long b_, 189 FT_Long c_ ) 190 { 191 FT_Int s = 1; 192 FT_UInt64 a, b, c, d; 193 FT_Long d_; 194 195 196 FT_MOVE_SIGN( FT_UInt64, a_, a, s ); 197 FT_MOVE_SIGN( FT_UInt64, b_, b, s ); 198 FT_MOVE_SIGN( FT_UInt64, c_, c, s ); 199 200 d = c > 0 ? a * b / c 201 : 0x7FFFFFFFUL; 202 203 d_ = (FT_Long)d; 204 205 return s < 0 ? NEG_LONG( d_ ) : d_; 206 } 207 208 209 /* documentation is in freetype.h */ 210 211 FT_EXPORT_DEF( FT_Long ) 212 FT_MulFix( FT_Long a_, 213 FT_Long b_ ) 214 { 215 #ifdef FT_CONFIG_OPTION_INLINE_MULFIX 216 217 return FT_MulFix_64( a_, b_ ); 218 219 #else 220 221 FT_Int64 ab = MUL_INT64( a_, b_ ); 222 223 /* this requires arithmetic right shift of signed numbers */ 224 return (FT_Long)( ( ab + 0x8000L + ( ab >> 63 ) ) >> 16 ); 225 226 #endif /* FT_CONFIG_OPTION_INLINE_MULFIX */ 227 } 228 229 230 /* documentation is in freetype.h */ 231 232 FT_EXPORT_DEF( FT_Long ) 233 FT_DivFix( FT_Long a_, 234 FT_Long b_ ) 235 { 236 FT_Int s = 1; 237 FT_UInt64 a, b, q; 238 FT_Long q_; 239 240 241 FT_MOVE_SIGN( FT_UInt64, a_, a, s ); 242 FT_MOVE_SIGN( FT_UInt64, b_, b, s ); 243 244 q = b > 0 ? ( ( a << 16 ) + ( b >> 1 ) ) / b 245 : 0x7FFFFFFFUL; 246 247 q_ = (FT_Long)q; 248 249 return s < 0 ? NEG_LONG( q_ ) : q_; 250 } 251 252 253 #else /* !FT_INT64 */ 254 255 256 static void 257 ft_multo64( FT_UInt32 x, 258 FT_UInt32 y, 259 FT_Int64 *z ) 260 { 261 FT_UInt32 lo1, hi1, lo2, hi2, lo, hi, i1, i2; 262 263 264 lo1 = x & 0x0000FFFFU; hi1 = x >> 16; 265 lo2 = y & 0x0000FFFFU; hi2 = y >> 16; 266 267 lo = lo1 * lo2; 268 i1 = lo1 * hi2; 269 i2 = lo2 * hi1; 270 hi = hi1 * hi2; 271 272 /* Check carry overflow of i1 + i2 */ 273 i1 += i2; 274 hi += (FT_UInt32)( i1 < i2 ) << 16; 275 276 hi += i1 >> 16; 277 i1 = i1 << 16; 278 279 /* Check carry overflow of i1 + lo */ 280 lo += i1; 281 hi += ( lo < i1 ); 282 283 z->lo = lo; 284 z->hi = hi; 285 } 286 287 288 static FT_UInt32 289 ft_div64by32( FT_UInt32 hi, 290 FT_UInt32 lo, 291 FT_UInt32 y ) 292 { 293 FT_UInt32 r, q; 294 FT_Int i; 295 296 297 if ( hi >= y ) 298 return (FT_UInt32)0x7FFFFFFFL; 299 300 /* We shift as many bits as we can into the high register, perform */ 301 /* 32-bit division with modulo there, then work through the remaining */ 302 /* bits with long division. This optimization is especially noticeable */ 303 /* for smaller dividends that barely use the high register. */ 304 305 i = 31 - FT_MSB( hi ); 306 r = ( hi << i ) | ( lo >> ( 32 - i ) ); lo <<= i; /* left 64-bit shift */ 307 q = r / y; 308 r -= q * y; /* remainder */ 309 310 i = 32 - i; /* bits remaining in low register */ 311 do 312 { 313 q <<= 1; 314 r = ( r << 1 ) | ( lo >> 31 ); lo <<= 1; 315 316 if ( r >= y ) 317 { 318 r -= y; 319 q |= 1; 320 } 321 } while ( --i ); 322 323 return q; 324 } 325 326 327 static void 328 FT_Add64( FT_Int64* x, 329 FT_Int64* y, 330 FT_Int64 *z ) 331 { 332 FT_UInt32 lo, hi; 333 334 335 lo = x->lo + y->lo; 336 hi = x->hi + y->hi + ( lo < x->lo ); 337 338 z->lo = lo; 339 z->hi = hi; 340 } 341 342 343 /* The FT_MulDiv function has been optimized thanks to ideas from */ 344 /* Graham Asher and Alexei Podtelezhnikov. The trick is to optimize */ 345 /* a rather common case when everything fits within 32-bits. */ 346 /* */ 347 /* We compute 'a*b+c/2', then divide it by 'c' (all positive values). */ 348 /* */ 349 /* The product of two positive numbers never exceeds the square of */ 350 /* its mean values. Therefore, we always avoid the overflow by */ 351 /* imposing */ 352 /* */ 353 /* (a + b) / 2 <= sqrt(X - c/2) , */ 354 /* */ 355 /* where X = 2^32 - 1, the maximum unsigned 32-bit value, and using */ 356 /* unsigned arithmetic. Now we replace `sqrt' with a linear function */ 357 /* that is smaller or equal for all values of c in the interval */ 358 /* [0;X/2]; it should be equal to sqrt(X) and sqrt(3X/4) at the */ 359 /* endpoints. Substituting the linear solution and explicit numbers */ 360 /* we get */ 361 /* */ 362 /* a + b <= 131071.99 - c / 122291.84 . */ 363 /* */ 364 /* In practice, we should use a faster and even stronger inequality */ 365 /* */ 366 /* a + b <= 131071 - (c >> 16) */ 367 /* */ 368 /* or, alternatively, */ 369 /* */ 370 /* a + b <= 129894 - (c >> 17) . */ 371 /* */ 372 /* FT_MulFix, on the other hand, is optimized for a small value of */ 373 /* the first argument, when the second argument can be much larger. */ 374 /* This can be achieved by scaling the second argument and the limit */ 375 /* in the above inequalities. For example, */ 376 /* */ 377 /* a + (b >> 8) <= (131071 >> 4) */ 378 /* */ 379 /* covers the practical range of use. The actual test below is a bit */ 380 /* tighter to avoid the border case overflows. */ 381 /* */ 382 /* In the case of FT_DivFix, the exact overflow check */ 383 /* */ 384 /* a << 16 <= X - c/2 */ 385 /* */ 386 /* is scaled down by 2^16 and we use */ 387 /* */ 388 /* a <= 65535 - (c >> 17) . */ 389 390 /* documentation is in freetype.h */ 391 392 FT_EXPORT_DEF( FT_Long ) 393 FT_MulDiv( FT_Long a_, 394 FT_Long b_, 395 FT_Long c_ ) 396 { 397 FT_Int s = 1; 398 FT_UInt32 a, b, c; 399 400 401 /* XXX: this function does not allow 64-bit arguments */ 402 403 FT_MOVE_SIGN( FT_UInt32, a_, a, s ); 404 FT_MOVE_SIGN( FT_UInt32, b_, b, s ); 405 FT_MOVE_SIGN( FT_UInt32, c_, c, s ); 406 407 if ( c == 0 ) 408 a = 0x7FFFFFFFUL; 409 410 else if ( a + b <= 129894UL - ( c >> 17 ) ) 411 a = ( a * b + ( c >> 1 ) ) / c; 412 413 else 414 { 415 FT_Int64 temp, temp2; 416 417 418 ft_multo64( a, b, &temp ); 419 420 temp2.hi = 0; 421 temp2.lo = c >> 1; 422 423 FT_Add64( &temp, &temp2, &temp ); 424 425 /* last attempt to ditch long division */ 426 a = ( temp.hi == 0 ) ? temp.lo / c 427 : ft_div64by32( temp.hi, temp.lo, c ); 428 } 429 430 a_ = (FT_Long)a; 431 432 return s < 0 ? NEG_LONG( a_ ) : a_; 433 } 434 435 436 FT_BASE_DEF( FT_Long ) 437 FT_MulDiv_No_Round( FT_Long a_, 438 FT_Long b_, 439 FT_Long c_ ) 440 { 441 FT_Int s = 1; 442 FT_UInt32 a, b, c; 443 444 445 /* XXX: this function does not allow 64-bit arguments */ 446 447 FT_MOVE_SIGN( FT_UInt32, a_, a, s ); 448 FT_MOVE_SIGN( FT_UInt32, b_, b, s ); 449 FT_MOVE_SIGN( FT_UInt32, c_, c, s ); 450 451 if ( c == 0 ) 452 a = 0x7FFFFFFFUL; 453 454 else if ( a + b <= 131071UL ) 455 a = a * b / c; 456 457 else 458 { 459 FT_Int64 temp; 460 461 462 ft_multo64( a, b, &temp ); 463 464 /* last attempt to ditch long division */ 465 a = ( temp.hi == 0 ) ? temp.lo / c 466 : ft_div64by32( temp.hi, temp.lo, c ); 467 } 468 469 a_ = (FT_Long)a; 470 471 return s < 0 ? NEG_LONG( a_ ) : a_; 472 } 473 474 475 /* documentation is in freetype.h */ 476 477 FT_EXPORT_DEF( FT_Long ) 478 FT_MulFix( FT_Long a_, 479 FT_Long b_ ) 480 { 481 #ifdef FT_MULFIX_ASSEMBLER 482 483 return FT_MULFIX_ASSEMBLER( a_, b_ ); 484 485 #elif 0 486 487 /* 488 * This code is nonportable. See comment below. 489 * 490 * However, on a platform where right-shift of a signed quantity fills 491 * the leftmost bits by copying the sign bit, it might be faster. 492 */ 493 494 FT_Long sa, sb; 495 FT_UInt32 a, b; 496 497 498 /* 499 * This is a clever way of converting a signed number `a' into its 500 * absolute value (stored back into `a') and its sign. The sign is 501 * stored in `sa'; 0 means `a' was positive or zero, and -1 means `a' 502 * was negative. (Similarly for `b' and `sb'). 503 * 504 * Unfortunately, it doesn't work (at least not portably). 505 * 506 * It makes the assumption that right-shift on a negative signed value 507 * fills the leftmost bits by copying the sign bit. This is wrong. 508 * According to K&R 2nd ed, section `A7.8 Shift Operators' on page 206, 509 * the result of right-shift of a negative signed value is 510 * implementation-defined. At least one implementation fills the 511 * leftmost bits with 0s (i.e., it is exactly the same as an unsigned 512 * right shift). This means that when `a' is negative, `sa' ends up 513 * with the value 1 rather than -1. After that, everything else goes 514 * wrong. 515 */ 516 sa = ( a_ >> ( sizeof ( a_ ) * 8 - 1 ) ); 517 a = ( a_ ^ sa ) - sa; 518 sb = ( b_ >> ( sizeof ( b_ ) * 8 - 1 ) ); 519 b = ( b_ ^ sb ) - sb; 520 521 a = (FT_UInt32)a_; 522 b = (FT_UInt32)b_; 523 524 if ( a + ( b >> 8 ) <= 8190UL ) 525 a = ( a * b + 0x8000U ) >> 16; 526 else 527 { 528 FT_UInt32 al = a & 0xFFFFUL; 529 530 531 a = ( a >> 16 ) * b + al * ( b >> 16 ) + 532 ( ( al * ( b & 0xFFFFUL ) + 0x8000UL ) >> 16 ); 533 } 534 535 sa ^= sb; 536 a = ( a ^ sa ) - sa; 537 538 return (FT_Long)a; 539 540 #else /* 0 */ 541 542 FT_Int s = 1; 543 FT_UInt32 a, b; 544 545 546 /* XXX: this function does not allow 64-bit arguments */ 547 548 FT_MOVE_SIGN( FT_UInt32, a_, a, s ); 549 FT_MOVE_SIGN( FT_UInt32, b_, b, s ); 550 551 if ( a + ( b >> 8 ) <= 8190UL ) 552 a = ( a * b + 0x8000UL ) >> 16; 553 else 554 { 555 FT_UInt32 al = a & 0xFFFFUL; 556 557 558 a = ( a >> 16 ) * b + al * ( b >> 16 ) + 559 ( ( al * ( b & 0xFFFFUL ) + 0x8000UL ) >> 16 ); 560 } 561 562 a_ = (FT_Long)a; 563 564 return s < 0 ? NEG_LONG( a_ ) : a_; 565 566 #endif /* 0 */ 567 568 } 569 570 571 /* documentation is in freetype.h */ 572 573 FT_EXPORT_DEF( FT_Long ) 574 FT_DivFix( FT_Long a_, 575 FT_Long b_ ) 576 { 577 FT_Int s = 1; 578 FT_UInt32 a, b, q; 579 FT_Long q_; 580 581 582 /* XXX: this function does not allow 64-bit arguments */ 583 584 FT_MOVE_SIGN( FT_UInt32, a_, a, s ); 585 FT_MOVE_SIGN( FT_UInt32, b_, b, s ); 586 587 if ( b == 0 ) 588 { 589 /* check for division by 0 */ 590 q = 0x7FFFFFFFUL; 591 } 592 else if ( a <= 65535UL - ( b >> 17 ) ) 593 { 594 /* compute result directly */ 595 q = ( ( a << 16 ) + ( b >> 1 ) ) / b; 596 } 597 else 598 { 599 /* we need more bits; we have to do it by hand */ 600 FT_Int64 temp, temp2; 601 602 603 temp.hi = a >> 16; 604 temp.lo = a << 16; 605 temp2.hi = 0; 606 temp2.lo = b >> 1; 607 608 FT_Add64( &temp, &temp2, &temp ); 609 q = ft_div64by32( temp.hi, temp.lo, b ); 610 } 611 612 q_ = (FT_Long)q; 613 614 return s < 0 ? NEG_LONG( q_ ) : q_; 615 } 616 617 618 #endif /* !FT_INT64 */ 619 620 621 /* documentation is in ftglyph.h */ 622 623 FT_EXPORT_DEF( void ) 624 FT_Matrix_Multiply( const FT_Matrix* a, 625 FT_Matrix *b ) 626 { 627 FT_Fixed xx, xy, yx, yy; 628 629 630 if ( !a || !b ) 631 return; 632 633 xx = ADD_LONG( FT_MulFix( a->xx, b->xx ), 634 FT_MulFix( a->xy, b->yx ) ); 635 xy = ADD_LONG( FT_MulFix( a->xx, b->xy ), 636 FT_MulFix( a->xy, b->yy ) ); 637 yx = ADD_LONG( FT_MulFix( a->yx, b->xx ), 638 FT_MulFix( a->yy, b->yx ) ); 639 yy = ADD_LONG( FT_MulFix( a->yx, b->xy ), 640 FT_MulFix( a->yy, b->yy ) ); 641 642 b->xx = xx; 643 b->xy = xy; 644 b->yx = yx; 645 b->yy = yy; 646 } 647 648 649 /* documentation is in ftglyph.h */ 650 651 FT_EXPORT_DEF( FT_Error ) 652 FT_Matrix_Invert( FT_Matrix* matrix ) 653 { 654 FT_Pos delta, xx, yy; 655 656 657 if ( !matrix ) 658 return FT_THROW( Invalid_Argument ); 659 660 /* compute discriminant */ 661 delta = FT_MulFix( matrix->xx, matrix->yy ) - 662 FT_MulFix( matrix->xy, matrix->yx ); 663 664 if ( !delta ) 665 return FT_THROW( Invalid_Argument ); /* matrix can't be inverted */ 666 667 matrix->xy = -FT_DivFix( matrix->xy, delta ); 668 matrix->yx = -FT_DivFix( matrix->yx, delta ); 669 670 xx = matrix->xx; 671 yy = matrix->yy; 672 673 matrix->xx = FT_DivFix( yy, delta ); 674 matrix->yy = FT_DivFix( xx, delta ); 675 676 return FT_Err_Ok; 677 } 678 679 680 /* documentation is in ftcalc.h */ 681 682 FT_BASE_DEF( void ) 683 FT_Matrix_Multiply_Scaled( const FT_Matrix* a, 684 FT_Matrix *b, 685 FT_Long scaling ) 686 { 687 FT_Fixed xx, xy, yx, yy; 688 689 FT_Long val = 0x10000L * scaling; 690 691 692 if ( !a || !b ) 693 return; 694 695 xx = ADD_LONG( FT_MulDiv( a->xx, b->xx, val ), 696 FT_MulDiv( a->xy, b->yx, val ) ); 697 xy = ADD_LONG( FT_MulDiv( a->xx, b->xy, val ), 698 FT_MulDiv( a->xy, b->yy, val ) ); 699 yx = ADD_LONG( FT_MulDiv( a->yx, b->xx, val ), 700 FT_MulDiv( a->yy, b->yx, val ) ); 701 yy = ADD_LONG( FT_MulDiv( a->yx, b->xy, val ), 702 FT_MulDiv( a->yy, b->yy, val ) ); 703 704 b->xx = xx; 705 b->xy = xy; 706 b->yx = yx; 707 b->yy = yy; 708 } 709 710 711 /* documentation is in ftcalc.h */ 712 713 FT_BASE_DEF( FT_Bool ) 714 FT_Matrix_Check( const FT_Matrix* matrix ) 715 { 716 FT_Fixed xx, xy, yx, yy; 717 FT_Fixed val; 718 FT_Int shift; 719 FT_ULong temp1, temp2; 720 721 722 if ( !matrix ) 723 return 0; 724 725 xx = matrix->xx; 726 xy = matrix->xy; 727 yx = matrix->yx; 728 yy = matrix->yy; 729 val = FT_ABS( xx ) | FT_ABS( xy ) | FT_ABS( yx ) | FT_ABS( yy ); 730 731 /* we only handle non-zero 32-bit values */ 732 if ( !val || val > 0x7FFFFFFFL ) 733 return 0; 734 735 /* Scale matrix to avoid the temp1 overflow, which is */ 736 /* more stringent than avoiding the temp2 overflow. */ 737 738 shift = FT_MSB( val ) - 12; 739 740 if ( shift > 0 ) 741 { 742 xx >>= shift; 743 xy >>= shift; 744 yx >>= shift; 745 yy >>= shift; 746 } 747 748 temp1 = 32U * (FT_ULong)FT_ABS( xx * yy - xy * yx ); 749 temp2 = (FT_ULong)( xx * xx ) + (FT_ULong)( xy * xy ) + 750 (FT_ULong)( yx * yx ) + (FT_ULong)( yy * yy ); 751 752 if ( temp1 <= temp2 ) 753 return 0; 754 755 return 1; 756 } 757 758 759 /* documentation is in ftcalc.h */ 760 761 FT_BASE_DEF( void ) 762 FT_Vector_Transform_Scaled( FT_Vector* vector, 763 const FT_Matrix* matrix, 764 FT_Long scaling ) 765 { 766 FT_Pos xz, yz; 767 768 FT_Long val = 0x10000L * scaling; 769 770 771 if ( !vector || !matrix ) 772 return; 773 774 xz = ADD_LONG( FT_MulDiv( vector->x, matrix->xx, val ), 775 FT_MulDiv( vector->y, matrix->xy, val ) ); 776 yz = ADD_LONG( FT_MulDiv( vector->x, matrix->yx, val ), 777 FT_MulDiv( vector->y, matrix->yy, val ) ); 778 779 vector->x = xz; 780 vector->y = yz; 781 } 782 783 784 /* documentation is in ftcalc.h */ 785 786 FT_BASE_DEF( FT_UInt32 ) 787 FT_Vector_NormLen( FT_Vector* vector ) 788 { 789 FT_Int32 x_ = vector->x; 790 FT_Int32 y_ = vector->y; 791 FT_Int32 b, z; 792 FT_UInt32 x, y, u, v, l; 793 FT_Int sx = 1, sy = 1, shift; 794 795 796 FT_MOVE_SIGN( FT_UInt32, x_, x, sx ); 797 FT_MOVE_SIGN( FT_UInt32, y_, y, sy ); 798 799 /* trivial cases */ 800 if ( x == 0 ) 801 { 802 if ( y > 0 ) 803 vector->y = sy * 0x10000; 804 return y; 805 } 806 else if ( y == 0 ) 807 { 808 if ( x > 0 ) 809 vector->x = sx * 0x10000; 810 return x; 811 } 812 813 /* Estimate length and prenormalize by shifting so that */ 814 /* the new approximate length is between 2/3 and 4/3. */ 815 /* The magic constant 0xAAAAAAAAUL (2/3 of 2^32) helps */ 816 /* achieve this in 16.16 fixed-point representation. */ 817 l = x > y ? x + ( y >> 1 ) 818 : y + ( x >> 1 ); 819 820 shift = 31 - FT_MSB( l ); 821 shift -= 15 + ( l >= ( 0xAAAAAAAAUL >> shift ) ); 822 823 if ( shift > 0 ) 824 { 825 x <<= shift; 826 y <<= shift; 827 828 /* re-estimate length for tiny vectors */ 829 l = x > y ? x + ( y >> 1 ) 830 : y + ( x >> 1 ); 831 } 832 else 833 { 834 x >>= -shift; 835 y >>= -shift; 836 l >>= -shift; 837 } 838 839 /* lower linear approximation for reciprocal length minus one */ 840 b = 0x10000 - (FT_Int32)l; 841 842 x_ = (FT_Int32)x; 843 y_ = (FT_Int32)y; 844 845 /* Newton's iterations */ 846 do 847 { 848 u = (FT_UInt32)( x_ + ( x_ * b >> 16 ) ); 849 v = (FT_UInt32)( y_ + ( y_ * b >> 16 ) ); 850 851 /* Normalized squared length in the parentheses approaches 2^32. */ 852 /* On two's complement systems, converting to signed gives the */ 853 /* difference with 2^32 even if the expression wraps around. */ 854 z = -(FT_Int32)( u * u + v * v ) / 0x200; 855 z = z * ( ( 0x10000 + b ) >> 8 ) / 0x10000; 856 857 b += z; 858 859 } while ( z > 0 ); 860 861 vector->x = sx < 0 ? -(FT_Pos)u : (FT_Pos)u; 862 vector->y = sy < 0 ? -(FT_Pos)v : (FT_Pos)v; 863 864 /* Conversion to signed helps to recover from likely wrap around */ 865 /* in calculating the prenormalized length, because it gives the */ 866 /* correct difference with 2^32 on two's complement systems. */ 867 l = (FT_UInt32)( 0x10000 + (FT_Int32)( u * x + v * y ) / 0x10000 ); 868 if ( shift > 0 ) 869 l = ( l + ( 1 << ( shift - 1 ) ) ) >> shift; 870 else 871 l <<= -shift; 872 873 return l; 874 } 875 876 877 /* documentation is in ftcalc.h */ 878 879 FT_BASE_DEF( FT_UInt32 ) 880 FT_SqrtFixed( FT_UInt32 v ) 881 { 882 if ( v == 0 ) 883 return 0; 884 885 #ifndef FT_INT64 886 887 /* Algorithm by Christophe Meessen (1993) with overflow fixed and */ 888 /* rounding added. Any unsigned fixed 16.16 argument is acceptable. */ 889 /* However, this algorithm is slower than the Babylonian method with */ 890 /* a good initial guess. We only use it for large 32-bit values when */ 891 /* 64-bit computations are not desirable. */ 892 else if ( v > 0x10000U ) 893 { 894 FT_UInt32 r = v >> 1; 895 FT_UInt32 q = ( v & 1 ) << 15; 896 FT_UInt32 b = 0x20000000; 897 FT_UInt32 t; 898 899 900 do 901 { 902 t = q + b; 903 if ( r >= t ) 904 { 905 r -= t; 906 q = t + b; /* equivalent to q += 2*b */ 907 } 908 r <<= 1; 909 b >>= 1; 910 911 } while ( b > 0x10 ); /* exactly 25 cycles */ 912 913 return ( q + 0x40 ) >> 7; 914 } 915 else 916 { 917 FT_UInt32 r = ( v << 16 ) - 1; 918 919 #else /* FT_INT64 */ 920 921 else 922 { 923 FT_UInt64 r = ( (FT_UInt64)v << 16 ) - 1; 924 925 #endif /* FT_INT64 */ 926 927 FT_UInt32 q = 1 << ( ( 17 + FT_MSB( v ) ) >> 1 ); 928 FT_UInt32 t; 929 930 931 /* Babylonian method with rounded-up division */ 932 do 933 { 934 t = q; 935 q = ( t + (FT_UInt32)( r / t ) + 1 ) >> 1; 936 937 } while ( q != t ); /* less than 6 cycles */ 938 939 return q; 940 } 941 } 942 943 944 /* documentation is in ftcalc.h */ 945 946 FT_BASE_DEF( FT_Int ) 947 ft_corner_orientation( FT_Pos in_x, 948 FT_Pos in_y, 949 FT_Pos out_x, 950 FT_Pos out_y ) 951 { 952 /* we silently ignore overflow errors since such large values */ 953 /* lead to even more (harmless) rendering errors later on */ 954 955 #ifdef FT_INT64 956 957 FT_Int64 delta = SUB_INT64( MUL_INT64( in_x, out_y ), 958 MUL_INT64( in_y, out_x ) ); 959 960 961 return ( delta > 0 ) - ( delta < 0 ); 962 963 #else 964 965 FT_Int64 z1, z2; 966 FT_Int result; 967 968 969 if ( (FT_ULong)FT_ABS( in_x ) + (FT_ULong)FT_ABS( out_y ) <= 92681UL ) 970 { 971 z1.lo = (FT_UInt32)in_x * (FT_UInt32)out_y; 972 z1.hi = (FT_UInt32)( (FT_Int32)z1.lo >> 31 ); /* sign-expansion */ 973 } 974 else 975 ft_multo64( (FT_UInt32)in_x, (FT_UInt32)out_y, &z1 ); 976 977 if ( (FT_ULong)FT_ABS( in_y ) + (FT_ULong)FT_ABS( out_x ) <= 92681UL ) 978 { 979 z2.lo = (FT_UInt32)in_y * (FT_UInt32)out_x; 980 z2.hi = (FT_UInt32)( (FT_Int32)z2.lo >> 31 ); /* sign-expansion */ 981 } 982 else 983 ft_multo64( (FT_UInt32)in_y, (FT_UInt32)out_x, &z2 ); 984 985 if ( (FT_Int32)z1.hi > (FT_Int32)z2.hi ) 986 result = +1; 987 else if ( (FT_Int32)z1.hi < (FT_Int32)z2.hi ) 988 result = -1; 989 else if ( z1.lo > z2.lo ) 990 result = +1; 991 else if ( z1.lo < z2.lo ) 992 result = -1; 993 else 994 result = 0; 995 996 /* XXX: only the sign of return value, +1/0/-1 must be used */ 997 return result; 998 999 #endif 1000 } 1001 1002 1003 /* documentation is in ftcalc.h */ 1004 1005 FT_BASE_DEF( FT_Int ) 1006 ft_corner_is_flat( FT_Pos in_x, 1007 FT_Pos in_y, 1008 FT_Pos out_x, 1009 FT_Pos out_y ) 1010 { 1011 FT_Pos ax = in_x + out_x; 1012 FT_Pos ay = in_y + out_y; 1013 1014 FT_Pos d_in, d_out, d_hypot; 1015 1016 1017 /* The idea of this function is to compare the length of the */ 1018 /* hypotenuse with the `in' and `out' length. The `corner' */ 1019 /* represented by `in' and `out' is flat if the hypotenuse's */ 1020 /* length isn't too large. */ 1021 /* */ 1022 /* This approach has the advantage that the angle between */ 1023 /* `in' and `out' is not checked. In case one of the two */ 1024 /* vectors is `dominant', that is, much larger than the */ 1025 /* other vector, we thus always have a flat corner. */ 1026 /* */ 1027 /* hypotenuse */ 1028 /* x---------------------------x */ 1029 /* \ / */ 1030 /* \ / */ 1031 /* in \ / out */ 1032 /* \ / */ 1033 /* o */ 1034 /* Point */ 1035 1036 d_in = FT_HYPOT( in_x, in_y ); 1037 d_out = FT_HYPOT( out_x, out_y ); 1038 d_hypot = FT_HYPOT( ax, ay ); 1039 1040 /* now do a simple length comparison: */ 1041 /* */ 1042 /* d_in + d_out < 17/16 d_hypot */ 1043 1044 return ( d_in + d_out - d_hypot ) < ( d_hypot >> 4 ); 1045 } 1046 1047 1048 /* END */