Menu
▾
▴
Re: finite summation - sum [k=1:4] k
|
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.
|
