[100]022 Chapter 22 User Function, Tutorial by
[100]022 revised on 2012.10.28 / cemmathThe Simple is the Best
Chapter 22 User Function
//------
// grammar of function
//------
Caution : double semi-colons (;;) activate screen-output inside a function, whereas a semi-colon (;) suppresses screen-output.
If not specified explicitly, a function argument is regarded as 'double'. In this regard, the following declarations
double f(x) = x ;
double f(double x) = x ;
double f(double x) { return x; } // standard C grammar
are exactly all the same in Cemmath. A use of '=' is just for simple-expression function. Other functions are almost the same as in C-language.
//------
// a simple example, double ftn(x) = double ftn(double x)
//------
#> double fn(x) = |x|*sin(x) ;
#> fn ;
double fn(x)
#> fn(3) ;
ans = 0.42336002
#> fn'(3) ; // derivative for double(double) function
ans = -2.8288575
#> fn''(3) ; // 2nd-order derivative for double(double) function
ans = -2.4033449
#> fn~(0,3); // integration over (0,3) for double(double) function
ans = 3.1110975
#> fn++( [1,2,3] ) ; // upgrade for matrix
ans = [ 0.841471 1.81859 0.42336 ]
A double function can be upgraded to accept vector argument by attaching the postfix '++'. This upgrade alleviates the use of do-loop.
//------
// return type
//------
Several types including 'void' (no type) can be returned by a function. An example is
matrix test(double n) { // matrix test(n) assumes n is also a matrix
(a,b,c) = (0,0,0);
switch(n) {
case 1: a = 10; break; // case must be followed by an integer only
default: c = 5;
case 2: b = 20; break;
case 3: c = 100;
}
return [a,b,c];
}
#> test(1);
ans = [ 10 0 0 ]
#> test(2);
ans = [ 0 20 0 ]
#> test(3);
ans = [ 0 0 100 ]
#> test(4);
ans = [ 0 20 5 ]
//------
// function upgrade, matrix mftn(double x)
//------
#> matrix ut(double x) = [ x; 1; x^2 ];
#> ut( 3 );
ans =
[ 3 ]
[ 1 ]
[ 9 ]
#> ut++( [ 3,4,5 ] );
ans =
[ 3 4 5 ]
[ 1 1 1 ]
[ 9 16 25 ]
#> ut++( [ 3,4; 5,6 ] );
ans =
[ 3 4 ]
[ 5 6 ]
[ 1 1 ]
[ 1 1 ]
[ 9 16 ]
[ 25 36 ]
//------
// piecewise functions, matrix mftn(double x)
//------
matrix g(double x) = // piecewise functions
[0,1, x ; // 0 <= x <= 1, y = x
1,2, x^2-1; // 1 <= x <= 2, y = x^2 - 1
2,3, x^3-8; // 2 <= x <= 3, y = x^3 - 8
3,4, x^4-81 ]; // 3 <= x <= 4, y = x^4 - 81
#> g.spl( 1.5 ); // evaluate at x
#> g.spl([0.5,1.5,2.5,3.5]); // evaluate at multipoints [x]
ans = 1.25
ans = [ 0.5 1.25 7.625 69.0625 ]
It is also possible to define a typical 'double(double)' function utilizing 'if-else' clause as follows.
double gg(double x) {
if ( 0 <= x & x <= 1 ) return x;
else if( 1 <= x & x <= 2 ) return x^2-1;
else if( 2 <= x & x <= 3 ) return x^3-8;
else if( 3 <= x & x <= 4 ) return x^4-81;
}
#> gg( 1.5 ); // evaluate at x
#> gg++([0.5,1.5,2.5,3.5]); // evaluate at multipoints [x]
ans = 1.25
ans = [ 0.5 1.25 7.625 69.0625 ]
//------
// copying argument
//------
Arguments are by default copied inside of a function, and thus original variables are unchanged.
void swap1(x,y) { t = x; x = y; y = t; }
#> x = 3; y = 5;
x = 3
y = 5
#> swap1(x,y); x; y;
x = 3
y = 5
Nothing happens to x and y.
//------
// reference
//------
A variable can be refered by a prefix '&', and the original variables are delivered inside a function.
void swap2(&x,&y) { t = x; x = y; y = t; }
#> x = 3; y = 5;
x = 3
y = 5
#> swap2(x,y); x; y;
x = 5
y = 3
Note that x and y are changed.
//------
// array argument
//------
A prefix '*' is attached to array arguments and data delivery is by reference.
void array(matrix *A){
A[0] = .zeros(3);
A[2] = [ 1,2,3,4; 5,6,7,8 ];
}
#> matrix An[4]; array(An); An;
An = matrix [4]
[0] =
[ 0 0 0 ]
[ 0 0 0 ]
[ 0 0 0 ]
[1] =
[ undefined ]
[2] =
[ 1 2 3 4 ]
[ 5 6 7 8 ]
[3] =
[ undefined ]
//------
// repeating types
//------
When the same data types are repeated, only the leading variable can be declared such that
(double a,b,c, vertex u,v)
which is the same as
(double a, double b, double c, vertex u, vertex v)
%> an example
complex leading(double a,b,c, vertex u,v) {
z = a+b+c + (u.x*v.x + u.y*v.y)! ;
return z;
}
#> leading(3,4,5, <1,2,-1>,<4,-1,5>);
ans = 12 + 2!
//------
// recurrence function
//------
A function can call itself inside its definition, and this is the case of recurrence function
%> factorial function
double fac(n) {
if( n == 1 ) return 1;
return n*fac(n-1);
}
#> fac(5);
ans = 120
//------
// dynamic assignment
//------
Arguments of a function can be expressed by the arguments defined before.
double area( a, b=a ) = a*b; // b is expressed by a
#> area ; double area( a, b=a )
#> area(3);
ans = 9
#> area(3,4);
ans = 12
//------
// function argument
//------
A special grammar is innovated for a function argument as follows.
double ffx(x) = x^2;
double ggx(x) = x*sin(x);
double pftn(f(x),x) = f(x)^2; // function argument f(x)
#> pftn(ffx,3); pftn(ggx,3);
ans = 81
ans = 0.17923371
//------
// auto tracing
//------
A history of calling functions is traced automatically.
void f1(&x1) { x1 = sin(x1);; z1=1;; }
void f2(&x2) { x2 = cos(x2);; f1(x2);; z2=2;; }
void f3(&x3) { x3 = exp(x3);; f2(x3);; z3=3;; }
void f4(&x4) { x4 = x4^2;; f3(x4);; z4=4;; }
#> y=3; f4(y);
y = 3
f4.y = 9
f4.f3.y = 8103.0839
f4.f3.f2.y = -0.6086217
f4.f3.f2.f1.y = -0.5717372
f4.f3.f2.f1.z1 = 1
f4.f3.f2.z2 = 2
f4.f3.z3 = 3
f4.z4 = 4
%> case of no reference
void g1(x1) { x1 = sin(x1);; z1=1;; }
void g2(x2) { x2 = cos(x2);; g1(x2);; z2=2;; }
void g3(x3) { x3 = exp(x3);; g2(x3);; z3=3;; }
void g4(x4) { x4 = x4^2;; g3(x4);; z4=4;; }
#> y=3; g4(y);
y = 3
g4.x4 = 9
g4.g3.x3 = 8103.0839
g4.g3.g2.x2 = -0.6086217
g4.g3.g2.g1.x1 = -0.5717372
g4.g3.g2.g1.z1 = 1
g4.g3.g2.z2 = 2
g4.g3.z3 = 3
g4.z4 = 4
//------
// Compatibility with C-Language
//------
Consider the following two subroutines developed in C language.
double x_fun(double);
double x_fun(double x) { return exp(x)-sin(x)-2; }
double x_bisect(double (*fun)(double),double,double);
double x_bisect(double (*fun)(double),double xa,double xb ) {
double iter,eps,fa,fb,x,f;
iter = 0; eps = 1.e-7;
fa = fun(xa); if( fabs(fa) < eps ) return xa;
fb = fun(xb); if( fabs(fb) < eps ) return xb;
if( fa*fb > 0 ) { return 0; }
while( ++iter < 200 ) {
x = 0.5*(xa+xb);
f = fun(x);
if( fabs(f) < eps) return x;
if( fa*f > 0 ) { xa = x; fa = f; }
else { xb = x; fb = f; }
}
return 0;
}
It is easy to find that the role of these two subroutines is to find a root by the bisection method in numerical analysis. The grammars with the function pointer, prototypes are all implemented in Cemmath too. Therefore, the above two subroutines can be compiled in Cemmath and thus the following command
#> bisect_x3 = x_bisect(x_fun,0,4);
results in
bisect_x3 = 1.0541271
1