主题:[讨论]任意函数表达式求值模块
weixing1531 [专家分:2580] 发布于 2008-05-05 19:03:00
不久之前,我发过一张贴——《我见过最牛B的程序》。讲的是函数表达式求值。很可惜是C
语言描述。
今天我特地贴一张fortran 版的!
这个模块是 www.fortran.com 这个网站 Dr Stuart Midgley 写的。
他的电子邮箱:sdm900@gmail.com,大家如有疑问可以找他(该老外有问必答,我试过)。
废话不多说,程序见附件。
program test
use params
use extrafunc
use interpreter
implicit none
character(len=255) :: s, f,error1,error2
real(fdp) :: r
real(fdp), dimension(:),allocatable :: varvals
integer :: funcnum
write(*,*)"依次输入函数的自变量:"
read(*,'(a)')s
allocate(varvals(f_numbervars(s))) ! f_numbervars(s)为自变量的个数
write(*,*)"输入关于自变量的函数:"
read(*,'(a)')f
write(*,*)"依次输入对应自变量的值:"
read(*,*)varvals
call s_createfn(f, s, funcnum,error1) ! 创建函数
write(*,'(a)')trim(error1) ! 提示出错信息
call s_evaluatefn(funcnum, varvals, r,error2) ! 计算函数值
write(*,'(a)')trim(error2) ! 提示出错信息
write(*,*)"函数值等于:",r
call s_destroyfn(funcnum) ! 销毁函数
stop
end program test
应用举例:
输入函数的自变量:
x
输入关于自变量的函数:
x**2+sqrt(x)/x
输入自变量的值:
2.0
输出:
OK
OK
函数值等于: 4.7071067811865475
另外模块还提供别外两个接口:s_errprecision 和 s_evaluateerr
前者好象是把实数的左数第二位起四舍五入(我个人理解,不常用)
后者如果自变量输入:x[2.0] y[3.0] 则计算二元函数 f(x,y) 在△x=2.0 △y=3.0 时的全微分。
此外模块还提供了数学上常用的函数如gamma函数等,详见 module extrafunc,输入表达式过程
中,**和^同为乘方,支持小、中、大括号,其中if函数非常类似于C语言中的三元操作符。
我个人认为,该模块结合dislin可以画出任意函数的图像,是不是有点matlab的意味?
最后更新于:2008-05-05 22:45:00
回复列表 (共11个回复)
沙发
weixing1531 [专家分:2580] 发布于 2008-05-05 20:07:00
*********************************************************************
This software is provided AS IS with NO WARRANTY. I also only provide
it for "not for profit" use. If you wish to use the software for any
other purpose, please contact me on sdm900@gmail.com .
*********************************************************************
Hi Fortran Programmers
OK, the function interpreter is now finished (at least to a level that
satisfies me).
PLEASE, this code comes AS IS without any warranty of any kind. I have
developed this code enough to satisfy my own needs and no body elses.
If you find a use for it, well and good. If you manage to improve on
it, please me know and I will put the improvements in my own code.
The interpreter does 57 functions
cos(x)
sin(x)
tan(x)
exp(x)
log(x)
log10(x)
logn(x,y)
sqrt(x)
cbrt(x)
acos(x)
asin(x)
atan(x)
cosh(x)
sinh(x)
tanh(x)
anint(x)
aint(x)
abs(x)
delta(x)
step(x)
hat(x)
min(x,y)
max(x,y)
besj0(x)
besj1(x)
besjn(n,x)
besy0(x)
besy1(x)
besyn(n,x)
besi0(x)
besi1(x)
besin(n,x)
besk0(x)
besk1(x)
beskn(n,x)
erf(x)
erfc(x)
ierf(x)
ierfc(x)
gamma(x)
lgamma(x)
csch(x)
sech(x)
coth(x)
if(conditional, true, false)
gauss(x)
sinc(x)
fresc(x)
fress(x)
expi(x)
sini(x)
cosi(x)
logi(x)
elle(x)
ellk(x)
ielle(x, phi)
iellf(x, phi)
15 operators
+
-
*
/
**
^
>
>=
=>
<
<=
=<
==
=
!=
brackets and correct order of operation.
Numbers are all assumed to be double precision real and are entered in
the usual fortran double precision way (or any way that read(*,*) num
will interpret)
# = integer
#
.#
#.#
#.#d#
#.#d-#
#.#d+#
#.#e#
#.#e-#
#.#e+#
#.#D#
#.#D-#
#.#D+#
#.#E#
#.#E-#
#.#E+#
etc...
The functions are simply strings with any number of possible variables
**** THE CODE NOW HAS ERROR DETECTION ****
Yes, I got around to adding it at the request of a user of my code.
Basically, the code will try to detect when an illegal function is
entered. It is quite basic, but seems to work quite well.
Basically, when calling the function, you have the option of
including an error field ( character(len=#)::error ). If no error
occurs, this error field will contain OK. If an error is detected
then error will contain the first # chacters of the error explanation.
It will still return a valid number (that obtained previously) but you
can at least have some primitive idea of what went wrong.
********************************************
**** THE CODE NOW DOES ERROR LIMIT COMPUTATIONS ****
The code now allows the computation of error calculations of the form
D = delta (error)
d = partial derivative
Df = df/dx Dx + df/dy Dy
You input the variables using a square bracket
x[0.1]
would give an error of 0.1 on the variable x.
******************************************************
Anyway, the code isn't documented, though it should be really easy to
follow.
KNOWN PROBLEMS
No known problems.
If you have any problems, let me know.
Stu.
http://stu.ods.org
sdm900@gmail.com
板凳
weixing1531 [专家分:2580] 发布于 2008-05-05 20:08:00
The code has been broken up into fairly simple pieces, which should allow
great flexibility and speed.
subroutine s_createfn
subroutine s_evaluatefn
subroutine s_evaluateerr
subroutine s_destroyfn
The Fortran90 Interface blocks would look like
interface
subroutine s_createfn(func, vars, funcnumber, error)
character(len = *), intent(in) :: func, vars
integer, intent(out) :: funcnumber
character(len=*), intent(out), optional :: error
end subroutine
subroutine s_evaluatefn(funcnumber, varvals, number, error)
integer, intent(in) :: funcnumber
real(fdp), dimension(:), intent(in) :: varvals
real(fdp), intent(out) :: number
character(len=*), intent(out), optional :: error
end subroutine
subroutine s_evaluateerr(funcnumber, varvals, number, error)
integer, intent(in) :: funcnumber
real(fdp), dimension(:), intent(in) :: varvals
real(fdp), intent(out) :: number
character(len=*), intent(out), optional :: error
end subroutine
subroutine s_destroyfn(funcnumber)
integer, intent(in) :: funcnumber
end subroutine
end interface
where
fdp = 8 Double precision float
So, as an example you could do
read(*,'(a)') f1
call s_createfn(f1, 'x[0.1] y[0.01] z[0.02]', funcnum1, error)
write(*,*) trim(error)
call s_evaluatefn(funcnum1, (/1.0d0, 2.0d0, 3.0d0/), r, error)
write(*,*) r
call s_evaluateerr(funcnum1, (/1.0d0, 2.0d0, 3.0d0/), e, error)
write(*,*) e
call s_destroyfn(funcnum1)
Once you have created a function and have a function number for it, that function will remained stored in memory in an optimised fashon until you destroy it or exit the program.
The basic idea is that you load in a function from some input (STDIN or text file) and call s_createfn which basically parses the text function into an internal representation for fast execution. Then, whenever you need that function computed for various values of the input variables, you call s_evaluatefn. Once you have finished, you call s_destroyfn and exit the program.
3 楼
f2003 [专家分:7960] 发布于 2008-05-05 21:17:00
楼主练练手,把它打磨打磨,弄个界面,搞个函数计算器怎么样?
4 楼
weixing1531 [专家分:2580] 发布于 2008-05-06 21:55:00
源程序下载链接:
http://stu.ods.org/fortran/parser_20070302.zip
为了能在CVF6.6c上运行,得改一个地方
在module params中
把 sqrt2 = sqrt(2.0d0) 改为 sqrt2 = 1.4142135623730950488d0
如果是 g95 或 gfortran 的话,就不用改
另外为了方便得到自变量的个数,我把 f_numbervars 设为公有了
5 楼
weixing1531 [专家分:2580] 发布于 2008-06-04 00:41:00
[quote]楼主练练手,把它打磨打磨,弄个界面,搞个函数计算器怎么样?[/quote]
这个模块的作者已经做了个画函数图像的界面
用Xeffort搞的
下载地址:
[url=http://www.xeffort.com/xeffort/samples/xgraph/xgraph.htm]http://www.xeffort.com/xeffort/samples/xgraph/xgraph.htm[/url]
6 楼
vqimwr [专家分:1320] 发布于 2008-06-04 09:21:00
谢楼上,好东西,还有源代码。
7 楼
臭石头雪球 [专家分:23030] 发布于 2008-08-26 12:20:00
唉,出了一次野外,这好帖子居然错过了....
一直很羡慕 Microsoft 的那个加强版 Calc.exe ,直接输入表达式求值,很好玩。
8 楼
jason388 [专家分:6150] 发布于 2008-08-26 12:53:00
谢谢分享!
c和c++表达式求值的项目有很多,但fortran版的还是第一次见到。受教了。
9 楼
weixing1531 [专家分:2580] 发布于 2008-09-07 21:21:00
应用再举一个很有用的例子
IMSL 库中有一个用二分法求函数根的函数 DZBREN
下面有一段程序引至彭国伦《FORTRAN 95》:
program main
use IMSL
implicit none
real(8), parameter :: ERRABS = 0.0
real(8), parameter :: ERRREL = 1e-14
integer :: MAXFN = 100
real(8) :: A,B ! 二分法求根区间
real(8), external :: F
A=0.4_8
B=0.5_8
call DZBREN (F, ERRABS, ERRREL, A, B, MAXFN)
write(*,*) B,MAXFN
end program
real(8) function F(X)
implicit none
real(8),intent(in)::X
F=tan(x)+x-1
return
end function
这种方法有一个很大的缺陷:函数F(X)必须在编译之前确定,即早绑定。如果你要计算另外一个F(X),
你不得不重新编译,真可谓“多次编译,多次运行”!
接下来介绍另外一种方法——运行时绑定:
program main
use IMSL
use params
use extrafunc
use interpreter
implicit none
real(8), parameter :: ERRABS = 0.0
real(8), parameter :: ERRREL = 1e-14
integer :: MAXFN = 100,funcnum
real(8) :: A,B
real(8), external :: F
character(len=255) ::hs,s='x',error1 ! 因为只一个自变量,故s='x'
common funcnum ! 与子程序共享数据
write(*,*)'输入关于自变量x的函数:'
read(*,'(a)')hs
write(*,*)'输入二分法求根区间:'
read(*,*)A,B
call s_createfn(hs, s, funcnum,error1) ! 创建函数
write(*,'(a)')trim(error1) ! 提示创建函数时的出错信息
call DZBREN (F, ERRABS, ERRREL, A, B, MAXFN)
write(*,*) B,MAXFN
call s_destroyfn(funcnum) ! 销毁函数
end program
real(8) function F(X)
use params
use extrafunc
use interpreter
implicit none
real(8),intent(in)::X
character(len=255) :: error2
integer :: funcnum
common funcnum ! 与主程序共享数据
call s_evaluatefn(funcnum, (/x/), F,error2) ! 计算函数值
write(*,'(a)')trim(error2) ! 提示计算函数值时的出错信息
return
end function
这种方法必须要用到前面给出下载链接的模块use params 、use extrafunc 、use interpreter
它能在程序运行时再确定函数,做到“一次编译,永久运用”!
是不是用点像 matlab ?!但效率应该比 matlab 要好一点,毕竟 matlab 是解释执行。
但这种方法和上一种方法比起来:灵活有余,效率不足!
10 楼
kaierme [专家分:970] 发布于 2009-09-17 22:50:00
感谢楼主,也感谢开发者,免费提供这样的代码
我一直有这样的想法,去做这样一个东西,都没能实现
我来回复