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