Mercurial > hg > octave-nkf > gnulib-hg

annotate lib/sqrtl.c @ 7292:17785d5bede0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix docstrings
author Sergey Poznyakoff <gray@gnu.org.ua>
date Sun, 10 Sep 2006 11:52:44 +0000
parents a48fb0e98c8c
children c115624bf99a
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
1 /* Emulation for sqrtl.
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
2 Contributed by Paolo Bonzini
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
3
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
4 Copyright 2002, 2003 Free Software Foundation, Inc.
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
5
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
6 This file is part of gnulib.
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
7
4435
d4efdb886b2a Switch from LGPL to GPL.
Paul Eggert <eggert@cs.ucla.edu>
parents: 4239
diff changeset
8 This program is free software; you can redistribute it and/or modify
d4efdb886b2a Switch from LGPL to GPL.
Paul Eggert <eggert@cs.ucla.edu>
parents: 4239
diff changeset
9 it under the terms of the GNU General Public License as published by
d4efdb886b2a Switch from LGPL to GPL.
Paul Eggert <eggert@cs.ucla.edu>
parents: 4239
diff changeset
10 the Free Software Foundation; either version 2, or (at your option)
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
11 any later version.
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
12
4435
d4efdb886b2a Switch from LGPL to GPL.
Paul Eggert <eggert@cs.ucla.edu>
parents: 4239
diff changeset
13 This program is distributed in the hope that it will be useful,
d4efdb886b2a Switch from LGPL to GPL.
Paul Eggert <eggert@cs.ucla.edu>
parents: 4239
diff changeset
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
d4efdb886b2a Switch from LGPL to GPL.
Paul Eggert <eggert@cs.ucla.edu>
parents: 4239
diff changeset
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
d4efdb886b2a Switch from LGPL to GPL.
Paul Eggert <eggert@cs.ucla.edu>
parents: 4239
diff changeset
16 GNU General Public License for more details.
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
17
4435
d4efdb886b2a Switch from LGPL to GPL.
Paul Eggert <eggert@cs.ucla.edu>
parents: 4239
diff changeset
18 You should have received a copy of the GNU General Public License along
d4efdb886b2a Switch from LGPL to GPL.
Paul Eggert <eggert@cs.ucla.edu>
parents: 4239
diff changeset
19 with this program; if not, write to the Free Software Foundation,
5848
a48fb0e98c8c *** empty log message ***
Paul Eggert <eggert@cs.ucla.edu>
parents: 4435
diff changeset
20 Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
21
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
22 #include <float.h>
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
23 #include <math.h>
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
24
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
25 #include "mathl.h"
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
26
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
27 /* A simple Newton-Raphson method. */
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
28 long double
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
29 sqrtl(long double x)
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
30 {
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
31 long double delta, y;
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
32 int exponent;
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
33
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
34 /* Check for negative numbers */
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
35 if (x < 0.0L)
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
36 return (long double) sqrt(-1);
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
37
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
38 /* Check for zero, NANs and infinites */
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
39 if (x + x == x || x != x)
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
40 return x;
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
41
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
42 frexpl (x, &exponent);
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
43 y = ldexpl (x, -exponent / 2);
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
44
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
45 do
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
46 {
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
47 delta = y;
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
48 y = (y + x / y) * 0.5L;
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
49 delta -= y;
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
50 }
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
51 while (delta != 0.0L);
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
52
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
53 return y;
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
54 }