2 * exp-deriv.c : Expression derivation
4 * Copyright (C) 2016 Morten Welinder (terra@gnome.org)
6 * This program is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU General Public License as
8 * published by the Free Software Foundation; either version 2 of the
9 * License, or (at your option) version 3.
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
21 #include <gnumeric-config.h>
22 #include <glib/gi18n-lib.h>
24 #include <expr-deriv.h>
25 #include <expr-impl.h>
33 /* ------------------------------------------------------------------------- */
35 struct GnmExprDeriv_
{
41 * gnm_expr_deriv_info_new:
43 * Returns: (transfer full): A new #GnmExprDeriv.
46 gnm_expr_deriv_info_new (void)
48 GnmExprDeriv
*res
= g_new0 (GnmExprDeriv
, 1);
54 * gnm_expr_deriv_info_unref:
55 * @deriv: (transfer full) (nullable): #GnmExprDeriv
58 gnm_expr_deriv_info_unref (GnmExprDeriv
*deriv
)
60 if (!deriv
|| deriv
->ref_count
-- > 1)
66 * gnm_expr_deriv_info_ref:
67 * @deriv: (transfer none) (nullable): #GnmExprDeriv
69 * Returns: (transfer full) (nullable): a new reference to @deriv.
72 gnm_expr_deriv_info_ref (GnmExprDeriv
*deriv
)
80 gnm_expr_deriv_info_get_type (void)
85 t
= g_boxed_type_register_static ("GnmExprDeriv",
86 (GBoxedCopyFunc
)gnm_expr_deriv_info_ref
,
87 (GBoxedFreeFunc
)gnm_expr_deriv_info_unref
);
93 * gnm_expr_deriv_info_set_var:
94 * @deriv: #GnmExprDeriv
95 * @var: (transfer none): location of variable
98 gnm_expr_deriv_info_set_var (GnmExprDeriv
*deriv
, GnmEvalPos
const *var
)
103 /* ------------------------------------------------------------------------- */
105 static GnmExpr
const *
106 gnm_value_deriv (GnmValue
const *v
)
108 if (VALUE_IS_NUMBER (v
))
109 return gnm_expr_new_constant (value_new_float (0));
114 static GnmExpr
const *madd (GnmExpr
const *l
, gboolean copyl
, GnmExpr
const *r
, gboolean copyr
);
115 static GnmExpr
const *msub (GnmExpr
const *l
, gboolean copyl
, GnmExpr
const *r
, gboolean copyr
);
116 static GnmExpr
const *mmul (GnmExpr
const *l
, gboolean copyl
, GnmExpr
const *r
, gboolean copyr
);
117 static GnmExpr
const *mdiv (GnmExpr
const *l
, gboolean copyl
, GnmExpr
const *r
, gboolean copyr
);
118 static GnmExpr
const *mneg (GnmExpr
const *l
, gboolean copyl
);
119 static GnmExpr
const *optimize_sum (GnmExpr
const *e
);
125 is_any_const (GnmExpr
const *e
, gnm_float
*c
)
127 GnmValue
const *v
= gnm_expr_get_constant (e
);
128 if (v
&& VALUE_IS_FLOAT (v
)) {
129 if (c
) *c
= value_get_as_float (v
);
136 is_const (GnmExpr
const *e
, gnm_float c
)
138 GnmValue
const *v
= gnm_expr_get_constant (e
);
139 return v
&& VALUE_IS_FLOAT (v
) && value_get_as_float (v
) == c
;
143 is_neg (GnmExpr
const *e
)
145 return (GNM_EXPR_GET_OPER (e
) == GNM_EXPR_OP_UNARY_NEG
);
149 is_lcmul (GnmExpr
const *e
, gnm_float
*c
)
151 return (GNM_EXPR_GET_OPER (e
) == GNM_EXPR_OP_MULT
&&
152 is_any_const (e
->binary
.value_a
, c
));
156 // Optimizing constructor for "+". Takes ownership of "l" and "r"
157 // if the corresponding "copy" argument is false.
159 // Note, that this plays fast-and-loose with semantics when errors are
161 static GnmExpr
const *
162 madd (GnmExpr
const *l
, gboolean copyl
, GnmExpr
const *r
, gboolean copyr
)
164 if (is_const (l
, 0)) {
165 if (!copyl
) gnm_expr_free (l
);
166 if (copyr
) r
= gnm_expr_copy (r
);
170 if (is_const (r
, 0)) {
171 if (!copyr
) gnm_expr_free (r
);
172 if (copyl
) l
= gnm_expr_copy (l
);
176 if (copyl
) l
= gnm_expr_copy (l
);
177 if (copyr
) r
= gnm_expr_copy (r
);
178 return gnm_expr_new_binary (l
, GNM_EXPR_OP_ADD
, r
);
181 // Optimizing constructor for unary "-". Takes ownership of "l"
182 // if the corresponding "copy" argument is false.
184 // Note, that this plays fast-and-loose with semantics when errors are
186 static GnmExpr
const *
187 mneg (GnmExpr
const *l
, gboolean copyl
)
190 if (is_any_const (l
, &x
)) {
191 if (!copyl
) gnm_expr_free (l
);
192 return gnm_expr_new_constant (value_new_float (-x
));
195 if (is_lcmul (l
, &x
)) {
196 GnmExpr
const *res
= mmul (gnm_expr_new_constant (value_new_float (-x
)), 0,
197 l
->binary
.value_b
, 1);
198 if (!copyl
) gnm_expr_free (l
);
202 if (copyl
) l
= gnm_expr_copy (l
);
203 return gnm_expr_new_unary (GNM_EXPR_OP_UNARY_NEG
, l
);
206 // Optimizing constructor for "-". Takes ownership of "l" and "r"
207 // if the corresponding "copy" argument is false.
209 // Note, that this plays fast-and-loose with semantics when errors are
211 static GnmExpr
const *
212 msub (GnmExpr
const *l
, gboolean copyl
, GnmExpr
const *r
, gboolean copyr
)
214 if (is_const (r
, 0)) {
215 if (!copyr
) gnm_expr_free (r
);
216 if (copyl
) l
= gnm_expr_copy (l
);
220 if (is_const (l
, 0)) {
221 if (!copyl
) gnm_expr_free (l
);
222 return mneg (r
, copyr
);
225 if (copyl
) l
= gnm_expr_copy (l
);
226 if (copyr
) r
= gnm_expr_copy (r
);
227 return gnm_expr_new_binary (l
, GNM_EXPR_OP_SUB
, r
);
230 // Optimizing constructor for "*". Takes ownership of "l" and "r"
231 // if the corresponding "copy" argument is false.
233 // Note, that this plays fast-and-loose with semantics when errors are
235 static GnmExpr
const *
236 mmul (GnmExpr
const *l
, gboolean copyl
, GnmExpr
const *r
, gboolean copyr
)
238 if (is_const (l
, 1) || is_const (r
, 0)) {
239 if (!copyl
) gnm_expr_free (l
);
240 if (copyr
) r
= gnm_expr_copy (r
);
244 if (is_const (l
, 0) || is_const (r
, 1)) {
245 if (!copyr
) gnm_expr_free (r
);
246 if (copyl
) l
= gnm_expr_copy (l
);
250 if (is_const (l
, -1)) {
251 if (!copyl
) gnm_expr_free (l
);
252 return mneg (r
, copyr
);
256 GnmExpr
const *res
= mneg (mmul (l
->unary
.value
, 1, r
, copyr
), 0);
257 if (!copyl
) gnm_expr_free (l
);
262 GnmExpr
const *res
= mneg (mmul (l
, copyl
, r
->unary
.value
, 1), 0);
263 if (!copyr
) gnm_expr_free (r
);
267 if (is_lcmul (l
, NULL
)) {
268 GnmExpr
const *res
= mmul (l
->binary
.value_a
, 1,
269 mmul (l
->binary
.value_b
, 1,
271 if (!copyl
) gnm_expr_free (l
);
275 if (copyl
) l
= gnm_expr_copy (l
);
276 if (copyr
) r
= gnm_expr_copy (r
);
277 return gnm_expr_new_binary (l
, GNM_EXPR_OP_MULT
, r
);
280 // Optimizing constructor for "/". Takes ownership of "l" and "r"
281 // if the corresponding "copy" argument is false.
283 // Note, that this plays fast-and-loose with semantics when errors are
285 static GnmExpr
const *
286 mdiv (GnmExpr
const *l
, gboolean copyl
, GnmExpr
const *r
, gboolean copyr
)
288 if (is_const (l
, 0) || is_const (r
, 1)) {
289 if (!copyr
) gnm_expr_free (r
);
290 if (copyl
) l
= gnm_expr_copy (l
);
294 if (copyl
) l
= gnm_expr_copy (l
);
295 if (copyr
) r
= gnm_expr_copy (r
);
296 return gnm_expr_new_binary (l
, GNM_EXPR_OP_DIV
, r
);
299 // Optimizing constructor for "^". Takes ownership of "l" and "r"
300 // if the corresponding "copy" argument is false.
302 // Note, that this plays fast-and-loose with semantics when errors are
304 static GnmExpr
const *
305 mexp (GnmExpr
const *l
, gboolean copyl
, GnmExpr
const *r
, gboolean copyr
)
307 if (is_const (r
, 1)) {
308 if (!copyr
) gnm_expr_free (r
);
309 if (copyl
) l
= gnm_expr_copy (l
);
313 if (copyl
) l
= gnm_expr_copy (l
);
314 if (copyr
) r
= gnm_expr_copy (r
);
315 return gnm_expr_new_binary (l
, GNM_EXPR_OP_EXP
, r
);
318 static GnmExpr
const *
319 msum (GnmExprList
*as
)
321 GnmFunc
*fsum
= gnm_func_lookup_or_add_placeholder ("SUM");
322 GnmExpr
const *res
= gnm_expr_new_funcall (fsum
, as
);
323 GnmExpr
const *opt
= optimize_sum (res
);
333 static GnmExpr
const *
334 optimize_sum (GnmExpr
const *e
)
336 int argc
= e
->func
.argc
;
337 GnmExprConstPtr
*argv
= e
->func
.argv
;
338 gboolean all_neg
= (argc
> 0);
339 gboolean all_lcmul
= (argc
> 0);
343 for (i
= 0; i
< argc
; i
++) {
344 GnmExpr
const *a
= argv
[i
];
347 all_neg
= all_neg
&& is_neg (a
);
349 all_lcmul
= all_lcmul
&&
350 is_lcmul (a
, i
? &x
: &cl
) &&
355 GnmExprList
*as
= NULL
;
356 for (i
= argc
; i
-- > 0;) {
357 GnmExpr
const *a
= argv
[i
];
358 as
= g_slist_prepend (as
, (gpointer
)gnm_expr_copy (a
->unary
.value
));
360 return mneg (msum (as
), 0);
364 GnmExprList
*as
= NULL
;
365 for (i
= argc
; i
-- > 0;) {
366 GnmExpr
const *a
= argv
[i
];
367 as
= g_slist_prepend (as
, (gpointer
)gnm_expr_copy (a
->binary
.value_b
));
369 return mmul (gnm_expr_new_constant (value_new_float (cl
)), 0,
376 static GnmExpr
const *
377 optimize (GnmExpr
const *e
)
379 GnmExprOp op
= GNM_EXPR_GET_OPER (e
);
382 case GNM_EXPR_OP_FUNCALL
: {
383 GnmFunc
*f
= gnm_expr_get_func_def (e
);
384 GnmFunc
*fsum
= gnm_func_lookup_or_add_placeholder ("SUM");
387 return optimize_sum (e
);
395 /* ------------------------------------------------------------------------- */
397 struct cb_arg_collect
{
399 GnmCellRef
const *cr0
;
400 GnmEvalPos
const *ep
;
404 cb_arg_collect (GnmCellIter
const *iter
, gpointer user_
)
406 struct cb_arg_collect
*user
= user_
;
407 GnmCell
const *cell
= iter
->cell
;
410 gnm_cellref_init (&cr
, user
->cr0
->sheet
,
411 cell
->pos
.col
, cell
->pos
.row
,
413 parse_pos_init_evalpos (&pp
, user
->ep
);
414 gnm_cellref_set_col_ar (&cr
, &pp
, user
->cr0
->col_relative
);
415 gnm_cellref_set_row_ar (&cr
, &pp
, user
->cr0
->row_relative
);
416 user
->args
= gnm_expr_list_prepend
418 (gpointer
)gnm_expr_new_cellref (&cr
));
423 * gnm_expr_deriv_collect:
425 * @ep: evaluation position
426 * @info: extra information, not currently used
428 * Returns: (type GSList) (transfer full) (element-type GnmExpr): list of
429 * expressions expanded from @expr
432 gnm_expr_deriv_collect (GnmExpr
const *expr
,
433 GnmEvalPos
const *ep
,
436 struct cb_arg_collect user
;
441 for (i
= 0; i
< expr
->func
.argc
; i
++) {
442 GnmExpr
const *e
= expr
->func
.argv
[i
];
443 GnmValue
const *v
= gnm_expr_get_constant (e
);
445 if (!v
|| !VALUE_IS_CELLRANGE (v
)) {
446 user
.args
= gnm_expr_list_prepend
447 (user
.args
, (gpointer
)gnm_expr_copy (e
));
451 user
.cr0
= &value_get_rangeref (v
)->a
;
452 workbook_foreach_cell_in_range (ep
, v
,
453 CELL_ITER_IGNORE_BLANK
,
458 return g_slist_reverse (user
.args
);
461 /* ------------------------------------------------------------------------- */
463 #define MAYBE_FREE(e) do { if (e) gnm_expr_free (e); } while (0)
465 #define COMMON_BINARY_START \
466 GnmExpr const *a = expr->binary.value_a; /* Not owned */ \
467 GnmExpr const *da = gnm_expr_deriv (a, ep, info); \
468 GnmExpr const *b = expr->binary.value_b; /* Not owned */ \
469 GnmExpr const *db = gnm_expr_deriv (b, ep, info); \
476 #define COMMON_BINARY_END }
481 * @ep: position for @expr
482 * @info: Derivative information
484 * Returns: (transfer full) (nullable): the derivative of @expr with respect
488 gnm_expr_deriv (GnmExpr
const *expr
,
489 GnmEvalPos
const *ep
,
492 GnmExprOp op
= GNM_EXPR_GET_OPER (expr
);
495 case GNM_EXPR_OP_RANGE_CTOR
:
496 case GNM_EXPR_OP_INTERSECT
:
497 case GNM_EXPR_OP_NAME
:
498 case GNM_EXPR_OP_ARRAY_CORNER
:
499 case GNM_EXPR_OP_ARRAY_ELEM
:
500 case GNM_EXPR_OP_SET
:
502 case GNM_EXPR_OP_EQUAL
:
505 case GNM_EXPR_OP_GTE
:
506 case GNM_EXPR_OP_LTE
:
507 case GNM_EXPR_OP_NOT_EQUAL
:
508 case GNM_EXPR_OP_CAT
:
509 case GNM_EXPR_OP_PERCENTAGE
:
513 case GNM_EXPR_OP_PAREN
:
514 case GNM_EXPR_OP_UNARY_PLUS
:
515 return gnm_expr_deriv (expr
->unary
.value
, ep
, info
);
517 case GNM_EXPR_OP_UNARY_NEG
: {
518 GnmExpr
const *d
= gnm_expr_deriv (expr
->unary
.value
, ep
, info
);
519 return d
? mneg (d
, 0) : NULL
;
522 case GNM_EXPR_OP_ADD
: {
524 return madd (da
, 0, db
, 0);
528 case GNM_EXPR_OP_SUB
: {
530 return msub (da
, 0, db
, 0);
534 case GNM_EXPR_OP_MULT
: {
536 GnmExpr
const *t1
= mmul (da
, 0, b
, 1);
537 GnmExpr
const *t2
= mmul (a
, 1, db
, 0);
538 return madd (t1
, 0, t2
, 0);
542 case GNM_EXPR_OP_DIV
: {
544 GnmExpr
const *t1
= mmul (da
, 0, b
, 1);
545 GnmExpr
const *t2
= mmul (a
, 1, db
, 0);
546 GnmExpr
const *d
= msub (t1
, 0, t2
, 0);
547 GnmExpr
const *n
= mmul (b
, 1, b
, 1);
548 return mdiv (d
, 0, n
, 0);
552 case GNM_EXPR_OP_EXP
: {
554 GnmFunc
*fln
= gnm_func_lookup ("ln", NULL
);
556 if (is_any_const (b
, &cb
)) {
557 GnmExpr
const *bm1
= gnm_expr_new_constant (value_new_float (cb
- 1));
558 GnmExpr
const *t1
= mexp (a
, 1, bm1
, 0);
560 return mmul (mmul (b
, 1, t1
, 0), 0, da
, 0);
562 // a^b = exp(b*log(a))
563 // (a^b)' = a^b * (a'*b/a + b'*ln(a))
564 GnmExpr
const *t1
= mdiv (mmul (da
, 0, b
, 1), 0, a
, 1);
565 GnmExpr
const *t2
= mmul
567 gnm_expr_new_funcall1 (fln
, gnm_expr_copy (a
)), 0);
568 GnmExpr
const *s
= madd (t1
, 0, t2
, 0);
569 return mmul (expr
, 1, s
, 0);
578 case GNM_EXPR_OP_FUNCALL
: {
579 GnmFunc
*f
= gnm_expr_get_func_def (expr
);
580 GnmExpr
const *res
= gnm_func_derivative (f
, expr
, ep
, info
);
581 GnmExpr
const *opt
= res
? optimize (res
) : NULL
;
589 case GNM_EXPR_OP_CONSTANT
:
590 return gnm_value_deriv (expr
->constant
.value
);
592 case GNM_EXPR_OP_CELLREF
: {
598 GnmExprTop
const *texpr
;
599 GnmExprTop
const *texpr2
;
600 GnmExprRelocateInfo rinfo
;
602 gnm_cellref_make_abs (&r
, &expr
->cellref
.ref
, ep
);
603 sheet
= eval_sheet (r
.sheet
, ep
->sheet
);
605 if (sheet
== info
->var
.sheet
&&
606 r
.col
== info
->var
.eval
.col
&&
607 r
.row
== info
->var
.eval
.row
)
608 return gnm_expr_new_constant (value_new_float (1));
610 cell
= sheet_cell_get (sheet
, r
.col
, r
.row
);
612 return gnm_expr_new_constant (value_new_float (0));
613 if (!gnm_cell_has_expr (cell
))
614 return gnm_value_deriv (cell
->value
);
616 eval_pos_init_cell (&ep2
, cell
);
617 res
= gnm_expr_deriv (cell
->base
.texpr
->expr
, &ep2
, info
);
621 // The just-computed derivative is relative to the wrong
624 texpr
= gnm_expr_top_new (res
);
625 parse_pos_init_evalpos (&rinfo
.pos
, &ep2
);
626 rinfo
.reloc_type
= GNM_EXPR_RELOCATE_MOVE_RANGE
;
627 rinfo
.origin
.start
= rinfo
.origin
.end
= ep2
.eval
;
628 rinfo
.origin_sheet
= ep2
.sheet
;
629 rinfo
.target_sheet
= ep
->sheet
;
630 rinfo
.col_offset
= ep
->eval
.col
- ep2
.eval
.col
;
631 rinfo
.row_offset
= ep
->eval
.row
- ep2
.eval
.row
;
632 texpr2
= gnm_expr_top_relocate (texpr
, &rinfo
, FALSE
);
635 res
= gnm_expr_copy (texpr2
->expr
);
636 gnm_expr_top_unref (texpr2
);
638 res
= gnm_expr_copy (texpr
->expr
);
640 gnm_expr_top_unref (texpr
);
645 #ifndef DEBUG_SWITCH_ENUM
647 g_assert_not_reached ();
653 /* ------------------------------------------------------------------------- */
656 * gnm_expr_top_deriv:
658 * @ep: Evaluation position
659 * @info: Derivative information
661 * Returns: (transfer full) (nullable): The derivative of @texpr with
665 gnm_expr_top_deriv (GnmExprTop
const *texpr
,
666 GnmEvalPos
const *ep
,
671 g_return_val_if_fail (GNM_IS_EXPR_TOP (texpr
), NULL
);
672 g_return_val_if_fail (ep
!= NULL
, NULL
);
673 g_return_val_if_fail (info
!= NULL
, NULL
);
675 expr
= gnm_expr_deriv (texpr
->expr
, ep
, info
);
676 if (gnm_debug_flag ("deriv")) {
677 GnmParsePos pp
, ppvar
;
679 Sheet
*sheet
= ep
->sheet
;
680 GnmConventions
const *convs
= sheet_get_conventions (sheet
);
682 parse_pos_init_evalpos (&ppvar
, &info
->var
);
683 parse_pos_init_evalpos (&pp
, ep
);
685 s
= gnm_expr_top_as_string (texpr
, &pp
, convs
);
686 g_printerr ("Derivative of %s with respect to %s:%s",
687 s
, parsepos_as_string (&ppvar
),
688 expr
? "\n" : " cannot compute.\n");
691 s
= gnm_expr_as_string (expr
, &pp
, convs
);
692 g_printerr ("%s\n\n", s
);
697 return gnm_expr_top_new (expr
);
701 * gnm_expr_cell_deriv:
705 * Returns: (transfer full) (nullable): The derivative of cell @y with
706 * respect to cell @x.
709 gnm_expr_cell_deriv (GnmCell
*y
, GnmCell
*x
)
711 GnmExprTop
const *res
;
715 g_return_val_if_fail (y
!= NULL
, NULL
);
716 g_return_val_if_fail (gnm_cell_has_expr (y
), NULL
);
717 g_return_val_if_fail (x
!= NULL
, NULL
);
719 eval_pos_init_cell (&ep
, y
);
721 info
= gnm_expr_deriv_info_new ();
722 eval_pos_init_cell (&var
, x
);
723 gnm_expr_deriv_info_set_var (info
, &var
);
725 res
= gnm_expr_top_deriv (y
->base
.texpr
, &ep
, info
);
727 gnm_expr_deriv_info_unref (info
);
733 * gnm_expr_cell_deriv_value:
737 * Returns: The derivative of cell @y with respect to cell @x at the
738 * current value of @x. Returns NaN on error.
741 gnm_expr_cell_deriv_value (GnmCell
*y
, GnmCell
*x
)
743 GnmExprTop
const *dydx
;
748 g_return_val_if_fail (y
!= NULL
, gnm_nan
);
749 g_return_val_if_fail (x
!= NULL
, gnm_nan
);
751 dydx
= gnm_expr_cell_deriv (y
, x
);
755 eval_pos_init_cell (&ep
, y
);
756 v
= gnm_expr_top_eval (dydx
, &ep
, GNM_EXPR_EVAL_SCALAR_NON_EMPTY
);
757 res
= VALUE_IS_NUMBER (v
) ? value_get_as_float (v
) : gnm_nan
;
760 gnm_expr_top_unref (dydx
);
765 /* ------------------------------------------------------------------------- */
768 * gnm_expr_deriv_chain:
769 * @expr: #GnmExpr for a function call with one argument
770 * @deriv: (transfer full) (nullable): Derivative of @expr's function.
771 * @ep: position for @expr
772 * @info: Derivative information
774 * Applies the chain rule to @expr.
776 * Returns: (transfer full) (nullable): the derivative of @expr with respect
780 gnm_expr_deriv_chain (GnmExpr
const *expr
,
781 GnmExpr
const *deriv
,
782 GnmEvalPos
const *ep
,
785 GnmExpr
const *deriv2
;
790 deriv2
= gnm_expr_deriv (gnm_expr_get_func_arg (expr
, 0), ep
, info
);
792 gnm_expr_free (deriv
);
796 return mmul (deriv
, 0, deriv2
, 0);
799 /* ------------------------------------------------------------------------- */
802 gnm_expr_deriv_shutdown_ (void)
806 /* ------------------------------------------------------------------------- */