src/base/fttrigon.c

   1 /***************************************************************************/
   2 /*                                                                         */
   3 /*  fttrigon.c                                                             */
   4 /*                                                                         */
   5 /*    FreeType trigonometric functions (body).                             */
   6 /*                                                                         */
   7 /*  Copyright 2001-2016 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   /* This is a fixed-point CORDIC implementation of trigonometric          */
  21   /* functions as well as transformations between Cartesian and polar      */
  22   /* coordinates.  The angles are represented as 16.16 fixed-point values  */
  23   /* in degrees, i.e., the angular resolution is 2^-16 degrees.  Note that */
  24   /* only vectors longer than 2^16*180/pi (or at least 22 bits) on a       */
  25   /* discrete Cartesian grid can have the same or better angular           */
  26   /* resolution.  Therefore, to maintain this precision, some functions    */
  27   /* require an interim upscaling of the vectors, whereas others operate   */
  28   /* with 24-bit long vectors directly.                                    */
  29   /*                                                                       */
  30   /*************************************************************************/
  31
  32 #include <ft2build.h>
  33 #include FT_INTERNAL_OBJECTS_H
  34 #include FT_INTERNAL_CALC_H
  35 #include FT_TRIGONOMETRY_H
  36
  37
  38   /* the Cordic shrink factor 0.858785336480436 * 2^32 */
  39 #define FT_TRIG_SCALE      0xDBD95B16UL
  40
  41   /* the highest bit in overflow-safe vector components, */
  42   /* MSB of 0.858785336480436 * sqrt(0.5) * 2^30         */
  43 #define FT_TRIG_SAFE_MSB   29
  44
  45   /* this table was generated for FT_PI = 180L << 16, i.e. degrees */
  46 #define FT_TRIG_MAX_ITERS  23
  47
  48   static const FT_Angle
  49   ft_trig_arctan_table[] =
  50   {
  51     1740967L, 919879L, 466945L, 234379L, 117304L, 58666L, 29335L,
  52     14668L, 7334L, 3667L, 1833L, 917L, 458L, 229L, 115L,
  53     57L, 29L, 14L, 7L, 4L, 2L, 1L
  54   };
  55
  56
  57 #ifdef FT_LONG64
  58
  59   /* multiply a given value by the CORDIC shrink factor */
  60   static FT_Fixed
  61   ft_trig_downscale( FT_Fixed  val )
  62   {
  63     FT_Int  s = 1;
  64
  65
  66     if ( val < 0 )
  67     {
  68        val = -val;
  69        s = -1;
  70     }
  71
  72     /* 0x40000000 comes from regression analysis between true */
  73     /* and CORDIC hypotenuse, so it minimizes the error       */
  74     val = (FT_Fixed)(
  75             ( (FT_UInt64)val * FT_TRIG_SCALE + 0x40000000UL ) >> 32 );
  76
  77     return s < 0 ? -val : val;
  78   }
  79
  80 #else /* !FT_LONG64 */
  81
  82   /* multiply a given value by the CORDIC shrink factor */
  83   static FT_Fixed
  84   ft_trig_downscale( FT_Fixed  val )
  85   {
  86     FT_Int     s = 1;
  87     FT_UInt32  lo1, hi1, lo2, hi2, lo, hi, i1, i2;
  88
  89
  90     if ( val < 0 )
  91     {
  92        val = -val;
  93        s = -1;
  94     }
  95
  96     lo1 = (FT_UInt32)val & 0x0000FFFFU;
  97     hi1 = (FT_UInt32)val >> 16;
  98     lo2 = FT_TRIG_SCALE & 0x0000FFFFU;
  99     hi2 = FT_TRIG_SCALE >> 16;
 100
 101     lo = lo1 * lo2;
 102     i1 = lo1 * hi2;
 103     i2 = lo2 * hi1;
 104     hi = hi1 * hi2;
 105
 106     /* Check carry overflow of i1 + i2 */
 107     i1 += i2;
 108     hi += (FT_UInt32)( i1 < i2 ) << 16;
 109
 110     hi += i1 >> 16;
 111     i1  = i1 << 16;
 112
 113     /* Check carry overflow of i1 + lo */
 114     lo += i1;
 115     hi += ( lo < i1 );
 116
 117     /* 0x40000000 comes from regression analysis between true */
 118     /* and CORDIC hypotenuse, so it minimizes the error       */
 119
 120     /* Check carry overflow of lo + 0x40000000 */
 121     lo += 0x40000000UL;
 122     hi += ( lo < 0x40000000UL );
 123
 124     val = (FT_Fixed)hi;
 125
 126     return s < 0 ? -val : val;
 127   }
 128
 129 #endif /* !FT_LONG64 */
 130
 131
 132   /* undefined and never called for zero vector */
 133   static FT_Int
 134   ft_trig_prenorm( FT_Vector*  vec )
 135   {
 136     FT_Pos  x, y;
 137     FT_Int  shift;
 138
 139
 140     x = vec->x;
 141     y = vec->y;
 142
 143     shift = FT_MSB( (FT_UInt32)( FT_ABS( x ) | FT_ABS( y ) ) );
 144
 145     if ( shift <= FT_TRIG_SAFE_MSB )
 146     {
 147       shift  = FT_TRIG_SAFE_MSB - shift;
 148       vec->x = (FT_Pos)( (FT_ULong)x << shift );
 149       vec->y = (FT_Pos)( (FT_ULong)y << shift );
 150     }
 151     else
 152     {
 153       shift -= FT_TRIG_SAFE_MSB;
 154       vec->x = x >> shift;
 155       vec->y = y >> shift;
 156       shift  = -shift;
 157     }
 158
 159     return shift;
 160   }
 161
 162
 163   static void
 164   ft_trig_pseudo_rotate( FT_Vector*  vec,
 165                          FT_Angle    theta )
 166   {
 167     FT_Int           i;
 168     FT_Fixed         x, y, xtemp, b;
 169     const FT_Angle  *arctanptr;
 170
 171
 172     x = vec->x;
 173     y = vec->y;
 174
 175     /* Rotate inside [-PI/4,PI/4] sector */
 176     while ( theta < -FT_ANGLE_PI4 )
 177     {
 178       xtemp  =  y;
 179       y      = -x;
 180       x      =  xtemp;
 181       theta +=  FT_ANGLE_PI2;
 182     }
 183
 184     while ( theta > FT_ANGLE_PI4 )
 185     {
 186       xtemp  = -y;
 187       y      =  x;
 188       x      =  xtemp;
 189       theta -=  FT_ANGLE_PI2;
 190     }
 191
 192     arctanptr = ft_trig_arctan_table;
 193
 194     /* Pseudorotations, with right shifts */
 195     for ( i = 1, b = 1; i < FT_TRIG_MAX_ITERS; b <<= 1, i++ )
 196     {
 197       if ( theta < 0 )
 198       {
 199         xtemp  = x + ( ( y + b ) >> i );
 200         y      = y - ( ( x + b ) >> i );
 201         x      = xtemp;
 202         theta += *arctanptr++;
 203       }
 204       else
 205       {
 206         xtemp  = x - ( ( y + b ) >> i );
 207         y      = y + ( ( x + b ) >> i );
 208         x      = xtemp;
 209         theta -= *arctanptr++;
 210       }
 211     }
 212
 213     vec->x = x;
 214     vec->y = y;
 215   }
 216
 217
 218   static void
 219   ft_trig_pseudo_polarize( FT_Vector*  vec )
 220   {
 221     FT_Angle         theta;
 222     FT_Int           i;
 223     FT_Fixed         x, y, xtemp, b;
 224     const FT_Angle  *arctanptr;
 225
 226
 227     x = vec->x;
 228     y = vec->y;
 229
 230     /* Get the vector into [-PI/4,PI/4] sector */
 231     if ( y > x )
 232     {
 233       if ( y > -x )
 234       {
 235         theta =  FT_ANGLE_PI2;
 236         xtemp =  y;
 237         y     = -x;
 238         x     =  xtemp;
 239       }
 240       else
 241       {
 242         theta =  y > 0 ? FT_ANGLE_PI : -FT_ANGLE_PI;
 243         x     = -x;
 244         y     = -y;
 245       }
 246     }
 247     else
 248     {
 249       if ( y < -x )
 250       {
 251         theta = -FT_ANGLE_PI2;
 252         xtemp = -y;
 253         y     =  x;
 254         x     =  xtemp;
 255       }
 256       else
 257       {
 258         theta = 0;
 259       }
 260     }
 261
 262     arctanptr = ft_trig_arctan_table;
 263
 264     /* Pseudorotations, with right shifts */
 265     for ( i = 1, b = 1; i < FT_TRIG_MAX_ITERS; b <<= 1, i++ )
 266     {
 267       if ( y > 0 )
 268       {
 269         xtemp  = x + ( ( y + b ) >> i );
 270         y      = y - ( ( x + b ) >> i );
 271         x      = xtemp;
 272         theta += *arctanptr++;
 273       }
 274       else
 275       {
 276         xtemp  = x - ( ( y + b ) >> i );
 277         y      = y + ( ( x + b ) >> i );
 278         x      = xtemp;
 279         theta -= *arctanptr++;
 280       }
 281     }
 282
 283     /* round theta to acknowledge its error that mostly comes */
 284     /* from accumulated rounding errors in the arctan table   */
 285     if ( theta >= 0 )
 286       theta = FT_PAD_ROUND( theta, 16 );
 287     else
 288       theta = -FT_PAD_ROUND( -theta, 16 );
 289
 290     vec->x = x;
 291     vec->y = theta;
 292   }
 293
 294
 295   /* documentation is in fttrigon.h */
 296
 297   FT_EXPORT_DEF( FT_Fixed )
 298   FT_Cos( FT_Angle  angle )
 299   {
 300     FT_Vector  v;
 301
 302
 303     FT_Vector_Unit( &v, angle );
 304
 305     return v.x;
 306   }
 307
 308
 309   /* documentation is in fttrigon.h */
 310
 311   FT_EXPORT_DEF( FT_Fixed )
 312   FT_Sin( FT_Angle  angle )
 313   {
 314     FT_Vector  v;
 315
 316
 317     FT_Vector_Unit( &v, angle );
 318
 319     return v.y;
 320   }
 321
 322
 323   /* documentation is in fttrigon.h */
 324
 325   FT_EXPORT_DEF( FT_Fixed )
 326   FT_Tan( FT_Angle  angle )
 327   {
 328     FT_Vector  v;
 329
 330
 331     FT_Vector_Unit( &v, angle );
 332
 333     return FT_DivFix( v.y, v.x );
 334   }
 335
 336
 337   /* documentation is in fttrigon.h */
 338
 339   FT_EXPORT_DEF( FT_Angle )
 340   FT_Atan2( FT_Fixed  dx,
 341             FT_Fixed  dy )
 342   {
 343     FT_Vector  v;
 344
 345
 346     if ( dx == 0 && dy == 0 )
 347       return 0;
 348
 349     v.x = dx;
 350     v.y = dy;
 351     ft_trig_prenorm( &v );
 352     ft_trig_pseudo_polarize( &v );
 353
 354     return v.y;
 355   }
 356
 357
 358   /* documentation is in fttrigon.h */
 359
 360   FT_EXPORT_DEF( void )
 361   FT_Vector_Unit( FT_Vector*  vec,
 362                   FT_Angle    angle )
 363   {
 364     if ( !vec )
 365       return;
 366
 367     vec->x = FT_TRIG_SCALE >> 8;
 368     vec->y = 0;
 369     ft_trig_pseudo_rotate( vec, angle );
 370     vec->x = ( vec->x + 0x80L ) >> 8;
 371     vec->y = ( vec->y + 0x80L ) >> 8;
 372   }
 373
 374
 375   /* these macros return 0 for positive numbers,
 376      and -1 for negative ones */
 377 #define FT_SIGN_LONG( x )   ( (x) >> ( FT_SIZEOF_LONG * 8 - 1 ) )
 378 #define FT_SIGN_INT( x )    ( (x) >> ( FT_SIZEOF_INT * 8 - 1 ) )
 379 #define FT_SIGN_INT32( x )  ( (x) >> 31 )
 380 #define FT_SIGN_INT16( x )  ( (x) >> 15 )
 381
 382
 383   /* documentation is in fttrigon.h */
 384
 385   FT_EXPORT_DEF( void )
 386   FT_Vector_Rotate( FT_Vector*  vec,
 387                     FT_Angle    angle )
 388   {
 389     FT_Int     shift;
 390     FT_Vector  v;
 391
 392
 393     if ( !vec || !angle )
 394       return;
 395
 396     v = *vec;
 397
 398     if ( v.x == 0 && v.y == 0 )
 399       return;
 400
 401     shift = ft_trig_prenorm( &v );
 402     ft_trig_pseudo_rotate( &v, angle );
 403     v.x = ft_trig_downscale( v.x );
 404     v.y = ft_trig_downscale( v.y );
 405
 406     if ( shift > 0 )
 407     {
 408       FT_Int32  half = (FT_Int32)1L << ( shift - 1 );
 409
 410
 411       vec->x = ( v.x + half + FT_SIGN_LONG( v.x ) ) >> shift;
 412       vec->y = ( v.y + half + FT_SIGN_LONG( v.y ) ) >> shift;
 413     }
 414     else
 415     {
 416       shift  = -shift;
 417       vec->x = (FT_Pos)( (FT_ULong)v.x << shift );
 418       vec->y = (FT_Pos)( (FT_ULong)v.y << shift );
 419     }
 420   }
 421
 422
 423   /* documentation is in fttrigon.h */
 424
 425   FT_EXPORT_DEF( FT_Fixed )
 426   FT_Vector_Length( FT_Vector*  vec )
 427   {
 428     FT_Int     shift;
 429     FT_Vector  v;
 430
 431
 432     if ( !vec )
 433       return 0;
 434
 435     v = *vec;
 436
 437     /* handle trivial cases */
 438     if ( v.x == 0 )
 439     {
 440       return FT_ABS( v.y );
 441     }
 442     else if ( v.y == 0 )
 443     {
 444       return FT_ABS( v.x );
 445     }
 446
 447     /* general case */
 448     shift = ft_trig_prenorm( &v );
 449     ft_trig_pseudo_polarize( &v );
 450
 451     v.x = ft_trig_downscale( v.x );
 452
 453     if ( shift > 0 )
 454       return ( v.x + ( 1L << ( shift - 1 ) ) ) >> shift;
 455
 456     return (FT_Fixed)( (FT_UInt32)v.x << -shift );
 457   }
 458
 459
 460   /* documentation is in fttrigon.h */
 461
 462   FT_EXPORT_DEF( void )
 463   FT_Vector_Polarize( FT_Vector*  vec,
 464                       FT_Fixed   *length,
 465                       FT_Angle   *angle )
 466   {
 467     FT_Int     shift;
 468     FT_Vector  v;
 469
 470
 471     if ( !vec || !length || !angle )
 472       return;
 473
 474     v = *vec;
 475
 476     if ( v.x == 0 && v.y == 0 )
 477       return;
 478
 479     shift = ft_trig_prenorm( &v );
 480     ft_trig_pseudo_polarize( &v );
 481
 482     v.x = ft_trig_downscale( v.x );
 483
 484     *length = shift >= 0 ?                      ( v.x >>  shift )
 485                          : (FT_Fixed)( (FT_UInt32)v.x << -shift );
 486     *angle  = v.y;
 487   }
 488
 489
 490   /* documentation is in fttrigon.h */
 491
 492   FT_EXPORT_DEF( void )
 493   FT_Vector_From_Polar( FT_Vector*  vec,
 494                         FT_Fixed    length,
 495                         FT_Angle    angle )
 496   {
 497     if ( !vec )
 498       return;
 499
 500     vec->x = length;
 501     vec->y = 0;
 502
 503     FT_Vector_Rotate( vec, angle );
 504   }
 505
 506
 507   /* documentation is in fttrigon.h */
 508
 509   FT_EXPORT_DEF( FT_Angle )
 510   FT_Angle_Diff( FT_Angle  angle1,
 511                  FT_Angle  angle2 )
 512   {
 513     FT_Angle  delta = angle2 - angle1;
 514
 515
 516     while ( delta <= -FT_ANGLE_PI )
 517       delta += FT_ANGLE_2PI;
 518
 519     while ( delta > FT_ANGLE_PI )
 520       delta -= FT_ANGLE_2PI;
 521
 522     return delta;
 523   }
 524
 525
 526 /* END */