Try the code below--it is very fast and should give you enough precision as long as you know it to enough precision for convergence the root you're looking for is not to large.
- Robert On Sep 4, 2007, at 9:35 AM, John Voight wrote: > > Hello Carl, > >> As Robert Bradshaw mentioned, I did spend a lot of time working on >> real root isolation; unfortunately, I have not yet found the time to >> polish it off and contribute it to SAGE (although I'm still planning >> to!). >> >> What do you need for real root finding for your problem? In >> particular, are your polynomials known exactly (integral or rational >> coefficients)? What degree are your polynomials? How big are the >> coefficients? Are your polynomials known to be squarefree (no >> multiple roots)? What accuracy do you need on the resulting roots? > > In my situation, I have the following absurdly easy case: I have a > real interval in which I know there is exactly one real root which I > need to know to maybe 6 or 10 digits of precision. The polynomials > are monic, have small integer coefficients (bounded in absolute value > by maybe 100 or so), and have small degree (<= 11). The problem is > that there are zillions of them--so I need this data very quickly! > > Any advice would be most appreciated. > > Yours, > > John Voight > Assistant Professor of Mathematics > University of Vermont > [EMAIL PROTECTED] > [EMAIL PROTECTED] > http://www.cems.uvm.edu/~voight/ {{{id=5| %cython cdef double_newton_raphson_c(double* coeffs, double x, Py_ssize_t n, int iter): cdef Py_ssize_t i, j cdef double eval_f, eval_df for j from 0 <= j < iter: eval_f = coeffs[n] for i from n > i >= 0: eval_f = eval_f*x + coeffs[i] eval_df = coeffs[n]*n for i from n > i > 0: eval_df = eval_df*x + coeffs[i]*i x -= eval_f/eval_df return x def double_newton_raphson(f, x, iter=8): f = list(f) cdef Py_ssize_t i, n = len(f)-1 cdef double* coeffs = <double*>sage_malloc(sizeof(double)*(n+1)) if coeffs == NULL: raise MemoryError for i from 0 <= i <= n: coeffs[i] = f[i] root = double_newton_raphson_c(coeffs, x, n, iter) sage_free(coeffs) return root }}} {{{id=11| x = ZZ['x'].0 }}} {{{id=29| f = prod([x-i for i in range(1,7)]); f /// x^6 - 21*x^5 + 175*x^4 - 735*x^3 + 1624*x^2 - 1764*x + 720 }}} {{{id=12| double_newton_raphson(f, 6.2) /// 5.9999999999999227 }}} {{{id=43| %time f = list(f) for _ in range(10^5): z = double_newton_raphson(f, 6.2r) /// CPU time: 0.25 s, Wall time: 0.26 s }}} --~--~---------~--~----~------------~-------~--~----~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~----------~----~----~----~------~----~------~--~---
