Return to support.c CVS log

Up to [CSRG BSD Unix] / 43BSDTahoe / usr.lib / libm / ieee

Annotation of 43BSDTahoe/usr.lib/libm/ieee/support.c, revision 1.1

1.1     ! root        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 are permitted
        !             6:  * provided that the above copyright notice and this paragraph are
        !             7:  * duplicated in all such forms and that any documentation,
        !             8:  * advertising materials, and other materials related to such
        !             9:  * distribution and use acknowledge that the software was developed
        !            10:  * by the University of California, Berkeley.  The name of the
        !            11:  * University may not be used to endorse or promote products derived
        !            12:  * from this software without specific prior written permission.
        !            13:  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
        !            14:  * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
        !            15:  * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
        !            16:  *
        !            17:  * All recipients should regard themselves as participants in an ongoing
        !            18:  * research project and hence should feel obligated to report their
        !            19:  * experiences (good or bad) with these elementary function codes, using
        !            20:  * the sendbug(8) program, to the authors.
        !            21:  */
        !            22: 
        !            23: #ifndef lint
        !            24: static char sccsid[] = "@(#)support.c  5.3 (Berkeley) 6/30/88";
        !            25: #endif /* not lint */
        !            26: 
        !            27: /* 
        !            28:  * Some IEEE standard 754 recommended functions and remainder and sqrt for 
        !            29:  * supporting the C elementary functions.
        !            30:  ******************************************************************************
        !            31:  * WARNING:
        !            32:  *      These codes are developed (in double) to support the C elementary
        !            33:  * functions temporarily. They are not universal, and some of them are very
        !            34:  * slow (in particular, drem and sqrt is extremely inefficient). Each 
        !            35:  * computer system should have its implementation of these functions using 
        !            36:  * its own assembler.
        !            37:  ******************************************************************************
        !            38:  *
        !            39:  * IEEE 754 required operations:
        !            40:  *     drem(x,p) 
        !            41:  *              returns  x REM y  =  x - [x/y]*y , where [x/y] is the integer
        !            42:  *              nearest x/y; in half way case, choose the even one.
        !            43:  *     sqrt(x) 
        !            44:  *              returns the square root of x correctly rounded according to 
        !            45:  *             the rounding mod.
        !            46:  *
        !            47:  * IEEE 754 recommended functions:
        !            48:  * (a) copysign(x,y) 
        !            49:  *              returns x with the sign of y. 
        !            50:  * (b) scalb(x,N) 
        !            51:  *              returns  x * (2**N), for integer values N.
        !            52:  * (c) logb(x) 
        !            53:  *              returns the unbiased exponent of x, a signed integer in 
        !            54:  *              double precision, except that logb(0) is -INF, logb(INF) 
        !            55:  *              is +INF, and logb(NAN) is that NAN.
        !            56:  * (d) finite(x) 
        !            57:  *              returns the value TRUE if -INF < x < +INF and returns 
        !            58:  *              FALSE otherwise.
        !            59:  *
        !            60:  *
        !            61:  * CODED IN C BY K.C. NG, 11/25/84;
        !            62:  * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
        !            63:  */
        !            64: 
        !            65: 
        !            66: #if defined(vax)||defined(tahoe)      /* VAX D format */
        !            67: #include <errno.h>
        !            68:     static unsigned short msign=0x7fff , mexp =0x7f80 ;
        !            69:     static short  prep1=57, gap=7, bias=129           ;   
        !            70:     static double novf=1.7E38, nunf=3.0E-39, zero=0.0 ;
        !            71: #else  /* defined(vax)||defined(tahoe) */
        !            72:     static unsigned short msign=0x7fff, mexp =0x7ff0  ;
        !            73:     static short prep1=54, gap=4, bias=1023           ;
        !            74:     static double novf=1.7E308, nunf=3.0E-308,zero=0.0;
        !            75: #endif /* defined(vax)||defined(tahoe) */
        !            76: 
        !            77: double scalb(x,N)
        !            78: double x; int N;
        !            79: {
        !            80:         int k;
        !            81:         double scalb();
        !            82: 
        !            83: #ifdef national
        !            84:         unsigned short *px=(unsigned short *) &x + 3;
        !            85: #else  /* national */
        !            86:         unsigned short *px=(unsigned short *) &x;
        !            87: #endif /* national */
        !            88: 
        !            89:         if( x == zero )  return(x); 
        !            90: 
        !            91: #if defined(vax)||defined(tahoe)
        !            92:         if( (k= *px & mexp ) != ~msign ) {
        !            93:             if (N < -260)
        !            94:                return(nunf*nunf);
        !            95:            else if (N > 260) {
        !            96:                extern double infnan(),copysign();
        !            97:                return(copysign(infnan(ERANGE),x));
        !            98:            }
        !            99: #else  /* defined(vax)||defined(tahoe) */
        !           100:         if( (k= *px & mexp ) != mexp ) {
        !           101:             if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf);
        !           102:             if( k == 0 ) {
        !           103:                  x *= scalb(1.0,(int)prep1);  N -= prep1; return(scalb(x,N));}
        !           104: #endif /* defined(vax)||defined(tahoe) */
        !           105: 
        !           106:             if((k = (k>>gap)+ N) > 0 )
        !           107:                 if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap);
        !           108:                 else x=novf+novf;               /* overflow */
        !           109:             else
        !           110:                 if( k > -prep1 ) 
        !           111:                                         /* gradual underflow */
        !           112:                     {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);}
        !           113:                 else
        !           114:                 return(nunf*nunf);
        !           115:             }
        !           116:         return(x);
        !           117: }
        !           118: 
        !           119: 
        !           120: double copysign(x,y)
        !           121: double x,y;
        !           122: {
        !           123: #ifdef national
        !           124:         unsigned short  *px=(unsigned short *) &x+3,
        !           125:                         *py=(unsigned short *) &y+3;
        !           126: #else  /* national */
        !           127:         unsigned short  *px=(unsigned short *) &x,
        !           128:                         *py=(unsigned short *) &y;
        !           129: #endif /* national */
        !           130: 
        !           131: #if defined(vax)||defined(tahoe)
        !           132:         if ( (*px & mexp) == 0 ) return(x);
        !           133: #endif /* defined(vax)||defined(tahoe) */
        !           134: 
        !           135:         *px = ( *px & msign ) | ( *py & ~msign );
        !           136:         return(x);
        !           137: }
        !           138: 
        !           139: double logb(x)
        !           140: double x; 
        !           141: {
        !           142: 
        !           143: #ifdef national
        !           144:         short *px=(short *) &x+3, k;
        !           145: #else  /* national */
        !           146:         short *px=(short *) &x, k;
        !           147: #endif /* national */
        !           148: 
        !           149: #if defined(vax)||defined(tahoe)
        !           150:         return (int)(((*px&mexp)>>gap)-bias);
        !           151: #else  /* defined(vax)||defined(tahoe) */
        !           152:         if( (k= *px & mexp ) != mexp )
        !           153:             if ( k != 0 )
        !           154:                 return ( (k>>gap) - bias );
        !           155:             else if( x != zero)
        !           156:                 return ( -1022.0 );
        !           157:             else        
        !           158:                 return(-(1.0/zero));    
        !           159:         else if(x != x)
        !           160:             return(x);
        !           161:         else
        !           162:             {*px &= msign; return(x);}
        !           163: #endif /* defined(vax)||defined(tahoe) */
        !           164: }
        !           165: 
        !           166: finite(x)
        !           167: double x;    
        !           168: {
        !           169: #if defined(vax)||defined(tahoe)
        !           170:         return(1);
        !           171: #else  /* defined(vax)||defined(tahoe) */
        !           172: #ifdef national
        !           173:         return( (*((short *) &x+3 ) & mexp ) != mexp );
        !           174: #else  /* national */
        !           175:         return( (*((short *) &x ) & mexp ) != mexp );
        !           176: #endif /* national */
        !           177: #endif /* defined(vax)||defined(tahoe) */
        !           178: }
        !           179: 
        !           180: double drem(x,p)
        !           181: double x,p;
        !           182: {
        !           183:         short sign;
        !           184:         double hp,dp,tmp,drem(),scalb();
        !           185:         unsigned short  k; 
        !           186: #ifdef national
        !           187:         unsigned short
        !           188:               *px=(unsigned short *) &x  +3, 
        !           189:               *pp=(unsigned short *) &p  +3,
        !           190:               *pd=(unsigned short *) &dp +3,
        !           191:               *pt=(unsigned short *) &tmp+3;
        !           192: #else  /* national */
        !           193:         unsigned short
        !           194:               *px=(unsigned short *) &x  , 
        !           195:               *pp=(unsigned short *) &p  ,
        !           196:               *pd=(unsigned short *) &dp ,
        !           197:               *pt=(unsigned short *) &tmp;
        !           198: #endif /* national */
        !           199: 
        !           200:         *pp &= msign ;
        !           201: 
        !           202: #if defined(vax)||defined(tahoe)
        !           203:         if( ( *px & mexp ) == ~msign ) /* is x a reserved operand? */
        !           204: #else  /* defined(vax)||defined(tahoe) */
        !           205:         if( ( *px & mexp ) == mexp )
        !           206: #endif /* defined(vax)||defined(tahoe) */
        !           207:                return  (x-p)-(x-p);    /* create nan if x is inf */
        !           208:        if (p == zero) {
        !           209: #if defined(vax)||defined(tahoe)
        !           210:                extern double infnan();
        !           211:                return(infnan(EDOM));
        !           212: #else  /* defined(vax)||defined(tahoe) */
        !           213:                return zero/zero;
        !           214: #endif /* defined(vax)||defined(tahoe) */
        !           215:        }
        !           216: 
        !           217: #if defined(vax)||defined(tahoe)
        !           218:         if( ( *pp & mexp ) == ~msign ) /* is p a reserved operand? */
        !           219: #else  /* defined(vax)||defined(tahoe) */
        !           220:         if( ( *pp & mexp ) == mexp )
        !           221: #endif /* defined(vax)||defined(tahoe) */
        !           222:                { if (p != p) return p; else return x;}
        !           223: 
        !           224:         else  if ( ((*pp & mexp)>>gap) <= 1 ) 
        !           225:                 /* subnormal p, or almost subnormal p */
        !           226:             { double b; b=scalb(1.0,(int)prep1);
        !           227:               p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);}
        !           228:         else  if ( p >= novf/2)
        !           229:             { p /= 2 ; x /= 2; return(drem(x,p)*2);}
        !           230:         else 
        !           231:             {
        !           232:                 dp=p+p; hp=p/2;
        !           233:                 sign= *px & ~msign ;
        !           234:                 *px &= msign       ;
        !           235:                 while ( x > dp )
        !           236:                     {
        !           237:                         k=(*px & mexp) - (*pd & mexp) ;
        !           238:                         tmp = dp ;
        !           239:                         *pt += k ;
        !           240: 
        !           241: #if defined(vax)||defined(tahoe)
        !           242:                         if( x < tmp ) *pt -= 128 ;
        !           243: #else  /* defined(vax)||defined(tahoe) */
        !           244:                         if( x < tmp ) *pt -= 16 ;
        !           245: #endif /* defined(vax)||defined(tahoe) */
        !           246: 
        !           247:                         x -= tmp ;
        !           248:                     }
        !           249:                 if ( x > hp )
        !           250:                     { x -= p ;  if ( x >= hp ) x -= p ; }
        !           251: 
        !           252: #if defined(vax)||defined(tahoe)
        !           253:                if (x)
        !           254: #endif /* defined(vax)||defined(tahoe) */
        !           255:                        *px ^= sign;
        !           256:                 return( x);
        !           257: 
        !           258:             }
        !           259: }
        !           260: double sqrt(x)
        !           261: double x;
        !           262: {
        !           263:         double q,s,b,r;
        !           264:         double logb(),scalb();
        !           265:         double t,zero=0.0;
        !           266:         int m,n,i,finite();
        !           267: #if defined(vax)||defined(tahoe)
        !           268:         int k=54;
        !           269: #else  /* defined(vax)||defined(tahoe) */
        !           270:         int k=51;
        !           271: #endif /* defined(vax)||defined(tahoe) */
        !           272: 
        !           273:     /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
        !           274:         if(x!=x||x==zero) return(x);
        !           275: 
        !           276:     /* sqrt(negative) is invalid */
        !           277:         if(x<zero) {
        !           278: #if defined(vax)||defined(tahoe)
        !           279:                extern double infnan();
        !           280:                return (infnan(EDOM));  /* NaN */
        !           281: #else  /* defined(vax)||defined(tahoe) */
        !           282:                return(zero/zero);
        !           283: #endif /* defined(vax)||defined(tahoe) */
        !           284:        }
        !           285: 
        !           286:     /* sqrt(INF) is INF */
        !           287:         if(!finite(x)) return(x);               
        !           288: 
        !           289:     /* scale x to [1,4) */
        !           290:         n=logb(x);
        !           291:         x=scalb(x,-n);
        !           292:         if((m=logb(x))!=0) x=scalb(x,-m);       /* subnormal number */
        !           293:         m += n; 
        !           294:         n = m/2;
        !           295:         if((n+n)!=m) {x *= 2; m -=1; n=m/2;}
        !           296: 
        !           297:     /* generate sqrt(x) bit by bit (accumulating in q) */
        !           298:             q=1.0; s=4.0; x -= 1.0; r=1;
        !           299:             for(i=1;i<=k;i++) {
        !           300:                 t=s+1; x *= 4; r /= 2;
        !           301:                 if(t<=x) {
        !           302:                     s=t+t+2, x -= t; q += r;}
        !           303:                 else
        !           304:                     s *= 2;
        !           305:                 }
        !           306:             
        !           307:     /* generate the last bit and determine the final rounding */
        !           308:             r/=2; x *= 4; 
        !           309:             if(x==zero) goto end; 100+r; /* trigger inexact flag */
        !           310:             if(s<x) {
        !           311:                 q+=r; x -=s; s += 2; s *= 2; x *= 4;
        !           312:                 t = (x-s)-5; 
        !           313:                 b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */
        !           314:                 b=1.0+r/4;   if(b>1.0) t=1;    /* b>1 : Round-to-(+INF) */
        !           315:                 if(t>=0) q+=r; }             /* else: Round-to-nearest */
        !           316:             else { 
        !           317:                 s *= 2; x *= 4; 
        !           318:                 t = (x-s)-1; 
        !           319:                 b=1.0+3*r/4; if(b==1.0) goto end;
        !           320:                 b=1.0+r/4;   if(b>1.0) t=1;
        !           321:                 if(t>=0) q+=r; }
        !           322:             
        !           323: end:        return(scalb(q,n));
        !           324: }
        !           325: 
        !           326: #if 0
        !           327: /* DREM(X,Y)
        !           328:  * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
        !           329:  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
        !           330:  * INTENDED FOR ASSEMBLY LANGUAGE
        !           331:  * CODED IN C BY K.C. NG, 3/23/85, 4/8/85.
        !           332:  *
        !           333:  * Warning: this code should not get compiled in unless ALL of
        !           334:  * the following machine-dependent routines are supplied.
        !           335:  * 
        !           336:  * Required machine dependent functions (not on a VAX):
        !           337:  *     swapINX(i): save inexact flag and reset it to "i"
        !           338:  *     swapENI(e): save inexact enable and reset it to "e"
        !           339:  */
        !           340: 
        !           341: double drem(x,y)       
        !           342: double x,y;
        !           343: {
        !           344: 
        !           345: #ifdef national                /* order of words in floating point number */
        !           346:        static n0=3,n1=2,n2=1,n3=0;
        !           347: #else /* VAX, SUN, ZILOG, TAHOE */
        !           348:        static n0=0,n1=1,n2=2,n3=3;
        !           349: #endif
        !           350: 
        !           351:        static unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390;
        !           352:        static double zero=0.0;
        !           353:        double hy,y1,t,t1;
        !           354:        short k;
        !           355:        long n;
        !           356:        int i,e; 
        !           357:        unsigned short xexp,yexp, *px  =(unsigned short *) &x  , 
        !           358:                        nx,nf,    *py  =(unsigned short *) &y  ,
        !           359:                        sign,     *pt  =(unsigned short *) &t  ,
        !           360:                                  *pt1 =(unsigned short *) &t1 ;
        !           361: 
        !           362:        xexp = px[n0] & mexp ;  /* exponent of x */
        !           363:        yexp = py[n0] & mexp ;  /* exponent of y */
        !           364:        sign = px[n0] &0x8000;  /* sign of x     */
        !           365: 
        !           366: /* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */
        !           367:        if(x!=x) return(x); if(y!=y) return(y);      /* x or y is NaN */
        !           368:        if( xexp == mexp )   return(zero/zero);      /* x is INF */
        !           369:        if(y==zero) return(y/y);
        !           370: 
        !           371: /* save the inexact flag and inexact enable in i and e respectively
        !           372:  * and reset them to zero
        !           373:  */
        !           374:        i=swapINX(0);   e=swapENI(0);   
        !           375: 
        !           376: /* subnormal number */
        !           377:        nx=0;
        !           378:        if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;}
        !           379: 
        !           380: /* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */
        !           381:        if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;}
        !           382: 
        !           383:        nf=nx;
        !           384:        py[n0] &= 0x7fff;       
        !           385:        px[n0] &= 0x7fff;
        !           386: 
        !           387: /* mask off the least significant 27 bits of y */
        !           388:        t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t;
        !           389: 
        !           390: /* LOOP: argument reduction on x whenever x > y */
        !           391: loop:
        !           392:        while ( x > y )
        !           393:        {
        !           394:            t=y;
        !           395:            t1=y1;
        !           396:            xexp=px[n0]&mexp;     /* exponent of x */
        !           397:            k=xexp-yexp-m25;
        !           398:            if(k>0)     /* if x/y >= 2**26, scale up y so that x/y < 2**26 */
        !           399:                {pt[n0]+=k;pt1[n0]+=k;}
        !           400:            n=x/t; x=(x-n*t1)-n*(t-t1);
        !           401:        }       
        !           402:     /* end while (x > y) */
        !           403: 
        !           404:        if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;}
        !           405: 
        !           406: /* final adjustment */
        !           407: 
        !           408:        hy=y/2.0;
        !           409:        if(x>hy||((x==hy)&&n%2==1)) x-=y; 
        !           410:        px[n0] ^= sign;
        !           411:        if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;}
        !           412: 
        !           413: /* restore inexact flag and inexact enable */
        !           414:        swapINX(i); swapENI(e); 
        !           415: 
        !           416:        return(x);      
        !           417: }
        !           418: #endif
        !           419: 
        !           420: #if 0
        !           421: /* SQRT
        !           422:  * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT
        !           423:  * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE
        !           424:  * CODED IN C BY K.C. NG, 3/22/85.
        !           425:  *
        !           426:  * Warning: this code should not get compiled in unless ALL of
        !           427:  * the following machine-dependent routines are supplied.
        !           428:  * 
        !           429:  * Required machine dependent functions:
        !           430:  *     swapINX(i)  ...return the status of INEXACT flag and reset it to "i"
        !           431:  *     swapRM(r)   ...return the current Rounding Mode and reset it to "r"
        !           432:  *     swapENI(e)  ...return the status of inexact enable and reset it to "e"
        !           433:  *     addc(t)     ...perform t=t+1 regarding t as a 64 bit unsigned integer
        !           434:  *     subc(t)     ...perform t=t-1 regarding t as a 64 bit unsigned integer
        !           435:  */
        !           436: 
        !           437: static unsigned long table[] = {
        !           438: 0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740,
        !           439: 58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478,
        !           440: 21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, };
        !           441: 
        !           442: double newsqrt(x)
        !           443: double x;
        !           444: {
        !           445:         double y,z,t,addc(),subc(),b54=134217728.*134217728.; /* b54=2**54 */
        !           446:         long mx,scalx,mexp=0x7ff00000;
        !           447:         int i,j,r,e,swapINX(),swapRM(),swapENI();       
        !           448:         unsigned long *py=(unsigned long *) &y   ,
        !           449:                       *pt=(unsigned long *) &t   ,
        !           450:                       *px=(unsigned long *) &x   ;
        !           451: #ifdef national         /* ordering of word in a floating point number */
        !           452:         int n0=1, n1=0; 
        !           453: #else
        !           454:         int n0=0, n1=1; 
        !           455: #endif
        !           456: /* Rounding Mode:  RN ...round-to-nearest 
        !           457:  *                 RZ ...round-towards 0
        !           458:  *                 RP ...round-towards +INF
        !           459:  *                RM ...round-towards -INF
        !           460:  */
        !           461:         int RN=0,RZ=1,RP=2,RM=3;/* machine dependent: work on a Zilog Z8070
        !           462:                                  * and a National 32081 & 16081
        !           463:                                  */
        !           464: 
        !           465: /* exceptions */
        !           466:        if(x!=x||x==0.0) return(x);  /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
        !           467:        if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */
        !           468:         if((mx=px[n0]&mexp)==mexp) return(x);  /* sqrt(+INF) is +INF */
        !           469: 
        !           470: /* save, reset, initialize */
        !           471:         e=swapENI(0);   /* ...save and reset the inexact enable */
        !           472:         i=swapINX(0);   /* ...save INEXACT flag */
        !           473:         r=swapRM(RN);   /* ...save and reset the Rounding Mode to RN */
        !           474:         scalx=0;
        !           475: 
        !           476: /* subnormal number, scale up x to x*2**54 */
        !           477:         if(mx==0) {x *= b54 ; scalx-=0x01b00000;}
        !           478: 
        !           479: /* scale x to avoid intermediate over/underflow:
        !           480:  * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */
        !           481:         if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;}
        !           482:         if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;}
        !           483: 
        !           484: /* magic initial approximation to almost 8 sig. bits */
        !           485:         py[n0]=(px[n0]>>1)+0x1ff80000;
        !           486:         py[n0]=py[n0]-table[(py[n0]>>15)&31];
        !           487: 
        !           488: /* Heron's rule once with correction to improve y to almost 18 sig. bits */
        !           489:         t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0;
        !           490: 
        !           491: /* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */
        !           492:         t=y*y; z=t;  pt[n0]+=0x00100000; t+=z; z=(x-z)*y; 
        !           493:         t=z/(t+x) ;  pt[n0]+=0x00100000; y+=t;
        !           494: 
        !           495: /* twiddle last bit to force y correctly rounded */ 
        !           496:         swapRM(RZ);     /* ...set Rounding Mode to round-toward-zero */
        !           497:         swapINX(0);     /* ...clear INEXACT flag */
        !           498:         swapENI(e);     /* ...restore inexact enable status */
        !           499:         t=x/y;          /* ...chopped quotient, possibly inexact */
        !           500:         j=swapINX(i);   /* ...read and restore inexact flag */
        !           501:         if(j==0) { if(t==y) goto end; else t=subc(t); }  /* ...t=t-ulp */
        !           502:         b54+0.1;        /* ..trigger inexact flag, sqrt(x) is inexact */
        !           503:         if(r==RN) t=addc(t);            /* ...t=t+ulp */
        !           504:         else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */
        !           505:         y=y+t;                          /* ...chopped sum */
        !           506:         py[n0]=py[n0]-0x00100000;       /* ...correctly rounded sqrt(x) */
        !           507: end:    py[n0]=py[n0]+scalx;            /* ...scale back y */
        !           508:         swapRM(r);                      /* ...restore Rounding Mode */
        !           509:         return(y);
        !           510: }
        !           511: #endif

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.