/* @(#)e_hypot.c 5.1 93/09/24 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include #if defined(LIBM_SCCS) && !defined(lint) __RCSID("$NetBSD: e_hypot.c,v 1.13 2008/04/25 22:21:53 christos Exp $"); #endif /* __ieee754_hypot(x,y) * * Method : * If (assume round-to-nearest) z=x*x+y*y * has error less than sqrt(2)/2 ulp, than * sqrt(z) has error less than 1 ulp (exercise). * * So, compute sqrt(x*x+y*y) with some care as * follows to get the error below 1 ulp: * * Assume x>y>0; * (if possible, set rounding to round-to-nearest) * 1. if x > 2y use * x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y * where x1 = x with lower 32 bits cleared, x2 = x-x1; else * 2. if x <= 2y use * t1*yy1+((x-y)*(x-y)+(t1*y2+t2*y)) * where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1, * yy1= y with lower 32 bits chopped, y2 = y-yy1. * * NOTE: scaling may be necessary if some argument is too * large or too tiny * * Special cases: * hypot(x,y) is INF if x or y is +INF or -INF; else * hypot(x,y) is NAN if x or y is NAN. * * Accuracy: * hypot(x,y) returns sqrt(x^2+y^2) with error less * than 1 ulps (units in the last place) */ #include "math.h" #include "math_private.h" double __ieee754_hypot(double x, double y) { double a=x,b=y,t1,t2,yy1,y2,w; int32_t j,k,ha,hb; GET_HIGH_WORD(ha,x); ha &= 0x7fffffff; GET_HIGH_WORD(hb,y); hb &= 0x7fffffff; if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;} SET_HIGH_WORD(a,ha); /* a <- |a| */ SET_HIGH_WORD(b,hb); /* b <- |b| */ if((ha-hb)>0x3c00000) {return a+b;} /* x/y > 2**60 */ k=0; if(ha > 0x5f300000) { /* a>2**500 */ if(ha >= 0x7ff00000) { /* Inf or NaN */ u_int32_t low; w = a+b; /* for sNaN */ GET_LOW_WORD(low,a); if(((ha&0xfffff)|low)==0) w = a; GET_LOW_WORD(low,b); if(((hb^0x7ff00000)|low)==0) w = b; return w; } /* scale a and b by 2**-600 */ ha -= 0x25800000; hb -= 0x25800000; k += 600; SET_HIGH_WORD(a,ha); SET_HIGH_WORD(b,hb); } if(hb < 0x20b00000) { /* b < 2**-500 */ if(hb <= 0x000fffff) { /* subnormal b or 0 */ u_int32_t low; GET_LOW_WORD(low,b); if((hb|low)==0) return a; t1=0; SET_HIGH_WORD(t1,0x7fd00000); /* t1=2^1022 */ b *= t1; a *= t1; k -= 1022; } else { /* scale a and b by 2^600 */ ha += 0x25800000; /* a *= 2^600 */ hb += 0x25800000; /* b *= 2^600 */ k -= 600; SET_HIGH_WORD(a,ha); SET_HIGH_WORD(b,hb); } } /* medium size a and b */ w = a-b; if (w>b) { t1 = 0; SET_HIGH_WORD(t1,ha); t2 = a-t1; w = __ieee754_sqrt(t1*t1-(b*(-b)-t2*(a+t1))); } else { a = a+a; yy1 = 0; SET_HIGH_WORD(yy1,hb); y2 = b - yy1; t1 = 0; SET_HIGH_WORD(t1,ha+0x00100000); t2 = a - t1; w = __ieee754_sqrt(t1*yy1-(w*(-w)-(t1*y2+t2*b))); } if(k!=0) { u_int32_t high; t1 = 1.0; GET_HIGH_WORD(high,t1); SET_HIGH_WORD(t1,high+(k<<20)); return t1*w; } else return w; } olspan='3' class='logmsg'> Since timestamp was introduced, it was too small to store full dmesg of one boot. And putting this behind of #include <m68k/param.h> had no effect. This was fixed on many arch in 1997 but x68k was not... http://mail-index.netbsd.org/source-changes/1997/09/20/0021.html 2012-02-10proper sizing of kmem_arena on different portspara PR port-i386/45946: Kernel locks up in VMEM system 2011-02-08Remove clause 3 (UCB advertising clause) from the University of Utahrmind copyright. Confirmed by Mike Hibler, mike at cs.utah.edu - thanks! Also, merge UCB and Utah copyright texts back into one, as they originally were. Extra verification by snj@. 2009-12-11Refactoring MD segment related definitions on m68k ports:tsutsui - move macro SEGSHIFT, NBSEG, and SEGOFSET from <machine/param.h> to <m68k/pmap_motorola.h> where they should belong - also simplify NBSEG macro for both 020/030 and 040/060 cases - also move m68k_trunc_seg(), m68k_round_seg(), and m68k_seg_offset() (that use SEGOFSET) from <m68k/param.h> to <m68k/pmap_motorola.h> - add comments about each segment size value used on pmap_motorola implementation - add TIA_SIZE, TIA_SHIFT, TIB_SIZE, and TIB_SHIFT macro which can be used for ste/pte index sizes - use proper vaddr_t cast instead of unsigned for SEGOFSET/PGOFSET macro Tested on atari, hp300, and news68k. 2009-12-09Fix wrong NBSEG values for all hp300 pmap derived m68k ports.tsutsui They were incorrect since 1997 on amiga and atari, and since 2002 on other ports, but fortunately they don't look so fatal. Anyway, these values will be moved into <m68k/pmap_motorola.h> soon since they are quite pmap implementation dependent. 2007-12-31Remove COMPAT_HPUX.ad 2007-03-11Merge #ifdef _KERNEL stuff into one place.isaki Suggested by tsutsui@. 2007-03-11Include <machine/intr.h> only in _KERNEL.isaki Suggested by tsutsui@. 2005-12-11merge ktrace-lwp.christos 2005-01-18de-__P, remove register, ansify.chs 2003-08-07Move UCB-licensed code from 4-clause to 3-clause licence.agc Patches provided by Joel Baker in PR 22364, verified by myself. 2003-01-18Merge the nathanw_sa branch.thorpej 2002-11-03merge pmap.h for all 4k-page motorola-MMU (or compatible) m68k platforms.chs 2001-05-30use _KERNEL_OPTmrg 2001-01-11Sync. w/ luna68k; eliminate irrelavant code imposed byminoura the restriction of hp300 hardware and HP-UX compatibility. 2000-02-11Update for the NKMEMPAGES changes.thorpej 1999-12-04CL* discarding.ragge 1999-05-05Use mvme68k-derived delay routine.minoura Now that we support various models, the old constant-loop delay routine may cause problems. 1998-12-13Splitted out the interrupt handling code to a separate file.minoura Disallow getting down the interrupt level by splXX. 1998-06-30sync with -current after a long silencemsaitoh 1998-06-25defopt COMPAT_HPUXthorpej