home *** CD-ROM | disk | FTP | other *** search
/ Fresh Fish 4 / FreshFish_May-June1994.bin / bbs / gnu / libm-src.lha / src / amiga / libm / mc68881 / atan2.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-09-23  |  5.3 KB  |  145 lines
  1. /*-
  2.  * Copyright (c) 1990 The Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * This code is derived from software contributed to Berkeley by
  6.  * the Systems Programming Group of the University of Utah Computer
  7.  * Science Department.
  8.  *
  9.  * Redistribution and use in source and binary forms, with or without
  10.  * modification, are permitted provided that the following conditions
  11.  * are met:
  12.  * 1. Redistributions of source code must retain the above copyright
  13.  *    notice, this list of conditions and the following disclaimer.
  14.  * 2. Redistributions in binary form must reproduce the above copyright
  15.  *    notice, this list of conditions and the following disclaimer in the
  16.  *    documentation and/or other materials provided with the distribution.
  17.  * 3. All advertising materials mentioning features or use of this software
  18.  *    must display the following acknowledgement:
  19.  *    This product includes software developed by the University of
  20.  *    California, Berkeley and its contributors.
  21.  * 4. Neither the name of the University nor the names of its contributors
  22.  *    may be used to endorse or promote products derived from this software
  23.  *    without specific prior written permission.
  24.  *
  25.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  26.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  27.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  28.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  29.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  30.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  31.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  32.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  33.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  34.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  35.  * SUCH DAMAGE.
  36.  */
  37.  
  38. #ifndef lint
  39. static char sccsid[] = "@(#)atan2.c    5.1 (Berkeley) 5/17/90";
  40. #endif /* not lint */
  41.  
  42. /*
  43.  * ATAN2(Y,X)
  44.  * RETURN ARG (X+iY)
  45.  * DOUBLE PRECISION (IEEE DOUBLE 53 BITS)
  46.  *
  47.  * Scaled down version to weed out special cases.  "Normal" cases are
  48.  * handled by calling atan2__A(), an assembly coded support routine in
  49.  * support.s.
  50.  *
  51.  * Required system supported functions :
  52.  *    copysign(x,y)
  53.  *    atan2__A(y,x)
  54.  *    
  55.  * Method :
  56.  *    1. Deal with special cases
  57.  *    2. Call atan2__A() to do the others
  58.  *
  59.  * Special cases:
  60.  * Notations: atan2(y,x) == ARG (x+iy) == ARG(x,y).
  61.  *
  62.  *    ARG( NAN , (anything) ) is NaN;
  63.  *    ARG( (anything), NaN ) is NaN;
  64.  *    ARG(+(anything but NaN), +-0) is +-0  ;
  65.  *    ARG(-(anything but NaN), +-0) is +-PI ;
  66.  *    ARG( 0, +-(anything but 0 and NaN) ) is +-PI/2;
  67.  *    ARG( +INF,+-(anything but INF and NaN) ) is +-0 ;
  68.  *    ARG( -INF,+-(anything but INF and NaN) ) is +-PI;
  69.  *    ARG( +INF,+-INF ) is +-PI/4 ;
  70.  *    ARG( -INF,+-INF ) is +-3PI/4;
  71.  *    ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2;
  72.  *
  73.  * Accuracy:
  74.  *    atan2(y,x) returns (PI/pi) * the exact ARG (x+iy) nearly rounded, 
  75.  *    where
  76.  *
  77.  *    in decimal:
  78.  *        pi = 3.141592653589793 23846264338327 ..... 
  79.  *    53 bits   PI = 3.141592653589793 115997963 ..... ,
  80.  *    56 bits   PI = 3.141592653589793 227020265 ..... ,  
  81.  *
  82.  *    in hexadecimal:
  83.  *        pi = 3.243F6A8885A308D313198A2E....
  84.  *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18    error=.276ulps
  85.  *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    error=.206ulps
  86.  *    
  87.  *    In a test run with 356,000 random argument on [-1,1] * [-1,1] on a
  88.  *    VAX, the maximum observed error was 1.41 ulps (units of the last place)
  89.  *    compared with (PI/pi)*(the exact ARG(x+iy)).
  90.  *
  91.  * Note:
  92.  *    We use machine PI (the true pi rounded) in place of the actual
  93.  *    value of pi for all the trig and inverse trig functions. In general, 
  94.  *    if trig is one of sin, cos, tan, then computed trig(y) returns the 
  95.  *    exact trig(y*pi/PI) nearly rounded; correspondingly, computed arctrig 
  96.  *    returns the exact arctrig(y)*PI/pi nearly rounded. These guarantee the 
  97.  *    trig functions have period PI, and trig(arctrig(x)) returns x for
  98.  *    all critical values x.
  99.  *    
  100.  * Constants:
  101.  * The hexadecimal values are the intended ones for the following constants.
  102.  * The decimal values may be used, provided that the compiler will convert
  103.  * from decimal to binary accurately enough to produce the hexadecimal values
  104.  * shown.
  105.  */
  106.  
  107. static double 
  108. PIo4   =  7.8539816339744827900E-1    , /*Hex  2^ -1   *  1.921FB54442D18 */
  109. PIo2   =  1.5707963267948965580E0     , /*Hex  2^  0   *  1.921FB54442D18 */
  110. PI     =  3.1415926535897931160E0     ; /*Hex  2^  1   *  1.921FB54442D18 */
  111.  
  112. double atan2(y,x)
  113. double  y,x;
  114. {  
  115.     static double zero=0, one=1;
  116.     double copysign(),atan2__A(),signy,signx;
  117.     int finite();
  118.  
  119.     /* if x or y is NAN */
  120.     if(x!=x) return(x); if(y!=y) return(y);
  121.  
  122.     /* copy down the sign of y and x */
  123.     signy = copysign(one,y);
  124.     signx = copysign(one,x);
  125.  
  126.     /* when y = 0 */
  127.     if(y==zero) return((signx==one)?y:copysign(PI,signy));
  128.  
  129.     /* when x = 0 */
  130.     if(x==zero) return(copysign(PIo2,signy));
  131.         
  132.     /* when x is INF */
  133.     if(!finite(x))
  134.         if(!finite(y)) 
  135.         return(copysign((signx==one)?PIo4:3*PIo4,signy));
  136.         else
  137.         return(copysign((signx==one)?zero:PI,signy));
  138.  
  139.     /* when y is INF */
  140.     if(!finite(y)) return(copysign(PIo2,signy));
  141.  
  142.     /* else let atan2__A do the work */
  143.     return(atan2__A(y,x));
  144. }
  145.