Menu
▾
▴
gnuplot-beta
|
From: <mw...@gm...> - 2011-01-11 17:00:29
|
Hi,
I want to add finite summation to gnuplot, that is a function sum(a, b, x, f), where a, b (int) are the summation boundary, x (double) is a varialbe, f (string) is the name of the function to eveluate a term. Example usage:
f(n, x) = sin(n*x)
plot sum(1, 10, x, "f")
I am stuck at how to call f from sum (see call(f) below) and would appreciate your help very much. (I assume f is the string representing the name of the function to call.)
- for (i=a; i<=b; ++i) {
- // calculate f_i(x);
- push(Gcomplex(&tmp, x.real, x.imag)); // todo check the order
- push(Ginteger(&tmp, i));
-
- //call(f); // assume f is the name (string) of the function to call
-
- pop(&tmp);
- ret.real += real(&tmp);
- ret.imag += imag(&tmp);
- }
Thank you very much for your time and effort.
--- specfun.c 2011-01-11 17:40:19.000000000 +0100
+++ specfun.c.orig 2011-01-11 17:37:27.000000000 +0100
@@ -1933,44 +1933,3 @@
push(Gcomplex(&a, x, 0.0));
}
-void
-f_sum(union argument *arg)
-{
- struct value tmp;
- int a, b;
- struct cmplx x, ret, result;
- int i;
-
- (void) arg; /* avoid -Wunused warning */
-
- // pop x
- pop(&tmp);
- x.real = real(&tmp);
- x.imag = imag(&tmp);
-
- // pop b
- if (pop(&tmp)->type != INTGR)
- int_warn(NO_CARET, "expecting an integer as second argument to sum");
- b = tmp.v.int_val;
-
- // pop a
- if (pop(&tmp)->type != INTGR)
- int_warn(NO_CARET, "expecting an integer as first argument to sum");
- a = tmp.v.int_val;
-
- ret.real = 0;
- ret.imag = 0;
- for (i=a; i<=b; ++i) {
- // calculate f_i(x);
- push(Gcomplex(&tmp, x.real, x.imag)); // todo check the order
- push(Ginteger(&tmp, i));
-
- //call(f);
-
- pop(&tmp);
- ret.real += real(&tmp);
- ret.imag += imag(&tmp);
- }
-
- push(Gcomplex(&tmp, ret.real, ret.imag));
-}
--- eval.c 2011-01-11 16:25:01.000000000 +0100
+++ eval.c.orig 2011-01-11 17:38:00.000000000 +0100
@@ -176,8 +176,6 @@
{"atanh", f_atanh},
{"lambertw", f_lambertw}, /* HBB, from G.Kuhnle 20001107 */
- {"sum", f_sum}, /* summation */
-
{"column", f_column}, /* for using */
{"valid", f_valid}, /* for using */
{"timecolumn", f_timecolumn}, /* for using */
--- specfun.h 2011-01-11 16:22:05.000000000 +0100
+++ specfun.h.orig 2011-01-11 17:37:21.000000000 +0100
@@ -61,6 +61,4 @@
void f_inverse_erf __PROTO((union argument *x));
void f_lambertw __PROTO((union argument *x));
-void f_sum __PROTO((union argument *x));
-
#endif /* GNUPLOT_SPECFUN_H */
--
GMX DSL Doppel-Flat ab 19,99 Euro/mtl.! Jetzt mit
gratis Handy-Flat! http://portal.gmx.net/de/go/dsl
|
|
From: Ethan M. <merritt@u.washington.edu> - 2011-01-11 17:22:59
|
On Tuesday, January 11, 2011 09:00:19 am mw...@gm... wrote:
> Hi,
>
> I want to add finite summation to gnuplot, that is a function sum(a, b, x, f), where a, b (int) are the summation boundary, x (double) is a varialbe, f (string) is the name of the function to eveluate a term. Example usage:
>
> f(n, x) = sin(n*x)
> plot sum(1, 10, x, "f")
That pair of commands does not seem to be compatible with your description.
How is the code to know there is a second parameter to f()? And where would
it get this second parameter from?
> I am stuck at how to call f from sum (see call(f) below) and would appreciate your help very much. (I assume f is the string representing the name of the function to call.)
>
> - for (i=a; i<=b; ++i) {
> - // calculate f_i(x);
> - push(Gcomplex(&tmp, x.real, x.imag)); // todo check the order
> - push(Ginteger(&tmp, i));
> -
> - //call(f); // assume f is the name (string) of the function to call
There is no existing routine that does either call-by-name or
retrieve-user-function-by-name; you would have to write one.
Something like (untested):
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;
}
It would belong in eval.c next to add_udv_by_name().
Then you would do something like:
udft_entry *udf = get_udf_by_name( f );
if (udf)
execute_at(udf->at);
else
int_error(NOCARET, "No user-defined function %s",f);
> -
> - pop(&tmp);
> - ret.real += real(&tmp);
> - ret.imag += imag(&tmp);
> - }
>
> Thank you very much for your time and effort.
>
> --- specfun.c 2011-01-11 17:40:19.000000000 +0100
> +++ specfun.c.orig 2011-01-11 17:37:27.000000000 +0100
> @@ -1933,44 +1933,3 @@
> push(Gcomplex(&a, x, 0.0));
> }
>
> -void
> -f_sum(union argument *arg)
> -{
> - struct value tmp;
> - int a, b;
> - struct cmplx x, ret, result;
> - int i;
> -
> - (void) arg; /* avoid -Wunused warning */
> -
> - // pop x
> - pop(&tmp);
> - x.real = real(&tmp);
> - x.imag = imag(&tmp);
> -
> - // pop b
> - if (pop(&tmp)->type != INTGR)
> - int_warn(NO_CARET, "expecting an integer as second argument to sum");
> - b = tmp.v.int_val;
> -
> - // pop a
> - if (pop(&tmp)->type != INTGR)
> - int_warn(NO_CARET, "expecting an integer as first argument to sum");
> - a = tmp.v.int_val;
> -
> - ret.real = 0;
> - ret.imag = 0;
> - for (i=a; i<=b; ++i) {
> - // calculate f_i(x);
> - push(Gcomplex(&tmp, x.real, x.imag)); // todo check the order
> - push(Ginteger(&tmp, i));
> -
> - //call(f);
> -
> - pop(&tmp);
> - ret.real += real(&tmp);
> - ret.imag += imag(&tmp);
> - }
> -
> - push(Gcomplex(&tmp, ret.real, ret.imag));
> -}
>
> --- eval.c 2011-01-11 16:25:01.000000000 +0100
> +++ eval.c.orig 2011-01-11 17:38:00.000000000 +0100
> @@ -176,8 +176,6 @@
> {"atanh", f_atanh},
> {"lambertw", f_lambertw}, /* HBB, from G.Kuhnle 20001107 */
>
> - {"sum", f_sum}, /* summation */
> -
> {"column", f_column}, /* for using */
> {"valid", f_valid}, /* for using */
> {"timecolumn", f_timecolumn}, /* for using */
> --- specfun.h 2011-01-11 16:22:05.000000000 +0100
> +++ specfun.h.orig 2011-01-11 17:37:21.000000000 +0100
> @@ -61,6 +61,4 @@
> void f_inverse_erf __PROTO((union argument *x));
> void f_lambertw __PROTO((union argument *x));
>
> -void f_sum __PROTO((union argument *x));
> -
> #endif /* GNUPLOT_SPECFUN_H */
>
>
|
|
From: Hans-Bernhard B. <HBB...@t-...> - 2011-01-11 20:07:09
|
On 11.01.2011 18:00, mw...@gm... wrote:
> I want to add finite summation to gnuplot, that is a function sum(a,
> b, x, f), where a, b (int) are the summation boundary, x (double) is
> a varialbe, f (string) is the name of the function to eveluate a
> term. Example usage:
Without changing gnuplot, one can already do this recursive function:
sum_f(a,b,delta) = (a < b) ? 0 : (f(a) + sum_f(a+delta,b,delta))
For more tricky examples, see bivariat.dem
The only aspect of your job this doesn't cover is passing in a function
as another function's argument. gnuplot is a plotting program after
all, not a full-fledged functional programming language.
I don't think adding 'sum' as a new builtin is worth the effort. If a
new builtin is to be added, then that should be an eval() function, i.e.
a function that takes a string and treats it as an expression to
evaluate, so the above could be expanded to
sum(f,a,b,delta) = (a<b)?0 : (eval("f(a)")+sum_f(a+delta,b,delta))
|
|
From: <mw...@gm...> - 2011-01-11 20:37:52
|
> Ethan Merritt <merritt@u.washington.edu> wrote:
> On Tuesday, January 11, 2011 09:00:19 am mw...@gm... wrote:
>> Hi,
>>
>> I want to add finite summation to gnuplot, that is a function sum(a, b,
>> x, f), where a, b (int) are the summation boundary, x (double) is a
>> varialbe, f (string) is the name of the function to eveluate a term.
>> Example usage:
>>
>> f(n, x) = sin(n*x)
>> plot sum(1, 10, x, "f")
>
> That pair of commands does not seem to be compatible with your
> description. How is the code to know there is a second parameter to
> f()? And where would it get this second parameter from?
sum requires its 4-th parameter to be the name of a 2 parametric function. Before the call to f in sum, two values are pushed. I think I missed the point of what you were trying to tell me.
>> I am stuck at how to call f from sum (see call(f) below) and would
>> appreciate your help very much. (I assume f is the string representing
>> the name of the function to call.)
>>
>> - for (i=a; i<=b; ++i) {
>> - // calculate f_i(x);
>> - push(Gcomplex(&tmp, x.real, x.imag)); // todo check the order
>> - push(Ginteger(&tmp, i));
>> -
>> - //call(f); // assume f is the name (string) of the function to call
>
> There is no existing routine that does either call-by-name or
> retrieve-user-function-by-name; you would have to write one. Something
> like (untested):
>
> 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;
> }
>
> It would belong in eval.c next to add_udv_by_name().
>
> Then you would do something like:
>
> udft_entry *udf = get_udf_by_name( f );
> if (udf)
> execute_at(udf->at);
> else
> int_error(NOCARET, "No user-defined function %s",f);
>
>
>> -
>> - pop(&tmp);
>> - ret.real += real(&tmp);
>> - ret.imag += imag(&tmp);
>> - }
>>
>> Thank you very much for your time and effort.
>>
With your generous help I was able to get stuck a little bit later. I assume its related to your initial comment about the parameters of f. I thought, that I have to push two values (int n, cmplx x) onto the stack and then execute_at(..). I could verify (by calling pop()/more_on_stack() repediatly) that execute_at(..) doesn't pop values from the stack. In fact I believe that initially the stack containes two values (n, x) and then execute_at(..) adds another. Did I do it wrong? How can I pass n, x to the function if not by pushing them?
+ for (i=a; i<=b; ++i) {
+
+ // calculate f_i(x);
+ push(Gcomplex(&tmp, x.real, x.imag)); // todo check the order
+ push(Ginteger(&tmp, i));
+
+ if (udf = get_udf_by_name(f))
+ execute_at(udf->at);
+ else
+ int_error(NO_CARET, "No user-defined function %s", f);
+
+ /*
+ int_warn(NO_CARET, "stack_size == 3 (%d)", more_on_stack()); pop(&tmp);
+ int_warn(NO_CARET, "stack_size == 2 (%d)", more_on_stack()); pop(&tmp);
+ int_warn(NO_CARET, "stack_size == 1 (%d)", more_on_stack()); pop(&tmp);
+ int_warn(NO_CARET, "stack_size == 0 (%d)", more_on_stack());
+ */
+
+ pop(&tmp);
+ ret.real += real(&tmp);
+ ret.imag += imag(&tmp);
+ }
Thank you Sir.
--- specfun.c.orig 2011-01-11 17:37:27.000000000 +0100
+++ specfun.c 2011-01-11 21:06:01.000000000 +0100
@@ -1933,3 +1933,57 @@
push(Gcomplex(&a, x, 0.0));
}
+void
+f_sum(union argument *arg)
+{
+ struct value tmp;
+ int a, b;
+ struct cmplx x, ret, result;
+ int i;
+ udft_entry *udf;
+ char *f = "f";
+
+ (void) arg; /* avoid -Wunused warning */
+
+ // pop x
+ pop(&tmp);
+ x.real = real(&tmp);
+ x.imag = imag(&tmp);
+
+ // pop b
+ if (pop(&tmp)->type != INTGR)
+ int_warn(NO_CARET, "expecting an integer as second argument to sum");
+ b = tmp.v.int_val;
+
+ // pop a
+ if (pop(&tmp)->type != INTGR)
+ int_warn(NO_CARET, "expecting an integer as first argument to sum");
+ a = tmp.v.int_val;
+
+ ret.real = 0;
+ ret.imag = 0;
+ for (i=a; i<=b; ++i) {
+
+ // calculate f_i(x);
+ push(Gcomplex(&tmp, x.real, x.imag)); // todo check the order
+ push(Ginteger(&tmp, i));
+
+ if (udf = get_udf_by_name(f))
+ execute_at(udf->at);
+ else
+ int_error(NO_CARET, "No user-defined function %s", f);
+
+ /*
+ int_warn(NO_CARET, "stack_size == 3 (%d)", more_on_stack()); pop(&tmp);
+ int_warn(NO_CARET, "stack_size == 2 (%d)", more_on_stack()); pop(&tmp);
+ int_warn(NO_CARET, "stack_size == 1 (%d)", more_on_stack()); pop(&tmp);
+ int_warn(NO_CARET, "stack_size == 0 (%d)", more_on_stack());
+ */
+
+ pop(&tmp);
+ ret.real += real(&tmp);
+ ret.imag += imag(&tmp);
+ }
+
+ push(Gcomplex(&tmp, ret.real, ret.imag));
+}
--- eval.c.orig 2011-01-11 17:38:00.000000000 +0100
+++ eval.c 2011-01-11 20:41:11.000000000 +0100
@@ -176,6 +176,8 @@
{"atanh", f_atanh},
{"lambertw", f_lambertw}, /* HBB, from G.Kuhnle 20001107 */
+ {"sum", f_sum}, /* summation */
+
{"column", f_column}, /* for using */
{"valid", f_valid}, /* for using */
{"timecolumn", f_timecolumn}, /* for using */
@@ -689,6 +691,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));
zsh: exit 1 diff -u eval.c.orig eval.c
--- specfun.h.orig 2011-01-11 17:37:21.000000000 +0100
+++ specfun.h 2011-01-11 16:22:05.000000000 +0100
@@ -61,4 +61,6 @@
void f_inverse_erf __PROTO((union argument *x));
void f_lambertw __PROTO((union argument *x));
+void f_sum __PROTO((union argument *x));
+
#endif /* GNUPLOT_SPECFUN_H */
--
NEU: FreePhone - kostenlos mobil telefonieren und surfen!
Jetzt informieren: http://www.gmx.net/de/go/freephone
|
|
From: Ethan M. <merritt@u.washington.edu> - 2011-01-11 21:32:11
|
On Tuesday, January 11, 2011 12:37:42 pm mw...@gm... wrote:
> > Ethan Merritt <merritt@u.washington.edu> wrote:
> > On Tuesday, January 11, 2011 09:00:19 am mw...@gm... wrote:
> >> Hi,
> >>
> >> I want to add finite summation to gnuplot, that is a function sum(a, b,
> >> x, f), where a, b (int) are the summation boundary, x (double) is a
> >> varialbe, f (string) is the name of the function to eveluate a term.
> >> Example usage:
> >>
> >> f(n, x) = sin(n*x)
> >> plot sum(1, 10, x, "f")
> >
> > That pair of commands does not seem to be compatible with your
> > description. How is the code to know there is a second parameter to
> > f()? And where would it get this second parameter from?
>
> sum requires its 4-th parameter to be the name of a 2 parametric function. Before the call to f in sum, two values are pushed. I think I missed the point of what you were trying to tell me.
I guess I don't understand what you are trying to do.
(1) I thought f was a string parameter passed to f_sum, containing the name
of a user-defined function. But your code never pops it as a parameter
and instead refers to it as a string constant "f". That may be just for
debugging but it contradicts your definition of the sum() function so I'm
not sure.
(2) I still don't understand what is "n", or why you want to pass a second
parameter to f.
Are you thinking that you want an infinite class of functions f_n(x)?
That would take a lot of work to support.
> With your generous help I was able to get stuck a little bit later. I assume its related to your initial comment about the parameters of f. I thought, that I have to push two values (int n, cmplx x) onto the stack and then execute_at(..). I could verify (by calling pop()/more_on_stack() repediatly) that execute_at(..) doesn't pop values from the stack. In fact I believe that initially the stack containes two values (n, x) and then execute_at(..) adds another. Did I do it wrong? How can I pass n, x to the function if not by pushing them?
You either have to duplicate the relevant parts of f_calln(),
or else you have to set up all the bookkeeping variables needed by
f_calln() and call f_calln() itself to do the work.
Ethan
|
|
From: <mw...@gm...> - 2011-01-12 15:50:08
|
Hi,
> Ethan Merritt <merritt@u.washington.edu> wrote:
> On Tuesday, January 11, 2011 12:37:42 pm mw...@gm... wrote:
>>> Ethan Merritt <merritt@u.washington.edu> wrote:
>>> On Tuesday, January 11, 2011 09:00:19 am mw...@gm... wrote:
>>>> Hi,
>>>>
>>>> I want to add finite summation to gnuplot, that is a function sum(a, b,
>>>> x, f), where a, b (int) are the summation boundary, x (double) is a
>>>> varialbe, f (string) is the name of the function to eveluate a term.
>>>> Example usage:
>>>>
>>>> f(n, x) = sin(n*x)
>>>> plot sum(1, 10, x, "f")
>>>
>>> That pair of commands does not seem to be compatible with your
>>> description. How is the code to know there is a second parameter to
>>> f()? And where would it get this second parameter from?
>>
>> sum requires its 4-th parameter to be the name of a 2 parametric
>> function. Before the call to f in sum, two values are pushed. I think I
>> missed the point of what you were trying to tell me.
>
> I guess I don't understand what you are trying to do.
>
> (1) I thought f was a string parameter passed to f_sum, containing the
> name of a user-defined function. But your code never pops it as a parameter
> and instead refers to it as a string constant "f". That may be just for
> debugging but it contradicts your definition of the sum() function so I'm
> not sure.
>
> (2) I still don't understand what is "n", or why you want to pass a second
> parameter to f.
> Are you thinking that you want an infinite class of functions f_n(x)?
> That would take a lot of work to support.
>
>> With your generous help I was able to get stuck a little bit later. I
> assume its related to your initial comment about the parameters of f. I
> thought, that I have to push two values (int n, cmplx x) onto the stack and
> then execute_at(..). I could verify (by calling pop()/more_on_stack()
> repediatly) that execute_at(..) doesn't pop values from the stack. In fact I
> believe that initially the stack containes two values (n, x) and then
> execute_at(..) adds another. Did I do it wrong? How can I pass n, x to the
> function if > not by pushing them?
>
> You either have to duplicate the relevant parts of f_calln(),
> or else you have to set up all the bookkeeping variables needed by
> f_calln() and call f_calln() itself to do the work.
Thank you again for pointing me into the right direction. I didn't realize
fast enought that I haven't explained enough, what I try to do, so let me
explain. 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.) 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 have added below a working patch but would like you to have a look at it. 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?
* 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).
* Do I have to alloc memeory for the values I push on the stack befor the call to f(.)?
* This is only curiosity is arg in foo(union argument *arg) every used and if, what for?
I would also like you to take a look at my addition to gnuplot.doc and proofread it.
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).
Thank you again for your generous help without I would never have gone so far.
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
--
NEU: FreePhone - kostenlos mobil telefonieren und surfen!
Jetzt informieren: http://www.gmx.net/de/go/freephone
|
|
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
|
|
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
|
|
From: <mw...@gm...> - 2011-01-19 18:16:08
|
Hi,
I have implemented 'sum' which can be used as folows:
a = sum [k=1:4] k
f(x) = sum [k=1:4] sin(k*x)
print sum [k=1:4] k
plot sum [k=1:4] sin(k*x)
Motivation: A friend asked me how to plot the sum of the first 10, 100, 1000
fourier coefficients of a (particular) square wave. This patch allows to do so:
fourier(k, x) = sin(3./2*k)/k * 2./3*cos(k*x)
plot 1./2 + sum [k=1:1000] fourier(k, x)
# alternatively
plot 1./2 + sum [k=1:1000] sin(3./2*k)/k * 2./3*cos(k*x)
Questions: I am not sure where to put the documentation in gnuplot.doc:
| The `sum` keyword allows the calulation of the finite sum
| \sum_{i=a}^b f_i.
|
| Example:
| # fourier coefficients of some square wave
| fourier(k, x) = sin(3./2*k)/k * 2./3*cos(k*x)
| plot 1./2 + sum [k=1:1000] fourier(k, x)
|
| plot 1./2 + sum [k=1:1000] sin(3./2*k)/k * 2./3*cos(k*x)
> Ethan Merritt <merritt@u.washington.edu> wrote:
> On Wednesday, January 12, 2011 07:49:57 am mw...@gm... wrote:
[...]
> My first thought would be to aim for
> fourier(x) = for [k=1:100] sin(3./2*k)/k * 2./3*cos(k*x)
I had to rename it to 'sum', since 'for' turned out to be ambiguous (for
example: "plot for [k=1:10] sin(k*x)")
>> 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.
I suppose an experienced gnuplot developer knows all the tricks by heart whereas
I had a much harder time and ultimately was not able to solve the original
problem of plotting the sum of 1000 terms. I tend to believe that '?:' results
in hard to read code and 'sum' is easier to use and understand by the "normal"
user. 'sum' also removes the recursion.
[...]
>> 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 not done so yet, but will do so if my patch is not complete garbage.
Thank you for your consideration,
Micha Wiedenmann
diff --git gnuplot-cvs/demo/bivariat.dem b/demo/bivariat.dem
index 1865a59..21c46be 100644
--- gnuplot-cvs/demo/bivariat.dem
+++ b/demo/bivariat.dem
@@ -115,5 +115,27 @@ set title "Greatest Common Divisor (for integers only)"
plot gcd(x, 60) with impulses
pause -1 "Hit return to continue"
+#
+# This definition computes the sum of the first 10, 100, 1000 fourier
+# coefficients of a (particular) square wave.
+
+set title "Finite summation of 10, 100, 1000 fourier coefficients"
+
+set samples 500
+set xrange [-10:10]
+set yrange [-0.4:1.2]
+set key bottom right
+
+fourier(k, x) = sin(3./2*k)/k * 2./3*cos(k*x)
+sum10(x) = 1./2 + sum [k=1:10] fourier(k, x)
+sum100(x) = 1./2 + sum [k=1:100] fourier(k, x)
+sum1000(x) = 1./2 + sum [k=1:1000] fourier(k, x)
+
+plot \
+ sum10(x) title "1./2 + sum [k=1:10] sin(3./2*k)/k * 2./3*cos(k*x)", \
+ sum100(x) title "1./2 + sum [k=1:100] sin(3./2*k)/k * 2./3*cos(k*x)", \
+ sum1000(x) title "1./2 + sum [k=1:1000] sin(3./2*k)/k * 2./3*cos(k*x)"
+pause -1 "Hit return to continue"
+
reset
diff --git gnuplot-cvs/src/eval.c b/src/eval.c
index 8b551d8..64517ab 100644
--- gnuplot-cvs/src/eval.c
+++ b/src/eval.c
@@ -89,6 +89,7 @@ const struct ft_entry GPFAR ft[] =
{"pop", f_pop},
{"call", f_call},
{"calln", f_calln},
+ {"sum", f_sum},
{"lnot", f_lnot},
{"bnot", f_bnot},
{"uminus", f_uminus},
@@ -669,6 +670,20 @@ add_udv_by_name(char *key)
return (*udv_ptr);
}
+struct udvt_entry *
+get_udv_by_name(char *key)
+{
+ struct udvt_entry *udv = first_udv;
+
+ while (udv) {
+ if (!strcmp(key, udv->udv_name))
+ return udv;
+
+ udv = udv->next_udv;
+ }
+
+ return NULL;
+}
static void update_plot_bounds __PROTO((void));
static void fill_gpval_axis __PROTO((AXIS_INDEX axis));
diff --git gnuplot-cvs/src/eval.h b/src/eval.h
index 79edb6e..ed02262 100644
--- gnuplot-cvs/src/eval.h
+++ b/src/eval.h
@@ -53,7 +53,7 @@
enum operators {
/* keep this in line with table in eval.c */
PUSH, PUSHC, PUSHD1, PUSHD2, PUSHD, POP,
- CALL, CALLN, LNOT, BNOT, UMINUS,
+ CALL, CALLN, SUM, LNOT, BNOT, UMINUS,
LOR, LAND, BOR, XOR, BAND, EQ, NE, GT, LT, GE, LE, PLUS, MINUS, MULT,
DIV, MOD, POWER, FACTORIAL, BOOLE,
DOLLARS, /* for using extension - div */
@@ -161,6 +161,7 @@ void execute_at __PROTO((struct at_type *at_ptr));
void evaluate_at __PROTO((struct at_type *at_ptr, struct value *val_ptr));
void free_at __PROTO((struct at_type *at_ptr));
struct udvt_entry * add_udv_by_name __PROTO((char *key));
+struct udvt_entry * get_udv_by_name __PROTO((char *key));
/* update GPVAL_ variables available to user */
void update_gpval_variables __PROTO((int from_plot_command));
diff --git gnuplot-cvs/src/internal.c b/src/internal.c
index b4cd25f..706df95 100644
--- gnuplot-cvs/src/internal.c
+++ b/src/internal.c
@@ -200,6 +200,54 @@ f_calln(union argument *x)
void
+f_sum(union argument *arg)
+{
+ struct value beg, end, varname; /* [<var> = <start>:<end>] */
+ udft_entry *udf; /* function to evaluate */
+ udvt_entry *udv; /* iteration variable */
+ struct value ret; /* result */
+ struct value z;
+ int i;
+
+ (void) pop(&end);
+ (void) pop(&beg);
+ (void) pop(&varname);
+
+ if (beg.type != INTGR || end.type != INTGR)
+ int_error(NO_CARET, "range specifiers of sum must have integer values");
+ if (varname.type != STRING)
+ int_error(NO_CARET, "internal error: f_sum expects argument (varname) of type string.");
+
+ udv = get_udv_by_name(varname.v.string_val);
+ if (!udv)
+ int_error(NO_CARET, "internal error: f_sum could not access iteration variable.");
+ udv->udv_undef = false;
+
+ udf = arg->udf_arg;
+ if (!udf)
+ int_error(NO_CARET, "internal error: f_sum could not access summation coefficient function");
+
+ Gcomplex(&ret, 0, 0);
+ for (i=beg.v.int_val; i<=end.v.int_val; ++i) {
+ double x, y;
+
+ /* calculate f_i = f() with user defined variable i */
+ Ginteger(&udv->udv_value, i);
+ execute_at(udf->at);
+
+ pop(&z);
+ x = real(&ret) + real(&z);
+ y = imag(&ret) + imag(&z);
+ Gcomplex(&ret, x, y);
+ }
+
+ gpfree_string(&varname);
+
+ push(Gcomplex(&z, real(&ret), imag(&ret)));
+}
+
+
+void
f_lnot(union argument *arg)
{
struct value a;
diff --git gnuplot-cvs/src/internal.h b/src/internal.h
index 9bb1a3a..7fbdbf8 100644
--- gnuplot-cvs/src/internal.h
+++ b/src/internal.h
@@ -55,6 +55,7 @@ void f_pushd __PROTO((union argument *x));
void f_pop __PROTO((union argument *x));
void f_call __PROTO((union argument *x));
void f_calln __PROTO((union argument *x));
+void f_sum __PROTO((union argument *x));
void f_lnot __PROTO((union argument *x));
void f_bnot __PROTO((union argument *x));
void f_lor __PROTO((union argument *x));
diff --git gnuplot-cvs/src/parse.c b/src/parse.c
index 2a665ab..954d9fe 100644
--- gnuplot-cvs/src/parse.c
+++ b/src/parse.c
@@ -89,6 +89,7 @@ static void parse_relational_expression __PROTO((void));
static void parse_additive_expression __PROTO((void));
static void parse_multiplicative_expression __PROTO((void));
static void parse_unary_expression __PROTO((void));
+static void parse_sum_expression __PROTO((void));
static int parse_assignment_expression __PROTO((void));
static int is_builtin_function __PROTO((int t_num));
@@ -190,7 +191,8 @@ string_or_express(struct at_type **atptr)
has_dummies = FALSE;
for (i = 0; i < at->a_count; i++) {
enum operators op_index = at->actions[i].index;
- if ( op_index == PUSHD1 || op_index == PUSHD2 || op_index == PUSHD ) {
+ if ( op_index == PUSHD1 || op_index == PUSHD2 || op_index == PUSHD
+ || op_index == SUM ) {
has_dummies = TRUE;
break;
}
@@ -488,6 +490,8 @@ parse_primary_expression()
c_token++;
add_action(call_type)->udf_arg = add_udf(tok);
}
+ } else if (equals(c_token, "sum")) {
+ parse_sum_expression();
/* dummy_func==NULL is a flag to say no dummy variables active */
} else if (dummy_func) {
if (equals(c_token, c_dummy_var[0])) {
@@ -834,6 +838,112 @@ parse_unary_expression()
parse_primary_expression();
}
+
+/* create action code for 'sum' expressions */
+static void
+parse_sum_expression()
+{
+ /* Design: Use a user defined variable (udv) as iterator variable (k). The
+ * original idea was to treat the expression after the range as a function
+ * f(k). Consider 'g(x) = sum [k=1:4] f(k)', there are two dummy variables
+ * 'x' and 'k' from different functions 'g' and 'f' which cannot be handled
+ * by the parser. */
+
+ char *errormsg = "Expecting 'sum [<var> = <start>:<end>]'\n";
+ char *varname = NULL;
+ union argument *arg;
+ struct udft_entry *udf;
+
+ struct at_type * save_at;
+ int save_at_size;
+ int i;
+
+ if (!equals(c_token, "sum"))
+ return;
+ c_token++;
+
+ if (!equals(c_token, "["))
+ int_error(c_token, errormsg);
+ c_token++;
+
+ /* <var> */
+ if (!isletter(c_token))
+ int_error(c_token, errormsg);
+ /* create a user defined variable and pass it to f_sum via the action
+ * table, since the argument of f_sum is already used by the udf */
+ m_capture(&varname, c_token, c_token);
+ add_udv(c_token);
+ arg = add_action(PUSHC);
+ Gstring(&(arg->v_arg), varname);
+ c_token++;
+
+ if (!equals(c_token, "="))
+ int_error(c_token, errormsg);
+ c_token++;
+
+ /* <start> */
+ if (!isanumber(c_token))
+ int_error(c_token, errormsg);
+ arg = add_action(PUSHC);
+ convert(&(arg->v_arg), c_token);
+ c_token++;
+
+ if (!equals(c_token, ":"))
+ int_error(c_token, errormsg);
+ c_token++;
+
+ /* <end> */
+ if (!isanumber(c_token))
+ int_error(c_token, errormsg);
+ arg = add_action(PUSHC);
+ convert(&(arg->v_arg), c_token);
+ c_token++;
+
+ /* TODO add increment */
+ if (!equals(c_token, "]"))
+ int_error(c_token, errormsg);
+ c_token++;
+
+ /* parse the next expression and convert it to an action table. */
+ /* save environment to restart parsing */
+ save_at = at;
+ save_at_size = at_size;
+
+ at = (struct at_type *) gp_alloc(sizeof(struct at_type), "action table");
+ at->a_count = 0;
+ /* taken from temp_at()
+ * XXX Why is it necessary to reset the action table? Shouldn't it be
+ * either sizeof(struct at_type) or a_count = 0? */
+ memset(at, 0, sizeof(*at));
+ at_size = MAX_AT_LEN;
+
+ /* Q: Do I have to save and restore parse_recursion_level?
+ * A: parse_recursion_level is used to abort parsing after the strings
+ * ('-\pi', '-\pi/2') in ('-\pi' -pi, '-\pi/2' -pi/2.), otherwise it would
+ * try to subtract pi from the string '-\pi'. This is only in effect for
+ * parse_recursion_level == 1 and string_result_only == true. My conclusion
+ * is thus to not touch parse_recursion_level. */
+ parse_expression();
+
+ /* save action table in a user defined function */
+ udf = (struct udft_entry *) gp_alloc(sizeof(struct udft_entry), "sum");
+ udf->next_udf = (struct udft_entry *) NULL;
+ udf->udf_name = NULL; /* TODO maybe add a name and definition */
+ udf->at = at;
+ udf->definition = NULL;
+ udf->dummy_num = 0;
+ for (i = 0; i < MAX_NUM_VAR; i++)
+ (void) Ginteger(&(udf->dummy_values[i]), 0);
+
+ /* restore environment */
+ at = save_at;
+ at_size = save_at_size;
+
+ /* pass the udf to f_sum using the argument */
+ add_action(SUM)->udf_arg = udf;
+}
+
+
/* find or add value and return pointer */
struct udvt_entry *
add_udv(int t_num)
--
Empfehlen Sie GMX DSL Ihren Freunden und Bekannten und wir
belohnen Sie mit bis zu 50,- Euro! https://freundschaftswerbung.gmx.de
|
|
From: <pl...@pi...> - 2011-01-19 20:49:40
|
On 01/19/11 19:15, mw...@gm... wrote:
> Hi,
>
> I have implemented 'sum' which can be used as folows:
>
> a = sum [k=1:4] k
> f(x) = sum [k=1:4] sin(k*x)
> print sum [k=1:4] k
> plot sum [k=1:4] sin(k*x)
>
> Motivation: A friend asked me how to plot the sum of the first 10, 100, 1000
> fourier coefficients of a (particular) square wave. This patch allows to do so:
>
> fourier(k, x) = sin(3./2*k)/k * 2./3*cos(k*x)
> plot 1./2 + sum [k=1:1000] fourier(k, x)
>
> # alternatively
> plot 1./2 + sum [k=1:1000] sin(3./2*k)/k * 2./3*cos(k*x)
>
>
> Questions: I am not sure where to put the documentation in gnuplot.doc:
> | The `sum` keyword allows the calulation of the finite sum
> | \sum_{i=a}^b f_i.
> |
> | Example:
> | # fourier coefficients of some square wave
> | fourier(k, x) = sin(3./2*k)/k * 2./3*cos(k*x)
> | plot 1./2 + sum [k=1:1000] fourier(k, x)
> |
> | plot 1./2 + sum [k=1:1000] sin(3./2*k)/k * 2./3*cos(k*x)
>
>
>> Ethan Merritt<merritt@u.washington.edu> wrote:
>> On Wednesday, January 12, 2011 07:49:57 am mw...@gm... wrote:
> [...]
>> My first thought would be to aim for
>> fourier(x) = for [k=1:100] sin(3./2*k)/k * 2./3*cos(k*x)
>
> I had to rename it to 'sum', since 'for' turned out to be ambiguous (for
> example: "plot for [k=1:10] sin(k*x)")
>
>>> 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.
>
> I suppose an experienced gnuplot developer knows all the tricks by heart whereas
> I had a much harder time and ultimately was not able to solve the original
> problem of plotting the sum of 1000 terms. I tend to believe that '?:' results
> in hard to read code and 'sum' is easier to use and understand by the "normal"
> user. 'sum' also removes the recursion.
>
> [...]
>
>>> 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 not done so yet, but will do so if my patch is not complete garbage.
>
> Thank you for your consideration,
> Micha Wiedenmann
>
>
> diff --git gnuplot-cvs/demo/bivariat.dem b/demo/bivariat.dem
> index 1865a59..21c46be 100644
> --- gnuplot-cvs/demo/bivariat.dem
> +++ b/demo/bivariat.dem
> @@ -115,5 +115,27 @@ set title "Greatest Common Divisor (for integers only)"
> plot gcd(x, 60) with impulses
> pause -1 "Hit return to continue"
>
> +#
> +# This definition computes the sum of the first 10, 100, 1000 fourier
> +# coefficients of a (particular) square wave.
> +
> +set title "Finite summation of 10, 100, 1000 fourier coefficients"
> +
> +set samples 500
> +set xrange [-10:10]
> +set yrange [-0.4:1.2]
> +set key bottom right
> +
> +fourier(k, x) = sin(3./2*k)/k * 2./3*cos(k*x)
> +sum10(x) = 1./2 + sum [k=1:10] fourier(k, x)
> +sum100(x) = 1./2 + sum [k=1:100] fourier(k, x)
> +sum1000(x) = 1./2 + sum [k=1:1000] fourier(k, x)
> +
> +plot \
> + sum10(x) title "1./2 + sum [k=1:10] sin(3./2*k)/k * 2./3*cos(k*x)", \
> + sum100(x) title "1./2 + sum [k=1:100] sin(3./2*k)/k * 2./3*cos(k*x)", \
> + sum1000(x) title "1./2 + sum [k=1:1000] sin(3./2*k)/k * 2./3*cos(k*x)"
> +pause -1 "Hit return to continue"
> +
> reset
>
> diff --git gnuplot-cvs/src/eval.c b/src/eval.c
> index 8b551d8..64517ab 100644
> --- gnuplot-cvs/src/eval.c
> +++ b/src/eval.c
> @@ -89,6 +89,7 @@ const struct ft_entry GPFAR ft[] =
> {"pop", f_pop},
> {"call", f_call},
> {"calln", f_calln},
> + {"sum", f_sum},
> {"lnot", f_lnot},
> {"bnot", f_bnot},
> {"uminus", f_uminus},
> @@ -669,6 +670,20 @@ add_udv_by_name(char *key)
> return (*udv_ptr);
> }
>
> +struct udvt_entry *
> +get_udv_by_name(char *key)
> +{
> + struct udvt_entry *udv = first_udv;
> +
> + while (udv) {
> + if (!strcmp(key, udv->udv_name))
> + return udv;
> +
> + udv = udv->next_udv;
> + }
> +
> + return NULL;
> +}
>
> static void update_plot_bounds __PROTO((void));
> static void fill_gpval_axis __PROTO((AXIS_INDEX axis));
> diff --git gnuplot-cvs/src/eval.h b/src/eval.h
> index 79edb6e..ed02262 100644
> --- gnuplot-cvs/src/eval.h
> +++ b/src/eval.h
> @@ -53,7 +53,7 @@
> enum operators {
> /* keep this in line with table in eval.c */
> PUSH, PUSHC, PUSHD1, PUSHD2, PUSHD, POP,
> - CALL, CALLN, LNOT, BNOT, UMINUS,
> + CALL, CALLN, SUM, LNOT, BNOT, UMINUS,
> LOR, LAND, BOR, XOR, BAND, EQ, NE, GT, LT, GE, LE, PLUS, MINUS, MULT,
> DIV, MOD, POWER, FACTORIAL, BOOLE,
> DOLLARS, /* for using extension - div */
> @@ -161,6 +161,7 @@ void execute_at __PROTO((struct at_type *at_ptr));
> void evaluate_at __PROTO((struct at_type *at_ptr, struct value *val_ptr));
> void free_at __PROTO((struct at_type *at_ptr));
> struct udvt_entry * add_udv_by_name __PROTO((char *key));
> +struct udvt_entry * get_udv_by_name __PROTO((char *key));
>
> /* update GPVAL_ variables available to user */
> void update_gpval_variables __PROTO((int from_plot_command));
> diff --git gnuplot-cvs/src/internal.c b/src/internal.c
> index b4cd25f..706df95 100644
> --- gnuplot-cvs/src/internal.c
> +++ b/src/internal.c
> @@ -200,6 +200,54 @@ f_calln(union argument *x)
>
>
> void
> +f_sum(union argument *arg)
> +{
> + struct value beg, end, varname; /* [<var> =<start>:<end>] */
> + udft_entry *udf; /* function to evaluate */
> + udvt_entry *udv; /* iteration variable */
> + struct value ret; /* result */
> + struct value z;
> + int i;
> +
> + (void) pop(&end);
> + (void) pop(&beg);
> + (void) pop(&varname);
> +
> + if (beg.type != INTGR || end.type != INTGR)
> + int_error(NO_CARET, "range specifiers of sum must have integer values");
> + if (varname.type != STRING)
> + int_error(NO_CARET, "internal error: f_sum expects argument (varname) of type string.");
> +
> + udv = get_udv_by_name(varname.v.string_val);
> + if (!udv)
> + int_error(NO_CARET, "internal error: f_sum could not access iteration variable.");
> + udv->udv_undef = false;
> +
> + udf = arg->udf_arg;
> + if (!udf)
> + int_error(NO_CARET, "internal error: f_sum could not access summation coefficient function");
> +
> + Gcomplex(&ret, 0, 0);
> + for (i=beg.v.int_val; i<=end.v.int_val; ++i) {
> + double x, y;
> +
> + /* calculate f_i = f() with user defined variable i */
> + Ginteger(&udv->udv_value, i);
> + execute_at(udf->at);
> +
> + pop(&z);
> + x = real(&ret) + real(&z);
> + y = imag(&ret) + imag(&z);
> + Gcomplex(&ret, x, y);
> + }
> +
> + gpfree_string(&varname);
> +
> + push(Gcomplex(&z, real(&ret), imag(&ret)));
> +}
> +
> +
> +void
> f_lnot(union argument *arg)
> {
> struct value a;
> diff --git gnuplot-cvs/src/internal.h b/src/internal.h
> index 9bb1a3a..7fbdbf8 100644
> --- gnuplot-cvs/src/internal.h
> +++ b/src/internal.h
> @@ -55,6 +55,7 @@ void f_pushd __PROTO((union argument *x));
> void f_pop __PROTO((union argument *x));
> void f_call __PROTO((union argument *x));
> void f_calln __PROTO((union argument *x));
> +void f_sum __PROTO((union argument *x));
> void f_lnot __PROTO((union argument *x));
> void f_bnot __PROTO((union argument *x));
> void f_lor __PROTO((union argument *x));
> diff --git gnuplot-cvs/src/parse.c b/src/parse.c
> index 2a665ab..954d9fe 100644
> --- gnuplot-cvs/src/parse.c
> +++ b/src/parse.c
> @@ -89,6 +89,7 @@ static void parse_relational_expression __PROTO((void));
> static void parse_additive_expression __PROTO((void));
> static void parse_multiplicative_expression __PROTO((void));
> static void parse_unary_expression __PROTO((void));
> +static void parse_sum_expression __PROTO((void));
> static int parse_assignment_expression __PROTO((void));
> static int is_builtin_function __PROTO((int t_num));
>
> @@ -190,7 +191,8 @@ string_or_express(struct at_type **atptr)
> has_dummies = FALSE;
> for (i = 0; i< at->a_count; i++) {
> enum operators op_index = at->actions[i].index;
> - if ( op_index == PUSHD1 || op_index == PUSHD2 || op_index == PUSHD ) {
> + if ( op_index == PUSHD1 || op_index == PUSHD2 || op_index == PUSHD
> + || op_index == SUM ) {
> has_dummies = TRUE;
> break;
> }
> @@ -488,6 +490,8 @@ parse_primary_expression()
> c_token++;
> add_action(call_type)->udf_arg = add_udf(tok);
> }
> + } else if (equals(c_token, "sum")) {
> + parse_sum_expression();
> /* dummy_func==NULL is a flag to say no dummy variables active */
> } else if (dummy_func) {
> if (equals(c_token, c_dummy_var[0])) {
> @@ -834,6 +838,112 @@ parse_unary_expression()
> parse_primary_expression();
> }
>
> +
> +/* create action code for 'sum' expressions */
> +static void
> +parse_sum_expression()
> +{
> + /* Design: Use a user defined variable (udv) as iterator variable (k). The
> + * original idea was to treat the expression after the range as a function
> + * f(k). Consider 'g(x) = sum [k=1:4] f(k)', there are two dummy variables
> + * 'x' and 'k' from different functions 'g' and 'f' which cannot be handled
> + * by the parser. */
> +
> + char *errormsg = "Expecting 'sum [<var> =<start>:<end>]'\n";
> + char *varname = NULL;
> + union argument *arg;
> + struct udft_entry *udf;
> +
> + struct at_type * save_at;
> + int save_at_size;
> + int i;
> +
> + if (!equals(c_token, "sum"))
> + return;
> + c_token++;
> +
> + if (!equals(c_token, "["))
> + int_error(c_token, errormsg);
> + c_token++;
> +
> + /*<var> */
> + if (!isletter(c_token))
> + int_error(c_token, errormsg);
> + /* create a user defined variable and pass it to f_sum via the action
> + * table, since the argument of f_sum is already used by the udf */
> + m_capture(&varname, c_token, c_token);
> + add_udv(c_token);
> + arg = add_action(PUSHC);
> + Gstring(&(arg->v_arg), varname);
> + c_token++;
> +
> + if (!equals(c_token, "="))
> + int_error(c_token, errormsg);
> + c_token++;
> +
> + /*<start> */
> + if (!isanumber(c_token))
> + int_error(c_token, errormsg);
> + arg = add_action(PUSHC);
> + convert(&(arg->v_arg), c_token);
> + c_token++;
> +
> + if (!equals(c_token, ":"))
> + int_error(c_token, errormsg);
> + c_token++;
> +
> + /*<end> */
> + if (!isanumber(c_token))
> + int_error(c_token, errormsg);
> + arg = add_action(PUSHC);
> + convert(&(arg->v_arg), c_token);
> + c_token++;
> +
> + /* TODO add increment */
> + if (!equals(c_token, "]"))
> + int_error(c_token, errormsg);
> + c_token++;
> +
> + /* parse the next expression and convert it to an action table. */
> + /* save environment to restart parsing */
> + save_at = at;
> + save_at_size = at_size;
> +
> + at = (struct at_type *) gp_alloc(sizeof(struct at_type), "action table");
> + at->a_count = 0;
> + /* taken from temp_at()
> + * XXX Why is it necessary to reset the action table? Shouldn't it be
> + * either sizeof(struct at_type) or a_count = 0? */
> + memset(at, 0, sizeof(*at));
> + at_size = MAX_AT_LEN;
> +
> + /* Q: Do I have to save and restore parse_recursion_level?
> + * A: parse_recursion_level is used to abort parsing after the strings
> + * ('-\pi', '-\pi/2') in ('-\pi' -pi, '-\pi/2' -pi/2.), otherwise it would
> + * try to subtract pi from the string '-\pi'. This is only in effect for
> + * parse_recursion_level == 1 and string_result_only == true. My conclusion
> + * is thus to not touch parse_recursion_level. */
> + parse_expression();
> +
> + /* save action table in a user defined function */
> + udf = (struct udft_entry *) gp_alloc(sizeof(struct udft_entry), "sum");
> + udf->next_udf = (struct udft_entry *) NULL;
> + udf->udf_name = NULL; /* TODO maybe add a name and definition */
> + udf->at = at;
> + udf->definition = NULL;
> + udf->dummy_num = 0;
> + for (i = 0; i< MAX_NUM_VAR; i++)
> + (void) Ginteger(&(udf->dummy_values[i]), 0);
> +
> + /* restore environment */
> + at = save_at;
> + at_size = save_at_size;
> +
> + /* pass the udf to f_sum using the argument */
> + add_action(SUM)->udf_arg = udf;
> +}
> +
> +
> /* find or add value and return pointer */
> struct udvt_entry *
> add_udv(int t_num)
>
The idea is interesting but this kind of ad hoc , special case addition
to syntax can only end up as a mess in the long run. As Hans pointed out
gnuplot is not supposed to be a full blown programming language. If that
ever becomes a desire it will require top-down design, it can't happen
by evolution.
If that is attempted it will reach a breaking point where everything
needs restructuring and the conviction gnuplot has towards backwards
compatibility will have to be broken.
This is really trying to add a structured programming feature without
structured programming. Gnuplot does not even have a proper if-then-else
syntax yet.
If the latter was to be implemented fully it would probably require
adding the structured programming infrastructure that would make adding
this "sum" idea in a clean and coherent way fairly easy.
But adding that infrastructure will need some careful design.
regards, Peter.
|
