jehanne/qa/lib/c/sqrt.c
Giacomo Tesio e70feee4a3 libc: introduce "jehanne_" namespace
With this commit all functions declared in libc.h have been renamed
with the "jehanne_" prefix. This is done for several reason:

- it removes conflicts during symbol resolution when linking
  standard C libraries like newlib or musl
- it allows programs depending on a standard C library to directly
  link to a library depending on our non standard libc (eg libsec).

To ease transiction two files are provided:

- sys/include/lib9.h that can be included instead of <libc.h> to use
  the old names (via a simple set of macros)
- sys/src/lib/c/lib9.c that can be compiled with a program where the
  macro provided by lib9.h are too dumb (see for example rc or grep).

In the kernel port/lib.h has been modified accordingly and some of
the functions it directly provides has been renamed too (eg malloc
in qmalloc.c and print in devcons.c).
2017-04-19 23:48:21 +02:00

84 lines
1.4 KiB
C

/*
* This file is part of the UCB release of Plan 9. It is subject to the license
* terms in the LICENSE file found in the top-level directory of this
* distribution and at http://akaros.cs.berkeley.edu/files/Plan9License. No
* part of the UCB release of Plan 9, including this file, may be copied,
* modified, propagated, or distributed except according to the terms contained
* in the LICENSE file.
*/
#include <u.h>
#include <lib9.h>
/*
sqrtC returns the square root of its floating
point argument. Newton's method.
calls frexp
*/
double
sqrtC(double arg)
{
double x, temp;
int exp, i;
if(arg <= 0) {
if(arg < 0)
return NaN();
return 0;
}
if(isInf(arg, 1))
return arg;
x = frexp(arg, &exp);
while(x < 0.5) {
x *= 2;
exp--;
}
/*
* NOTE
* this wont work on 1's comp
*/
if(exp & 1) {
x *= 2;
exp--;
}
temp = 0.5 * (1.0+x);
while(exp > 60) {
temp *= (1L<<30);
exp -= 60;
}
while(exp < -60) {
temp /= (1L<<30);
exp += 60;
}
if(exp >= 0)
temp *= 1L << (exp/2);
else
temp /= 1L << (-exp/2);
for(i=0; i<=4; i++)
temp = 0.5*(temp + arg/temp);
return temp;
}
void
main(void)
{
double v;
char *a, *b;
for(v = 4; v < 65536; v += v + 1) {
a = smprint("%f", sqrtC(v));
b = smprint("%f", sqrt(v));
if(strcmp(a, b)){
print("FAIL\n");
exits("FAIL");
}
free(a);
free(b);
}
print("PASS\n");
exits("PASS");
}