home *** CD-ROM | disk | FTP | other *** search
/ Fresh Fish 4 / FreshFish_May-June1994.bin / bsd / src / libm / libm-amiga / ieee / support.c < prev   
C/C++ Source or Header  |  1993-09-23  |  18KB  |  525 lines
  1. /*
  2.  * Copyright (c) 1985 Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms, with or without
  6.  * modification, are permitted provided that the following conditions
  7.  * are met:
  8.  * 1. Redistributions of source code must retain the above copyright
  9.  *    notice, this list of conditions and the following disclaimer.
  10.  * 2. Redistributions in binary form must reproduce the above copyright
  11.  *    notice, this list of conditions and the following disclaimer in the
  12.  *    documentation and/or other materials provided with the distribution.
  13.  * 3. All advertising materials mentioning features or use of this software
  14.  *    must display the following acknowledgement:
  15.  *    This product includes software developed by the University of
  16.  *    California, Berkeley and its contributors.
  17.  * 4. Neither the name of the University nor the names of its contributors
  18.  *    may be used to endorse or promote products derived from this software
  19.  *    without specific prior written permission.
  20.  *
  21.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  22.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  23.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  24.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  25.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  26.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  27.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  28.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  29.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  30.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  31.  * SUCH DAMAGE.
  32.  */
  33.  
  34. #ifndef lint
  35. static char sccsid[] = "@(#)support.c    5.6 (Berkeley) 10/9/90";
  36. #endif /* not lint */
  37.  
  38. /* 
  39.  * Some IEEE standard 754 recommended functions and remainder and sqrt for 
  40.  * supporting the C elementary functions.
  41.  ******************************************************************************
  42.  * WARNING:
  43.  *      These codes are developed (in double) to support the C elementary
  44.  * functions temporarily. They are not universal, and some of them are very
  45.  * slow (in particular, drem and sqrt is extremely inefficient). Each 
  46.  * computer system should have its implementation of these functions using 
  47.  * its own assembler.
  48.  ******************************************************************************
  49.  *
  50.  * IEEE 754 required operations:
  51.  *     drem(x,p) 
  52.  *              returns  x REM y  =  x - [x/y]*y , where [x/y] is the integer
  53.  *              nearest x/y; in half way case, choose the even one.
  54.  *     sqrt(x) 
  55.  *              returns the square root of x correctly rounded according to 
  56.  *        the rounding mod.
  57.  *
  58.  * IEEE 754 recommended functions:
  59.  * (a) copysign(x,y) 
  60.  *              returns x with the sign of y. 
  61.  * (b) scalb(x,N) 
  62.  *              returns  x * (2**N), for integer values N.
  63.  * (c) logb(x) 
  64.  *              returns the unbiased exponent of x, a signed integer in 
  65.  *              double precision, except that logb(0) is -INF, logb(INF) 
  66.  *              is +INF, and logb(NAN) is that NAN.
  67.  * (d) finite(x) 
  68.  *              returns the value TRUE if -INF < x < +INF and returns 
  69.  *              FALSE otherwise.
  70.  *
  71.  *
  72.  * CODED IN C BY K.C. NG, 11/25/84;
  73.  * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
  74.  */
  75.  
  76. #include "mathimpl.h"
  77.  
  78. #if defined(vax)||defined(tahoe)      /* VAX D format */
  79. #include <errno.h>
  80.     static const unsigned short msign=0x7fff , mexp =0x7f80 ;
  81.     static const short  prep1=57, gap=7, bias=129           ;   
  82.     static const double novf=1.7E38, nunf=3.0E-39, zero=0.0 ;
  83. #else    /* defined(vax)||defined(tahoe) */
  84.     static const unsigned short msign=0x7fff, mexp =0x7ff0  ;
  85.     static const short prep1=54, gap=4, bias=1023           ;
  86.     static const double novf=1.7E308, nunf=3.0E-308,zero=0.0;
  87. #endif    /* defined(vax)||defined(tahoe) */
  88.  
  89. double scalb(x,N)
  90. double x; int N;
  91. {
  92.         int k;
  93.  
  94. #ifdef national
  95.         unsigned short *px=(unsigned short *) &x + 3;
  96. #else    /* national */
  97.         unsigned short *px=(unsigned short *) &x;
  98. #endif    /* national */
  99.  
  100.         if( x == zero )  return(x); 
  101.  
  102. #if defined(vax)||defined(tahoe)
  103.         if( (k= *px & mexp ) != ~msign ) {
  104.             if (N < -260)
  105.         return(nunf*nunf);
  106.         else if (N > 260) {
  107.         return(copysign(infnan(ERANGE),x));
  108.         }
  109. #else    /* defined(vax)||defined(tahoe) */
  110.         if( (k= *px & mexp ) != mexp ) {
  111.             if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf);
  112.             if( k == 0 ) {
  113.                  x *= scalb(1.0,(int)prep1);  N -= prep1; return(scalb(x,N));}
  114. #endif    /* defined(vax)||defined(tahoe) */
  115.  
  116.             if((k = (k>>gap)+ N) > 0 )
  117.                 if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap);
  118.                 else x=novf+novf;               /* overflow */
  119.             else
  120.                 if( k > -prep1 ) 
  121.                                         /* gradual underflow */
  122.                     {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);}
  123.                 else
  124.                 return(nunf*nunf);
  125.             }
  126.         return(x);
  127. }
  128.  
  129.  
  130. double copysign(x,y)
  131. double x,y;
  132. {
  133. #ifdef national
  134.         unsigned short  *px=(unsigned short *) &x+3,
  135.                         *py=(unsigned short *) &y+3;
  136. #else    /* national */
  137.         unsigned short  *px=(unsigned short *) &x,
  138.                         *py=(unsigned short *) &y;
  139. #endif    /* national */
  140.  
  141. #if defined(vax)||defined(tahoe)
  142.         if ( (*px & mexp) == 0 ) return(x);
  143. #endif    /* defined(vax)||defined(tahoe) */
  144.  
  145.         *px = ( *px & msign ) | ( *py & ~msign );
  146.         return(x);
  147. }
  148.  
  149. double logb(x)
  150. double x; 
  151. {
  152.  
  153. #ifdef national
  154.         short *px=(short *) &x+3, k;
  155. #else    /* national */
  156.         short *px=(short *) &x, k;
  157. #endif    /* national */
  158.  
  159. #if defined(vax)||defined(tahoe)
  160.         return (int)(((*px&mexp)>>gap)-bias);
  161. #else    /* defined(vax)||defined(tahoe) */
  162.         if( (k= *px & mexp ) != mexp )
  163.             if ( k != 0 )
  164.                 return ( (k>>gap) - bias );
  165.             else if( x != zero)
  166.                 return ( -1022.0 );
  167.             else        
  168.                 return(-(1.0/zero));    
  169.         else if(x != x)
  170.             return(x);
  171.         else
  172.             {*px &= msign; return(x);}
  173. #endif    /* defined(vax)||defined(tahoe) */
  174. }
  175.  
  176. finite(x)
  177. double x;    
  178. {
  179. #if defined(vax)||defined(tahoe)
  180.         return(1);
  181. #else    /* defined(vax)||defined(tahoe) */
  182. #ifdef national
  183.         return( (*((short *) &x+3 ) & mexp ) != mexp );
  184. #else    /* national */
  185.         return( (*((short *) &x ) & mexp ) != mexp );
  186. #endif    /* national */
  187. #endif    /* defined(vax)||defined(tahoe) */
  188. }
  189.  
  190. double drem(x,p)
  191. double x,p;
  192. {
  193.         short sign;
  194.         double hp,dp,tmp;
  195.         unsigned short  k; 
  196. #ifdef national
  197.         unsigned short
  198.               *px=(unsigned short *) &x  +3, 
  199.               *pp=(unsigned short *) &p  +3,
  200.               *pd=(unsigned short *) &dp +3,
  201.               *pt=(unsigned short *) &tmp+3;
  202. #else    /* national */
  203.         unsigned short
  204.               *px=(unsigned short *) &x  , 
  205.               *pp=(unsigned short *) &p  ,
  206.               *pd=(unsigned short *) &dp ,
  207.               *pt=(unsigned short *) &tmp;
  208. #endif    /* national */
  209.  
  210.         *pp &= msign ;
  211.  
  212. #if defined(vax)||defined(tahoe)
  213.         if( ( *px & mexp ) == ~msign )    /* is x a reserved operand? */
  214. #else    /* defined(vax)||defined(tahoe) */
  215.         if( ( *px & mexp ) == mexp )
  216. #endif    /* defined(vax)||defined(tahoe) */
  217.         return  (x-p)-(x-p);    /* create nan if x is inf */
  218.     if (p == zero) {
  219. #if defined(vax)||defined(tahoe)
  220.         return(infnan(EDOM));
  221. #else    /* defined(vax)||defined(tahoe) */
  222.         return zero/zero;
  223. #endif    /* defined(vax)||defined(tahoe) */
  224.     }
  225.  
  226. #if defined(vax)||defined(tahoe)
  227.         if( ( *pp & mexp ) == ~msign )    /* is p a reserved operand? */
  228. #else    /* defined(vax)||defined(tahoe) */
  229.         if( ( *pp & mexp ) == mexp )
  230. #endif    /* defined(vax)||defined(tahoe) */
  231.         { if (p != p) return p; else return x;}
  232.  
  233.         else  if ( ((*pp & mexp)>>gap) <= 1 ) 
  234.                 /* subnormal p, or almost subnormal p */
  235.             { double b; b=scalb(1.0,(int)