diff options
Diffstat (limited to 'libm/ldouble/arcdotl.c')
-rw-r--r-- | libm/ldouble/arcdotl.c | 108 |
1 files changed, 0 insertions, 108 deletions
diff --git a/libm/ldouble/arcdotl.c b/libm/ldouble/arcdotl.c deleted file mode 100644 index 952f027c6..000000000 --- a/libm/ldouble/arcdotl.c +++ /dev/null @@ -1,108 +0,0 @@ -/* arcdot.c - * - * Angle between two vectors - * - * - * - * - * SYNOPSIS: - * - * long double p[3], q[3], arcdotl(); - * - * y = arcdotl( p, q ); - * - * - * - * DESCRIPTION: - * - * For two vectors p, q, the angle A between them is given by - * - * p.q / (|p| |q|) = cos A . - * - * where "." represents inner product, "|x|" the length of vector x. - * If the angle is small, an expression in sin A is preferred. - * Set r = q - p. Then - * - * p.q = p.p + p.r , - * - * |p|^2 = p.p , - * - * |q|^2 = p.p + 2 p.r + r.r , - * - * p.p^2 + 2 p.p p.r + p.r^2 - * cos^2 A = ---------------------------- - * p.p (p.p + 2 p.r + r.r) - * - * p.p + 2 p.r + p.r^2 / p.p - * = --------------------------- , - * p.p + 2 p.r + r.r - * - * sin^2 A = 1 - cos^2 A - * - * r.r - p.r^2 / p.p - * = -------------------- - * p.p + 2 p.r + r.r - * - * = (r.r - p.r^2 / p.p) / q.q . - * - * ACCURACY: - * - * About 1 ULP. See arcdot.c. - * - */ - -/* -Cephes Math Library Release 2.3: November, 1995 -Copyright 1995 by Stephen L. Moshier -*/ - -#include <math.h> -#ifdef ANSIPROT -extern long double sqrtl ( long double ); -extern long double acosl ( long double ); -extern long double asinl ( long double ); -extern long double atanl ( long double ); -#else -long double sqrtl(), acosl(), asinl(), atanl(); -#endif -extern long double PIL; - -long double arcdotl(p,q) -long double p[], q[]; -{ -long double pp, pr, qq, rr, rt, pt, qt, pq; -int i; - -pq = 0.0L; -qq = 0.0L; -pp = 0.0L; -pr = 0.0L; -rr = 0.0L; -for (i=0; i<3; i++) - { - pt = p[i]; - qt = q[i]; - pq += pt * qt; - qq += qt * qt; - pp += pt * pt; - rt = qt - pt; - pr += pt * rt; - rr += rt * rt; - } -if (rr == 0.0L || pp == 0.0L || qq == 0.0L) - return 0.0L; -rt = (rr - (pr * pr) / pp) / qq; -if (rt <= 0.75L) - { - rt = sqrtl(rt); - qt = asinl(rt); - if (pq < 0.0L) - qt = PIL - qt; - } -else - { - pt = pq / sqrtl(pp*qq); - qt = acosl(pt); - } -return qt; -} |