[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Toon-members] tag Makefile.am Makefile.in libtool src/five_po...
|
From: |
Gerhard Reitmayr |
|
Subject: |
[Toon-members] tag Makefile.am Makefile.in libtool src/five_po... |
|
Date: |
Wed, 22 Apr 2009 23:23:13 +0000 |
CVSROOT: /cvsroot/toon
Module name: tag
Changes by: Gerhard Reitmayr <gerhard> 09/04/22 23:23:13
Modified files:
. : Makefile.am Makefile.in libtool
src : five_point.cpp
Added files:
src : solve.h sturm.c util.c
Removed files:
src : ccmcomplex.h plrt.c
Log message:
replaced polynomial root solver with sturm chain-based one from Graphic
Gems 1 -> works
CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/tag/Makefile.am?cvsroot=toon&r1=1.12&r2=1.13
http://cvs.savannah.gnu.org/viewcvs/tag/Makefile.in?cvsroot=toon&r1=1.12&r2=1.13
http://cvs.savannah.gnu.org/viewcvs/tag/libtool?cvsroot=toon&r1=1.16&r2=1.17
http://cvs.savannah.gnu.org/viewcvs/tag/src/five_point.cpp?cvsroot=toon&r1=1.3&r2=1.4
http://cvs.savannah.gnu.org/viewcvs/tag/src/solve.h?cvsroot=toon&rev=1.1
http://cvs.savannah.gnu.org/viewcvs/tag/src/sturm.c?cvsroot=toon&rev=1.1
http://cvs.savannah.gnu.org/viewcvs/tag/src/util.c?cvsroot=toon&rev=1.1
http://cvs.savannah.gnu.org/viewcvs/tag/src/ccmcomplex.h?cvsroot=toon&r1=1.1&r2=0
http://cvs.savannah.gnu.org/viewcvs/tag/src/plrt.c?cvsroot=toon&r1=1.1&r2=0
Patches:
Index: Makefile.am
===================================================================
RCS file: /cvsroot/toon/tag/Makefile.am,v
retrieving revision 1.12
retrieving revision 1.13
diff -u -b -r1.12 -r1.13
--- Makefile.am 22 Apr 2009 19:59:36 -0000 1.12
+++ Makefile.am 22 Apr 2009 23:23:11 -0000 1.13
@@ -2,7 +2,7 @@
lib_LTLIBRARIES = libtoontag.la
libtoontag_la_SOURCES = src/quartic.cpp src/threepointpose.cpp
src/absorient.cpp \
- src/handeye.cpp src/five_point.cpp
src/five_point_buildmatrix.cpp src/plrt.c
+ src/handeye.cpp src/five_point.cpp
src/five_point_buildmatrix.cpp src/sturm.c src/util.c
# the header files to be installed
pkginclude_HEADERS = \
Index: Makefile.in
===================================================================
RCS file: /cvsroot/toon/tag/Makefile.in,v
retrieving revision 1.12
retrieving revision 1.13
diff -u -b -r1.12 -r1.13
--- Makefile.in 22 Apr 2009 19:59:36 -0000 1.12
+++ Makefile.in 22 Apr 2009 23:23:11 -0000 1.13
@@ -59,7 +59,8 @@
LTLIBRARIES = $(lib_LTLIBRARIES)
libtoontag_la_LIBADD =
am_libtoontag_la_OBJECTS = quartic.lo threepointpose.lo absorient.lo \
- handeye.lo five_point.lo five_point_buildmatrix.lo plrt.lo
+ handeye.lo five_point.lo five_point_buildmatrix.lo sturm.lo \
+ util.lo
libtoontag_la_OBJECTS = $(am_libtoontag_la_OBJECTS)
DEFAULT_INCLUDES = address@hidden@
depcomp = $(SHELL) $(top_srcdir)/depcomp
@@ -207,7 +208,7 @@
top_srcdir = @top_srcdir@
lib_LTLIBRARIES = libtoontag.la
libtoontag_la_SOURCES = src/quartic.cpp src/threepointpose.cpp
src/absorient.cpp \
- src/handeye.cpp src/five_point.cpp
src/five_point_buildmatrix.cpp src/plrt.c
+ src/handeye.cpp src/five_point.cpp
src/five_point_buildmatrix.cpp src/sturm.c src/util.c
# the header files to be installed
@@ -309,9 +310,10 @@
@AMDEP_TRUE@@am__include@ @address@hidden/$(DEPDIR)/address@hidden@
@AMDEP_TRUE@@am__include@ @address@hidden/$(DEPDIR)/address@hidden@
@AMDEP_TRUE@@am__include@ @address@hidden/$(DEPDIR)/address@hidden@
address@hidden@@am__include@ @address@hidden/$(DEPDIR)/address@hidden@
@AMDEP_TRUE@@am__include@ @address@hidden/$(DEPDIR)/address@hidden@
address@hidden@@am__include@ @address@hidden/$(DEPDIR)/address@hidden@
@AMDEP_TRUE@@am__include@ @address@hidden/$(DEPDIR)/address@hidden@
address@hidden@@am__include@ @address@hidden/$(DEPDIR)/address@hidden@
.c.o:
@am__fastdepCC_TRUE@ $(COMPILE) -MT $@ -MD -MP -MF $(DEPDIR)/$*.Tpo -c -o $@
$<
@@ -334,12 +336,19 @@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE)
$(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(LTCOMPILE) -c -o $@ $<
-plrt.lo: src/plrt.c
address@hidden@ $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS)
--mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS)
$(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT plrt.lo -MD -MP -MF $(DEPDIR)/plrt.Tpo
-c -o plrt.lo `test -f 'src/plrt.c' || echo '$(srcdir)/'`src/plrt.c
address@hidden@ mv -f $(DEPDIR)/plrt.Tpo $(DEPDIR)/plrt.Plo
address@hidden@@am__fastdepCC_FALSE@ source='src/plrt.c' object='plrt.lo'
libtool=yes @AMDEPBACKSLASH@
+sturm.lo: src/sturm.c
address@hidden@ $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS)
--mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS)
$(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT sturm.lo -MD -MP -MF $(DEPDIR)/sturm.Tpo
-c -o sturm.lo `test -f 'src/sturm.c' || echo '$(srcdir)/'`src/sturm.c
address@hidden@ mv -f $(DEPDIR)/sturm.Tpo $(DEPDIR)/sturm.Plo
address@hidden@@am__fastdepCC_FALSE@ source='src/sturm.c' object='sturm.lo'
libtool=yes @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE)
$(depcomp) @AMDEPBACKSLASH@
address@hidden@ $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS)
--mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS)
$(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o plrt.lo `test -f 'src/plrt.c' || echo
'$(srcdir)/'`src/plrt.c
address@hidden@ $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS)
--mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS)
$(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o sturm.lo `test -f 'src/sturm.c' ||
echo '$(srcdir)/'`src/sturm.c
+
+util.lo: src/util.c
address@hidden@ $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS)
--mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS)
$(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT util.lo -MD -MP -MF $(DEPDIR)/util.Tpo
-c -o util.lo `test -f 'src/util.c' || echo '$(srcdir)/'`src/util.c
address@hidden@ mv -f $(DEPDIR)/util.Tpo $(DEPDIR)/util.Plo
address@hidden@@am__fastdepCC_FALSE@ source='src/util.c' object='util.lo'
libtool=yes @AMDEPBACKSLASH@
address@hidden@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE)
$(depcomp) @AMDEPBACKSLASH@
address@hidden@ $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS)
--mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS)
$(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o util.lo `test -f 'src/util.c' || echo
'$(srcdir)/'`src/util.c
.cpp.o:
@am__fastdepCXX_TRUE@ $(CXXCOMPILE) -MT $@ -MD -MP -MF $(DEPDIR)/$*.Tpo -c -o
$@ $<
Index: libtool
===================================================================
RCS file: /cvsroot/toon/tag/libtool,v
retrieving revision 1.16
retrieving revision 1.17
diff -u -b -r1.16 -r1.17
--- libtool 22 Apr 2009 19:59:36 -0000 1.16
+++ libtool 22 Apr 2009 23:23:11 -0000 1.17
@@ -44,7 +44,7 @@
# ### BEGIN LIBTOOL CONFIG
-# Libtool was configured on host pirx.eng.cam.ac.uk:
+# Libtool was configured on host cpc2-cmbg6-0-0-cust797.cmbg.cable.ntl.com:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -6809,7 +6809,7 @@
# End:
# ### BEGIN LIBTOOL TAG CONFIG: CXX
-# Libtool was configured on host pirx.eng.cam.ac.uk:
+# Libtool was configured on host cpc2-cmbg6-0-0-cust797.cmbg.cable.ntl.com:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
Index: src/five_point.cpp
===================================================================
RCS file: /cvsroot/toon/tag/src/five_point.cpp,v
retrieving revision 1.3
retrieving revision 1.4
diff -u -b -r1.3 -r1.4
--- src/five_point.cpp 22 Apr 2009 20:47:47 -0000 1.3
+++ src/five_point.cpp 22 Apr 2009 23:23:12 -0000 1.4
@@ -1,20 +1,69 @@
#include <tag/five_point.h>
+#include <tag/stdpp.h>
#include <TooN/helpers.h>
#include <TooN/gauss_jordan.h>
#include <TooN/SVD.h>
-#include "ccmcomplex.h"
-
using namespace TooN;
using namespace std;
using namespace std::tr1;
-//CCMath's polyroot
-extern "C" int plrt(double *cof,int n,Cpx *root,double ra,double rb);
+extern "C" {
+#include "solve.h"
+}
namespace tag {
+vector<double> get_roots(const Vector<11> & p){
+ // uses Graphics Gem I code
+ // for details see here:
+ // http://tog.acm.org/GraphicsGems/
+
+ poly sseq[10];
+ Vector<11,double,Reference> v(sseq[0].coef);
+ v = p;
+
+ int num_poly = buildsturm(10, sseq);
+ int atmin, atmax;
+ int nroots = numroots(num_poly, sseq, &atmin, &atmax);
+
+ vector<double> roots(nroots);
+ if(nroots == 0)
+ return roots;
+
+ double min = -1.0;
+ int nchanges = numchanges(num_poly, sseq, min);
+ for (int i = 0; nchanges != atmin && i != MAXPOW; i++) {
+ min *= 10.0;
+ nchanges = numchanges(num_poly, sseq, min);
+ }
+
+ if (nchanges != atmin) {
+ cerr << "solve: unable to bracket all negative roots"
<< endl;
+ atmin = nchanges;
+ }
+
+ double max = 1.0;
+ nchanges = numchanges(num_poly, sseq, max);
+ for (int i = 0; nchanges != atmax && i != MAXPOW; i++) {
+ max *= 10.0;
+ nchanges = numchanges(num_poly, sseq, max);
+ }
+
+ if (nchanges != atmax) {
+ cerr << "solve: unable to bracket all positive roots"
<< endl;
+ atmax = nchanges;
+ }
+
+ nroots = atmin - atmax;
+
+ /*
+ * perform the bisection.
+ */
+ sbisect(num_poly, sseq, min, max, atmin, atmax, &roots[0]);
+}
+
void build_matrix(const Vector<9>& X, const Vector<9>& Y, const Vector<9>& Z,
const Vector<9>& W, Matrix<10,20>& R);
Vector<9> stack_points(const Vector<3>&p, const Vector<3>& q)
@@ -60,6 +109,23 @@
return val;
}
+template<int N> Vector<N-1> poly_diff(const Vector<N> & v)
+{
+ Vector<N-1> ret;
+ for(int i = 1; i < N; ++i){
+ ret[i-1] = v[i]*i;
+ }
+ return ret;
+}
+
+template<int N> pair<Vector<2>, Vector<N-1> > poly_div(const Vector<N> & num,
const Vector<N-1> & denom)
+{
+ Vector<2> f;
+ Vector<N-1> r;
+
+
+}
+
Matrix<3, 3, double, Reference::RowMajor> as_matrix(Vector<9>& v)
{
return Matrix<3, 3, double, Reference::RowMajor>(&v[0]);
@@ -165,23 +231,12 @@
//The polynomial is...
Vector<11> n = poly_mul(p1, b_31) + poly_mul(p2, b_32) + poly_mul(p3,
b_33);
-
- //Use the CCMath root finder.
- //-1, 1 seem to work OK.
- //
- // Return value means:
- // 0 -> normal exit
- // m>0 -> convergence failure for roots k<m
- Cpx roots[10];
- int num = plrt(&n[0], 10, roots, 1, -1);
-
+ vector<double> roots = get_roots(n);
vector<Matrix<3> > Es;
- for(int i=num; i < 10; i++)
- {
- if(abs(roots[i].im) < 1e-10)
+ for(int i=0; i <roots.size(); i++)
{
- double z = roots[i].re;
+ double z = roots[i];
//Solve the linear system in x, y for the forst two
rows of b
// [b11 b12 b13] [x] [ 0 ]
@@ -199,7 +254,6 @@
Es.push_back(x * as_matrix(X) + y*as_matrix(Y) +
z*as_matrix(Z) + as_matrix(W));
}
- }
return Es;
}
Index: src/solve.h
===================================================================
RCS file: src/solve.h
diff -N src/solve.h
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ src/solve.h 22 Apr 2009 23:23:12 -0000 1.1
@@ -0,0 +1,41 @@
+
+/*
+ * solve.h
+ *
+ * some useful constants and types.
+ */
+#define MAX_ORDER 12
+/* maximum order for a polynomial */
+
+#define RELERROR 1.0e-14
+/* smallest relative error we want */
+
+#define MAXPOW 32
+/* max power of 10 we wish to search to */
+
+#define MAXIT 800
+/* max number of iterations */
+
+/* a coefficient smaller than SMALL_ENOUGH is considered to
+ be zero (0.0). */
+
+#define SMALL_ENOUGH 1.0e-12
+
+
+/*
+ * structure type for representing a polynomial
+ */
+typedef struct p {
+ int ord;
+ double coef[MAX_ORDER+1];
+} poly;
+
+int modrf();
+int numroots(int np, poly * sseq, int * atneg, int * atpos);
+int numchanges(int np, poly * sseq, double a);
+int buildsturm(int ord, poly *sseq);
+
+double evalpoly();
+
+void sbisect(int np, poly * sseq, double min, double max, int atmin, int
atmax,double * roots);
+
Index: src/sturm.c
===================================================================
RCS file: src/sturm.c
diff -N src/sturm.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ src/sturm.c 22 Apr 2009 23:23:12 -0000 1.1
@@ -0,0 +1,287 @@
+
+/*
+ * sturm.c
+ *
+ * the functions to build and evaluate the Sturm sequence
+ */
+#include <math.h>
+#include <stdio.h>
+#include "solve.h"
+
+/*
+ * modp
+ *
+ * calculates the modulus of u(x) / v(x) leaving it in r, it
+ * returns 0 if r(x) is a constant.
+ * note: this function assumes the leading coefficient of v
+ * is 1 or -1
+ */
+static int
+modp(u, v, r)
+ poly *u, *v, *r;
+{
+ int k, j;
+ double *nr, *end, *uc;
+
+ nr = r->coef;
+ end = &u->coef[u->ord];
+
+ uc = u->coef;
+ while (uc <= end)
+ *nr++ = *uc++;
+
+ if (v->coef[v->ord] < 0.0) {
+
+
+ for (k = u->ord - v->ord - 1; k >= 0; k -= 2)
+ r->coef[k] = -r->coef[k];
+
+ for (k = u->ord - v->ord; k >= 0; k--)
+ for (j = v->ord + k - 1; j >= k; j--)
+ r->coef[j] = -r->coef[j] -
r->coef[v->ord + k]
+ * v->coef[j - k];
+ } else {
+ for (k = u->ord - v->ord; k >= 0; k--)
+ for (j = v->ord + k - 1; j >= k; j--)
+ r->coef[j] -= r->coef[v->ord + k] * v->coef[j -
k];
+ }
+
+ k = v->ord - 1;
+ while (k >= 0 && fabs(r->coef[k]) < SMALL_ENOUGH) {
+ r->coef[k] = 0.0;
+ k--;
+ }
+
+ r->ord = (k < 0) ? 0 : k;
+
+ return(r->ord);
+}
+
+/*
+ * buildsturm
+ *
+ * build up a sturm sequence for a polynomial in smat, returning
+ * the number of polynomials in the sequence
+ */
+int
+buildsturm(ord, sseq)
+ int ord;
+ poly *sseq;
+{
+ int i;
+ double f, *fp, *fc;
+ poly *sp;
+
+ sseq[0].ord = ord;
+ sseq[1].ord = ord - 1;
+
+
+ /*
+ * calculate the derivative and normalise the leading
+ * coefficient.
+ */
+ f = fabs(sseq[0].coef[ord] * ord);
+ fp = sseq[1].coef;
+ fc = sseq[0].coef + 1;
+ for (i = 1; i <= ord; i++)
+ *fp++ = *fc++ * i / f;
+
+ /*
+ * construct the rest of the Sturm sequence
+ */
+ for (sp = sseq + 2; modp(sp - 2, sp - 1, sp); sp++) {
+
+ /*
+ * reverse the sign and normalise
+ */
+ f = -fabs(sp->coef[sp->ord]);
+ for (fp = &sp->coef[sp->ord]; fp >= sp->coef; fp--)
+ *fp /= f;
+ }
+
+ sp->coef[0] = -sp->coef[0]; /* reverse the sign */
+
+ return(sp - sseq);
+}
+
+/*
+ * numroots
+ *
+ * return the number of distinct real roots of the polynomial
+ * described in sseq.
+ */
+int
+numroots(np, sseq, atneg, atpos)
+ int np;
+ poly *sseq;
+ int *atneg, *atpos;
+{
+ int atposinf, atneginf;
+ poly *s;
+ double f, lf;
+
+ atposinf = atneginf = 0;
+
+
+ /*
+ * changes at positive infinity
+ */
+ lf = sseq[0].coef[sseq[0].ord];
+
+ for (s = sseq + 1; s <= sseq + np; s++) {
+ f = s->coef[s->ord];
+ if (lf == 0.0 || lf * f < 0)
+ atposinf++;
+ lf = f;
+ }
+
+ /*
+ * changes at negative infinity
+ */
+ if (sseq[0].ord & 1)
+ lf = -sseq[0].coef[sseq[0].ord];
+ else
+ lf = sseq[0].coef[sseq[0].ord];
+
+ for (s = sseq + 1; s <= sseq + np; s++) {
+ if (s->ord & 1)
+ f = -s->coef[s->ord];
+ else
+ f = s->coef[s->ord];
+ if (lf == 0.0 || lf * f < 0)
+ atneginf++;
+ lf = f;
+ }
+
+ *atneg = atneginf;
+ *atpos = atposinf;
+
+ return(atneginf - atposinf);
+}
+
+/*
+ * numchanges
+ *
+ * return the number of sign changes in the Sturm sequence in
+ * sseq at the value a.
+ */
+int
+numchanges(np, sseq, a)
+ int np;
+ poly *sseq;
+ double a;
+
+{
+ int changes;
+ double f, lf;
+ poly *s;
+
+ changes = 0;
+
+ lf = evalpoly(sseq[0].ord, sseq[0].coef, a);
+
+ for (s = sseq + 1; s <= sseq + np; s++) {
+ f = evalpoly(s->ord, s->coef, a);
+ if (lf == 0.0 || lf * f < 0)
+ changes++;
+ lf = f;
+ }
+
+ return(changes);
+}
+
+/*
+ * sbisect
+ *
+ * uses a bisection based on the sturm sequence for the polynomial
+ * described in sseq to isolate intervals in which roots occur,
+ * the roots are returned in the roots array in order of magnitude.
+ */
+sbisect(np, sseq, min, max, atmin, atmax, roots)
+ int np;
+ poly *sseq;
+ double min, max;
+ int atmin, atmax;
+ double *roots;
+{
+ double mid;
+ int n1 = 0, n2 = 0, its, atmid, nroot;
+
+ if ((nroot = atmin - atmax) == 1) {
+
+ /*
+ * first try a less expensive technique.
+ */
+ if (modrf(sseq->ord, sseq->coef, min, max, &roots[0]))
+ return;
+
+
+ /*
+ * if we get here we have to evaluate the root the hard
+ * way by using the Sturm sequence.
+ */
+ for (its = 0; its < MAXIT; its++) {
+ mid = (min + max) / 2;
+
+ atmid = numchanges(np, sseq, mid);
+
+ if (fabs(mid) > RELERROR) {
+ if (fabs((max - min) / mid) < RELERROR)
{
+ roots[0] = mid;
+ return;
+ }
+ } else if (fabs(max - min) < RELERROR) {
+ roots[0] = mid;
+ return;
+ }
+
+ if ((atmin - atmid) == 0)
+ min = mid;
+ else
+ max = mid;
+ }
+
+ if (its == MAXIT) {
+ fprintf(stderr, "sbisect: overflow min %f max
%f\
+ diff %e nroot %d n1 %d n2 %d\n",
+ min, max, max - min, nroot, n1, n2);
+ roots[0] = mid;
+ }
+
+ return;
+ }
+
+ /*
+ * more than one root in the interval, we have to bisect...
+ */
+ for (its = 0; its < MAXIT; its++) {
+
+ mid = (min + max) / 2;
+
+ atmid = numchanges(np, sseq, mid);
+
+ n1 = atmin - atmid;
+ n2 = atmid - atmax;
+
+
+ if (n1 != 0 && n2 != 0) {
+ sbisect(np, sseq, min, mid, atmin, atmid,
roots);
+ sbisect(np, sseq, mid, max, atmid, atmax,
&roots[n1]);
+ break;
+ }
+
+ if (n1 == 0)
+ min = mid;
+ else
+ max = mid;
+ }
+
+ if (its == MAXIT) {
+ fprintf(stderr, "sbisect: roots too close together\n");
+ fprintf(stderr, "sbisect: overflow min %f max %f diff
%e\
+ nroot %d n1 %d n2 %d\n",
+ min, max, max - min, nroot, n1, n2);
+ for (n1 = atmax; n1 < atmin; n1++)
+ roots[n1 - atmax] = mid;
+ }
+}
Index: src/util.c
===================================================================
RCS file: src/util.c
diff -N src/util.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ src/util.c 22 Apr 2009 23:23:12 -0000 1.1
@@ -0,0 +1,119 @@
+
+/*
+ * util.c
+ *
+ * some utlity functions for root polishing and evaluating
+ * polynomials.
+ */
+#include <math.h>
+#include <stdio.h>
+#include "solve.h"
+
+/*
+ * evalpoly
+ *
+ * evaluate polynomial defined in coef returning its value.
+ */
+double
+evalpoly (ord, coef, x)
+ int ord;
+ double *coef, x;
+{
+ double *fp, f;
+
+ fp = &coef[ord];
+ f = *fp;
+
+ for (fp--; fp >= coef; fp--)
+ f = x * f + *fp;
+
+ return(f);
+}
+
+
+/*
+ * modrf
+ *
+ * uses the modified regula-falsi method to evaluate the root
+ * in interval [a,b] of the polynomial described in coef. The
+ * root is returned is returned in *val. The routine returns zero
+ * if it can't converge.
+ */
+int
+modrf(ord, coef, a, b, val)
+ int ord;
+ double *coef;
+ double a, b, *val;
+{
+ int its;
+ double fa, fb, x, fx, lfx;
+ double *fp, *scoef, *ecoef;
+
+ scoef = coef;
+ ecoef = &coef[ord];
+
+ fb = fa = *ecoef;
+ for (fp = ecoef - 1; fp >= scoef; fp--) {
+ fa = a * fa + *fp;
+ fb = b * fb + *fp;
+ }
+
+ /*
+ * if there is no sign difference the method won't work
+ */
+ if (fa * fb > 0.0)
+ return(0);
+
+ if (fabs(fa) < RELERROR) {
+ *val = a;
+ return(1);
+ }
+
+ if (fabs(fb) < RELERROR) {
+ *val = b;
+ return(1);
+ }
+
+ lfx = fa;
+
+
+ for (its = 0; its < MAXIT; its++) {
+
+ x = (fb * a - fa * b) / (fb - fa);
+
+ fx = *ecoef;
+ for (fp = ecoef - 1; fp >= scoef; fp--)
+ fx = x * fx + *fp;
+
+ if (fabs(x) > RELERROR) {
+ if (fabs(fx / x) < RELERROR) {
+ *val = x;
+ return(1);
+ }
+ } else if (fabs(fx) < RELERROR) {
+ *val = x;
+ return(1);
+ }
+
+ if ((fa * fx) < 0) {
+ b = x;
+ fb = fx;
+ if ((lfx * fx) > 0)
+ fa /= 2;
+ } else {
+ a = x;
+ fa = fx;
+ if ((lfx * fx) > 0)
+ fb /= 2;
+ }
+
+ lfx = fx;
+ }
+
+ // fprintf(stderr, "modrf overflow %f %f %f\n", a, b, fx);
+
+ return(0);
+}
+
+
+
Index: src/ccmcomplex.h
===================================================================
RCS file: src/ccmcomplex.h
diff -N src/ccmcomplex.h
--- src/ccmcomplex.h 22 Apr 2009 17:05:25 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,12 +0,0 @@
-/* complex.h CCMATH mathematics library source code.
- *
- * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
- * This code may be redistributed under the terms of the GNU library
- * public license (LGPL). ( See the lgpl.license file for details.)
- * ------------------------------------------------------------------------
- */
-#ifndef CPX
-struct complex {double re,im;};
-typedef struct complex Cpx;
-#define CPX 1
-#endif
Index: src/plrt.c
===================================================================
RCS file: src/plrt.c
diff -N src/plrt.c
--- src/plrt.c 22 Apr 2009 17:05:25 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,47 +0,0 @@
-/* plrt.c CCMATH mathematics library source code.
- *
- * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
- * This code may be redistributed under the terms of the GNU library
- * public license (LGPL). ( See the lgpl.license file for details.)
- * ------------------------------------------------------------------------
- */
-#include <stdlib.h>
-#include <math.h>
-#include "ccmcomplex.h"
-int plrt(double *cof,int n,Cpx *root,double ra,double rb)
-{ double a,b,s,t,w; int itr,pat; struct complex *pr;
- double *cs,*cf,*hf,*p,*q,*ul,test=1.e-28;
- cs=cf=(double *)calloc(2*n,sizeof(double)); hf=cf+n;
- pr=root+n-1; ul=hf+n-1;
- if(rb<0.) rb=ra*ra-rb*rb; else rb=ra*ra+rb*rb; ra*= -2.;
- q=cof+n; s= *q--;
- for(p=cf; p<hf ;) *p++ = *q-- /s;
- for(itr=pat=0;;){
- if(itr==0){
- if(n>2){ a=ra; b=rb;}
- else if(n==2){ a= *cf++; b= *cf;}
- else if(n==1){ pr->re= -(*cf); pr->im=0.; free(cs); return 0;}
- }
- s= -a/2.; t=s*s-b;
- if(t>=0.){ t=sqrt(t); pr->re=s+t; (pr--)->im=0.;
- pr->re=s-t; (pr++)->im=0.; }
- else{ t=sqrt(-t); pr->re=s; (pr--)->im=t;
- pr->re=s; (pr++)->im= -t; }
- if(n==2){ free(cs); return 0;}
- for(p=hf,q=cf; q<hf ;) *p++ = *q++;
- for(p=hf,w=1.; p<ul ;){
- *p++ -=a*w; *p-=b*w; w=*(p-1); }
- t= -(*p--); t+=pr->re* *p; s=pr->im* *p;
- if(t*t+s*s<test){ pr-=2; ul-=2; n-=2; itr=pat=0;
- for(p=cf,q=hf; p<ul ;) *p++ = *q++;
- }
- else if(++itr<30){
- for(p=hf,w=1.; p<ul-2 ;){
- *p++ -=w*a; *p-=w*b; w= *(p-1); }
- t= *p++; q=p+1; s=t*t+w*(b*w-a*t);
- b+=(w*(*p*b- *q*a)+ *q*t)/s; a+=(*p*t- *q*w)/s;
- }
- else{ if(pat==3){ free(cs); return n;}
- itr=0; if(pat++ %2) ra= -ra; else rb= -rb;}
- }
-}