home *** CD-ROM | disk | FTP | other *** search
/ Amiga ISO Collection / AmigaUtilCD2.iso / Programming / GCC / GERLIB_DEV08B.LHA / gerlib / Bonus / normal / support.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-12  |  17.5 KB  |  535 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. #if defined(LIBC_SCCS) && !defined(lint)
  35. static char sccsid[] = "@(#)support.c    5.6 (Berkeley) 10/9/90";
  36. #endif /* LIBC_SCCS and 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. /* extracted copysign and finite, made some #amiga-specific things : Sat Oct  2 19:58:56 1993 GM */
  77.  
  78. #ifndef amiga
  79. #include "mathimpl.h"
  80. #else
  81.  
  82. #endif
  83.  
  84. #if defined(vax)||defined(tahoe)      /* VAX D format */
  85. #include <errno.h>
  86.     static const unsigned short msign=0x7fff , mexp =0x7f80 ;
  87.     static const short  prep1=57, gap=7, bias=129           ;   
  88.     static const double novf=1.7E38, nunf=3.0E-39, zero=0.0 ;
  89. #else    /* defined(vax)||defined(tahoe) */
  90.     static const unsigned short msign=0x7fff, mexp =0x7ff0  ;
  91.     static const short prep1=54, gap=4, bias=1023           ;
  92.     static const double novf=1.7E308, nunf=3.0E-308,zero=0.0;
  93. #endif    /* defined(vax)||defined(tahoe) */
  94.  
  95. double scalb(x,N)
  96. double x; int N;
  97. {
  98.         int k;
  99.  
  100. #ifdef national
  101.         unsigned short *px=(unsigned short *) &x + 3;
  102. #else    /* national */
  103.         unsigned short *px=(unsigned short *) &x;
  104. #endif    /* national */
  105.  
  106.         if( x == zero )  return(x);
  107.  
  108. #if defined(vax)||defined(tahoe)
  109.         if( (k= *px & mexp ) != ~msign ) {
  110.             if (N < -260)
  111.         return(nunf*nunf);
  112.         else if (N > 260) {
  113.         return(copysign(infnan(ERANGE),x));
  114.         }
  115. #else    /* defined(vax)||defined(tahoe) */
  116.         if( (k= *px & mexp ) != mexp ) {
  117.             if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf);
  118.             if( k == 0 ) {
  119.                  x *= scalb(1.0,(int)prep1);  N -= prep1; return(scalb(x,N));}
  120. #endif    /* defined(vax)||defined(tahoe) */
  121.  
  122.             if((k = (k>>gap)+ N) > 0 )
  123.                 if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap);
  124.                 else x=novf+novf;               /* overflow */
  125.             else
  126.                 if( k > -prep1 )
  127.                                         /* gradual underflow */
  128.                     {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);}
  129.                 else
  130.                 return(nunf*nunf);
  131.             }
  132.         return(x);
  133. }
  134.  
  135. #ifndef amiga
  136. double copysign(x,y)
  137. double x,y;
  138. {
  139. #ifdef national
  140.         unsigned short  *px=(unsigned short *) &x+3,
  141.                         *py=(unsigned short *) &y+3;
  142. #else    /* national */
  143.         unsigned short  *px=(unsigned short *) &x,
  144.                         *py=(unsigned short *) &y;
  145. #endif    /* national */
  146.  
  147. #if defined(vax)||defined(tahoe)
  148.         if ( (*px & mexp) == 0 ) return(x);
  149. #endif    /* defined(vax)||defined(tahoe) */
  150.  
  151.         *px = ( *px & msign ) | ( *py & ~msign );
  152.         return(x);
  153. }
  154. #endif
  155.  
  156. double logb(x)
  157. double x;
  158. {
  159.  
  160. #ifdef national
  161.         short *px=(short *) &x+3, k;
  162. #else    /* national */
  163.         short *px=(short *) &x, k;
  164. #endif    /* national */
  165.  
  166. #if defined(vax)||defined(tahoe)
  167.         return (int)(((*px&mexp)>>gap)-bias);
  168. #else    /* defined(vax)||defined(tahoe) */
  169.         if( (k= *px & mexp ) != mexp )
  170.             if ( k != 0 )
  171.                 return ( (k>>gap) - bias );
  172.             else if( x != zero)
  173.                 return ( -1022.0 );
  174.             else        
  175.                 return(-(1.0/zero));
  176.         else if(x != x)
  177.             return(x);
  178.         else
  179.             {*px &= msign; return(x);}
  180. #endif    /* defined(vax)||defined(tahoe) */
  181. }
  182.  
  183. #ifndef amiga
  184. finite(x)
  185. double x;
  186. {
  187. #if defined(vax)||defined(tahoe)
  188.         return(1);
  189. #else    /* defined(vax)||defined(tahoe) */
  190. #ifdef national
  191.         return( (*((short *) &x+3 ) & mexp ) != mexp );
  192. #else    /* national */
  193.         return( (*((short *) &x ) & mexp ) != mexp );
  194. #endif    /* national */
  195. #endif    /* defined(vax)||defined(tahoe) */
  196. }
  197. #endif
  198.  
  199. double drem(x,p)
  200. double x,p;
  201. {
  202.         short sign;
  203.         double hp,dp,tmp;
  204.         unsigned short  k; 
  205. #ifdef national
  206.         unsigned short
  207.               *px=(unsigned short *) &x  +3, 
  208.               *pp=(unsigned short *) &p  +3,
  209.               *pd=(unsigned short *) &dp +3,
  210.               *pt=(unsigned short *) &tmp+3;
  211. #else    /* national */
  212.         unsigned short
  213.               *px=(unsigned short *) &x  , 
  214.               *pp=(unsigned short *) &p  ,
  215.               *pd=(unsigned short *) &dp ,
  216.               *pt=(unsigned short *) &tmp;
  217. #endif    /* national */
  218.  
  219.         *pp &= msign ;
  220.  
  221. #if defined(vax)||defined(tahoe)
  222.         if( ( *px & mexp ) == ~msign )    /* is x a reserved operand? */
  223. #else    /* defined(vax)||defined(tahoe) */
  224.         if( ( *px & mexp ) == mexp )
  225. #endif    /* defined(vax)||defined(tahoe) */
  226.         return  (x-p)-(x-p);    /* create nan if x is inf */
  227.     if (p == zero) {
  228. #if defined(vax)||defined(tahoe)
  229.         return(infnan(EDOM));
  230. #else    /* defined(vax)||defined(tahoe) */
  231.         return zero/zero;
  232. #endif    /* defined(vax)||defined(tahoe) */
  233.     }
  234.  
  235. #if defined(vax)||defined(tahoe)
  236.         if( ( *pp & mexp ) == ~msign )    /* is p a reserved operand? */
  237. #else    /* defined(vax)||defined(tahoe) */
  238.         if( ( *pp & mexp ) == mexp )
  239. #endif    /* defined(vax)||defined(tahoe) */
  240.         { if (p != p) return p; else return x;}
  241.  
  242.         else  if ( ((*pp & mexp)>>gap) <= 1 ) 
  243.                 /* subnormal p, or almost subnormal p */
  244.             { double b; b=scalb(1.0,(int)prep1);
  245.               p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);}
  246.         else  if ( p >= novf/2)
  247.             { p /= 2 ; x /= 2; return(drem(x,p)*2);}
  248.         else 
  249.             {
  250.                 dp=p+p; hp=p/2;
  251.                 sign= *px & ~msign ;
  252.                 *px &= msign       ;
  253.                 while ( x > dp )
  254.                     {
  255.                         k=(*px & mexp) - (*pd & mexp) ;
  256.                         tmp = dp ;
  257.                         *pt += k ;
  258.  
  259. #if defined(vax)||defined(tahoe)
  260.                         if( x < tmp ) *pt -= 128 ;
  261. #else    /* defined(vax)||defined(tahoe) */
  262.                         if( x < tmp ) *pt -= 16 ;
  263. #endif    /* defined(vax)||defined(tahoe) */
  264.  
  265.                         x -= tmp ;
  266.                     }
  267.                 if ( x > hp )
  268.                     { x -= p ;  if ( x >= hp ) x -= p ; }
  269.  
  270. #if defined(vax)||defined(tahoe)
  271.         if (x)
  272. #endif    /* defined(vax)||defined(tahoe) */
  273.             *px ^= sign;
  274.                 return( x);
  275.  
  276.             }
  277. }
  278.  
  279. #ifndef amiga
  280. double sqrt(x)
  281. double x;
  282. {
  283.         double q,s,b,r;
  284.         double t;
  285.     double const zero=0.0;
  286.         int m,n,i;
  287. #if defined(vax)||defined(tahoe)
  288.         int k=54;
  289. #else    /* defined(vax)||defined(tahoe) */
  290.         int k=51;
  291. #endif    /* defined(vax)||defined(tahoe) */
  292.  
  293.     /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
  294.         if(x!=x||x==zero) return(x);
  295.  
  296.     /* sqrt(negative) is invalid */
  297.         if(x<zero) {
  298. #if defined(vax)||defined(tahoe)
  299.         return (infnan(EDOM));    /* NaN */
  300. #else    /* defined(vax)||defined(tahoe) */
  301.         return(zero/zero);
  302. #endif    /* defined(vax)||defined(tahoe) */
  303.     }
  304.  
  305.     /* sqrt(INF) is INF */
  306.         if(!finite(x)) return(x);
  307.  
  308.     /* scale x to [1,4) */
  309.         n=logb(x);
  310.         x=scalb(x,-n);
  311.         if((m=logb(x))!=0) x=scalb(x,-m);       /* subnormal number */
  312.         m += n; 
  313.         n = m/2;
  314.         if((n+n)!=m) {x *= 2; m -=1; n=m/2;}
  315.  
  316.     /* generate sqrt(x) bit by bit (accumulating in q) */
  317.             q=1.0; s=4.0; x -= 1.0; r=1;
  318.             for(i=1;i<=k;i++) {
  319.                 t=s+1; x *= 4; r /= 2;
  320.                 if(t<=x) {
  321.                     s=t+t+2, x -= t; q += r;}
  322.                 else
  323.                     s *= 2;
  324.                 }
  325.             
  326.     /* generate the last bit and determine the final rounding */
  327.             r/=2; x *= 4; 
  328.             if(x==zero) goto end; 100+r; /* trigger inexact flag */
  329.             if(s<x) {
  330.                 q+=r; x -=s; s += 2; s *= 2; x *= 4;
  331.                 t = (x-s)-5;
  332.                 b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */
  333.                 b=1.0+r/4;   if(b>1.0) t=1;    /* b>1 : Round-to-(+INF) */
  334.                 if(t>=0) q+=r; }          /* else: Round-to-nearest */
  335.             else { 
  336.                 s *= 2; x *= 4;
  337.                 t = (x-s)-1;
  338.                 b=1.0+3*r/4; if(b==1.0) goto end;
  339.                 b=1.0+r/4;   if(b>1.0) t=1;
  340.                 if(t>=0) q+=r; }
  341.  
  342. end:        return(scalb(q,n));
  343. }
  344. #endif
  345.  
  346. #if 0
  347. /* DREM(X,Y)
  348.  * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
  349.  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
  350.  * INTENDED FOR ASSEMBLY LANGUAGE
  351.  * CODED IN C BY K.C. NG, 3/23/85, 4/8/85.
  352.  *
  353.  * Warning: this code should not get compiled in unless ALL of
  354.  * the following machine-dependent routines are supplied.
  355.  * 
  356.  * Required machine dependent functions (not on a VAX):
  357.  *     swapINX(i): save inexact flag and reset it to "i"
  358.  *     swapENI(e): save inexact enable and reset it to "e"
  359.  */
  360.  
  361. double drem(x,y)    
  362. double x,y;
  363. {
  364.  
  365. #ifdef national        /* order of words in floating point number */
  366.     static const n0=3,n1=2,n2=1,n3=0;
  367. #else /* VAX, SUN, ZILOG, TAHOE */
  368.     static const n0=0,n1=1,n2=2,n3=3;
  369. #endif
  370.  
  371.         static const unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390;
  372.     static const double zero=0.0;
  373.     double hy,y1,t,t1;
  374.     short k;
  375.     long n;
  376.     int i,e; 
  377.     unsigned short xexp,yexp, *px  =(unsigned short *) &x  , 
  378.                   nx,nf,      *py  =(unsigned short *) &y  ,
  379.                   sign,      *pt  =(unsigned short *) &t  ,
  380.                         *pt1 =(unsigned short *) &t1 ;
  381.  
  382.     xexp = px[n0] & mexp ;    /* exponent of x */
  383.     yexp = py[n0] & mexp ;    /* exponent of y */
  384.     sign = px[n0] &0x8000;    /* sign of x     */
  385.  
  386. /* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */
  387.     if(x!=x) return(x); if(y!=y) return(y);         /* x or y is NaN */
  388.     if( xexp == mexp )   return(zero/zero);      /* x is INF */
  389.     if(y==zero) return(y/y);
  390.  
  391. /* save the inexact flag and inexact enable in i and e respectively
  392.  * and reset them to zero
  393.  */
  394.     i=swapINX(0);    e=swapENI(0);    
  395.  
  396. /* subnormal number */
  397.     nx=0;
  398.     if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;}
  399.  
  400. /* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */
  401.     if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;}
  402.  
  403.     nf=nx;
  404.     py[n0] &= 0x7fff;    
  405.     px[n0] &= 0x7fff;
  406.  
  407. /* mask off the least significant 27 bits of y */
  408.     t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t;
  409.  
  410. /* LOOP: argument reduction on x whenever x > y */
  411. loop:
  412.     while ( x > y )
  413.     {
  414.         t=y;
  415.         t1=y1;
  416.         xexp=px[n0]&mexp;      /* exponent of x */
  417.         k=xexp-yexp-m25;
  418.         if(k>0)     /* if x/y >= 2**26, scale up y so that x/y < 2**26 */
  419.         {pt[n0]+=k;pt1[n0]+=k;}
  420.         n=x/t; x=(x-n*t1)-n*(t-t1);
  421.     }    
  422.     /* end while (x > y) */
  423.  
  424.     if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;}
  425.  
  426. /* final adjustment */
  427.  
  428.     hy=y/2.0;
  429.     if(x>hy||((x==hy)&&n%2==1)) x-=y; 
  430.     px[n0] ^= sign;
  431.     if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;}
  432.  
  433. /* restore inexact flag and inexact enable */
  434.     swapINX(i); swapENI(e);    
  435.  
  436.     return(x);    
  437. }
  438. #endif
  439.  
  440. #if 0
  441. /* SQRT
  442.  * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT
  443.  * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE
  444.  * CODED IN C BY K.C. NG, 3/22/85.
  445.  *
  446.  * Warning: this code should not get compiled in unless ALL of
  447.  * the following machine-dependent routines are supplied.
  448.  * 
  449.  * Required machine dependent functions:
  450.  *     swapINX(i)  ...return the status of INEXACT flag and reset it to "i"
  451.  *     swapRM(r)   ...return the current Rounding Mode and reset it to "r"
  452.  *     swapENI(e)  ...return the status of inexact enable and reset it to "e"
  453.  *     addc(t)     ...perform t=t+1 regarding t as a 64 bit unsigned integer
  454.  *     subc(t)     ...perform t=t-1 regarding t as a 64 bit unsigned integer
  455.  */
  456.  
  457. static const unsigned long table[] = {
  458. 0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740,
  459. 58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478,
  460. 21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, };
  461.  
  462. double newsqrt(x)
  463. double x;
  464. {
  465.         double y,z,t,addc(),subc()
  466.     double const b54=134217728.*134217728.; /* b54=2**54 */
  467.         long mx,scalx;
  468.     long const mexp=0x7ff00000;
  469.         int i,j,r,e,swapINX(),swapRM(),swapENI();       
  470.         unsigned long *py=(unsigned long *) &y   ,
  471.                       *pt=(unsigned long *) &t   ,
  472.                       *px=(unsigned long *) &x   ;
  473. #ifdef national         /* ordering of word in a floating point number */
  474.         const int n0=1, n1=0; 
  475. #else
  476.         const int n0=0, n1=1; 
  477. #endif
  478. /* Rounding Mode:  RN ...round-to-nearest 
  479.  *                 RZ ...round-towards 0
  480.  *                 RP ...round-towards +INF
  481.  *           RM ...round-towards -INF
  482.  */
  483.         const int RN=0,RZ=1,RP=2,RM=3;
  484.                 /* machine dependent: work on a Zilog Z8070
  485.                                  * and a National 32081 & 16081
  486.                                  */
  487.  
  488. /* exceptions */
  489.     if(x!=x||x==0.0) return(x);  /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
  490.     if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */
  491.         if((mx=px[n0]&mexp)==mexp) return(x);  /* sqrt(+INF) is +INF */
  492.  
  493. /* save, reset, initialize */
  494.         e=swapENI(0);   /* ...save and reset the inexact enable */
  495.         i=swapINX(0);   /* ...save INEXACT flag */
  496.         r=swapRM(RN);   /* ...save and reset the Rounding Mode to RN */
  497.         scalx=0;
  498.  
  499. /* subnormal number, scale up x to x*2**54 */
  500.         if(mx==0) {x *= b54 ; scalx-=0x01b00000;}
  501.  
  502. /* scale x to avoid intermediate over/underflow:
  503.  * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */
  504.         if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;}
  505.         if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;}
  506.  
  507. /* magic initial approximation to almost 8 sig. bits */
  508.         py[n0]=(px[n0]>>1)+0x1ff80000;
  509.         py[n0]=py[n0]-table[(py[n0]>>15)&31];
  510.  
  511. /* Heron's rule once with correction to improve y to almost 18 sig. bits */
  512.         t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0;
  513.  
  514. /* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */
  515.         t=y*y; z=t;  pt[n0]+=0x00100000; t+=z; z=(x-z)*y; 
  516.         t=z/(t+x) ;  pt[n0]+=0x00100000; y+=t;
  517.  
  518. /* twiddle last bit to force y correctly rounded */ 
  519.         swapRM(RZ);     /* ...set Rounding Mode to round-toward-zero */
  520.         swapINX(0);     /* ...clear INEXACT flag */
  521.         swapENI(e);     /* ...restore inexact enable status */
  522.         t=x/y;          /* ...chopped quotient, possibly inexact */
  523.         j=swapINX(i);   /* ...read and restore inexact flag */
  524.         if(j==0) { if(t==y) goto end; else t=subc(t); }  /* ...t=t-ulp */
  525.         b54+0.1;        /* ..trigger inexact flag, sqrt(x) is inexact */
  526.         if(r==RN) t=addc(t);            /* ...t=t+ulp */
  527.         else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */
  528.         y=y+t;                          /* ...chopped sum */
  529.         py[n0]=py[n0]-0x00100000;       /* ...correctly rounded sqrt(x) */
  530. end:    py[n0]=py[n0]+scalx;            /* ...scale back y */
  531.         swapRM(r);                      /* ...restore Rounding Mode */
  532.         return(y);
  533. }
  534. #endif
  535.