Menu
▾
▴
Re: finite summation - how to call a function from parameter stack
|
From: Ethan M. <merritt@u.washington.edu> - 2011-01-12 19:44:36
|
On Wednesday, January 12, 2011 07:49:57 am mw...@gm... wrote:
[snip]
> I have implemented sum(a, b, "f") and sum(a, b, x, "g") in order to calculate finite sums:
>
> # calculate a finite sum
> f(n) = n
> print sum(1, 10, "f") # -> 55
>
> # calculate a finite sum with parameter x
>
> # fourier coefficients of some square wave
> fourier(k, x) = sin(3./2*k)/k * 2./3*cos(k*x)
> s(x) = 1./2 + sumx(1, 100, x, "fourier")
> plot s(x)
>
> I hope this is sufficient to understand why I want a parameter x (so one can plot some square waves and other things.)
I find that syntax rather confusing, but maybe it's just me.
My first thought would be to aim for
fourier(x) = for [k=1:100] sin(3./2*k)/k * 2./3*cos(k*x)
since that matches the iteration syntax already used elsewhere
in the set of gnuplot commands. I suspect the internal implementation
could still work more or less the way you have it, but the user
interface would be (for me anyhow) more obvious.
> As Mr. Bröker pointed out, there is a trick involving the tertinary operator '?' and recursion to sorta workaround. Originally I tried to plot the sum with '?' but I could only do so for 1..10, in case of 1..100 it said something about stack overflow.
I am inclined to agree with Hans-Bernhard that such a summation
meta-function is not required.
But feel free to provide a contrary argument.
The summation can easily be expressed as a simple recursive function.
The problem you saw with stack overflow probably resulted from failing to
order the terms so that the calculation involved tail-recursion.
E.g.
gnuplot> tailsum(N) = (N>0) ? tailsum(N-1) + N : 0
gnuplot> headsum(N) = (N>0) ? N + headsum(N-1) : 0
gnuplot> print tailsum(250)
31375
gnuplot> print headsum(250)
stack overflow
Now it is true that the current version of gnuplot imposes a further
limit on the recursion depth:
gnuplot> print tailsum(251)
recursion depth limit exceeded
It is currently a hard limit determined by
eval.h:#define STACK_DEPTH 250 /* maximum size of the execution stack */
But it's largely an arbitrary number to prevent run-away function evaluation.
It could be made a user-accessible parameter if there's a need for it.
> I have added below a working patch but would like you to have a look at it.
I don't really have much time right now.
Could you please upload it to the patch tracker on Sourceforge?
That way it won't get lost, you can update it whenever you want,
and testers can post comments and provide feedback.
> I have some things that are not clear to me, in particular:
>
> * Do I have to free the values that I have gotten using pop?
No.
> * In what order do I have to push i and x befor the call to f(i, x). I thought that I have to push them from right to left but that seems to contradicts my testing, now I am very unsure about it (even though it seems to work correctly).
If you are calling your own function then you need only be
consistent at both ends of the call. But in general the stack
is LIFO. pop() returns the value most recently push()ed
> * Do I have to alloc memeory for the values I push on the stack befor the call to f(.)?
The push() operation makes a copy.
It's up to you to manage allocate/free of the original.
> * This is only curiosity is arg in foo(union argument *arg) every used and if, what for?
I have wondered that myself :-)
> I would also like you to take a look at my addition to gnuplot.doc and proofread it.
Please put it on the patch tracker along with the patch itself.
> Finally it would be nice but not necessary if the function could be passed as 'f' instead of '"f", that is sum(1, 10, f), and even better as sum(n=[1:10], n) but I consider this as optional and very hard to do (at least with my knowledge).
That would not be a good idea, as it would prevent the use of string variables.
Consider:
ident(N) = N
f = "ident"
print sum(1, 10, f)
> Thank you again for your generous help without I would never have gone so far.
If you have this working now, I'd be interested to see some benchmarks.
Is there a performance benefit over using the existing mechanism of
defining a recursive function using tail-recursion?
Is there a real-world use that requires summing over more than 250 terms?
Ethan
> Thank you,
> Micha Wiedenmann
>
> --- eval.c.orig 2011-01-11 17:38:00.000000000 +0100
> +++ eval.c 2011-01-12 16:19:05.000000000 +0100
> @@ -176,6 +176,9 @@
> {"atanh", f_atanh},
> {"lambertw", f_lambertw}, /* HBB, from G.Kuhnle 20001107 */
>
> + {"sum", f_sum}, /* summation sum(a, b, "f") = \sum_{i=a}^b f(i) */
> + {"sumx", f_sumx}, /* summation sum(a, b, x, "f") = \sum_{i=a}^b f(i, x) (with parameter x) */
> +
> {"column", f_column}, /* for using */
> {"valid", f_valid}, /* for using */
> {"timecolumn", f_timecolumn}, /* for using */
> @@ -689,6 +692,20 @@
> return (*udv_ptr);
> }
>
> +struct udft_entry *
> +get_udf_by_name(char *key)
> +{
> + struct udft_entry *udf_ptr = first_udf;
> +
> + while (udf_ptr) {
> + if (!strcmp(key, udf_ptr->udf_name))
> + return udf_ptr;
> +
> + udf_ptr = udf_ptr->next_udf;
> + }
> +
> + return NULL;
> +}
>
> static void update_plot_bounds __PROTO((void));
> static void fill_gpval_axis __PROTO((AXIS_INDEX axis));
> --- eval.h.orig 2011-01-11 17:38:21.000000000 +0100
> +++ eval.h 2011-01-12 15:54:30.000000000 +0100
> @@ -164,6 +164,7 @@
> void apollo_pfm_catch __PROTO((void));
> #endif
> struct udvt_entry * add_udv_by_name __PROTO((char *key));
> +struct udft_entry * get_udf_by_name __PROTO((char *key));
>
> /* update GPVAL_ variables available to user */
> void update_gpval_variables __PROTO((int from_plot_command));
> --- specfun.c.orig 2011-01-11 17:37:27.000000000 +0100
> +++ specfun.c 2011-01-12 16:14:28.000000000 +0100
> @@ -49,6 +49,7 @@
> #include "specfun.h"
> #include "stdfn.h"
> #include "util.h"
> +#include "internal.h"
>
> #define ITMAX 200
>
> @@ -1933,3 +1934,130 @@
> push(Gcomplex(&a, x, 0.0));
> }
>
> +/* Finite summation
> + *
> + * Calculate the finite sum(a, b, "f") = \sum_{i=a}^b f(i), where "f" is the
> + * name (a string) of f(i), a function taking one (integer) parameter i.
> + */
> +void
> +f_sum(union argument *arg)
> +{
> + struct value tmp;
> + int a, b; // bounds [a, b]
> + struct cmplx ret;
> + int i;
> + struct udft_entry *udf;
> + union argument call_arg;
> +
> + (void) arg; /* avoid -Wunused warning */
> +
> + // pop function name
> + if (pop(&tmp)->type != STRING)
> + int_error(NO_CARET, "expecting a name of a function with 1 parameters as 3rd argument to sum");
> + if (!(udf = get_udf_by_name(tmp.v.string_val)))
> + int_error(NO_CARET, "No user-defined function %s", tmp.v.string_val);
> +
> + // pop upper bound (b)
> + if (pop(&tmp)->type != INTGR)
> + int_warn(NO_CARET, "expecting an integer as 2nd argument to sum");
> + b = tmp.v.int_val;
> +
> + // pop lower bound (a)
> + if (pop(&tmp)->type != INTGR)
> + int_warn(NO_CARET, "expecting an integer as 1st argument to sum");
> + a = tmp.v.int_val;
> +
> + ret.real = 0;
> + ret.imag = 0;
> + for (i=a; i<=b; ++i) {
> +
> + /* calculate f_i = f(i); */
> +
> + push(Ginteger(&tmp, i));
> +
> + call_arg.udf_arg = udf;
> + f_call(&call_arg);
> +
> + pop(&tmp);
> + ret.real += real(&tmp);
> + ret.imag += imag(&tmp);
> + }
> +
> + push(Gcomplex(&tmp, ret.real, ret.imag));
> +
> +}
> +
> +/* Finite summation (with free parameter x)
> + *
> + * Calculate the finite sum(a, b, x, "f") = \sum_{i=a}^b f(i, x), where "f" is the
> + * name (a string) of f(i, x), a function taking two parameters i (integer) and x
> + * (complex). Example usage:
> + *
> + * # fourier coefficients of some square wave
> + * fourier(k, x) = sin(3./2*k)/k * 2./3*cos(k*x)
> + * s(x) = 1./2 + sum(1, 10, x, "fourier")
> + * plot s(x)
> + */
> +void
> +f_sumx(union argument *arg)
> +{
> + struct value tmp;
> + int a, b; // bounds [a, b]
> + struct cmplx x, ret;
> + int i;
> + udft_entry *udf;
> + union argument calln_arg;
> +
> + // XXX is arg every used in any other function if so, what is it used for?
> + (void) arg; /* avoid -Wunused warning */
> +
> + // XXX do I have to free the poped values?
> + // pop function name
> + if (pop(&tmp)->type != STRING)
> + int_error(NO_CARET, "expecting a name of a function with 2 parameters as 4th argument to sum");
> + if (!(udf = get_udf_by_name(tmp.v.string_val)))
> + int_error(NO_CARET, "No user-defined function %s", tmp.v.string_val);
> +
> + // pop x
> + pop(&tmp);
> + x.real = real(&tmp);
> + x.imag = imag(&tmp);
> +
> + // pop upper bound (b)
> + if (pop(&tmp)->type != INTGR)
> + int_warn(NO_CARET, "expecting an integer as 2nd argument to sum");
> + b = tmp.v.int_val;
> +
> + // pop lower bound (a)
> + if (pop(&tmp)->type != INTGR)
> + int_warn(NO_CARET, "expecting an integer as 1st argument to sum");
> + a = tmp.v.int_val;
> +
> + ret.real = 0;
> + ret.imag = 0;
> + for (i=a; i<=b; ++i) {
> +
> + /* calculate f_i(x) = f(i, x); */
> +
> + /* XXX in what order do I have to push?
> + * I thought that I would have to push from right to left but my
> + * testing contradicts that finding. I am very unsure about the order
> + * and need to understand this
> + */
> + // XXX do I have to malloc instead of pushing a temporary?
> + push(Ginteger(&tmp, i));
> + push(Gcomplex(&tmp, x.real, x.imag));
> + push(Ginteger(&tmp, 2)); // push argument count
> +
> + calln_arg.udf_arg = udf;
> + f_calln(&calln_arg);
> +
> + pop(&tmp);
> + ret.real += real(&tmp);
> + ret.imag += imag(&tmp);
> + }
> +
> + push(Gcomplex(&tmp, ret.real, ret.imag));
> +
> +}
> +
> --- specfun.h.orig 2011-01-11 17:37:21.000000000 +0100
> +++ specfun.h 2011-01-12 15:36:01.000000000 +0100
> @@ -61,4 +61,7 @@
> void f_inverse_erf __PROTO((union argument *x));
> void f_lambertw __PROTO((union argument *x));
>
> +void f_sum __PROTO((union argument *x));
> +void f_sumx __PROTO((union argument *x));
> +
> #endif /* GNUPLOT_SPECFUN_H */
> --- ../docs/gnuplot.doc.orig 2011-01-12 15:09:42.000000000 +0100
> +++ ../docs/gnuplot.doc 2011-01-12 16:22:47.000000000 +0100
> @@ -1091,6 +1091,30 @@
> %sinh(x)@any@$sinh~x$, hyperbolic sine of $x$ in radians
> The `sinh(x)` function returns the hyperbolic sine of its argument. `sinh`
> expects its argument to be in radians.
> +4 sum
> +?expressions functions sum
> +?functions sum
> +?sum
> +#sum(a, b, "f") & any & $\sum_{i=a}^b f(i)$, finite sum \\
> +%sum(a, b, "f")@any@$sum(i=a..b, f(i)$, finite sum
> + The `sum(a, b, "f")` function calculates the finite sum \sum_{i=a}^b f(i),
> + where "f" is the name (a string) of f(i), a function taking one (integer)
> + parameter i. For example, `f(n)=n; print sum(1, 10, "f")` prints 55.
> +4 sumx
> +?expressions functions sumx
> +?functions sumx
> +?sumx
> +#sumx(a, b, x, "f") & any & $\sum_{i=a}^b f(i, x)$, finite sum (with parameter $x$)\\
> +%sumx(a, b, x, "f")@any@$sumx(i=a..b, f(i, x)$, finite sum (with parameter $x$)
> + The `sumx(a, b, x, "f")` function calculates the finite sum
> + \sum_{i=a}^b f(i, x), where "f" is the name (a string) of f(i, x), a function
> + taking two parameters i (integer) and x (complex).
> +
> + Example:
> + # fourier coefficients of some square wave
> + fourier(k, x) = sin(3./2*k)/k * 2./3*cos(k*x)
> + s(x) = 1./2 + sumx(1, 10, x, "fourier")
> + plot s(x)
> 4 sqrt
> ?expressions functions sqrt
> ?functions sqrt
>
>
--
Ethan A Merritt
Biomolecular Structure Center, K-428 Health Sciences Bldg
University of Washington, Seattle 98195-7742
|
