Menu
▾
▴
Re: finite summation - how to call a function from parameter stack
|
From: <pl...@pi...> - 2011-01-12 23:29:59
|
On 12/01/11 20:40, Ethan Merritt wrote:
> 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
>>
>>
>
If I was hitting 251 recursive calls of a function to do a sum
calculation I would be inclined to ask myself if it was the best way of
doing things.
/P
|
