Optimized integrate function in module poly
Rewrote the integrate function since the old one was quite hacky. The new one is about 7 times faster in release and 3 times faster in debug
This commit is contained in:
parent
ce1399db2d
commit
0d2f8d5079
1 changed files with 16 additions and 10 deletions
|
|
@ -11,9 +11,6 @@ import math
|
||||||
import strutils
|
import strutils
|
||||||
import numeric
|
import numeric
|
||||||
|
|
||||||
|
|
||||||
import times #todo: remove
|
|
||||||
|
|
||||||
type
|
type
|
||||||
TPoly* = object
|
TPoly* = object
|
||||||
cofs:seq[float]
|
cofs:seq[float]
|
||||||
|
|
@ -141,10 +138,21 @@ proc integral*(p:TPoly):TPoly=
|
||||||
|
|
||||||
proc integrate*(p:TPoly;xmin,xmax:float):float=
|
proc integrate*(p:TPoly;xmin,xmax:float):float=
|
||||||
## Computes the definite integral of `p` between `xmin` and `xmax`
|
## Computes the definite integral of `p` between `xmin` and `xmax`
|
||||||
# TODO: this can be done faster using a modified horners method,
|
## quickly using a modified version of Horners method
|
||||||
# see 'diff' function above.
|
var
|
||||||
var igr=p.integral
|
n=p.degree
|
||||||
result=igr.eval(xmax)-igr.eval(xmin)
|
s1=p[n]/float(n+1)
|
||||||
|
s2=s1
|
||||||
|
fac:float
|
||||||
|
|
||||||
|
dec n
|
||||||
|
while n>=0:
|
||||||
|
fac=p[n]/float(n+1)
|
||||||
|
s1 = s1*xmin+fac
|
||||||
|
s2 = s2*xmax+fac
|
||||||
|
dec n
|
||||||
|
|
||||||
|
result=s2*xmax-s1*xmin
|
||||||
|
|
||||||
proc initPoly*(cofs:varargs[float]):TPoly=
|
proc initPoly*(cofs:varargs[float]):TPoly=
|
||||||
## Initializes a polynomial with given coefficients.
|
## Initializes a polynomial with given coefficients.
|
||||||
|
|
@ -258,7 +266,7 @@ proc `/` *(p,q:TPoly):TPoly=
|
||||||
|
|
||||||
proc `mod` *(p,q:TPoly):TPoly=
|
proc `mod` *(p,q:TPoly):TPoly=
|
||||||
## Computes the polynomial modulo operation,
|
## Computes the polynomial modulo operation,
|
||||||
## that is the remainder op `p`/`q`
|
## that is the remainder of `p`/`q`
|
||||||
var dummy:TPoly
|
var dummy:TPoly
|
||||||
p.divMod(q,dummy,result)
|
p.divMod(q,dummy,result)
|
||||||
|
|
||||||
|
|
@ -363,5 +371,3 @@ proc roots*(p:TPoly,tol=1.0e-9,zerotol=1.0e-6,mergetol=1.0e-12,maxiter=1000):seq
|
||||||
addRoot(p,result,range.xmin,x,tol,zerotol,mergetol,maxiter)
|
addRoot(p,result,range.xmin,x,tol,zerotol,mergetol,maxiter)
|
||||||
range.xmin=x
|
range.xmin=x
|
||||||
addRoot(p,result,range.xmin,range.xmax,tol,zerotol,mergetol,maxiter)
|
addRoot(p,result,range.xmin,range.xmax,tol,zerotol,mergetol,maxiter)
|
||||||
|
|
||||||
|
|
||||||
Loading…
Add table
Add a link
Reference in a new issue