BSD 4_3_Reno release
[unix-history] / usr / src / lib / libm / common_source / pow.c
index e4b073a..d89348b 100644 (file)
@@ -1,20 +1,30 @@
-/* 
+/*
  * Copyright (c) 1985 Regents of the University of California.
  * Copyright (c) 1985 Regents of the University of California.
- * 
- * Use and reproduction of this software are granted  in  accordance  with
- * the terms and conditions specified in  the  Berkeley  Software  License
- * Agreement (in particular, this entails acknowledgement of the programs'
- * source, and inclusion of this notice) with the additional understanding
- * that  all  recipients  should regard themselves as participants  in  an
- * ongoing  research  project and hence should  feel  obligated  to report
- * their  experiences (good or bad) with these elementary function  codes,
- * using "sendbug 4bsd-bugs@BERKELEY", to the authors.
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms are permitted
+ * provided that: (1) source distributions retain this entire copyright
+ * notice and comment, and (2) distributions including binaries display
+ * the following acknowledgement:  ``This product includes software
+ * developed by the University of California, Berkeley and its contributors''
+ * in the documentation or other materials provided with the distribution
+ * and in all advertising materials mentioning features or use of this
+ * software. Neither the name of the University nor the names of its
+ * contributors may be used to endorse or promote products derived
+ * from this software without specific prior written permission.
+ * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ * All recipients should regard themselves as participants in an ongoing
+ * research project and hence should feel obligated to report their
+ * experiences (good or bad) with these elementary function codes, using
+ * the sendbug(8) program, to the authors.
  */
 
 #ifndef lint
  */
 
 #ifndef lint
-static char sccsid[] =
-"@(#)pow.c     4.5 (Berkeley) 8/21/85; 1.3 (ucb.elefunt) %G%";
-#endif not lint
+static char sccsid[] = "@(#)pow.c      5.6 (Berkeley) 6/1/90";
+#endif /* not lint */
 
 /* POW(X,Y)  
  * RETURN X**Y 
 
 /* POW(X,Y)  
  * RETURN X**Y 
@@ -84,48 +94,44 @@ static char sccsid[] =
  * shown.
  */
 
  * shown.
  */
 
-#ifdef VAX     /* VAX D format */
 #include <errno.h>
 #include <errno.h>
-extern double infnan();
-
-/* static double */
-/* ln2hi  =  6.9314718055829871446E-1    , Hex  2^  0   *  .B17217F7D00000 */
-/* ln2lo  =  1.6465949582897081279E-12   , Hex  2^-39   *  .E7BCD5E4F1D9CC */
-/* invln2 =  1.4426950408889634148E0     , Hex  2^  1   *  .B8AA3B295C17F1 */
-/* sqrt2  =  1.4142135623730950622E0     ; Hex  2^  1   *  .B504F333F9DE65 */
-static long     ln2hix[] = { 0x72174031, 0x0000f7d0};
-static long     ln2lox[] = { 0xbcd52ce7, 0xd9cce4f1};
-static long    invln2x[] = { 0xaa3b40b8, 0x17f1295c};
-static long     sqrt2x[] = { 0x04f340b5, 0xde6533f9};
-#define    ln2hi    (*(double*)ln2hix)
-#define    ln2lo    (*(double*)ln2lox)
-#define   invln2    (*(double*)invln2x)
-#define    sqrt2    (*(double*)sqrt2x)
-#else  /* IEEE double */
-static double
-ln2hi  =  6.9314718036912381649E-1    , /*Hex  2^ -1   *  1.62E42FEE00000 */
-ln2lo  =  1.9082149292705877000E-10   , /*Hex  2^-33   *  1.A39EF35793C76 */
-invln2 =  1.4426950408889633870E0     , /*Hex  2^  0   *  1.71547652B82FE */
-sqrt2  =  1.4142135623730951455E0     ; /*Hex  2^  0   *  1.6A09E667F3BCD */
+#include "mathimpl.h"
+
+vc(ln2hi,  6.9314718055829871446E-1  ,7217,4031,0000,f7d0,   0, .B17217F7D00000)
+vc(ln2lo,  1.6465949582897081279E-12 ,bcd5,2ce7,d9cc,e4f1, -39, .E7BCD5E4F1D9CC)
+vc(invln2, 1.4426950408889634148E0   ,aa3b,40b8,17f1,295c,   1, .B8AA3B295C17F1)
+vc(sqrt2,  1.4142135623730950622E0   ,04f3,40b5,de65,33f9,   1, .B504F333F9DE65)
+
+ic(ln2hi,  6.9314718036912381649E-1,   -1, 1.62E42FEE00000)
+ic(ln2lo,  1.9082149292705877000E-10, -33, 1.A39EF35793C76)
+ic(invln2, 1.4426950408889633870E0,     0, 1.71547652B82FE)
+ic(sqrt2,  1.4142135623730951455E0,     0, 1.6A09E667F3BCD)
+
+#ifdef vccast
+#define        ln2hi   vccast(ln2hi)
+#define        ln2lo   vccast(ln2lo)
+#define        invln2  vccast(invln2)
+#define        sqrt2   vccast(sqrt2)
 #endif
 
 #endif
 
-static double zero=0.0, half=1.0/2.0, one=1.0, two=2.0, negone= -1.0;
+const static double zero=0.0, half=1.0/2.0, one=1.0, two=2.0, negone= -1.0;
+
+static double pow_p();
 
 double pow(x,y)        
 double x,y;
 {
 
 double pow(x,y)        
 double x,y;
 {
-       double drem(),pow_p(),copysign(),t;
-       int finite();
+       double t;
 
        if     (y==zero)      return(one);
        else if(y==one
 
        if     (y==zero)      return(one);
        else if(y==one
-#ifndef VAX
+#if !defined(vax)&&!defined(tahoe)
                ||x!=x
                ||x!=x
-#endif
+#endif /* !defined(vax)&&!defined(tahoe) */
                ) return( x );      /* if x is NaN or y=1 */
                ) return( x );      /* if x is NaN or y=1 */
-#ifndef VAX
+#if !defined(vax)&&!defined(tahoe)
        else if(y!=y)         return( y );      /* if y is NaN */
        else if(y!=y)         return( y );      /* if y is NaN */
-#endif
+#endif /* !defined(vax)&&!defined(tahoe) */
        else if(!finite(y))                     /* if y is INF */
             if((t=copysign(x,one))==one) return(zero/zero);
             else if(t>one) return((y>zero)?y:zero);
        else if(!finite(y))                     /* if y is INF */
             if((t=copysign(x,one))==one) return(zero/zero);
             else if(t>one) return((y>zero)?y:zero);
@@ -147,38 +153,41 @@ double x,y;
        else if(x==zero)        /* x is -0 */
            return((y>zero)?-x:one/(-x));
        else {                  /* return NaN */
        else if(x==zero)        /* x is -0 */
            return((y>zero)?-x:one/(-x));
        else {                  /* return NaN */
-#ifdef VAX
+#if defined(vax)||defined(tahoe)
            return (infnan(EDOM));      /* NaN */
            return (infnan(EDOM));      /* NaN */
-#else  /* IEEE double */
+#else  /* defined(vax)||defined(tahoe) */
            return(zero/zero);
            return(zero/zero);
-#endif
+#endif /* defined(vax)||defined(tahoe) */
        }
 }
 
        }
 }
 
+#ifndef mc68881
 /* pow_p(x,y) return x**y for x with sign=1 and finite y */
 static double pow_p(x,y)       
 double x,y;
 {
 /* pow_p(x,y) return x**y for x with sign=1 and finite y */
 static double pow_p(x,y)       
 double x,y;
 {
-        double logb(),scalb(),copysign(),log__L(),exp__E();
         double c,s,t,z,tx,ty;
         double c,s,t,z,tx,ty;
+#ifdef tahoe
+       double tahoe_tmp;
+#endif /* tahoe */
         float sx,sy;
        long k=0;
         int n,m;
 
        if(x==zero||!finite(x)) {           /* if x is +INF or +0 */
         float sx,sy;
        long k=0;
         int n,m;
 
        if(x==zero||!finite(x)) {           /* if x is +INF or +0 */
-#ifdef VAX
+#if defined(vax)||defined(tahoe)
             return((y>zero)?x:infnan(ERANGE)); /* if y<zero, return +INF */
             return((y>zero)?x:infnan(ERANGE)); /* if y<zero, return +INF */
-#else
+#else  /* defined(vax)||defined(tahoe) */
             return((y>zero)?x:one/x);
             return((y>zero)?x:one/x);
-#endif
+#endif /* defined(vax)||defined(tahoe) */
        }
        if(x==1.0) return(x);   /* if x=1.0, return 1 since y is finite */
 
     /* reduce x to z in [sqrt(1/2)-1, sqrt(2)-1] */
         z=scalb(x,-(n=logb(x)));  
        }
        if(x==1.0) return(x);   /* if x=1.0, return 1 since y is finite */
 
     /* reduce x to z in [sqrt(1/2)-1, sqrt(2)-1] */
         z=scalb(x,-(n=logb(x)));  
-#ifndef VAX    /* IEEE double */       /* subnormal number */
+#if !defined(vax)&&!defined(tahoe)     /* IEEE double; subnormal number */
         if(n <= -1022) {n += (m=logb(z)); z=scalb(z,-m);} 
         if(n <= -1022) {n += (m=logb(z)); z=scalb(z,-m);} 
-#endif
+#endif /* !defined(vax)&&!defined(tahoe) */
         if(z >= sqrt2 ) {n += 1; z *= half;}  z -= one ;
 
     /* log(x) = nlog2+log(1+z) ~ nlog2 + t + tx */
         if(z >= sqrt2 ) {n += 1; z *= half;}  z -= one ;
 
     /* log(x) = nlog2+log(1+z) ~ nlog2 + t + tx */
@@ -203,7 +212,11 @@ double x,y;
           /* end of checking whether k==y */
 
                 sy=y; ty=y-sy;          /* y ~ sy + ty */
           /* end of checking whether k==y */
 
                 sy=y; ty=y-sy;          /* y ~ sy + ty */
+#ifdef tahoe
+               s = (tahoe_tmp = sx)*sy-k*ln2hi;
+#else  /* tahoe */
                s=(double)sx*sy-k*ln2hi;        /* (sy+ty)*(sx+tx)-kln2 */
                s=(double)sx*sy-k*ln2hi;        /* (sy+ty)*(sx+tx)-kln2 */
+#endif /* tahoe */
                z=(tx*ty-k*ln2lo);
                tx=tx*sy; ty=sx*ty;
                t=ty+z; t+=tx; t+=s;
                z=(tx*ty-k*ln2lo);
                tx=tx*sy; ty=sx*ty;
                t=ty+z; t+=tx; t+=s;
@@ -225,3 +238,4 @@ double x,y;
                        return(scalb(one, 5000)); 
 
 }
                        return(scalb(one, 5000)); 
 
 }
+#endif /* mc68881 */