Graph Framework
Loading...
Searching...
No Matches
arithmetic.hpp
Go to the documentation of this file.
1//------------------------------------------------------------------------------
6//------------------------------------------------------------------------------
7
8#ifndef arithmetic_h
9#define arithmetic_h
10
11#include "node.hpp"
12
13namespace graph {
14//------------------------------------------------------------------------------
23//------------------------------------------------------------------------------
24 template<jit::float_scalar T, bool SAFE_MATH=false>
25 bool is_constant_combinable(shared_leaf<T, SAFE_MATH> a,
26 shared_leaf<T, SAFE_MATH> b) {
27 if (a->is_constant() && b->is_constant()) {
28 auto a1 = piecewise_1D_cast(a);
29 auto a2 = piecewise_2D_cast(a);
30 auto b2 = piecewise_2D_cast(b);
31
32 return constant_cast(a).get() ||
33 constant_cast(b).get() ||
34 (a1.get() && a1->is_arg_match(b)) ||
35 (a2.get() && a2->is_arg_match(b)) ||
36 (a2.get() && (a2->is_row_match(b) || a2->is_col_match(b))) ||
37 (b2.get() && (b2->is_row_match(a) || b2->is_col_match(a)));
38 }
39 return false;
40 }
41
42//------------------------------------------------------------------------------
51//------------------------------------------------------------------------------
52 template<jit::float_scalar T, bool SAFE_MATH=false>
53 bool is_constant_promotable(shared_leaf<T, SAFE_MATH> a,
54 shared_leaf<T, SAFE_MATH> b) {
55 auto b1 = piecewise_1D_cast(b);
56 auto b2 = piecewise_2D_cast(b);
57
58 return a->is_constant() &&
59 (!b->is_constant() ||
60 (constant_cast(a).get() && (b1.get() || b2.get())) ||
61 (piecewise_1D_cast(a).get() && b2.get()));
62 }
63
64//------------------------------------------------------------------------------
73//------------------------------------------------------------------------------
74 template<jit::float_scalar T, bool SAFE_MATH=false>
75 bool is_variable_combinable(shared_leaf<T, SAFE_MATH> a,
76 shared_leaf<T, SAFE_MATH> b) {
77 return a->is_power_base_match(b);
78 }
79
80//------------------------------------------------------------------------------
89//------------------------------------------------------------------------------
90 template<jit::float_scalar T, bool SAFE_MATH=false>
91 bool is_variable_promotable(shared_leaf<T, SAFE_MATH> a,
92 shared_leaf<T, SAFE_MATH> b) {
93 return !b->is_constant() &&
94 (a->is_all_variables() &&
95 (!b->is_all_variables() ||
96 (b->is_all_variables() &&
97 a->get_complexity() < b->get_complexity())));
98 }
99
100//------------------------------------------------------------------------------
109//------------------------------------------------------------------------------
110 template<jit::float_scalar T, bool SAFE_MATH=false>
111 bool is_greater_exponent(shared_leaf<T, SAFE_MATH> a,
112 shared_leaf<T, SAFE_MATH> b) {
113 auto ae = constant_cast(a->get_power_exponent());
114 auto be = constant_cast(b->get_power_exponent());
115
116 return ae.get() && be.get() &&
117 std::abs(ae->evaluate().at(0)) > std::abs(be->evaluate().at(0));
118 }
119
120//******************************************************************************
121// Add node.
122//******************************************************************************
123//------------------------------------------------------------------------------
130//------------------------------------------------------------------------------
131 template<jit::float_scalar T, bool SAFE_MATH=false>
132 class add_node final : public branch_node<T, SAFE_MATH> {
133 private:
134//------------------------------------------------------------------------------
140//------------------------------------------------------------------------------
141 static std::string to_string(leaf_node<T, SAFE_MATH> *l,
142 leaf_node<T, SAFE_MATH> *r) {
143 return jit::format_to_string(reinterpret_cast<size_t> (l)) + "+" +
144 jit::format_to_string(reinterpret_cast<size_t> (r));
145 }
146
147 public:
148//------------------------------------------------------------------------------
153//------------------------------------------------------------------------------
158
159//------------------------------------------------------------------------------
165//------------------------------------------------------------------------------
166 virtual backend::buffer<T> evaluate() {
167 backend::buffer<T> l_result = this->left->evaluate();
168 backend::buffer<T> r_result = this->right->evaluate();
169 return l_result + r_result;
170 }
171
172//------------------------------------------------------------------------------
176//------------------------------------------------------------------------------
177 virtual shared_leaf<T, SAFE_MATH> reduce() {
178// Constant reductions.
179 auto l = constant_cast(this->left);
180 auto r = constant_cast(this->right);
181
182 if (l.get() && l->is(0)) {
183 return this->right;
184 } else if (r.get() && r->is(0)) {
185 return this->left;
186 } else if (l.get() && r.get()) {
187 return constant<T, SAFE_MATH> (this->evaluate());
188 } else if (r.get() && !l.get()) {
189 return this->right + this->left;
190 }
191
192 auto pl1 = piecewise_1D_cast(this->left);
193 auto pr1 = piecewise_1D_cast(this->right);
194
195 if (pl1.get() && (r.get() || pl1->is_arg_match(this->right))) {
196 return piecewise_1D(this->evaluate(), pl1->get_arg());
197 } else if (pr1.get() && (l.get() || pr1->is_arg_match(this->left))) {
198 return piecewise_1D(this->evaluate(), pr1->get_arg());
199 }
200
201 auto pl2 = piecewise_2D_cast(this->left);
202 auto pr2 = piecewise_2D_cast(this->right);
203
204 if (pl2.get() && (r.get() || pl2->is_arg_match(this->right))) {
205 return piecewise_2D(this->evaluate(),
206 pl2->get_num_columns(),
207 pl2->get_left(),
208 pl2->get_right());
209 } else if (pr2.get() && (l.get() || pr2->is_arg_match(this->left))) {
210 return piecewise_2D(this->evaluate(),
211 pr2->get_num_columns(),
212 pr2->get_left(),
213 pr2->get_right());
214 }
215
216// Combine 2D and 1D piecewise constants if a row or column matches.
217 if (pr2.get() && pr2->is_row_match(this->left)) {
218 backend::buffer<T> result = pl1->evaluate();
219 result.add_row(pr2->evaluate());
220 return piecewise_2D(result,
221 pr2->get_num_columns(),
222 pr2->get_left(),
223 pr2->get_right());
224 } else if (pr2.get() && pr2->is_col_match(this->left)) {
225 backend::buffer<T> result = pl1->evaluate();
226 result.add_col(pr2->evaluate());
227 return piecewise_2D(result,
228 pr2->get_num_columns(),
229 pr2->get_left(),
230 pr2->get_right());
231 } else if (pl2.get() && pl2->is_row_match(this->right)) {
232 backend::buffer<T> result = pl2->evaluate();
233 result.add_row(pr1->evaluate());
234 return piecewise_2D(result,
235 pl2->get_num_columns(),
236 pl2->get_left(),
237 pl2->get_right());
238 } else if (pl2.get() && pl2->is_col_match(this->right)) {
239 backend::buffer<T> result = pl2->evaluate();
240 result.add_col(pr1->evaluate());
241 return piecewise_2D(result,
242 pl2->get_num_columns(),
243 pl2->get_left(),
244 pl2->get_right());
245 }
246
247// Identity reductions.
248 if (this->left->is_match(this->right)) {
249 return 2.0*this->left;
250 }
251
252// Common factor reduction. If the left and right are both multiply nodes check
253// for a common factor. So you can change a*b + a*c -> a*(b + c).
254 auto lm = multiply_cast(this->left);
255 auto rm = multiply_cast(this->right);
256
257// v1 + -c*v2 -> v1 - c*v2
258// -c*v1 + v2 -> v2 - c*v1
259 if (rm.get() &&
260 rm->get_left()->is_constant() &&
261 rm->get_left()->evaluate().is_negative()) {
262 return this->left - (-this->right);
263 } else if (rm.get() &&
264 rm->get_left()->is_constant() &&
265 rm->get_left()->evaluate().is_negative()) {
266 return this->right - (-this->left);
267 }
268
269// a*b + c -> fma(a,b,c)
270// a + b*c -> fma(b,c,a)
271 if (lm.get()) {
272 return fma(lm->get_left(), lm->get_right(), this->right);
273 } else if (rm.get()) {
274 return fma(rm->get_left(), rm->get_right(), this->left);
275 }
276
277// Common denominator reduction. If the left and right are both divide nodes
278// for a common denominator. So you can change a/b + c/b -> (a + c)/d.
279 auto ld = divide_cast(this->left);
280 auto rd = divide_cast(this->right);
281
282// c is a constant.
283// a + -c/b -> a - c/b
284// a + (-c*d)/b -> a - (c*d)/b
285// -c/a + b -> b - c/a
286// (-c*d)/a + b -> b - (c*d)/a
287 if (rd.get()) {
288 auto rdlm = multiply_cast(rd->get_left());
289 if ((rd->get_left()->is_constant() &&
290 rd->get_left()->evaluate().is_negative()) ||
291 (rdlm.get() &&
292 (rdlm->get_left()->is_constant() &&
293 rdlm->get_left()->evaluate().is_negative()))) {
294 return this->left - (-rd->get_left())/rd->get_right();
295 }
296 } else if (ld.get()) {
297 auto ldlm = multiply_cast(ld->get_left());
298 if ((ld->get_left()->is_constant() &&
299 ld->get_left()->evaluate().is_negative()) ||
300 (ldlm.get() &&
301 (ldlm->get_left()->is_constant() &&
302 ldlm->get_left()->evaluate().is_negative()))) {
303 return this->right - (-ld->get_left())/ld->get_right();
304 }
305 }
306
307 if (ld.get() && rd.get()) {
308 if (ld->get_right()->is_match(rd->get_right())) {
309 return (ld->get_left() + rd->get_left())/ld->get_right();
310 }
311
312 auto ldlm = multiply_cast(ld->get_left());
313 auto rdlm = multiply_cast(rd->get_left());
314// a/b - c*a/d -> (1/b - c/d)*a
315// a/b - a*c/d -> (1/b - c/d)*a
316// c*a/b - a/d -> (c/b - 1/d)*a
317// a*c/b - a/d -> (c/b - 1/d)*a
318 if (rdlm.get()) {
319 if (ld->get_left()->is_match(rdlm->get_left())) {
320 return (1.0/ld->get_right() +
321 rdlm->get_right()/rd->get_right())*rdlm->get_left();
322 } else if (ld->get_left()->is_match(rdlm->get_right())) {
323 return (1.0/ld->get_right() +
324 rdlm->get_left()/rd->get_right())*rdlm->get_right();
325 }
326 } else if (ldlm.get()) {
327 if (rd->get_left()->is_match(ldlm->get_left())) {
328 return (ldlm->get_right()/ld->get_right() +
329 1.0/rd->get_right())*ldlm->get_left();
330 } else if (rd->get_left()->is_match(ldlm->get_right())) {
331 return (ldlm->get_left()/ld->get_right() +
332 1.0/rd->get_right())*ldlm->get_right();
333 }
334 }
335
336// c1*a/b + c2*a/d = c3*(a/b + c4*a/d)
337// a*b/c + d*b/e -> (a/c + d/e)*b
338// Make sure we prevent combining constants when we just need to factor out a
339// common term.
340// c1*a/b + c2*a/d -> (c1/b + c2/d)*a
341 if (ldlm.get() && rdlm.get()) {
342 if (is_constant_combinable(ldlm->get_left(),
343 rdlm->get_left()) &&
344 !ldlm->get_right()->is_match(rdlm->get_right())) {
345 return (ldlm->get_right()/ld->get_right() +
346 rdlm->get_left()/ldlm->get_left() *
347 rdlm->get_right()/rd->get_right())*ldlm->get_left();
348 }
349
350 if (ldlm->get_right()->is_match(rdlm->get_right())) {
351 return (ldlm->get_left()/ld->get_right() +
352 rdlm->get_left()/rd->get_right())*ldlm->get_right();
353 } else if (ldlm->get_right()->is_match(rdlm->get_left())) {
354 return (ldlm->get_left()/ld->get_right() +
355 rdlm->get_right()/rd->get_right())*ldlm->get_right();
356 } else if (ldlm->get_left()->is_match(rdlm->get_right())) {
357 return (ldlm->get_right()/ld->get_right() +
358 rdlm->get_left()/rd->get_right())*ldlm->get_left();
359 } else if (ldlm->get_left()->is_match(rdlm->get_left())) {
360 return (ldlm->get_right()/ld->get_right() +
361 rdlm->get_right()/rd->get_right())*ldlm->get_left();
362 }
363 }
364
365// (a/(c*b) + d/(e*c)) -> (a/b + d/e)/c
366// (a/(b*c) + d/(e*c)) -> (a/b + d/e)/c
367// (a/(c*b) + d/(c*e)) -> (a/b + d/e)/c
368// (a/(b*c) + d/(c*e)) -> (a/b + d/e)/c
369 auto ldrm = multiply_cast(ld->get_right());
370 auto rdrm = multiply_cast(rd->get_right());
371 if (ldrm.get() && rdrm.get()) {
372 if (ldrm->get_right()->is_match(rdrm->get_right())) {
373 return (ld->get_left()/ldrm->get_left() +
374 rd->get_left()/rdrm->get_left())/ldrm->get_right();
375 } else if (ldrm->get_right()->is_match(rdrm->get_left())) {
376 return (ld->get_left()/ldrm->get_left() +
377 rd->get_left()/rdrm->get_right())/ldrm->get_right();
378 } else if (ldrm->get_left()->is_match(rdrm->get_right())) {
379 return (ld->get_left()/ldrm->get_right() +
380 rd->get_left()/rdrm->get_left())/ldrm->get_left();
381 } else if (ldrm->get_left()->is_match(rdrm->get_left())) {
382 return (ld->get_left()/ldrm->get_right() +
383 rd->get_left()/rdrm->get_right())/ldrm->get_left();
384 }
385 }
386
387// a/b + c/(b*d) -> (a*b + c)/(b*d)
388// a/b + c/(d*b) -> (a*b + c)/(b*d)
389// a/(b*d) + c/b -> (c*b + a)/(b*d)
390// a/(d*b) + c/b -> (c*b + a)/(b*d)
391 if (rdrm.get()) {
392 if (ld->get_right()->is_match(rdrm->get_left())) {
393 return fma(ld->get_left(),
394 rdrm->get_right(),
395 rd->get_left()) /
396 rd->get_right();
397 } else if (ld->get_right()->is_match(rdrm->get_right())) {
398 return fma(ld->get_left(),
399 rdrm->get_left(),
400 rd->get_left()) /
401 rd->get_right();
402 }
403 } else if (ldrm.get()) {
404 if (rd->get_right()->is_match(ldrm->get_left())) {
405 return fma(rd->get_left(),
406 ldrm->get_right(),
407 ld->get_left()) /
408 ld->get_right();
409 } else if (rd->get_right()->is_match(ldrm->get_right())) {
410 return fma(rd->get_left(),
411 ldrm->get_left(),
412 ld->get_left()) /
413 ld->get_right();
414 }
415 }
416 }
417
418// Chained addition reductions.
419// a + (a + b) = fma(2,a,b)
420// a + (b + a) = fma(2,a,b)
421// (a + b) + a = fma(2,a,b)
422// (b + a) + a = fma(2,a,b)
423 auto la = add_cast(this->left);
424 if (la.get()) {
425 if (this->right->is_match(la->get_left())) {
426 return fma(2.0, this->right, la->get_right());
427 } else if (this->right->is_match(la->get_right())) {
428 return fma(2.0, this->right, la->get_left());
429 }
430 }
431 auto ra = add_cast(this->right);
432 if (ra.get()) {
433 if (this->left->is_match(ra->get_left())) {
434 return fma(2.0, this->left, ra->get_right());
435 } else if (this->left->is_match(ra->get_right())) {
436 return fma(2.0, this->left, ra->get_left());
437 }
438 }
439
440// Add subtraction reduction
441// (c1 - a) + c2 -> c3 - a
442// (a - c1) + c2 -> a + c2
443// These reductions are handled by moving constants to the right.
444// (a - b) + a -> 2a - b
445// (b - a) + a -> b
446 auto ls = subtract_cast(this->left);
447 if (ls.get()) {
448 if (ls->get_left()->is_match(this->right)) {
449 return static_cast<T> (2.0)*this->right - ls->get_right();
450 } else if (ls->get_right()->is_match(this->right)) {
451 return ls->get_left();
452 }
453 }
454
455// c1 + (c2 - a) -> c3 - a
456// c1 + (a - c2) -> c3 + a
457// a + (a - b) -> 2a - b
458// a + (b - a) -> b
459 auto rs = subtract_cast(this->right);
460 if (rs.get()) {
461 if (is_constant_combinable(this->left, rs->get_left())) {
462 return (this->left + rs->get_left()) - rs->get_right();
463 } else if (is_constant_combinable(this->left, rs->get_right())) {
464 return (this->left - rs->get_right()) + rs->get_left();
465 } else if (this->left->is_match(rs->get_left())) {
466 return static_cast<T> (2.0)*this->left - rs->get_right();
467 } else if (this->left->is_match(rs->get_right())) {
468 return rs->get_left();
469 }
470 }
471
472// Move cases like
473// (c1 + c2/x) + c3/y -> c1 + (c2/x + c3/y)
474// (c1 - c2/x) + c3/y -> c1 + (c3/y - c2/x)
475// in case of common denominators.
476 if (rd.get()) {
477 if (la.get() && divide_cast(la->get_right()).get()) {
478 return la->get_left() + (la->get_right() + this->right);
479 }
480
481 auto ls = subtract_cast(this->left);
482 if (ls.get() && divide_cast(ls->get_right()).get()) {
483 return ls->get_left() + (this->right - ls->get_right());
484 }
485 }
486
487 auto lfma = fma_cast(this->left);
488 auto rfma = fma_cast(this->right);
489 if (lfma.get()) {
490// fma(c,d,e) + a -> fma(c,d,e + a)
491 return fma(lfma->get_left(),
492 lfma->get_middle(),
493 lfma->get_right() + this->right);
494 } else if (rfma.get()) {
495// a + fma(c,d,e) -> fma(c,d,a + e)
496 return fma(rfma->get_left(),
497 rfma->get_middle(),
498 this->left + rfma->get_right());
499 }
500
501// fma(b,a,d) + fma(c,a,e) -> fma(a,b + c, d + e)
502// fma(a,b,d) + fma(c,a,e) -> fma(a,b + c, d + e)
503// fma(b,a,d) + fma(a,c,e) -> fma(a,b + c, d + e)
504// fma(a,b,d) + fma(a,c,e) -> fma(a,b + c, d + e)
505 if (lfma.get() && rfma.get()) {
506 if (lfma->get_middle()->is_match(rfma->get_middle())) {
507 return fma(lfma->get_middle(),
508 lfma->get_left() + rfma->get_left(),
509 lfma->get_right() + rfma->get_right());
510 } else if (lfma->get_left()->is_match(rfma->get_middle())) {
511 return fma(lfma->get_left(),
512 lfma->get_middle() + rfma->get_left(),
513 lfma->get_right() + rfma->get_right());
514 } else if (lfma->get_middle()->is_match(rfma->get_left())) {
515 return fma(lfma->get_middle(),
516 lfma->get_left() + rfma->get_middle(),
517 lfma->get_right() + rfma->get_right());
518 } else if (lfma->get_left()->is_match(rfma->get_left())) {
519 return fma(lfma->get_left(),
520 lfma->get_middle() + rfma->get_middle(),
521 lfma->get_right() + rfma->get_right());
522 }
523 }
524
525 auto pl = pow_cast(this->left);
526 auto pr = pow_cast(this->right);
527
528// (a*b)^c + (a*d)^c -> a^c*(b^c + d^c)
529// (b*a)^c + (a*d)^c -> a^c*(b^c + d^c)
530// (a*b)^c + (d*a)^c -> a^c*(b^c + d^c)
531// (b*a)^c + (d*a)^c -> a^c*(b^c + d^c)
532 if (pl.get() && pr.get() &&
533 pl->get_right()->is_match(pr->get_right())) {
534 auto plm = multiply_cast(pl->get_left());
535 auto prm = multiply_cast(pr->get_left());
536 if (plm.get() && prm.get()) {
537 if (plm->get_left()->is_match(prm->get_left())) {
538 return pow(plm->get_left(), pl->get_right())*
539 (pow(plm->get_right(), pl->get_right()) +
540 pow(prm->get_right(), pl->get_right()));
541 } else if (plm->get_left()->is_match(prm->get_right())) {
542 return pow(plm->get_left(), pl->get_right())*
543 (pow(plm->get_right(), pl->get_right()) +
544 pow(prm->get_left(), pl->get_right()));
545 } else if (plm->get_right()->is_match(prm->get_left())) {
546 return pow(plm->get_right(), pl->get_right())*
547 (pow(plm->get_left(), pl->get_right()) +
548 pow(prm->get_right(), pl->get_right()));
549 } else if (plm->get_right()->is_match(prm->get_right())) {
550 return pow(plm->get_right(), pl->get_right())*
551 (pow(plm->get_left(), pl->get_right()) +
552 pow(prm->get_left(), pl->get_right()));
553 }
554 }
555
556// cos(x)^2 + sin(x)^2 -> 1
557// sin(x)^2 + cos(x)^2 -> 1
558 auto plrc = constant_cast(pl->get_right());
559 if (plrc.get() && plrc->is(static_cast<T> (2.0))) {
560 auto pls = sin_cast(pl->get_left());
561 auto prc = cos_cast(pr->get_left());
562 auto plc = cos_cast(pl->get_left());
563 auto prs = sin_cast(pr->get_left());
564 if ((pls.get() && prc.get() && pls->get_arg()->is_match(prc->get_arg())) ||
565 (plc.get() && prs.get() && plc->get_arg()->is_match(prs->get_arg()))) {
566 return one<T, SAFE_MATH> ();
567 }
568 }
569 }
570
571// (a/y)^e + b/y^e -> (a^2 + b)/(y^e)
572// b/y^e + (a/y)^e -> (b + a^2)/(y^e)
573// (a/y)^e + (b/y)^e -> (a^2 + b^2)/(y^e)
574 if (pl.get() && rd.get()) {
575 auto rdp = pow_cast(rd->get_right());
576 if (rdp.get() && pl->get_right()->is_match(rdp->get_right())) {
577 auto plld = divide_cast(pl->get_left());
578 if (plld.get() &&
579 rdp->get_left()->is_match(plld->get_right())) {
580 return (pow(plld->get_left(), pl->get_right()) +
581 rd->get_left()) /
582 pow(rdp->get_left(), pl->get_right());
583 }
584 }
585 } else if (pr.get() && ld.get()) {
586 auto ldp = pow_cast(ld->get_right());
587 if (ldp.get() && pr->get_right()->is_match(ldp->get_right())) {
588 auto prld = divide_cast(pr->get_left());
589 if (prld.get() &&
590 ldp->get_left()->is_match(prld->get_right())) {
591 return (pow(prld->get_left(), pr->get_right()) +
592 ld->get_left()) /
593 pow(ldp->get_left(), pr->get_right());
594 }
595 }
596 } else if (pl.get() && pr.get()) {
597 if (pl->get_right()->is_match(pr->get_right())) {
598 auto pld = divide_cast(pl->get_left());
599 auto prd = divide_cast(pr->get_left());
600 if (pld.get() && prd.get() &&
601 pld->get_right()->is_match(prd->get_right())) {
602 return (pow(pld->get_left(), pl->get_right()) +
603 pow(prd->get_left(), pl->get_right())) /
604 pow(pld->get_right(), pl->get_right());
605 }
606 }
607 }
608
609 return this->shared_from_this();
610 }
611
612//------------------------------------------------------------------------------
619//------------------------------------------------------------------------------
620 virtual shared_leaf<T, SAFE_MATH>
621 df(shared_leaf<T, SAFE_MATH> x) {
622 if (this->is_match(x)) {
623 return one<T, SAFE_MATH> ();
624 }
625
626 const size_t hash = reinterpret_cast<size_t> (x.get());
627 if (this->df_cache.find(hash) == this->df_cache.end()) {
628 this->df_cache[hash] = this->left->df(x) + this->right->df(x);
629 }
630 return this->df_cache[hash];
631 }
632
633//------------------------------------------------------------------------------
641//------------------------------------------------------------------------------
642 virtual shared_leaf<T, SAFE_MATH>
643 compile(std::ostringstream &stream,
644 jit::register_map &registers,
645 const jit::register_map &thread_mem,
646 const jit::register_usage &usage) {
647 if (registers.find(this) == registers.end()) {
648 auto l = this->left->compile(stream, registers,
649 thread_mem, usage);
650 auto r = this->right->compile(stream, registers,
651 thread_mem, usage);
652
653 registers[this] = jit::to_string('r', this);
654 stream << " const ";
655 jit::add_type<T> (stream);
656 stream << " " << registers[this] << " = "
657 << registers[l.get()] << " + "
658 << registers[r.get()];
659 this->endline(stream, usage);
660 }
661
662 return this->shared_from_this();
663 }
664
665//------------------------------------------------------------------------------
670//------------------------------------------------------------------------------
671 virtual bool is_match(shared_leaf<T, SAFE_MATH> x) {
672 if (this == x.get()) {
673 return true;
674 }
675
676 auto x_cast = add_cast(x);
677 if (x_cast.get()) {
678// Addition is commutative.
679 return (this->left->is_match(x_cast->get_left()) &&
680 this->right->is_match(x_cast->get_right())) ||
681 (this->right->is_match(x_cast->get_left()) &&
682 this->left->is_match(x_cast->get_right()));
683 }
684
685 return false;
686 }
687
688//------------------------------------------------------------------------------
690//------------------------------------------------------------------------------
691 virtual void to_latex() const {
692 bool l_brackets = add_cast(this->left).get() ||
693 subtract_cast(this->left).get();
694 bool r_brackets = add_cast(this->right).get() ||
695 subtract_cast(this->right).get();
696 if (l_brackets) {
697 std::cout << "\\left(";
698 }
699 this->left->to_latex();
700 if (l_brackets) {
701 std::cout << "\\right)";
702 }
703 std::cout << "+";
704 if (r_brackets) {
705 std::cout << "\\left(";
706 }
707 this->right->to_latex();
708 if (r_brackets) {
709 std::cout << "\\right)";
710 }
711 }
712
713//------------------------------------------------------------------------------
717//------------------------------------------------------------------------------
718 virtual shared_leaf<T, SAFE_MATH> remove_pseudo() {
719 if (this->has_pseudo()) {
720 return this->left->remove_pseudo() +
721 this->right->remove_pseudo();
722 }
723 return this->shared_from_this();
724 }
725
726//------------------------------------------------------------------------------
732//------------------------------------------------------------------------------
733 virtual shared_leaf<T, SAFE_MATH> to_vizgraph(std::stringstream &stream,
734 jit::register_map &registers) {
735 if (registers.find(this) == registers.end()) {
736 const std::string name = jit::to_string('r', this);
737 registers[this] = name;
738 stream << " " << name
739 << " [label = \"+\", shape = oval, style = filled, fillcolor = blue, fontcolor = white];" << std::endl;
740
741 auto l = this->left->to_vizgraph(stream, registers);
742 stream << " " << name << " -- " << registers[l.get()] << ";" << std::endl;
743 auto r = this->right->to_vizgraph(stream, registers);
744 stream << " " << name << " -- " << registers[r.get()] << ";" << std::endl;
745 }
746
747 return this->shared_from_this();
748 }
749 };
750
751//------------------------------------------------------------------------------
762//------------------------------------------------------------------------------
763 template<jit::float_scalar T, bool SAFE_MATH=false>
764 shared_leaf<T, SAFE_MATH> add(shared_leaf<T, SAFE_MATH> l,
765 shared_leaf<T, SAFE_MATH> r) {
766 auto temp = std::make_shared<add_node<T, SAFE_MATH>> (l, r)->reduce();
767// Test for hash collisions.
768 for (size_t i = temp->get_hash();
769 i < std::numeric_limits<size_t>::max(); i++) {
770 if (leaf_node<T, SAFE_MATH>::caches.nodes.find(i) ==
771 leaf_node<T, SAFE_MATH>::caches.nodes.end()) {
772 leaf_node<T, SAFE_MATH>::caches.nodes[i] = temp;
773 return temp;
774 } else if (temp->is_match(leaf_node<T, SAFE_MATH>::caches.nodes[i])) {
775 return leaf_node<T, SAFE_MATH>::caches.nodes[i];
776 }
777 }
778#if defined(__clang__) || defined(__GNUC__)
779 __builtin_unreachable();
780#else
781 assert(false && "Should never reach.");
782#endif
783 }
784
785//------------------------------------------------------------------------------
796//------------------------------------------------------------------------------
797 template<jit::float_scalar T, bool SAFE_MATH=false>
802
803//------------------------------------------------------------------------------
815//------------------------------------------------------------------------------
816 template<jit::float_scalar T, jit::float_scalar L, bool SAFE_MATH=false>
821
822//------------------------------------------------------------------------------
834//------------------------------------------------------------------------------
835 template<jit::float_scalar T, jit::float_scalar R, bool SAFE_MATH=false>
840
842 template<jit::float_scalar T, bool SAFE_MATH=false>
843 using shared_add = std::shared_ptr<add_node<T, SAFE_MATH>>;
844
845//------------------------------------------------------------------------------
853//------------------------------------------------------------------------------
854 template<jit::float_scalar T, bool SAFE_MATH=false>
855 shared_add<T, SAFE_MATH> add_cast(shared_leaf<T, SAFE_MATH> x) {
856 return std::dynamic_pointer_cast<add_node<T, SAFE_MATH>> (x);
857 }
858
859//******************************************************************************
860// Subtract node.
861//******************************************************************************
862//------------------------------------------------------------------------------
869//------------------------------------------------------------------------------
870 template<jit::float_scalar T, bool SAFE_MATH=false>
871 class subtract_node final : public branch_node<T, SAFE_MATH> {
872 private:
873//------------------------------------------------------------------------------
879//------------------------------------------------------------------------------
880 static std::string to_string(leaf_node<T, SAFE_MATH> *l,
881 leaf_node<T, SAFE_MATH> *r) {
882 return jit::format_to_string(reinterpret_cast<size_t> (l)) + "-" +
883 jit::format_to_string(reinterpret_cast<size_t> (r));
884 }
885
886 public:
887//------------------------------------------------------------------------------
892//------------------------------------------------------------------------------
897
898//------------------------------------------------------------------------------
904//------------------------------------------------------------------------------
905 virtual backend::buffer<T> evaluate() {
906 backend::buffer<T> l_result = this->left->evaluate();
907 backend::buffer<T> r_result = this->right->evaluate();
908 return l_result - r_result;
909 }
910
911//------------------------------------------------------------------------------
915//------------------------------------------------------------------------------
916 virtual shared_leaf<T, SAFE_MATH> reduce() {
917// Identity reductions.
918 auto l = constant_cast(this->left);
919 if (this->left->is_match(this->right)) {
920 if (l.get() && l->is(0)) {
921 return this->left;
922 }
923
924 return zero<T, SAFE_MATH> ();
925 }
926
927// Constant reductions.
928 auto r = constant_cast(this->right);
929
930 if (l.get() && l->is(0)) {
931 return -this->right;
932 } else if (r.get() && r->is(0)) {
933 return this->left;
934 } else if (l.get() && r.get()) {
935 return constant<T, SAFE_MATH> (this->evaluate());
936 } else if (r.get() && r->evaluate().is_negative()) {
937 return this->left + -this->right;
938 }
939
940 auto pl1 = piecewise_1D_cast(this->left);
941 auto pr1 = piecewise_1D_cast(this->right);
942
943 if (pl1.get() && (r.get() || pl1->is_arg_match(this->right))) {
944 return piecewise_1D(this->evaluate(), pl1->get_arg());
945 } else if (pr1.get() && (l.get() || pr1->is_arg_match(this->left))) {
946 return piecewise_1D(this->evaluate(), pr1->get_arg());
947 }
948
949 auto pl2 = piecewise_2D_cast(this->left);
950 auto pr2 = piecewise_2D_cast(this->right);
951
952 if (pl2.get() && (r.get() || pl2->is_arg_match(this->right))) {
953 return piecewise_2D(this->evaluate(),
954 pl2->get_num_columns(),
955 pl2->get_left(),
956 pl2->get_right());
957 } else if (pr2.get() && (l.get() || pr2->is_arg_match(this->left))) {
958 return piecewise_2D(this->evaluate(),
959 pl2->get_num_columns(),
960 pr2->get_left(),
961 pr2->get_right());
962 }
963
964// Combine 2D and 1D piecewise constants if a row or column matches.
965 if (pr2.get() && pr2->is_row_match(this->left)) {
966 backend::buffer<T> result = pl1->evaluate();
967 result.subtract_row(pr2->evaluate());
968 return piecewise_2D(result,
969 pl2->get_num_columns(),
970 pr2->get_left(),
971 pr2->get_right());
972 } else if (pr2.get() && pr2->is_col_match(this->left)) {
973 backend::buffer<T> result = pl1->evaluate();
974 result.subtract_col(pr2->evaluate());
975 return piecewise_2D(result,
976 pl2->get_num_columns(),
977 pr2->get_left(),
978 pr2->get_right());
979 } else if (pl2.get() && pl2->is_row_match(this->right)) {
980 backend::buffer<T> result = pl2->evaluate();
981 result.subtract_row(pr1->evaluate());
982 return piecewise_2D(result,
983 pl2->get_num_columns(),
984 pl2->get_left(),
985 pl2->get_right());
986 } else if (pl2.get() && pl2->is_col_match(this->right)) {
987 backend::buffer<T> result = pl2->evaluate();
988 result.subtract_col(pr1->evaluate());
989 return piecewise_2D(result,
990 pl2->get_num_columns(),
991 pl2->get_left(),
992 pl2->get_right());
993 }
994// (c1 + a) - c2 -> c3 + a
995// c1 - (c2 + a) -> c3 - a
996 auto la = add_cast(this->left);
997 if (la.get()) {
998 if (is_constant_combinable(la->get_left(), this->right)) {
999 return (la->get_left() - this->right) + la->get_right();
1000 }
1001 }
1002 auto ra = add_cast(this->right);
1003 if (ra.get()) {
1004 if (is_constant_combinable(this->left, ra->get_left())) {
1005 return (this->left - ra->get_left()) - ra->get_right();
1006 }
1007 }
1008
1009// (c1 - a) - c2 -> c3 - a
1010// (a - c3) - c2 -> a + c3
1011 auto ls = subtract_cast(this->left);
1012 if (ls.get()) {
1013 if (is_constant_combinable(ls->get_left(), this->right)) {
1014 return (ls->get_left() - this->right) - ls->get_right();
1015 } else if (is_constant_combinable(ls->get_right(),
1016 this->right)) {
1017 return -(ls->get_right() + this->right) - ls->get_left();
1018 }
1019 }
1020// c1 - (c2 - a) -> c3 + a
1021// c1 - (a - c2) -> c3 - a
1022 auto rs = subtract_cast(this->right);
1023 if (rs.get()) {
1024 if (is_constant_combinable(this->left, rs->get_left())) {
1025 return (this->left - rs->get_left()) + rs->get_right();
1026 } else if (is_constant_combinable(this->left, rs->get_right())) {
1027 return (this->left + rs->get_right()) - rs->get_left();
1028 }
1029 }
1030
1031// Common factor reduction. If the left and right are both multiply nodes check
1032// for a common factor. So you can change a*b - a*c -> a*(b - c).
1033 auto lm = multiply_cast(this->left);
1034 auto rm = multiply_cast(this->right);
1035
1036// c1*(c2 + a) - c3 -> fma(c1,a,c4)
1037 if (lm.get()) {
1038 auto lmra = add_cast(lm->get_right());
1039 if (lmra.get()) {
1040 if (is_constant_combinable(lm->get_left(),
1041 lmra->get_left()) &&
1042 is_constant_combinable(lm->get_left(),
1043 this->right)) {
1044 return fma(lm->get_left(),
1045 lmra->get_right(),
1046 lm->get_left()*lmra->get_left() - this->right);
1047 }
1048 }
1049// c1*(c2 - a) - c3 -> c4 - c1*a
1050 auto lmrs = subtract_cast(lm->get_right());
1051 if (lmrs.get()) {
1052 if (is_constant_combinable(lm->get_left(),
1053 lmrs->get_left()) &&
1054 is_constant_combinable(lm->get_left(),
1055 this->right)) {
1056 return lm->get_left()*lmrs->get_left() - this->right -
1057 lm->get_left()*lmrs->get_right();
1058 }
1059 }
1060 }
1061
1062// (a + b) - a -> b
1063// (a + b) - b -> a
1064// a - (a + b) -> -b
1065// b - (a + b) -> -a
1066 if (la.get()) {
1067 if (la->get_left()->is_match(this->right)) {
1068 return la->get_right();
1069 } else if (la->get_right()->is_match(this->right)) {
1070 return la->get_left();
1071 }
1072 } else if (ra.get()) {
1073 if (ra->get_left()->is_match(this->left)) {
1074 return none<T, SAFE_MATH> ()*ra->get_right();
1075 } else if (ra->get_right()->is_match(this->left)) {
1076 return none<T, SAFE_MATH> ()*ra->get_left();
1077 }
1078 }
1079
1080// Assume constants are on the left.
1081// v1 - -c*v2 -> v1 + c*v2
1082 if (rm.get() &&
1083 rm->get_left()->is_constant() &&
1084 rm->get_left()->evaluate().is_negative()) {
1085 return this->left + (-this->right);
1086 }
1087
1088 if (lm.get()) {
1089// Assume constants are on the left.
1090// -a - b -> -(a + b)
1091 auto lmc = constant_cast(lm->get_left());
1092 if (lmc.get() && lmc->is(-1)) {
1093 return lm->get_left()*(lm->get_right() + this->right);
1094 }
1095
1096// a*v - v = (a - 1)*v
1097// v*a - v = (a - 1)*v
1098 if (this->right->is_match(lm->get_right())) {
1099 return (lm->get_left() - 1.0)*this->right;
1100 } else if (this->right->is_match(lm->get_left())) {
1101 return (lm->get_right() - 1.0)*this->right;
1102 }
1103 }
1104// v - a*v = (1 - a)*v
1105// v - v*a = (1 - a)*v
1106 if (rm.get()) {
1107 if (this->left->is_match(rm->get_right())) {
1108 return (1.0 - rm->get_left())*this->left;
1109 } else if (this->left->is_match(rm->get_left())) {
1110 return (1.0 - rm->get_right())*this->left;
1111 }
1112 }
1113
1114 if (lm.get() && rm.get()) {
1115 if (lm->get_left()->is_match(rm->get_left())) {
1116// a*b - a*c -> a*(b - c)
1117 return lm->get_left()*(lm->get_right() - rm->get_right());
1118 } else if (lm->get_left()->is_match(rm->get_right())) {
1119// a*b - c*a -> a*(b - c)
1120 return lm->get_left()*(lm->get_right() - rm->get_left());
1121 } else if (lm->get_right()->is_match(rm->get_left())) {
1122// b*a - a*c -> a*(b - c)
1123 return lm->get_right()*(lm->get_left() - rm->get_right());
1124 } else if (lm->get_right()->is_match(rm->get_right())) {
1125// b*a - c*a -> a*(b - c)
1126 return lm->get_right()*(lm->get_left() - rm->get_left());
1127 }
1128
1129// Change cases like c1*a - c2*b -> c1*(a - c2/c1*b)
1130// Note need to make sure c1 doesn't contain any zeros.
1131 if (lm->get_left()->is_constant() &&
1132 rm->get_left()->is_constant() &&
1133 !lm->get_left()->has_constant_zero()) {
1134 return lm->get_left()*(lm->get_right() -
1135 (rm->get_left()/lm->get_left())*rm->get_right());
1136 }
1137
1138// Handle case
1139 auto rmrm = multiply_cast(rm->get_right());
1140 if (rmrm.get()) {
1141// a*b - c*(d*b) -> (a - c*d)*b
1142 if (lm->get_right()->is_match(rmrm->get_right())) {
1143 return (lm->get_left() - rm->get_left()*rmrm->get_left())*lm->get_right();
1144 }
1145// a*b - c*(b*d) -> (a - c*d)*b
1146 if (lm->get_right()->is_match(rmrm->get_left())) {
1147 return (lm->get_left() - rm->get_left()*rmrm->get_right())*lm->get_right();
1148 }
1149// b*a - c*(d*b) -> (a - c*d)*b
1150 if (lm->get_left()->is_match(rmrm->get_right())) {
1151 return (lm->get_right() - rm->get_left()*rmrm->get_left())*lm->get_left();
1152 }
1153// b*a - c*(b*d) -> (a - c*d)*b
1154 if (lm->get_left()->is_match(rmrm->get_left())) {
1155 return (lm->get_right() - rm->get_left()*rmrm->get_right())*lm->get_left();
1156 }
1157 }
1158 auto lmrm = multiply_cast(lm->get_right());
1159 if (lmrm.get()) {
1160// c*(d*b) - a*b -> (c*d - a)*b
1161 if (rm->get_right()->is_match(lmrm->get_right())) {
1162 return (lm->get_left()*lmrm->get_left() - rm->get_left())*rm->get_right();
1163 }
1164// c*(b*d) - a*b -> (c*d - a)*b
1165 if (rm->get_right()->is_match(lmrm->get_left())) {
1166 return (lm->get_left()*lmrm->get_right() - rm->get_left())*rm->get_right();
1167 }
1168// c*(d*b) - b*a -> (c*d - a)*b
1169 if (rm->get_left()->is_match(lmrm->get_right())) {
1170 return (lm->get_left()*lmrm->get_left() - rm->get_right())*rm->get_left();
1171 }
1172// c*(b*d) - b*a -> (c*d - a)*b
1173 if (rm->get_left()->is_match(lmrm->get_left())) {
1174 return (lm->get_left()*lmrm->get_right() - rm->get_right())*rm->get_left();
1175 }
1176 }
1177
1178// a/b*c - d/b*e -> (a*b - d*e)/b
1179// a/b*c - d*e/b -> (a*b - d*e)/b
1180// a*c/b - d/b*e -> (a*b - d*e)/b
1181// a*c/b - d*e/b -> (a*b - d*e)/b
1182 auto lmld = divide_cast(lm->get_left());
1183 auto rmld = divide_cast(rm->get_left());
1184 auto lmrd = divide_cast(lm->get_right());
1185 auto rmrd = divide_cast(rm->get_right());
1186 if (lmld.get() && rmld.get() &&
1187 lmld->get_right()->is_match(rmld->get_right())) {
1188 return (lmld->get_left()*lm->get_right() -
1189 rmld->get_left()*rm->get_right())/lmld->get_right();
1190 } else if (lmld.get() && rmrd.get() &&
1191 lmld->get_right()->is_match(rmrd->get_right())) {
1192 return (lmld->get_left()*lm->get_right() -
1193 rmrd->get_left()*rm->get_left())/lmld->get_right();
1194 } else if (lmrd.get() && rmld.get() &&
1195 lmrd->get_right()->is_match(rmld->get_right())) {
1196 return (lmrd->get_left()*lm->get_left() -
1197 rmld->get_left()*rm->get_right())/lmrd->get_right();
1198 } else if (lmrd.get() && rmrd.get() &&
1199 lmrd->get_right()->is_match(rmrd->get_right())) {
1200 return (lmrd->get_left()*lm->get_left() -
1201 rmrd->get_left()*rm->get_left())/lmrd->get_right();
1202 }
1203 }
1204
1205// Chained subtraction reductions.
1206 if (ls.get()) {
1207 auto lrm = multiply_cast(ls->get_right());
1208 if (lrm.get() && rm.get()) {
1209 if (lrm->get_left()->is_match(rm->get_left())) {
1210// (a - c*b) - c*d -> a - (b + d)*c
1211 return ls->get_left() -
1212 (lrm->get_right() +
1213 rm->get_right())*rm->get_left();
1214 } else if (lrm->get_left()->is_match(rm->get_right())) {
1215// (a - c*b) - d*c -> a - (b + d)*c
1216 return ls->get_left() -
1217 (lrm->get_right() +
1218 rm->get_left())*rm->get_right();
1219 } else if (lrm->get_right()->is_match(rm->get_left())) {
1220// (a - c*b) - c*d -> a - (b + d)*c
1221 return ls->get_left() -
1222 (lrm->get_left() +
1223 rm->get_right())*rm->get_left();
1224 } else if (lrm->get_right()->is_match(rm->get_right())) {
1225// (a - c*b) - d*c -> a - (b + d)*c
1226 return ls->get_left() -
1227 (lrm->get_left() +
1228 rm->get_left())*rm->get_right();
1229 }
1230 }
1231 }
1232
1233// Common denominator reduction. If the left and right are both divide nodes
1234// for a common denominator. So you can change a/b - c/b -> (a - c)/d.
1235 auto ld = divide_cast(this->left);
1236 auto rd = divide_cast(this->right);
1237
1238// c is a constant.
1239// a - -c/b -> a + c/b
1240// a - (-c*d)/b -> a + (c*d)/b
1241// -c/a - b -> -(b + c/a)
1242// (-c*d)/a - b -> -(b + (c*d)/a)
1243 if (rd.get()) {
1244 auto rdlm = multiply_cast(rd->get_left());
1245 if ((rd->get_left()->is_constant() &&
1246 rd->get_left()->evaluate().is_negative()) ||
1247 (rdlm.get() &&
1248 (rdlm->get_left()->is_constant() &&
1249 rdlm->get_left()->evaluate().is_negative()))) {
1250 return this->left + -this->right;
1251 }
1252 } else if (ld.get()) {
1253 auto ldlm = multiply_cast(ld->get_left());
1254 if ((ld->get_left()->is_constant() &&
1255 ld->get_left()->evaluate().is_negative()) ||
1256 (ldlm.get() &&
1257 (ldlm->get_left()->is_constant() &&
1258 ldlm->get_left()->evaluate().is_negative()))) {
1259 return -(-this->left + this->right);
1260 }
1261 }
1262
1263 if (ld.get() && rd.get()) {
1264 if (ld->get_right()->is_match(rd->get_right())) {
1265 return (ld->get_left() - rd->get_left())/ld->get_right();
1266 }
1267
1268 auto ldlm = multiply_cast(ld->get_left());
1269 auto rdlm = multiply_cast(rd->get_left());
1270// a/b - c*a/d -> (1/b - c/d)*a
1271// a/b - a*c/d -> (1/b - c/d)*a
1272// c*a/b - a/d -> (c/b - 1/d)*a
1273// a*c/b - a/d -> (c/b - 1/d)*a
1274 if (rdlm.get()) {
1275 if (ld->get_left()->is_match(rdlm->get_left())) {
1276 return (1.0/ld->get_right() -
1277 rdlm->get_right()/rd->get_right())*rdlm->get_left();
1278 } else if (ld->get_left()->is_match(rdlm->get_right())) {
1279 return (1.0/ld->get_right() -
1280 rdlm->get_left()/rd->get_right())*rdlm->get_right();
1281 }
1282 } else if (ldlm.get()) {
1283 if (rd->get_left()->is_match(ldlm->get_left())) {
1284 return (ldlm->get_right()/ld->get_right() -
1285 1.0/rd->get_right())*ldlm->get_left();
1286 } else if (rd->get_left()->is_match(ldlm->get_right())) {
1287 return (ldlm->get_left()/ld->get_right() -
1288 1.0/rd->get_right())*ldlm->get_right();
1289 }
1290 }
1291
1292// c1*a/b - c2*e/d = c3*(a/b - c4*e/d)
1293// a*b/c - d*b/e -> (a/c - d/e)*b
1294// Make sure we prevent combining constants when we just need to factor out a
1295// common term.
1296// c1*a/b - c2*a/d -> (c1/b - c2/d)*a
1297 if (ldlm.get() && rdlm.get()) {
1298 if (is_constant_combinable(ldlm->get_left(),
1299 rdlm->get_left()) &&
1300 !ldlm->get_right()->is_match(rdlm->get_right())) {
1301 return (ldlm->get_right()/ld->get_right() -
1302 rdlm->get_left()/ldlm->get_left() *
1303 rdlm->get_right()/rd->get_right())*ldlm->get_left();
1304 }
1305
1306 if (ldlm->get_right()->is_match(rdlm->get_right())) {
1307 return (ldlm->get_left()/ld->get_right() -
1308 rdlm->get_left()/rd->get_right())*ldlm->get_right();
1309 } else if (ldlm->get_right()->is_match(rdlm->get_left())) {
1310 return (ldlm->get_left()/ld->get_right() -
1311 rdlm->get_right()/rd->get_right())*ldlm->get_right();
1312 } else if (ldlm->get_left()->is_match(rdlm->get_right())) {
1313 return (ldlm->get_right()/ld->get_right() -
1314 rdlm->get_left()/rd->get_right())*ldlm->get_left();
1315 } else if (ldlm->get_left()->is_match(rdlm->get_left())) {
1316 return (ldlm->get_right()/ld->get_right() -
1317 rdlm->get_right()/rd->get_right())*ldlm->get_left();
1318 }
1319 }
1320
1321// (a/(c*b) - d/(e*c)) -> (a/b - d/e)/c
1322// (a/(b*c) - d/(e*c)) -> (a/b - d/e)/c
1323// (a/(c*b) - d/(c*e)) -> (a/b - d/e)/c
1324// (a/(b*c) - d/(c*e)) -> (a/b - d/e)/c
1325 auto ldrm = multiply_cast(ld->get_right());
1326 auto rdrm = multiply_cast(rd->get_right());
1327 if (ldrm.get() && rdrm.get()) {
1328 if (ldrm->get_right()->is_match(rdrm->get_right())) {
1329 return (ld->get_left()/ldrm->get_left() -
1330 rd->get_left()/rdrm->get_left())/ldrm->get_right();
1331 } else if (ldrm->get_right()->is_match(rdrm->get_left())) {
1332 return (ld->get_left()/ldrm->get_left() -
1333 rd->get_left()/rdrm->get_right())/ldrm->get_right();
1334 } else if (ldrm->get_left()->is_match(rdrm->get_right())) {
1335 return (ld->get_left()/ldrm->get_right() -
1336 rd->get_left()/rdrm->get_left())/ldrm->get_left();
1337 } else if (ldrm->get_left()->is_match(rdrm->get_left())) {
1338 return (ld->get_left()/ldrm->get_right() -
1339 rd->get_left()/rdrm->get_right())/ldrm->get_left();
1340 }
1341 }
1342
1343// a/b - c/(b*d) -> (a*d - c)/(b*d)
1344// a/b - c/(d*b) -> (a*d - c)/(b*d)
1345// a/(b*d) - c/b -> (a - c*d)/(b*d)
1346// a/(d*b) - c/b -> (a - c*d)/(b*d)
1347 if (rdrm.get()) {
1348 if (ld->get_right()->is_match(rdrm->get_left())) {
1349 return (ld->get_left()*rdrm->get_right() - rd->get_left()) /
1350 rd->get_right();
1351 } else if (ld->get_right()->is_match(rdrm->get_right())) {
1352 return (ld->get_left()*rdrm->get_left() - rd->get_left()) /
1353 rd->get_right();
1354 }
1355 } else if (ldrm.get()) {
1356 if (rd->get_right()->is_match(ldrm->get_left())) {
1357 return (ld->get_left() - rd->get_left()*ldrm->get_right()) /
1358 ld->get_right();
1359 } else if (rd->get_right()->is_match(ldrm->get_right())) {
1360 return (ld->get_left() - rd->get_left()*ldrm->get_left()) /
1361 ld->get_right();
1362 }
1363 }
1364 }
1365
1366// Move cases like
1367// (c1 + c2/x) - c3/y -> c1 + (c2/x - c3/y)
1368// (c1 - c2/x) - c3/y -> c1 - (c2/x + c3/y)
1369// in case of common denominators.
1370 if (rd.get()) {
1371 auto la = add_cast(this->left);
1372 if (la.get() && divide_cast(la->get_right()).get()) {
1373 return la->get_left() + (la->get_right() - this->right);
1374 } else if (ls.get() && divide_cast(ls->get_right()).get()) {
1375 return ls->get_left() - (this->right + ls->get_right());
1376 }
1377 }
1378
1379// Handle cases like:
1380// (a/y)^e - b/y^e -> (a^2 - b)/(y^e)
1381// b/y^e - (a/y)^e -> (b - a^2)/(y^e)
1382// (a/y)^e - (b/y)^e -> (a^2 - b^2)/(y^e)
1383 auto pl = pow_cast(this->left);
1384 auto pr = pow_cast(this->right);
1385 if (pl.get() && rd.get()) {
1386 auto rdp = pow_cast(rd->get_right());
1387 if (rdp.get() && pl->get_right()->is_match(rdp->get_right())) {
1388 auto plld = divide_cast(pl->get_left());
1389 if (plld.get() &&
1390 rdp->get_left()->is_match(plld->get_right())) {
1391 return (pow(plld->get_left(), pl->get_right()) -
1392 rd->get_left()) /
1393 pow(rdp->get_left(), pl->get_right());
1394 }
1395 }
1396 } else if (pr.get() && ld.get()) {
1397 auto ldp = pow_cast(ld->get_right());
1398 if (ldp.get() && pr->get_right()->is_match(ldp->get_right())) {
1399 auto prld = divide_cast(pr->get_left());
1400 if (prld.get() &&
1401 ldp->get_left()->is_match(prld->get_right())) {
1402 return (pow(prld->get_left(), pr->get_right()) -
1403 ld->get_left()) /
1404 pow(ldp->get_left(), pr->get_right());
1405 }
1406 }
1407 } else if (pl.get() && pr.get()) {
1408 if (pl->get_right()->is_match(pr->get_right())) {
1409 auto pld = divide_cast(pl->get_left());
1410 auto prd = divide_cast(pr->get_left());
1411 if (pld.get() && prd.get() &&
1412 pld->get_right()->is_match(prd->get_right())) {
1413 return (pow(pld->get_left(), pl->get_right()) -
1414 pow(prd->get_left(), pl->get_right())) /
1415 pow(pld->get_right(), pl->get_right());
1416 }
1417 }
1418 }
1419
1420 auto lfma = fma_cast(this->left);
1421 auto rfma = fma_cast(this->right);
1422
1423 if (lfma.get() && rfma.get()) {
1424 if (lfma->get_middle()->is_match(rfma->get_middle())) {
1425 return fma(lfma->get_left() - rfma->get_left(),
1426 lfma->get_middle(),
1427 lfma->get_right() - rfma->get_right());
1428 }
1429 }
1430
1431// fma(c,d,e) - a -> fma(c,d,e - a)
1432 if (lfma.get() && !this->right->is_all_variables()) {
1433 return fma(lfma->get_left(),
1434 lfma->get_middle(),
1435 lfma->get_right() - this->right);
1436 }
1437
1438// Reduce cases chained subtract multiply divide.
1439 if (ls.get()) {
1440// (a - b*c) - d*e -> a - (b*c + d*e)
1441// (a - b/c) - d/e -> a - (b/c + d/e)
1442 auto lsrd = divide_cast(ls->get_right());
1443 if ((multiply_cast(ls->get_right()).get() && (rm.get() || rd.get())) ||
1444 (divide_cast(ls->get_right()).get() && (rm.get() || rd.get()))) {
1445 return ls->get_left() - (ls->get_right() + this->right);
1446 }
1447 }
1448
1449 return this->shared_from_this();
1450 }
1451
1452//------------------------------------------------------------------------------
1459//------------------------------------------------------------------------------
1460 virtual shared_leaf<T, SAFE_MATH>
1461 df(shared_leaf<T, SAFE_MATH> x) {
1462 if (this->is_match(x)) {
1463 return one<T, SAFE_MATH> ();
1464 }
1465
1466 const size_t hash = reinterpret_cast<size_t> (x.get());
1467 if (this->df_cache.find(hash) == this->df_cache.end()) {
1468 this->df_cache[hash] = this->left->df(x) - this->right->df(x);
1469 }
1470 return this->df_cache[hash];
1471 }
1472
1473//------------------------------------------------------------------------------
1481//------------------------------------------------------------------------------
1482 virtual shared_leaf<T, SAFE_MATH>
1483 compile(std::ostringstream &stream,
1484 jit::register_map &registers,
1485 const jit::register_map &thread_mem,
1486 const jit::register_usage &usage) {
1487 if (registers.find(this) == registers.end()) {
1488 auto l = this->left->compile(stream, registers,
1489 thread_mem, usage);
1490 auto r = this->right->compile(stream, registers,
1491 thread_mem, usage);
1492
1493 registers[this] = jit::to_string('r', this);
1494 stream << " const ";
1495 jit::add_type<T> (stream);
1496 stream << " " << registers[this] << " = "
1497 << registers[l.get()] << " - "
1498 << registers[r.get()];
1499 this->endline(stream, usage);
1500 }
1501
1502 return this->shared_from_this();
1503 }
1504
1505//------------------------------------------------------------------------------
1510//------------------------------------------------------------------------------
1511 virtual bool is_match(shared_leaf<T, SAFE_MATH> x) {
1512 if (this == x.get()) {
1513 return true;
1514 }
1515
1516 auto x_cast = subtract_cast(x);
1517 if (x_cast.get()) {
1518 return this->left->is_match(x_cast->get_left()) &&
1519 this->right->is_match(x_cast->get_right());
1520 }
1521
1522 return false;
1523 }
1524
1525//------------------------------------------------------------------------------
1527//------------------------------------------------------------------------------
1528 virtual void to_latex() const {
1529 bool l_brackets = add_cast(this->left).get() ||
1530 subtract_cast(this->left).get();
1531 bool r_brackets = add_cast(this->right).get() ||
1532 subtract_cast(this->right).get();
1533 if (l_brackets) {
1534 std::cout << "\\left(";
1535 }
1536 this->left->to_latex();
1537 if (l_brackets) {
1538 std::cout << "\\right)";
1539 }
1540 std::cout << "-";
1541 if (r_brackets) {
1542 std::cout << "\\left(";
1543 }
1544 this->right->to_latex();
1545 if (r_brackets) {
1546 std::cout << "\\right)";
1547 }
1548 }
1549
1550//------------------------------------------------------------------------------
1554//------------------------------------------------------------------------------
1555 virtual shared_leaf<T, SAFE_MATH> remove_pseudo() {
1556 if (this->has_pseudo()) {
1557 return this->left->remove_pseudo() -
1558 this->right->remove_pseudo();
1559 }
1560 return this->shared_from_this();
1561 }
1562
1563//------------------------------------------------------------------------------
1569//------------------------------------------------------------------------------
1570 virtual shared_leaf<T, SAFE_MATH> to_vizgraph(std::stringstream &stream,
1571 jit::register_map &registers) {
1572 if (registers.find(this) == registers.end()) {
1573 const std::string name = jit::to_string('r', this);
1574 registers[this] = name;
1575 stream << " " << name
1576 << " [label = \"-\", shape = oval, style = filled, fillcolor = blue, fontcolor = white];" << std::endl;
1577
1578 auto l = this->left->to_vizgraph(stream, registers);
1579 stream << " " << name << " -- " << registers[l.get()] << ";" << std::endl;
1580 auto r = this->right->to_vizgraph(stream, registers);
1581 stream << " " << name << " -- " << registers[r.get()] << ";" << std::endl;
1582 }
1583
1584 return this->shared_from_this();
1585 }
1586 };
1587
1588//------------------------------------------------------------------------------
1597//------------------------------------------------------------------------------
1598 template<jit::float_scalar T, bool SAFE_MATH=false>
1599 shared_leaf<T, SAFE_MATH> subtract(shared_leaf<T, SAFE_MATH> l,
1600 shared_leaf<T, SAFE_MATH> r) {
1601 auto temp = std::make_shared<subtract_node<T, SAFE_MATH>> (l, r)->reduce();
1602// Test for hash collisions.
1603 for (size_t i = temp->get_hash();
1604 i < std::numeric_limits<size_t>::max(); i++) {
1605 if (leaf_node<T, SAFE_MATH>::caches.nodes.find(i) ==
1606 leaf_node<T, SAFE_MATH>::caches.nodes.end()) {
1607 leaf_node<T, SAFE_MATH>::caches.nodes[i] = temp;
1608 return temp;
1609 } else if (temp->is_match(leaf_node<T, SAFE_MATH>::caches.nodes[i])) {
1610 return leaf_node<T, SAFE_MATH>::caches.nodes[i];
1611 }
1612 }
1613#if defined(__clang__) || defined(__GNUC__)
1614 __builtin_unreachable();
1615#else
1616 assert(false && "Should never reach.");
1617#endif
1618 }
1619
1620//------------------------------------------------------------------------------
1632//------------------------------------------------------------------------------
1633 template<jit::float_scalar T, bool SAFE_MATH=false>
1638
1639//------------------------------------------------------------------------------
1652//------------------------------------------------------------------------------
1653 template<jit::float_scalar T, jit::float_scalar L, bool SAFE_MATH=false>
1658
1659//------------------------------------------------------------------------------
1672//------------------------------------------------------------------------------
1673 template<jit::float_scalar T, jit::float_scalar R, bool SAFE_MATH=false>
1678
1679//------------------------------------------------------------------------------
1690//------------------------------------------------------------------------------
1691 template<jit::float_scalar T, bool SAFE_MATH=false>
1695
1697 template<jit::float_scalar T, bool SAFE_MATH=false>
1698 using shared_subtract = std::shared_ptr<subtract_node<T, SAFE_MATH>>;
1699
1700//------------------------------------------------------------------------------
1708//------------------------------------------------------------------------------
1709 template<jit::float_scalar T, bool SAFE_MATH=false>
1710 shared_subtract<T, SAFE_MATH> subtract_cast(shared_leaf<T, SAFE_MATH> x) {
1711 return std::dynamic_pointer_cast<subtract_node<T, SAFE_MATH>> (x);
1712 }
1713
1714//******************************************************************************
1715// Multiply node.
1716//******************************************************************************
1717//------------------------------------------------------------------------------
1722//------------------------------------------------------------------------------
1723 template<jit::float_scalar T, bool SAFE_MATH=false>
1724 class multiply_node final : public branch_node<T, SAFE_MATH> {
1725 private:
1726//------------------------------------------------------------------------------
1733//------------------------------------------------------------------------------
1734 shared_leaf<T, SAFE_MATH>
1735 reduce_nested_fma_times_constant(shared_leaf<T, SAFE_MATH> trial) {
1736 auto temp = fma_cast(trial);
1737 if (temp.get()) {
1738 if (is_constant_combinable(this->left, temp->get_left()) &&
1739 is_constant_combinable(this->left, temp->get_right())) {
1740 return fma(this->left*temp->get_left(),
1741 temp->get_middle(),
1742 this->left*temp->get_right());
1743 } else {
1744 auto temp2 = reduce_nested_fma_times_constant(temp->get_left());
1745 if (temp2.get()) {
1746 return fma(temp2,
1747 temp->get_middle(),
1748 this->left*temp->get_right());
1749 }
1750 }
1751 }
1752 return null_leaf<T, SAFE_MATH> ();
1753 }
1754
1755//------------------------------------------------------------------------------
1763//------------------------------------------------------------------------------
1764 shared_leaf<T, SAFE_MATH>
1765 expand_nested_fma_times_add(shared_leaf<T, SAFE_MATH> trial,
1766 shared_add<T, SAFE_MATH> add) {
1767 auto temp = fma_cast(trial);
1768 if (temp.get()) {
1769 if (add->get_right()->is_match(temp->get_middle()) &&
1770 is_constant_combinable(add->get_left(), temp->get_right())) {
1771 auto temp2 = expand_nested_fma_times_add2(temp->get_left(),
1772 temp, add);
1773 if (temp2.get()) {
1774 return fma(temp2,
1775 add->get_right(),
1776 temp->get_right()*add->get_left());
1777 } else if (is_constant_combinable(add->get_left(), temp->get_left())) {
1778 return fma(fma(temp->get_left(),
1779 add->get_right(),
1780 add->get_left()*temp->get_left() + temp->get_right()),
1781 add->get_right(),
1782 temp->get_right()*add->get_left());
1783 }
1784 }
1785 }
1786 return null_leaf<T, SAFE_MATH> ();
1787 }
1788
1789//------------------------------------------------------------------------------
1798//------------------------------------------------------------------------------
1799 shared_leaf<T, SAFE_MATH>
1800 expand_nested_fma_times_add2(shared_leaf<T, SAFE_MATH> trial,
1801 shared_leaf<T, SAFE_MATH> last,
1802 shared_add<T, SAFE_MATH> add) {
1803 auto temp = fma_cast(trial);
1804 auto temp2 = fma_cast(last);
1805 assert(temp2.get() && "Assumed a fma node.");
1806 if (temp.get()) {
1807 if (add->get_right()->is_match(temp->get_middle()) &&
1808 is_constant_combinable(add->get_left(), temp->get_left()) &&
1809 is_constant_combinable(add->get_left(), temp->get_right())) {
1810
1811 return fma(fma(temp->get_left(),
1812 add->get_right(),
1813 add->get_left()*temp->get_left() +
1814 temp->get_right()),
1815 add->get_right(),
1816 add->get_left()*temp->get_right() +
1817 temp2->get_right());
1818 } else {
1819 auto temp3 = expand_nested_fma_times_add2(temp->get_left(),
1820 temp, add);
1821 if (temp3.get()) {
1822 return fma(temp3,
1823 add->get_right(),
1824 add->get_left()*temp->get_right() +
1825 temp2->get_right());
1826 }
1827 }
1828 }
1829 return null_leaf<T, SAFE_MATH> ();
1830 }
1831
1832//------------------------------------------------------------------------------
1838//------------------------------------------------------------------------------
1839 static std::string to_string(leaf_node<T, SAFE_MATH> *l,
1840 leaf_node<T, SAFE_MATH> *r) {
1841 return jit::format_to_string(reinterpret_cast<size_t> (l)) + "*" +
1842 jit::format_to_string(reinterpret_cast<size_t> (r));
1843 }
1844
1845 public:
1846//------------------------------------------------------------------------------
1851//------------------------------------------------------------------------------
1855
1856//------------------------------------------------------------------------------
1862//------------------------------------------------------------------------------
1863 virtual backend::buffer<T> evaluate() {
1864 backend::buffer<T> l_result = this->left->evaluate();
1865
1866// If the left are right are same don't evaluate the right.
1867// NOTE: Do not use is_match here. Remove once power is implemented.
1868 if (this->left.get() == this->right.get()) {
1869 return l_result*l_result;
1870 }
1871
1872// If all the elements on the left are zero, return the left side without
1873// reevaluating the right side. Stop this loop early once the first non zero
1874// element is encountered.
1875 if (l_result.is_zero()) {
1876 return l_result;
1877 }
1878
1879 backend::buffer<T> r_result = this->right->evaluate();
1880 return l_result*r_result;
1881 }
1882
1883//------------------------------------------------------------------------------
1887//------------------------------------------------------------------------------
1888 virtual shared_leaf<T, SAFE_MATH> reduce() {
1889 auto l = constant_cast(this->left);
1890 auto r = constant_cast(this->right);
1891
1892 if (l.get() && l->is(1)) {
1893 return this->right;
1894 } else if (l.get() && l->is(0)) {
1895 return this->left;
1896 } else if (r.get() && r->is(1)) {
1897 return this->left;
1898 } else if (r.get() && r->is(0)) {
1899 return this->right;
1900 } else if (l.get() && r.get()) {
1901 return constant<T, SAFE_MATH> (this->evaluate());
1902 }
1903
1904 auto pl1 = piecewise_1D_cast(this->left);
1905 auto pr1 = piecewise_1D_cast(this->right);
1906
1907 if (pl1.get() && (r.get() || pl1->is_arg_match(this->right))) {
1908 return piecewise_1D(this->evaluate(), pl1->get_arg());
1909 } else if (pr1.get() && (l.get() || pr1->is_arg_match(this->left))) {
1910 return piecewise_1D(this->evaluate(), pr1->get_arg());
1911 }
1912
1913 auto pl2 = piecewise_2D_cast(this->left);
1914 auto pr2 = piecewise_2D_cast(this->right);
1915
1916 if (pl2.get() && (r.get() || pl2->is_arg_match(this->right))) {
1917 return piecewise_2D(this->evaluate(),
1918 pl2->get_num_columns(),
1919 pl2->get_left(),
1920 pl2->get_right());
1921 } else if (pr2.get() && (l.get() || pr2->is_arg_match(this->left))) {
1922 return piecewise_2D(this->evaluate(),
1923 pr2->get_num_columns(),
1924 pr2->get_left(),
1925 pr2->get_right());
1926 }
1927
1928// Combine 2D and 1D piecewise constants if a row or column matches.
1929 if (pr2.get() && pr2->is_row_match(this->left)) {
1930 backend::buffer<T> result = pl1->evaluate();
1931 result.multiply_row(pr2->evaluate());
1932 return piecewise_2D(result,
1933 pr2->get_num_columns(),
1934 pr2->get_left(),
1935 pr2->get_right());
1936 } else if (pr2.get() && pr2->is_col_match(this->left)) {
1937 backend::buffer<T> result = pl1->evaluate();
1938 result.multiply_col(pr2->evaluate());
1939 return piecewise_2D(result,
1940 pr2->get_num_columns(),
1941 pr2->get_left(),
1942 pr2->get_right());
1943 } else if (pl2.get() && pl2->is_row_match(this->right)) {
1944 backend::buffer<T> result = pl2->evaluate();
1945 result.multiply_row(pr1->evaluate());
1946 return piecewise_2D(result,
1947 pl2->get_num_columns(),
1948 pl2->get_left(),
1949 pl2->get_right());
1950 } else if (pl2.get() && pl2->is_col_match(this->right)) {
1951 backend::buffer<T> result = pl2->evaluate();
1952 result.multiply_col(pr1->evaluate());
1953 return piecewise_2D(result,
1954 pl2->get_num_columns(),
1955 pl2->get_left(),
1956 pl2->get_right());
1957 }
1958
1959// Move constants to the left.
1960 if (is_constant_promotable(this->right, this->left)) {
1961 return this->right*this->left;
1962 }
1963
1964// Disable if the right is power like to avoid infinite loop.
1965 if (is_variable_promotable(this->left, this->right)) {
1966 return this->right*this->left;
1967 }
1968
1969// Move trig to the right.
1970 auto cl = cos_cast(this->left);
1971 auto sl = sin_cast(this->left);
1972 if ((cl.get() && !this->right->is_power_like() &&
1973 !this->right->is_all_variables() &&
1974 !sin_cast(this->right).get()) ||
1975 (sl.get() && !this->right->is_power_like() &&
1976 !this->right->is_all_variables()) ||
1977 (sl.get() && cos_cast(this->right).get())) {
1978 return this->right*this->left;
1979 }
1980
1981// Reduce x*x to x^2
1982 if (this->left->is_match(this->right)) {
1983 return pow(this->left, 2.0);
1984 }
1985
1986// Gather common terms.
1987 auto lm = multiply_cast(this->left);
1988 if (lm.get()) {
1989// Promote constants before variables.
1990// (c*v1)*v2 -> c*(v1*v2)
1991 if (is_constant_promotable(lm->get_left(),
1992 lm->get_right())) {
1993 return lm->get_left()*(lm->get_right()*this->right);
1994 }
1995
1996// (a^c*b)*a^d -> a^(c+d)*b
1997// (b*a^c)*a^d -> a^(c+d)*b
1998 if (is_variable_combinable(this->right, lm->get_left())) {
1999 return (this->right*lm->get_left())*lm->get_right();
2000 } else if (is_variable_combinable(this->right, lm->get_right())) {
2001 return (this->right*lm->get_right())*lm->get_left();
2002 }
2003
2004// Assume variables, sqrt of variables, and powers of variables are on the
2005// right.
2006// (a*v)*b -> (a*b)*v
2007 if (is_variable_promotable(lm->get_right(), this->right)) {
2008 return (lm->get_left()*this->right)*lm->get_right();
2009 }
2010
2011// (a*(b*c)^e)*c^f -> a*b^e*c^(e+f)
2012 auto lmrp = pow_cast(lm->get_right());
2013 if (lmrp.get()) {
2014 auto lmrplm = multiply_cast(lmrp->get_left());
2015 if (lmrplm.get() &&
2016 is_variable_combinable(lmrplm->get_right(),
2017 this->right)) {
2018 return (lm->get_left()*pow(lmrplm->get_left(),
2019 lmrp->get_right()))*pow(this->right->get_power_base(),
2020 lmrp->get_right() +
2021 this->right->get_power_exponent());
2022 }
2023 }
2024 }
2025
2026 auto rm = multiply_cast(this->right);
2027 if (rm.get()) {
2028// Assume constants are on the left.
2029// c1*(c2*v) -> c3*v
2030 if (is_constant_combinable(this->left,
2031 rm->get_left())) {
2032 auto temp = this->left*rm->get_left();
2033 if (temp->is_normal()) {
2034 return temp*rm->get_right();
2035 }
2036 }
2037
2038// a*(a*b) -> a^2*b
2039// a*(b*a) -> a^2*b
2040 if (is_variable_combinable(this->left, rm->get_left())) {
2041 return (this->left*rm->get_left())*rm->get_right();
2042 } else if (is_variable_combinable(this->left, rm->get_right())) {
2043 return (this->left*rm->get_right())*rm->get_left();
2044 }
2045
2046// Assume variables are on the left.
2047// a*(b*v) -> (a*b)*v
2048 if (is_variable_promotable(rm->get_right(), this->left)) {
2049 return (this->left*rm->get_left())*rm->get_right();
2050 }
2051
2052// c1*(fma(c2,x,c3)*y)-> fma(c4,x,c5)*y
2053// c1*(fma(fma(c2,x,c3),x,c4)*y)-> fma(fma(c5,x,c6),x,c7)*y
2054// c1*(fma(fma(fma(c2,x,c3),x,c4),x,c5)*y)-> fma(fma(fma(c6,x,c7),x,c8),x,c9)*y
2055// etc...
2056 auto temp = this->reduce_nested_fma_times_constant(rm->get_left());
2057 if (temp.get()) {
2058 return temp*rm->get_right();
2059 }
2060 }
2061
2062// v1*(c*v2) -> c*(v1*v2)
2063 if (rm.get() &&
2064 is_constant_promotable(rm->get_left(), this->left)) {
2065 return rm->get_left()*(this->left*rm->get_right());
2066 }
2067
2068// Assume trig on the right.
2069// a*(b*sin) -> (a*b)*sin
2070// a*(b*cos) -> (a*b)*cos
2071// (a*sin)*b -> (a*b)*sin
2072// (a*cos)*b -> (a*b)*cos
2073 if (lm.get() &&
2074 (sin_cast(lm->get_right()).get() ||
2075 cos_cast(lm->get_right()).get()) &&
2076 !sin_cast(this->right).get() &&
2077 !this->right->is_power_like()) {
2078 return (lm->get_left()*this->right)*lm->get_right();
2079 } else if (rm.get() &&
2080 (sin_cast(rm->get_right()).get() ||
2081 cos_cast(rm->get_right()).get()) &&
2082 !this->left->is_constant()) {
2083 return (this->left*rm->get_left())*rm->get_right();
2084 }
2085
2086// Factor out common constants c*b*c*d -> c*c*b*d. c*c will get reduced to c on
2087// the second pass.
2088 if (lm.get() && rm.get()) {
2089 if (is_constant_combinable(lm->get_left(),
2090 rm->get_left())) {
2091 auto temp = lm->get_left()*rm->get_left();
2092 if (temp->is_normal()) {
2093 return temp*(lm->get_right()*rm->get_right());
2094 }
2095 } else if (is_constant_combinable(lm->get_left(),
2096 rm->get_right())) {
2097 auto temp = lm->get_left()*rm->get_right();
2098 if (temp->is_normal()) {
2099 return temp*(lm->get_right()*rm->get_left());
2100 }
2101 } else if (is_constant_combinable(lm->get_right(),
2102 rm->get_left())) {
2103 auto temp = lm->get_right()*rm->get_left();
2104 if (temp->is_normal()) {
2105 return temp*(lm->get_left()*rm->get_right());
2106 }
2107 } else if (is_constant_combinable(lm->get_right(),
2108 rm->get_right())) {
2109 auto temp = lm->get_right()*rm->get_right();
2110 if (temp->is_normal()) {
2111 return temp*(lm->get_left()*rm->get_left());
2112 }
2113 }
2114
2115// Gather common terms. This will help reduce sqrt(a)*sqrt(a).
2116 if (lm->get_left()->is_match(rm->get_left())) {
2117 return (lm->get_left()*rm->get_left()) *
2118 (lm->get_right()*rm->get_right());
2119 } else if (lm->get_right()->is_match(rm->get_left())) {
2120 return (lm->get_right()*rm->get_left()) *
2121 (lm->get_left()*rm->get_right());
2122 } else if (lm->get_left()->is_match(rm->get_right())) {
2123 return (lm->get_left()*rm->get_right()) *
2124 (lm->get_right()*rm->get_left());
2125 } else if (lm->get_right()->is_match(rm->get_right())) {
2126 return (lm->get_right()*rm->get_right()) *
2127 (lm->get_left()*rm->get_left());
2128 }
2129 }
2130
2131// Common factor reduction. (a/b)*(c/a) = c/b.
2132 auto ld = divide_cast(this->left);
2133 auto rd = divide_cast(this->right);
2134
2135// a*(b/c) -> (a*b)/c
2136// (a/c)*b -> (a*b)/c
2137 if (rd.get()) {
2138 return (this->left*rd->get_left())/rd->get_right();
2139 } else if (ld.get()) {
2140 return (ld->get_left()*this->right)/ld->get_right();
2141 }
2142
2143// (a/b)*(c/a) -> c/b
2144// (b/a)*(a/c) -> c/b
2145 if (ld.get() && rd.get()) {
2146 if (ld->get_left()->is_match(rd->get_right())) {
2147 return rd->get_left()/ld->get_right();
2148 } else if (ld->get_right()->is_match(rd->get_left())) {
2149 return ld->get_left()/rd->get_right();
2150 }
2151
2152// Convert (a/b)*(c/d) -> (a*c)/(b*d). This should help reduce cases like.
2153// (a/b)*(a/b) + (c/b)*(c/b).
2154 return (ld->get_left()*rd->get_left()) /
2155 (ld->get_right()*rd->get_right());
2156 }
2157
2158// Power reductions.
2159 if (is_variable_combinable(this->left, this->right)) {
2160 return pow(this->left->get_power_base(),
2161 this->left->get_power_exponent() +
2162 this->right->get_power_exponent());
2163 }
2164
2165// a*b^-c -> a/b^c
2166 auto rp = pow_cast(this->right);
2167 if (rp.get()) {
2168 auto exponent = constant_cast(rp->get_right());
2169 if (exponent.get() && exponent->evaluate().is_negative()) {
2170 return this->left/pow(rp->get_left(), -rp->get_right());
2171 }
2172 }
2173// b^-c*a -> a/b^c
2174 auto lp = pow_cast(this->left);
2175 if (lp.get()) {
2176 auto exponent = constant_cast(lp->get_right());
2177 if (exponent.get() && exponent->evaluate().is_negative()) {
2178 return this->right/pow(lp->get_left(), -lp->get_right());
2179 }
2180 }
2181// a^b*c^b -> (a*c)^b
2182 if (lp.get() && rp.get()) {
2183 if (lp->get_right()->is_match(rp->get_right())) {
2184 return pow(lp->get_left()*rp->get_left(), lp->get_right());
2185 }
2186 }
2187// (a*b^c)*d^c -> a*(b*d)^c
2188// (a^c*b)*d^c -> b*(a*d)^c
2189// a^c*(b*d^c) -> b*(a*d)^c
2190// a^c*(b^c*d) -> d*(a*b)^c
2191 if (lm.get() && rp.get()) {
2192 auto lmlp = pow_cast(lm->get_left());
2193 auto lmrp = pow_cast(lm->get_right());
2194 if (lmrp.get()) {
2195 if (lmrp->get_right()->is_match(rp->get_right())) {
2196 return lm->get_left()*pow(lmrp->get_left()*rp->get_left(),
2197 rp->get_right());
2198 }
2199 } else if (lmlp.get()) {
2200 if (lmlp->get_right()->is_match(rp->get_right())) {
2201 return lm->get_right()*pow(lmlp->get_left()*rp->get_left(),
2202 rp->get_right());
2203 }
2204 }
2205 } else if (rm.get() && lp.get()) {
2206 auto rmlp = pow_cast(rm->get_left());
2207 auto rmrp = pow_cast(rm->get_right());
2208 if (rmrp.get()) {
2209 if (rmrp->get_right()->is_match(lp->get_right())) {
2210 return rm->get_left()*pow(lp->get_left()*rmrp->get_left(),
2211 lp->get_right());
2212 }
2213 } else if (rmlp.get()) {
2214 if (rmlp->get_right()->is_match(lp->get_right())) {
2215 return rm->get_right()*pow(lp->get_left()*rmlp->get_left(),
2216 lp->get_right());
2217 }
2218 }
2219 }
2220
2221// (b*a)^c*a^d -> b^c*a^(c + d)
2222// (a*b)^c*a^d -> b^c*a^(c + d)
2223// a^d*(b*a)^c -> b^c*a^(c + d)
2224// a^d*(a*b)^c -> b^c*a^(c + d)
2225 if (lp.get() && rp.get()) {
2226 auto lplm = multiply_cast(lp->get_left());
2227 auto rplm = multiply_cast(rp->get_left());
2228 if (lplm.get()) {
2229 if (is_variable_combinable(lplm->get_right(),
2230 this->right)) {
2231 return pow(lplm->get_left()->get_power_base(),
2232 this->left->get_power_exponent())*
2233 pow(this->right->get_power_base(),
2234 this->left->get_power_exponent() +
2235 this->right->get_power_exponent());
2236 } else if (is_variable_combinable(lplm->get_left(),
2237 this->right)) {
2238 return pow(lplm->get_right()->get_power_base(),
2239 this->left->get_power_exponent())*
2240 pow(this->right->get_power_base(),
2241 this->left->get_power_exponent() +
2242 this->right->get_power_exponent());
2243 }
2244 }
2245
2246 if (rplm.get()) {
2247 if (is_variable_combinable(rplm->get_right(),
2248 this->left)) {
2249 return pow(rplm->get_left()->get_power_base(),
2250 this->right->get_power_exponent())*
2251 pow(this->left->get_power_base(),
2252 this->left->get_power_exponent() +
2253 this->right->get_power_exponent());
2254 } else if (is_variable_combinable(rplm->get_left(),
2255 this->left)) {
2256 return pow(rplm->get_right()->get_power_base(),
2257 this->right->get_power_exponent())*
2258 pow(this->left->get_power_base(),
2259 this->left->get_power_exponent() +
2260 this->right->get_power_exponent());
2261 }
2262 }
2263 }
2264
2265 auto lpd = divide_cast(this->left->get_power_base());
2266 if (lpd.get()) {
2267// (a/b)^c*b^d -> a^c*b^(c-d)
2268 if (is_variable_combinable(lpd->get_right(),
2269 this->right)) {
2270 return pow(lpd->get_left(), this->left->get_power_exponent()) *
2271 pow(this->right->get_power_base(),
2272 this->right->get_power_exponent() -
2273 this->left->get_power_exponent()*lpd->get_right()->get_power_exponent());
2274 }
2275// (b/a)^c*b^d -> b^(c+d)/a^c
2276 if (is_variable_combinable(lpd->get_left(), this->right)) {
2277 return pow(this->right->get_power_base(),
2278 this->right->get_power_exponent() +
2279 this->left->get_power_exponent()*lpd->get_left()->get_power_exponent()) /
2280 pow(lpd->get_right(), this->left->get_power_exponent());
2281 }
2282 }
2283 auto rpd = divide_cast(this->right->get_power_base());
2284 if (rpd.get()) {
2285// b^d*(a/b)^c -> a^c*b^(c-d)
2286 if (is_variable_combinable(rpd->get_right(),
2287 this->left)) {
2288 return pow(rpd->get_left(), this->right->get_power_exponent()) *
2289 pow(this->left->get_power_base(),
2290 this->left->get_power_exponent() -
2291 this->right->get_power_exponent()*rpd->get_right()->get_power_exponent());
2292 }
2293// b^d*(b/a)^c -> b^(c+d)/a^c
2294 if (is_variable_combinable(rpd->get_left(),
2295 this->left)) {
2296 return pow(this->right->get_power_base(),
2297 this->right->get_power_exponent() +
2298 this->right->get_power_exponent()*rpd->get_left()->get_power_exponent()) /
2299 pow(rpd->get_right(), this->right->get_power_exponent());
2300 }
2301 }
2302
2303// exp(a)*exp(b) -> exp(a + b)
2304 auto le = exp_cast(this->left);
2305 auto re = exp_cast(this->right);
2306 if (le.get() && re.get()) {
2307 return exp(le->get_arg() + re->get_arg());
2308 }
2309
2310// exp(a)*(exp(b)*c) -> c*(exp(a)*exp(b))
2311// exp(a)*(c*exp(b)) -> c*(exp(a)*exp(b))
2312 if (le.get() && rm.get()) {
2313 auto rmle = exp_cast(rm->get_left());
2314 if (rmle.get()) {
2315 return rm->get_right()*(this->left*rm->get_left());
2316 }
2317 auto rmre = exp_cast(rm->get_right());
2318 if (rmre.get()) {
2319 return rm->get_left()*(this->left*rm->get_right());
2320 }
2321 }
2322// (exp(a)*c)*exp(b) -> c*(exp(a)*exp(b))
2323// (c*exp(a))*exp(b) -> c*(exp(a)*exp(b))
2324 if (re.get() && lm.get()) {
2325 auto lmle = exp_cast(lm->get_left());
2326 if (lmle.get()) {
2327 return lm->get_right()*(this->right*lm->get_left());
2328 }
2329 auto lmre = exp_cast(lm->get_right());
2330 if (lmre.get()) {
2331 return lm->get_left()*(this->right*lm->get_right());
2332 }
2333 }
2334// (exp(a)*c)*(exp(b)*d) -> (c*d)*(exp(a)*exp(b))
2335// (exp(a)*c)*(d*exp(b)) -> (c*d)*(exp(a)*exp(b))
2336// (c*exp(a))*(exp(b)*d) -> (c*d)*(exp(a)*exp(b))
2337// (c*exp(a))*(d*exp(b)) -> (c*d)*(exp(a)*exp(b))
2338 if (lm.get() && rm.get()) {
2339 auto lmle = exp_cast(lm->get_left());
2340 if (lmle.get()) {
2341 auto rmle = exp_cast(rm->get_left());
2342 if (rmle.get()) {
2343 return (lm->get_right()*rm->get_right()) *
2344 (lm->get_left()*rm->get_left());
2345 }
2346 auto rmre = exp_cast(rm->get_right());
2347 if (rmre.get()) {
2348 return (lm->get_right()*rm->get_left()) *
2349 (lm->get_left()*rm->get_right());
2350 }
2351 }
2352 auto lmre = exp_cast(lm->get_right());
2353 if (lmre.get()) {
2354 auto rmle = exp_cast(rm->get_left());
2355 if (rmle.get()) {
2356 return (lm->get_left()*rm->get_right()) *
2357 (lm->get_right()*rm->get_left());
2358 }
2359 auto rmre = exp_cast(rm->get_right());
2360 if (rmre.get()) {
2361 return (lm->get_left()*rm->get_left()) *
2362 (lm->get_right()*rm->get_right());
2363 }
2364 }
2365 }
2366
2367 if (ld.get() && re.get()) {
2368// (c/exp(a))*exp(b) -> c*(exp(b)/exp(a))
2369 auto ldre = exp_cast(ld->get_right());
2370 if (ldre.get()) {
2371 return ld->get_left()*(this->right/ld->get_right());
2372 }
2373// (exp(a)/c)*exp(b) -> (exp(a)*exp(b))/c
2374 auto ldle = exp_cast(ld->get_left());
2375 if (ldle.get()) {
2376 return (ld->get_left()*this->right)/ld->get_right();
2377 }
2378 }
2379 if (rd.get() && le.get()) {
2380// exp(a)*(c/exp(a)) -> c*(exp(a)/exp(b))
2381 auto rdre = exp_cast(rd->get_right());
2382 if (rdre.get()) {
2383 return rd->get_left()*(this->left/rd->get_right());
2384 }
2385// exp(a)*(exp(b)/c) -> (exp(a)*exp(b))/c
2386 auto rdle = exp_cast(rd->get_left());
2387 if (rdle.get()) {
2388 return (this->left*rd->get_left())/rd->get_right();
2389 }
2390 }
2391
2392 if (ld.get() && rm.get()) {
2393 auto rmle = exp_cast(rm->get_left());
2394 if (rmle.get()) {
2395// (c/exp(a))*(exp(b)*d) -> (c*d)*(exp(b)/exp(a))
2396 auto ldre = exp_cast(ld->get_right());
2397 if (ldre.get()) {
2398 return (ld->get_left()*rm->get_right()) *
2399 (rm->get_left()/ld->get_right());
2400 }
2401// (exp(a)/c)*(exp(b)*d) -> (d/c)*(exp(a)*exp(b))
2402 auto ldle = exp_cast(ld->get_left());
2403 if (ldle.get()) {
2404 return (rm->get_right()/ld->get_right()) *
2405 (ld->get_left()*rm->get_left());
2406 }
2407 }
2408 auto rmre = exp_cast(rm->get_right());
2409 if (rmre.get()) {
2410// (c/exp(a))*(d*exp(b)) -> (c*d)*(exp(b)/exp(a))
2411 auto ldre = exp_cast(ld->get_right());
2412 if (ldre.get()) {
2413 return (ld->get_left()*rm->get_left()) *
2414 (rm->get_right()/ld->get_right());
2415 }
2416// (exp(a)/c)*(d*exp(b)) -> (d/c)*(exp(a)*exp(b))
2417 auto ldle = exp_cast(ld->get_left());
2418 if (ldle.get()) {
2419 return (rm->get_left()/ld->get_right()) *
2420 (ld->get_left()*rm->get_right());
2421 }
2422 }
2423 } else if (rd.get() && lm.get()) {
2424 auto lmre = exp_cast(lm->get_right());
2425 if (lmre.get()) {
2426// (c*exp(a))*(exp(b)/d) -> (c/d)*(exp(a)*exp(b))
2427 auto rdre = exp_cast(rd->get_left());
2428 if (rdre.get()) {
2429 return (lm->get_left()/rd->get_right()) *
2430 (lm->get_right()*rd->get_left());
2431 }
2432// (c*exp(a))*(d/exp(b)) -> (c*d)*(exp(a)/exp(b))
2433 auto rdle = exp_cast(rd->get_right());
2434 if (rdle.get()) {
2435 return (lm->get_left()*rd->get_left()) *
2436 (lm->get_right()/rd->get_right());
2437 }
2438 }
2439 auto lmle = exp_cast(lm->get_left());
2440 if (lmle.get()) {
2441// (exp(a)*c)*(d/exp(b)) -> (c*d)*(exp(a)/exp(b))
2442 auto rdle = exp_cast(rd->get_right());
2443 if (rdle.get()) {
2444 return (lm->get_right()*rd->get_left()) *
2445 (lm->get_left()/rd->get_right());
2446 }
2447// (exp(a)*c)*(exp(b)/d) -> (c/d)*(exp(a)*exp(b))
2448 auto rdre = exp_cast(rd->get_left());
2449 if (rdre.get()) {
2450 return (lm->get_right()/rd->get_right()) *
2451 (lm->get_left()*rd->get_left());
2452 }
2453 }
2454 }
2455
2456// c1*fma(c2,x,c3) -> fma(c4,x,c5)
2457// c1*fma(fma(c2,x,c3),x,c4) -> fma(fma(c5,x,c6),x,c7)
2458// c1*fma(fma(fma(c2,x,c3),x,c4),x,c5) -> fma(fma(fma(c6,x,c7),x,c8),x,c9)
2459// etc...
2460 auto fma_reduce = this->reduce_nested_fma_times_constant(this->right);
2461 if (fma_reduce.get()) {
2462 return fma_reduce;
2463 }
2464
2465// fma(c1,x,c2)*(c3 + x) -> fma(fma(c1,x,c4),x,c5)
2466// fma(fma(c1,x,c2),x,c3)*(c4 + x) -> fma(fma(fma(c1,x,c5),x,c6),x,c7)
2467// etc...
2468 auto ra = add_cast(this->right);
2469 if (ra.get()) {
2470 auto fma_expand = this->expand_nested_fma_times_add(this->left,
2471 ra);
2472 if (fma_expand.get()) {
2473 return fma_expand;
2474 }
2475 }
2476
2477// Sqrt(a)*Sqrt(b) -> Sqrt(a*b)
2478 auto lsqr = sqrt_cast(this->left);
2479 auto rsqr = sqrt_cast(this->right);
2480 if (lsqr.get() && rsqr.get()) {
2481 return sqrt(lsqr->get_arg()*rsqr->get_arg());
2482 }
2483
2484 if constexpr (std::floating_point<T>) {
2485// hypot(b,c)*Sqrt(a) -> Sqrt((b^2 + c^2)*a)
2486// Sqrt(a)*hypot(b,c) -> Sqrt(a*(b^2 + c^2))
2487 auto lhypot = hypot_cast(this->left);
2488 auto rhypot = hypot_cast(this->right);
2489 if (lhypot.get() && rsqr.get()) {
2490 return sqrt((pow(lhypot->get_left(), static_cast<T> (2)) +
2491 pow(lhypot->get_right(), static_cast<T> (2))) *
2492 rsqr->get_arg());
2493 } else if (rhypot.get() && lsqr.get()) {
2494 return sqrt((pow(rhypot->get_left(), static_cast<T> (2)) +
2495 pow(rhypot->get_right(), static_cast<T> (2))) *
2496 lsqr->get_arg());
2497 }
2498
2499// Sqrt(x^2)*copysign(1,x) -> x
2500// copysign(1,x)*Sqrt(x^2) -> x
2501 auto lcsc = copysign_cast(this->left);
2502 auto rcsc = copysign_cast(this->right);
2503 if (rcsc.get() && lsqr.get()) {
2504 auto rcsclc = constant_cast(rcsc->get_left());
2505 auto lsqrpc = pow_cast(lsqr->get_arg());
2506 if (lsqrpc.get() && rcsclc.get() && rcsclc->is(1) &&
2507 lsqrpc->get_left()->is_match(rcsc->get_right())) {
2508 return rcsc->get_right();
2509 }
2510 } else if (lcsc.get() && rsqr.get()) {
2511 auto lcsclc = constant_cast(lcsc->get_left());
2512 auto rsqrpc = pow_cast(rsqr->get_arg());
2513 if (rsqrpc.get() && lcsclc.get() && lcsclc->is(1) &&
2514 rsqrpc->get_left()->is_match(lcsc->get_right())) {
2515 return lcsc->get_right();
2516 }
2517 }
2518 }
2519
2520// (a + b/c)*c -> fma(a,c,b)
2521// c*(a + b/c) -> fma(a,c,b)
2522// (b/c + a)*c -> fma(a,c,b)
2523// c*(b/c + a) -> fma(a,c,b)
2524 auto la = add_cast(this->left);
2525 if (la.get()) {
2526 auto lald = divide_cast(la->get_left());
2527 auto lard = divide_cast(la->get_right());
2528 if (lald.get() && lald->get_right()->is_match(this->right)) {
2529 return fma(la->get_right(), this->right, lald->get_left());
2530 } else if (lard.get() && lard->get_right()->is_match(this->right)) {
2531 return fma(la->get_left(), this->right, lard->get_left());
2532 }
2533 } else if (ra.get()) {
2534 auto rald = divide_cast(ra->get_left());
2535 auto rard = divide_cast(ra->get_right());
2536 if (rald.get() && rald->get_right()->is_match(this->left)) {
2537 return fma(ra->get_right(), this->left, rald->get_left());
2538 } else if (rard.get() && rard->get_right()->is_match(this->left)) {
2539 return fma(ra->get_left(), this->left, rard->get_left());
2540 }
2541 }
2542
2543// (a - b/c)*c -> a*c - b
2544// c*(a - b/c) -> a*c - b
2545// (b/c - a)*c -> b - a*c
2546// c*(b/c - a) -> b - a*c
2547 auto ls = subtract_cast(this->left);
2548 auto rs = subtract_cast(this->right);
2549 if (ls.get()) {
2550 auto lsld = divide_cast(ls->get_left());
2551 auto lsrd = divide_cast(ls->get_right());
2552 if (lsld.get() && lsld->get_right()->is_match(this->right)) {
2553 return lsld->get_left() - ls->get_right()*this->right;
2554 } else if (lsrd.get() && lsrd->get_right()->is_match(this->right)) {
2555 return ls->get_left()*this->right - lsrd->get_left();
2556 }
2557 } else if (rs.get()) {
2558 auto rsld = divide_cast(rs->get_left());
2559 auto rsrd = divide_cast(rs->get_right());
2560 if (rsld.get() && rsld->get_right()->is_match(this->left)) {
2561 return rsld->get_right() - rs->get_right()*this->left;
2562 } else if (rsrd.get() && rsrd->get_right()->is_match(this->left)) {
2563 return rs->get_left()*this->left - rsld->get_left();
2564 }
2565 }
2566
2567// Cases like
2568// (c/exp(a))*(exp(b)/d) -> (c/d)*(exp(b)/exp(a))
2569// (c/exp(a))*(d/exp(b)) -> (c*e)/(exp(b)*exp(a))
2570// (exp(a)/c)*(d/exp(b)) -> (d/c)*(exp(a)/exp(b))
2571// (exp(a)/c)*(exp(b)/d) -> (exp(a)*exp(b))/(c*d)
2572// Are taken care of by (a/b)*(c/d) -> (a*c)/(b*d) conversion above.
2573
2574 return this->shared_from_this();
2575 }
2576
2577//------------------------------------------------------------------------------
2584//------------------------------------------------------------------------------
2585 virtual shared_leaf<T, SAFE_MATH> df(shared_leaf<T, SAFE_MATH> x) {
2586 if (this->is_match(x)) {
2587 return one<T, SAFE_MATH> ();
2588 }
2589
2590 const size_t hash = reinterpret_cast<size_t> (x.get());
2591 if (this->df_cache.find(hash) == this->df_cache.end()) {
2592 this->df_cache[hash] = this->left->df(x)*this->right
2593 + this->left*this->right->df(x);
2594 }
2595 return this->df_cache[hash];
2596 }
2597
2598//------------------------------------------------------------------------------
2606//------------------------------------------------------------------------------
2607 virtual shared_leaf<T, SAFE_MATH>
2608 compile(std::ostringstream &stream,
2609 jit::register_map &registers,
2610 const jit::register_map &thread_mem,
2611 const jit::register_usage &usage) {
2612 if (registers.find(this) == registers.end()) {
2613 auto l = this->left->compile(stream, registers,
2614 thread_mem, usage);
2615 auto r = this->right->compile(stream, registers,
2616 thread_mem, usage);
2617
2618 registers[this] = jit::to_string('r', this);
2619 stream << " const ";
2620 jit::add_type<T> (stream);
2621 stream << " " << registers[this] << " = ";
2622 if constexpr (SAFE_MATH) {
2623 stream << "(" << registers[l.get()] << " == ";
2624 if constexpr (jit::complex_scalar<T>) {
2625 jit::add_type<T> (stream);
2626 stream << "(0, 0)";
2627 } else {
2628 stream << "0";
2629 }
2630 stream << " || " << registers[r.get()] << " == ";
2631 if constexpr (jit::complex_scalar<T>) {
2632 jit::add_type<T> (stream);
2633 stream << "(0, 0)";
2634 } else {
2635 stream << "0";
2636 }
2637 stream << ") ? ";
2638 if constexpr (jit::complex_scalar<T>) {
2639 jit::add_type<T> (stream);
2640 stream << "(0, 0)";
2641 } else {
2642 stream << "0";
2643 }
2644 stream << " : ";
2645 }
2646 stream << registers[l.get()] << "*"
2647 << registers[r.get()];
2648 this->endline(stream, usage);
2649 }
2650
2651 return this->shared_from_this();
2652 }
2653
2654//------------------------------------------------------------------------------
2659//------------------------------------------------------------------------------
2660 virtual bool is_match(shared_leaf<T, SAFE_MATH> x) {
2661 if (this == x.get()) {
2662 return true;
2663 }
2664
2665 auto x_cast = multiply_cast(x);
2666 if (x_cast.get()) {
2667// Multiplication is commutative.
2668 if ((this->left->is_match(x_cast->get_left()) &&
2669 this->right->is_match(x_cast->get_right())) ||
2670 (this->right->is_match(x_cast->get_left()) &&
2671 this->left->is_match(x_cast->get_right()))) {
2672 return true;
2673 }
2674 }
2675
2676 return false;
2677 }
2678
2679//------------------------------------------------------------------------------
2681//------------------------------------------------------------------------------
2682 virtual void to_latex() const {
2683 if (constant_cast(this->left).get() ||
2684 add_cast(this->left).get() ||
2685 subtract_cast(this->left).get()) {
2686 std::cout << "\\left(";
2687 this->left->to_latex();
2688 std::cout << "\\right)";
2689 } else {
2690 this->left->to_latex();
2691 }
2692 std::cout << " ";
2693 if (constant_cast(this->right).get() ||
2694 add_cast(this->right).get() ||
2695 subtract_cast(this->right).get()) {
2696 std::cout << "\\left(";
2697 this->right->to_latex();
2698 std::cout << "\\right)";
2699 } else {
2700 this->right->to_latex();
2701 }
2702 }
2703
2704//------------------------------------------------------------------------------
2708//------------------------------------------------------------------------------
2709 virtual shared_leaf<T, SAFE_MATH> remove_pseudo() {
2710 if (this->has_pseudo()) {
2711 return this->left->remove_pseudo() *
2712 this->right->remove_pseudo();
2713 }
2714 return this->shared_from_this();
2715 }
2716
2717//------------------------------------------------------------------------------
2723//------------------------------------------------------------------------------
2724 virtual shared_leaf<T, SAFE_MATH> to_vizgraph(std::stringstream &stream,
2725 jit::register_map &registers) {
2726 if (registers.find(this) == registers.end()) {
2727 const std::string name = jit::to_string('r', this);
2728 registers[this] = name;
2729 stream << " " << name
2730 << " [label = \"⨉\", shape = oval, style = filled, fillcolor = blue, fontcolor = white];" << std::endl;
2731
2732 auto l = this->left->to_vizgraph(stream, registers);
2733 stream << " " << name << " -- " << registers[l.get()] << ";" << std::endl;
2734 auto r = this->right->to_vizgraph(stream, registers);
2735 stream << " " << name << " -- " << registers[r.get()] << ";" << std::endl;
2736 }
2737
2738 return this->shared_from_this();
2739 }
2740 };
2741
2742//------------------------------------------------------------------------------
2750//------------------------------------------------------------------------------
2751 template<jit::float_scalar T, bool SAFE_MATH=false>
2752 shared_leaf<T, SAFE_MATH> multiply(shared_leaf<T, SAFE_MATH> l,
2753 shared_leaf<T, SAFE_MATH> r) {
2754 auto temp = std::make_shared<multiply_node<T, SAFE_MATH>> (l, r)->reduce();
2755// Test for hash collisions.
2756 for (size_t i = temp->get_hash();
2757 i < std::numeric_limits<size_t>::max(); i++) {
2758 if (leaf_node<T, SAFE_MATH>::caches.nodes.find(i) ==
2759 leaf_node<T, SAFE_MATH>::caches.nodes.end()) {
2760 leaf_node<T, SAFE_MATH>::caches.nodes[i] = temp;
2761 return temp;
2762 } else if (temp->is_match(leaf_node<T, SAFE_MATH>::caches.nodes[i])) {
2763 return leaf_node<T, SAFE_MATH>::caches.nodes[i];
2764 }
2765 }
2766#if defined(__clang__) || defined(__GNUC__)
2767 __builtin_unreachable();
2768#else
2769 assert(false && "Should never reach.");
2770#endif
2771 }
2772
2773//------------------------------------------------------------------------------
2784//------------------------------------------------------------------------------
2785 template<jit::float_scalar T, bool SAFE_MATH=false>
2790
2791//------------------------------------------------------------------------------
2803//------------------------------------------------------------------------------
2804 template<jit::float_scalar T, jit::float_scalar L, bool SAFE_MATH=false>
2809
2810//------------------------------------------------------------------------------
2822//------------------------------------------------------------------------------
2823 template<jit::float_scalar T, jit::float_scalar R, bool SAFE_MATH=false>
2828
2830 template<jit::float_scalar T, bool SAFE_MATH=false>
2831 using shared_multiply = std::shared_ptr<multiply_node<T, SAFE_MATH>>;
2832
2833//------------------------------------------------------------------------------
2841//------------------------------------------------------------------------------
2842 template<jit::float_scalar T, bool SAFE_MATH=false>
2843 shared_multiply<T, SAFE_MATH> multiply_cast(shared_leaf<T, SAFE_MATH> x) {
2844 return std::dynamic_pointer_cast<multiply_node<T, SAFE_MATH>> (x);
2845 }
2846
2847//******************************************************************************
2848// Divide node.
2849//******************************************************************************
2850//------------------------------------------------------------------------------
2855//------------------------------------------------------------------------------
2856 template<jit::float_scalar T, bool SAFE_MATH=false>
2857 class divide_node final : public branch_node<T, SAFE_MATH> {
2858 private:
2859//------------------------------------------------------------------------------
2865//------------------------------------------------------------------------------
2866 static std::string to_string(leaf_node<T, SAFE_MATH> *l,
2867 leaf_node<T, SAFE_MATH> *r) {
2868 return jit::format_to_string(reinterpret_cast<size_t> (l)) + "/" +
2869 jit::format_to_string(reinterpret_cast<size_t> (r));
2870 }
2871
2872 public:
2873//------------------------------------------------------------------------------
2878//------------------------------------------------------------------------------
2883
2884//------------------------------------------------------------------------------
2890//------------------------------------------------------------------------------
2891 virtual backend::buffer<T> evaluate() {
2892 backend::buffer<T> l_result = this->left->evaluate();
2893
2894// If all the elements on the left are zero, return the left side without
2895// reevaluating the right side. Stop this loop early once the first non zero
2896// element is encountered.
2897 if (l_result.is_zero()) {
2898 return l_result;
2899 }
2900
2901 backend::buffer<T> r_result = this->right->evaluate();
2902 return l_result/r_result;
2903 }
2904
2905//------------------------------------------------------------------------------
2909//------------------------------------------------------------------------------
2910 virtual shared_leaf<T, SAFE_MATH> reduce() {
2911// Constant Reductions.
2912 auto l = constant_cast(this->left);
2913 auto r = constant_cast(this->right);
2914
2915 if ((l.get() && l->is(0)) ||
2916 (r.get() && r->is(1))) {
2917 return this->left;
2918 } else if (l.get() && r.get()) {
2919 return constant<T, SAFE_MATH> (this->evaluate());
2920 }
2921
2922 auto pl1 = piecewise_1D_cast(this->left);
2923 auto pr1 = piecewise_1D_cast(this->right);
2924
2925 if (pl1.get() && (r.get() || pl1->is_arg_match(this->right))) {
2926 return piecewise_1D(this->evaluate(), pl1->get_arg());
2927 } else if (pr1.get() && (l.get() || pr1->is_arg_match(this->left))) {
2928 return piecewise_1D(this->evaluate(), pr1->get_arg());
2929 }
2930
2931 auto pl2 = piecewise_2D_cast(this->left);
2932 auto pr2 = piecewise_2D_cast(this->right);
2933
2934 if (pl2.get() && (r.get() || pl2->is_arg_match(this->right))) {
2935 return piecewise_2D(this->evaluate(),
2936 pl2->get_num_columns(),
2937 pl2->get_left(),
2938 pl2->get_right());
2939 } else if (pr2.get() && (l.get() || pr2->is_arg_match(this->left))) {
2940 return piecewise_2D(this->evaluate(),
2941 pr2->get_num_columns(),
2942 pr2->get_left(),
2943 pr2->get_right());
2944 }
2945
2946// Combine 2D and 1D piecewise constants if a row or column matches.
2947 if (pr2.get() && pr2->is_row_match(this->left)) {
2948 backend::buffer<T> result = pl1->evaluate();
2949 result.divide_row(pr2->evaluate());
2950 return piecewise_2D(result,
2951 pr2->get_num_columns(),
2952 pr2->get_left(),
2953 pr2->get_right());
2954 } else if (pr2.get() && pr2->is_col_match(this->left)) {
2955 backend::buffer<T> result = pl1->evaluate();
2956 result.divide_col(pr2->evaluate());
2957 return piecewise_2D(result,
2958 pr2->get_num_columns(),
2959 pr2->get_left(),
2960 pr2->get_right());
2961 } else if (pl2.get() && pl2->is_row_match(this->right)) {
2962 backend::buffer<T> result = pl2->evaluate();
2963 result.divide_row(pr1->evaluate());
2964 return piecewise_2D(result,
2965 pl2->get_num_columns(),
2966 pl2->get_left(),
2967 pl2->get_right());
2968 } else if (pl2.get() && pl2->is_col_match(this->right)) {
2969 backend::buffer<T> result = pl2->evaluate();
2970 result.divide_col(pr1->evaluate());
2971 return piecewise_2D(result,
2972 pl2->get_num_columns(),
2973 pl2->get_left(),
2974 pl2->get_right());
2975 }
2976
2977 if (this->left->is_match(this->right)) {
2978 return one<T, SAFE_MATH> ();
2979 }
2980
2981// Reduce cases of a/c1 -> c2*a
2982 if (this->right->is_constant()) {
2983 return (1.0/this->right)*this->left;
2984 }
2985
2986// a/(b/c + d) -> a*c/(c*d + b)
2987// a/(d + b/c) -> a*c/(c*d + b)
2988 auto ra = add_cast(this->right);
2989 if (ra.get()) {
2990 auto rald = divide_cast(ra->get_left());
2991 auto rard = divide_cast(ra->get_right());
2992 if (rald.get()) {
2993 return this->left*rald->get_right() /
2994 fma(rald->get_right(),
2995 ra->get_right(),
2996 rald->get_left());
2997 } else if (rard.get()) {
2998 return this->left*rard->get_right() /
2999 fma(rard->get_right(),
3000 ra->get_left(),
3001 rard->get_left());
3002 }
3003 }
3004
3005// a/(b/c - d) -> a*c/(b - c*d)
3006// a/(d - b/c) -> a*c/(c*d - b)
3007 auto rs = subtract_cast(this->right);
3008 if (rs.get()) {
3009 auto rsld = divide_cast(rs->get_left());
3010 auto rsrd = divide_cast(rs->get_right());
3011 if (rsld.get()) {
3012 return this->left*rsld->get_right() /
3013 (rsld->get_left() -
3014 rsld->get_right()*rs->get_right());
3015 } else if (rsrd.get()) {
3016 return this->left*rsrd->get_right() /
3017 (rsrd->get_right()*rs->get_left() -
3018 rsrd->get_left());
3019 }
3020 }
3021
3022// fma(a,d,c*d)/d -> a + c
3023// fma(a,d,d*c)/d -> a + c
3024// fma(d,a,c*d)/d -> a + c
3025// fma(d,a,d*c)/d -> a + c
3026 auto lfma = fma_cast(this->left);
3027 if (lfma.get()) {
3028 auto fmarm = multiply_cast(lfma->get_right());
3029 if (fmarm.get()) {
3030 if (lfma->get_middle()->is_match(this->right) &&
3031 fmarm->get_right()->is_match(this->right)) {
3032 return lfma->get_left() + fmarm->get_left();
3033 } else if (lfma->get_middle()->is_match(this->right) &&
3034 fmarm->get_left()->is_match(this->right)) {
3035 return lfma->get_left() + fmarm->get_right();
3036 } else if (lfma->get_left()->is_match(this->right) &&
3037 fmarm->get_right()->is_match(this->right)) {
3038 return lfma->get_middle() + fmarm->get_left();
3039 } else if (lfma->get_left()->is_match(this->right) &&
3040 fmarm->get_left()->is_match(this->right)) {
3041 return lfma->get_middle() + fmarm->get_right();
3042 }
3043 }
3044 }
3045
3046// Common factor reduction. (a*b)/(a*c) = b/c.
3047 auto lm = multiply_cast(this->left);
3048 auto rm = multiply_cast(this->right);
3049
3050 if (lm.get() && rm.get()) {
3051 if (is_variable_combinable(lm->get_left(),
3052 rm->get_left()) ||
3053 is_variable_combinable(lm->get_right(),
3054 rm->get_right())) {
3055 return (lm->get_left()/rm->get_left()) *
3056 (lm->get_right()/rm->get_right());
3057 } else if (is_variable_combinable(lm->get_left(),
3058 rm->get_right()) ||
3059 is_variable_combinable(lm->get_right(),
3060 rm->get_left())) {
3061 return (lm->get_left()/rm->get_right()) *
3062 (lm->get_right()/rm->get_left());
3063 }
3064 }
3065
3066// Move constants to the numerator.
3067// a/(c1*b) -> (c2*a)/b
3068// a/(b*c1) -> (c2*a)/b
3069 if (rm.get()) {
3070 if (rm->get_left()->is_constant() &&
3071 rm->get_left()->is_normal()) {
3072 return ((1.0/rm->get_left())*this->left)/rm->get_right();
3073 } else if (rm->get_right()->is_constant() &&
3074 rm->get_right()->is_normal()) {
3075 return ((1.0/rm->get_right())*this->left)/rm->get_left();
3076 }
3077
3078// a/((b/c + d)*e) -> a*c/((c*d + b)*e)
3079// a/((d + b/c)*e) -> a*c/((c*d + b)*e)
3080// a/(e*(b/c + d)) -> a*c/((c*d + b)*e)
3081// a/(e*(d + b/c)) -> a*c/((c*d + b)*e)
3082 auto rmla = add_cast(rm->get_left());
3083 auto rmra = add_cast(rm->get_right());
3084 if (rmla.get()) {
3085 auto rmlald = divide_cast(rmla->get_left());
3086 auto rmlard = divide_cast(rmla->get_right());
3087 if (rmlald.get()) {
3088 return this->left*rmlald->get_right() /
3089 (fma(rmlald->get_right(),
3090 rmla->get_right(),
3091 rmlald->get_left())*rm->get_right());
3092 } else if (rmlard.get()) {
3093 return this->left*rmlard->get_right() /
3094 (fma(rmlard->get_right(),
3095 rmla->get_left(),
3096 rmlard->get_left())*rm->get_right());
3097 }
3098 }
3099 if (rmra.get()) {
3100 auto rmrald = divide_cast(rmra->get_left());
3101 auto rmrard = divide_cast(rmra->get_right());
3102 if (rmrald.get()) {
3103 return this->left*rmrald->get_right() /
3104 (fma(rmrald->get_right(),
3105 rmra->get_right(),
3106 rmrald->get_left())*rm->get_left());
3107 } else if (rmrard.get()) {
3108 return this->left*rmrard->get_right() /
3109 (fma(rmrard->get_right(),
3110 rmra->get_left(),
3111 rmrard->get_left())*rm->get_left());
3112 }
3113 }
3114
3115// a/((b/c - d)*e) -> a*c/((b - c*d)*e)
3116// a/(e*(b/c - d)) -> a*c/((b - c*d)*e)
3117// a/((d - b/c)*e) -> a*c/((c*d - b)*e)
3118// a/(e*(d - b/c)) -> a*c/((c*d - b)*e)
3119 auto rmls = subtract_cast(rm->get_left());
3120 auto rmrs = subtract_cast(rm->get_right());
3121 if (rmls.get()) {
3122 auto rmlsld = divide_cast(rmls->get_left());
3123 auto rmlsrd = divide_cast(rmls->get_right());
3124 if (rmlsld.get()) {
3125 return this->left*rmlsld->get_right() /
3126 ((rmlsld->get_left() -
3127 rmlsld->get_right()*rmls->get_right())*rm->get_right());
3128 } else if (rmlsrd.get()) {
3129 return this->left*rmlsrd->get_right() /
3130 ((rmlsrd->get_right()*rmls->get_left() -
3131 rmlsrd->get_left())*rm->get_right());
3132 }
3133 }
3134 if (rmrs.get()) {
3135 auto rmrsld = divide_cast(rmrs->get_left());
3136 auto rmrsrd = divide_cast(rmrs->get_right());
3137 if (rmrsld.get()) {
3138 return this->left*rmrsld->get_right() /
3139 ((rmrsld->get_left() -
3140 rmrsld->get_right()*rmrs->get_right())*rm->get_left());
3141 } else if (rmrsrd.get()) {
3142 return this->left*rmrsrd->get_right() /
3143 ((rmrsrd->get_right()*rmrs->get_left() -
3144 rmrsrd->get_left())*rm->get_left());
3145 }
3146 }
3147 }
3148
3149 if (lm.get() && rm.get()) {
3150// (a*b)/(a*c) -> b/c
3151// (b*a)/(a*c) -> b/c
3152// (a*b)/(c*a) -> b/c
3153// (b*a)/(c*a) -> b/c
3154 if (lm->get_left()->is_match(rm->get_left())) {
3155 return lm->get_right()/rm->get_right();
3156 } else if (lm->get_left()->is_match(rm->get_right())) {
3157 return lm->get_right()/rm->get_left();
3158 } else if (lm->get_right()->is_match(rm->get_left())) {
3159 return lm->get_left()/rm->get_right();
3160 } else if (lm->get_right()->is_match(rm->get_right())) {
3161 return lm->get_left()/rm->get_left();
3162 }
3163 }
3164
3165 if (lm.get()) {
3166// (v1*v2)/v1 -> v2
3167// (v2*v1)/v1 -> v2
3168 if (lm->get_left()->is_match(this->right)) {
3169 return lm->get_right();
3170 } else if (lm->get_right()->is_match(this->right)) {
3171 return lm->get_left();
3172 }
3173
3174// (v1^a*v2)/v1^b -> v2*(v1^a/v1^b)
3175// (v2*v1^a)/v1^b -> v2*(v1^a/v1^b)
3176 if (is_variable_combinable(lm->get_left(),
3177 this->right)) {
3178 return lm->get_right()*(lm->get_left()/this->right);
3179 } else if (is_variable_combinable(lm->get_right(),
3180 this->right)) {
3181 return lm->get_left()*(lm->get_right()/this->right);
3182 }
3183 }
3184
3185// (a/b)/c -> a/(b*c)
3186// a/(b/c) -> a*c/b
3187 auto ld = divide_cast(this->left);
3188 auto rd = divide_cast(this->right);
3189 if (ld.get()) {
3190 return ld->get_left()/(ld->get_right()*this->right);
3191 }
3192 if (rd.get()) {
3193 return this->left*rd->get_right()/rd->get_left();
3194 }
3195
3196// Power reductions.
3197 if (is_variable_combinable(this->left,
3198 this->right)) {
3199 return pow(this->left->get_power_base(),
3200 this->left->get_power_exponent() -
3201 this->right->get_power_exponent());
3202 }
3203
3204// a/b^-c -> a*b^c
3205 auto rp = pow_cast(this->right);
3206 if (rp.get()) {
3207 auto exponent = constant_cast(rp->get_right());
3208 if (exponent.get() && exponent->evaluate().is_negative()) {
3209 return this->left*pow(rp->get_left(), -rp->get_right());
3210 }
3211 }
3212
3213// (a*b)^c/(a^d) = a^(c - d)*b^c
3214// (b*a)^c/(a^d) = a^(c - d)*b^c
3215 auto lp = pow_cast(this->left);
3216 if (lp.get()) {
3217 auto lpm = multiply_cast(this->left->get_power_base());
3218 if (lpm.get()) {
3219 if (lpm->get_left()->is_match(this->right->get_power_base())) {
3220 return pow(this->right->get_power_base(),
3221 this->left->get_power_exponent() -
3222 this->right->get_power_exponent()) *
3223 pow(lpm->get_right(),
3224 this->left->get_power_exponent());
3225 } else if (lpm->get_right()->is_match(this->right->get_power_base())) {
3226 return pow(this->right->get_power_base(),
3227 this->left->get_power_exponent() -
3228 this->right->get_power_exponent()) *
3229 pow(lpm->get_left(),
3230 this->left->get_power_exponent());
3231 }
3232 }
3233 }
3234
3235// a^b/c^b -> (a/c)^b
3236 if (lp.get() && rp.get()) {
3237 if (lp->get_right()->is_match(rp->get_right())) {
3238 return pow(lp->get_left()/rp->get_left(), lp->get_right());
3239 }
3240 }
3241
3242// (a*b)^c/((a^d)*e) = a^(c - d)*b^c/e
3243// (b*a)^c/((a^d)*e) = a^(c - d)*b^c/e
3244// (a*b)^c/(e*(a^d)) = a^(c - d)*b^c/e
3245// (b*a)^c/(e*(a^d)) = a^(c - d)*b^c/e
3246 if (lp.get() && rm.get()) {
3247 auto lpm = multiply_cast(this->left->get_power_base());
3248 if (lpm.get()) {
3249 if (lpm->get_left()->is_match(rm->get_left()->get_power_base())) {
3250 return (pow(rm->get_left()->get_power_base(),
3251 this->left->get_power_exponent() -
3252 rm->get_left()->get_power_exponent()) *
3253 pow(lpm->get_right(),
3254 this->left->get_power_exponent())) /
3255 rm->get_right();
3256 } else if (lpm->get_right()->is_match(rm->get_left()->get_power_base())) {
3257 return (pow(rm->get_left()->get_power_base(),
3258 this->left->get_power_exponent() -
3259 rm->get_left()->get_power_exponent()) *
3260 pow(lpm->get_left(),
3261 this->left->get_power_exponent())) /
3262 rm->get_right();
3263 } else if (lpm->get_left()->is_match(rm->get_right()->get_power_base())) {
3264 return (pow(rm->get_right()->get_power_base(),
3265 this->left->get_power_exponent() -
3266 rm->get_right()->get_power_exponent()) *
3267 pow(lpm->get_right(),
3268 this->left->get_power_exponent())) /
3269 rm->get_left();
3270 } else if (lpm->get_right()->is_match(rm->get_right()->get_power_base())) {
3271 return (pow(rm->get_right()->get_power_base(),
3272 this->left->get_power_exponent() -
3273 rm->get_right()->get_power_exponent()) *
3274 pow(lpm->get_left(),
3275 this->left->get_power_exponent())) /
3276 rm->get_left();
3277 }
3278 }
3279 }
3280
3281 if (lm.get()) {
3282// a*(b*c)/c -> a*b
3283// a*(c*b)/c -> a*b
3284// (a*c)*b/c -> a*b
3285// (c*a)*b/c -> a*b
3286 auto lmrm = multiply_cast(lm->get_right());
3287 auto lmlm = multiply_cast(lm->get_left());
3288 if (lmrm.get()) {
3289 if (is_variable_combinable(lmrm->get_right(),
3290 this->right)) {
3291 return lm->get_left()*lmrm->get_left() *
3292 (lmrm->get_right()/this->right);
3293 } else if (is_variable_combinable(lmrm->get_left(),
3294 this->right)) {
3295 return lm->get_left()*lmrm->get_right() *
3296 (lmrm->get_left()/this->right);
3297 }
3298 } else if (lmlm.get()) {
3299 if (is_variable_combinable(lmlm->get_right(),
3300 this->right)) {
3301 return lm->get_right()*lmlm->get_left() *
3302 (lmlm->get_right()/this->right);
3303 } else if (is_variable_combinable(lmlm->get_left(),
3304 this->right)) {
3305 return lm->get_right()*lmlm->get_right() *
3306 (lmlm->get_left()/this->right);
3307 }
3308 }
3309
3310// (f*(a*b)^c)/(a^d) = f*a^(c - d)*b^c
3311// (f*(b*a)^c)/(a^d) = f*a^(c - d)*b^c
3312// (((a*b)^c)*f)/(a^d) = f*a^(c - d)*b^c
3313// (((b*a)^c)*f)/(a^d) = f*a^(c - d)*b^c
3314 auto lmlp = pow_cast(lm->get_left());
3315 auto lmrp = pow_cast(lm->get_right());
3316 if (lmlp.get()) {
3317 auto lmlpm = multiply_cast(lmlp->get_power_base());
3318 if (lmlpm.get()) {
3319 if (lmlpm->get_left()->is_match(this->right->get_power_base())) {
3320 return lm->get_right() *
3321 pow(this->right->get_power_base(),
3322 lmlp->get_power_exponent() -
3323 this->right->get_power_exponent()) *
3324 pow(lmlpm->get_right(),
3325 lmlp->get_power_exponent());
3326 } else if (lmlpm->get_right()->is_match(this->right->get_power_base())) {
3327 return lm->get_right() *
3328 pow(this->right->get_power_base(),
3329 lmlp->get_power_exponent() -
3330 this->right->get_power_exponent()) *
3331 pow(lmlpm->get_left(),
3332 lmlp->get_power_exponent());
3333 }
3334 }
3335 } else if (lmrp.get()) {
3336 auto lmrpm = multiply_cast(lmrp->get_power_base());
3337 if (lmrpm.get()) {
3338 if (lmrpm->get_left()->is_match(this->right->get_power_base())) {
3339 return lm->get_left() *
3340 pow(this->right->get_power_base(),
3341 lmrp->get_power_exponent() -
3342 this->right->get_power_exponent()) *
3343 pow(lmrpm->get_right(),
3344 lmrp->get_power_exponent());
3345 } else if (lmrpm->get_right()->is_match(this->right->get_power_base())) {
3346 return lm->get_left() *
3347 pow(this->right->get_power_base(),
3348 lmrp->get_power_exponent() -
3349 this->right->get_power_exponent()) *
3350 pow(lmrpm->get_left(),
3351 lmrp->get_power_exponent());
3352 }
3353 }
3354 }
3355 }
3356
3357// f*(a*b)^c/((a^d)*e) = a^(c - d)*b^c/e
3358// f*(b*a)^c/((a^d)*e) = a^(c - d)*b^c/e
3359// f*(a*b)^c/(e*(a^d)) = a^(c - d)*b^c/e
3360// f*(b*a)^c/(e*(a^d)) = a^(c - d)*b^c/e
3361// (a*b)^c*f/((a^d)*e) = a^(c - d)*b^c/e
3362// (b*a)^c*f/((a^d)*e) = a^(c - d)*b^c/e
3363// (a*b)^c*f/(e*(a^d)) = a^(c - d)*b^c/e
3364// (b*a)^c*f/(e*(a^d)) = a^(c - d)*b^c/e
3365 if (lm.get() && rm.get()) {
3366 auto lmlp = pow_cast(lm->get_left());
3367 auto lmrp = pow_cast(lm->get_right());
3368 if (lmlp.get()) {
3369 auto lmlpm = multiply_cast(lmlp->get_power_base());
3370 if (lmlpm.get()) {
3371 if (lmlpm->get_left()->is_match(rm->get_left()->get_power_base())) {
3372 return lm->get_right() *
3373 (pow(rm->get_left()->get_power_base(),
3374 lmlp->get_power_exponent() -
3375 rm->get_left()->get_power_exponent())) *
3376 pow(lmlpm->get_right(),
3377 lmlp->get_power_exponent()) /
3378 rm->get_right();
3379 } else if (lmlpm->get_right()->is_match(rm->get_left()->get_power_base())) {
3380 return lm->get_right() *
3381 (pow(rm->get_left()->get_power_base(),
3382 lmlp->get_power_exponent() -
3383 rm->get_left()->get_power_exponent())) *
3384 pow(lmlpm->get_left(),
3385 lmlp->get_power_exponent()) /
3386 rm->get_right();
3387 } else if (lmlpm->get_left()->is_match(rm->get_right()->get_power_base())) {
3388 return lm->get_right() *
3389 (pow(rm->get_left()->get_power_base(),
3390 lmlp->get_power_exponent() -
3391 rm->get_right()->get_power_exponent())) *
3392 pow(lmlpm->get_right(),
3393 lmlp->get_power_exponent()) /
3394 rm->get_left();
3395 } else if (lmlpm->get_right()->is_match(rm->get_right()->get_power_base())) {
3396 return lm->get_right() *
3397 (pow(rm->get_left()->get_power_base(),
3398 lmlp->get_power_exponent() -
3399 rm->get_right()->get_power_exponent())) *
3400 pow(lmlpm->get_left(),
3401 lmlp->get_power_exponent()) /
3402 rm->get_left();
3403 }
3404 }
3405 } else if (lmrp.get()) {
3406 auto lmrpm = multiply_cast(lmrp->get_power_base());
3407 if (lmrpm.get()) {
3408 if (lmrpm->get_left()->is_match(rm->get_left()->get_power_base())) {
3409 return lm->get_left() *
3410 (pow(rm->get_left()->get_power_base(),
3411 lmrp->get_power_exponent() -
3412 rm->get_left()->get_power_exponent())) *
3413 pow(lmrpm->get_right(),
3414 lmrp->get_power_exponent()) /
3415 rm->get_right();
3416 } else if (lmrpm->get_right()->is_match(rm->get_left()->get_power_base())) {
3417 return lm->get_left() *
3418 (pow(rm->get_left()->get_power_base(),
3419 lmrp->get_power_exponent() -
3420 rm->get_left()->get_power_exponent())) *
3421 pow(lmrpm->get_left(),
3422 lmrp->get_power_exponent()) /
3423 rm->get_right();
3424 } else if (lmrpm->get_left()->is_match(rm->get_right()->get_power_base())) {
3425 return lm->get_left() *
3426 (pow(rm->get_left()->get_power_base(),
3427 lmrp->get_power_exponent() -
3428 rm->get_right()->get_power_exponent())) *
3429 pow(lmrpm->get_right(),
3430 lmrp->get_power_exponent()) /
3431 rm->get_left();
3432 } else if (lmrpm->get_right()->is_match(rm->get_right()->get_power_base())) {
3433 return lm->get_left() *
3434 (pow(rm->get_left()->get_power_base(),
3435 lmrp->get_power_exponent() -
3436 rm->get_right()->get_power_exponent())) *
3437 pow(lmrpm->get_left(),
3438 lmrp->get_power_exponent()) /
3439 rm->get_left();
3440 }
3441 }
3442 }
3443 }
3444
3445// exp(a)/exp(b) -> exp(a - b)
3446 auto lexp = exp_cast(this->left);
3447 auto rexp = exp_cast(this->right);
3448 if (lexp.get() && rexp.get()) {
3449 return exp(lexp->get_arg() - rexp->get_arg());
3450 }
3451
3452// (c*exp(a))/exp(b) -> c*(exp(a)/exp(b))
3453// (exp(a)*c)/exp(b) -> c*(exp(a)/exp(b))
3454 if (rexp.get() && lm.get()) {
3455 auto lmre = exp_cast(lm->get_right());
3456 if (lmre.get()) {
3457 return lm->get_left()*(lm->get_right()/this->right);
3458 }
3459 auto lmle = exp_cast(lm->get_left());
3460 if (lmle.get()) {
3461 return lm->get_right()*(lm->get_left()/this->right);
3462 }
3463 }
3464// ((c*exp(a))*d)/exp(b)
3465// ((exp(a)*c)*d)/exp(b)
3466// (c*(exp(a)*d))/exp(b)
3467// (c*(d*exp(a)))/exp(b)
3468 if (rexp.get() && lm.get()) {
3469 auto lmlm = multiply_cast(lm->get_left());
3470 auto lmrm = multiply_cast(lm->get_right());
3471
3472 if (lmlm.get()) {
3473 if (exp_cast(lmlm->get_right()).get()) {
3474 return lmlm->get_left()*lm->get_right() *
3475 (lmlm->get_right()/this->right);
3476 } else if (exp_cast(lmlm->get_left()).get()) {
3477 return lmlm->get_right()*lm->get_right() *
3478 (lmlm->get_left()/this->right);
3479 }
3480 } else if (lmrm.get()) {
3481 if (exp_cast(lmrm->get_right()).get()) {
3482 return lmrm->get_left()*lm->get_left() *
3483 (lmrm->get_right()/this->right);
3484 } else if (exp_cast(lmrm->get_left()).get()) {
3485 return lmrm->get_right()*lm->get_left() *
3486 (lmrm->get_left()/this->right);
3487 }
3488 }
3489 }
3490
3491// exp(a)/(c*exp(b)) -> (exp(a)/exp(b))/c
3492// exp(a)/(exp(b)*c) -> (exp(a)/exp(b))/c
3493 if (lexp.get() && rm.get()) {
3494 auto rmre = exp_cast(rm->get_right());
3495 if (rmre.get()) {
3496 return (this->left/rm->get_right())/rm->get_left();
3497 }
3498 auto rmle = exp_cast(rm->get_left());
3499 if (rmle.get()) {
3500 return (this->left/rm->get_left())/rm->get_right();
3501 }
3502 }
3503
3504// (c*exp(a))/(d*exp(b)) -> (c/d)*(exp(a)/exp(b))
3505// (c*exp(a))/(exp(b)*d) -> (c/d)*(exp(a)/exp(b))
3506// (exp(a)*c)/(d*exp(b)) -> (c/d)*(exp(a)/exp(b))
3507// (exp(a)*c)/(exp(b)*d) -> (c/d)*(exp(a)/exp(b))
3508 if (lm.get() && rm.get()) {
3509 auto lmre = exp_cast(lm->get_right());
3510 if (lmre.get()) {
3511 auto rmre = exp_cast(rm->get_right());
3512 if (rmre.get()) {
3513 return (lm->get_left()/rm->get_left()) *
3514 (lm->get_right()/rm->get_right());
3515 }
3516 auto rmle = exp_cast(rm->get_left());
3517 if (rmle.get()) {
3518 return (lm->get_left()/rm->get_right()) *
3519 (lm->get_right()/rm->get_left());
3520 }
3521 }
3522 auto lmle = exp_cast(lm->get_left());
3523 if (lmle.get()) {
3524 auto rmre = exp_cast(rm->get_right());
3525 if (rmre.get()) {
3526 return (lm->get_right()/rm->get_left()) *
3527 (lm->get_left()/rm->get_right());
3528 }
3529 auto rmle = exp_cast(rm->get_left());
3530 if (rmle.get()) {
3531 return (lm->get_right()/rm->get_right()) *
3532 (lm->get_left()/rm->get_left());
3533 }
3534 }
3535 }
3536
3537// exp(a)/(c/exp(b)) -> (exp(a)*exp(b))/c
3538// exp(a)/(exp(b)/c) -> c*(exp(a)/exp(b))
3539 if (rd.get() && lexp.get()) {
3540 auto rdre = exp_cast(rd->get_right());
3541 if (rdre.get()) {
3542 return (this->left*rd->get_right())/rd->get_left();
3543 }
3544 auto rdle = exp_cast(rd->get_left());
3545 if (rdle.get()) {
3546 return rd->get_right()*(this->left/rd->get_left());
3547 }
3548 }
3549
3550// (c/exp(a))/exp(b) -> c/(exp(a)*exp(b))
3551// (exp(a)/c)/exp(b) -> exp(a)/(c*exp(b))
3552// (c/exp(a))/(d/exp(b)) -> (c*exp(b))/(d*exp(a))
3553// (c/exp(a))/(exp(b)/d) -> (c*d)/(exp(b)*exp(a))
3554// (exp(a)/c)/(d/exp(b)) -> (exp(a)*exp(b))/(d*c)
3555// (exp(a)/c)/(exp(b)/d) -> (exp(a)*d)/(exp(b)*c)
3556// Note cases like this are already transformed by the (a/b)/c -> a/(b*c)
3557// above.
3558
3559 return this->shared_from_this();
3560 }
3561
3562//------------------------------------------------------------------------------
3569//------------------------------------------------------------------------------
3570 virtual shared_leaf<T, SAFE_MATH>
3571 df(shared_leaf<T, SAFE_MATH> x) {
3572 if (this->is_match(x)) {
3573 return one<T, SAFE_MATH> ();
3574 }
3575
3576 const size_t hash = reinterpret_cast<size_t> (x.get());
3577 if (this->df_cache.find(hash) == this->df_cache.end()) {
3578 this->df_cache[hash] = this->left->df(x)/this->right
3579 - this->left*this->right->df(x)/(this->right*this->right);
3580 }
3581 return this->df_cache[hash];
3582 }
3583
3584//------------------------------------------------------------------------------
3592//------------------------------------------------------------------------------
3593 virtual shared_leaf<T, SAFE_MATH>
3594 compile(std::ostringstream &stream,
3595 jit::register_map &registers,
3596 const jit::register_map &thread_mem,
3597 const jit::register_usage &usage) {
3598 if (registers.find(this) == registers.end()) {
3599 auto l = this->left->compile(stream, registers,
3600 thread_mem, usage);
3601 auto r = this->right->compile(stream, registers,
3602 thread_mem, usage);
3603
3604 registers[this] = jit::to_string('r', this);
3605 stream << " const ";
3606 jit::add_type<T> (stream);
3607 stream << " " << registers[this] << " = ";
3608 if constexpr (SAFE_MATH) {
3609 stream << registers[l.get()] << " == ";
3610 if constexpr (jit::complex_scalar<T>) {
3611 jit::add_type<T> (stream);
3612 stream << "(0, 0)";
3613 } else {
3614 stream << "0";
3615 }
3616 stream << " ? ";
3617 if constexpr (jit::complex_scalar<T>) {
3618 jit::add_type<T> (stream);
3619 stream << "(0, 0)";
3620 } else {
3621 stream << "0";
3622 }
3623 stream << " : ";
3624 }
3625 stream << registers[l.get()] << "/"
3626 << registers[r.get()];
3627 this->endline(stream, usage);
3628 }
3629 return this->shared_from_this();
3630 }
3631
3632//------------------------------------------------------------------------------
3637//------------------------------------------------------------------------------
3638 virtual bool is_match(shared_leaf<T, SAFE_MATH> x) {
3639 if (this == x.get()) {
3640 return true;
3641 }
3642
3643 auto x_cast = divide_cast(x);
3644 if (x_cast.get()) {
3645 return this->left->is_match(x_cast->get_left()) &&
3646 this->right->is_match(x_cast->get_right());
3647 }
3648
3649 return false;
3650 }
3651
3652//------------------------------------------------------------------------------
3654//------------------------------------------------------------------------------
3655 virtual void to_latex() const {
3656 std::cout << "\\frac{";
3657 this->left->to_latex();
3658 std::cout << "}{";
3659 this->right->to_latex();
3660 std::cout << "}";
3661 }
3662
3663//------------------------------------------------------------------------------
3667//------------------------------------------------------------------------------
3668 virtual shared_leaf<T, SAFE_MATH> remove_pseudo() {
3669 if (this->has_pseudo()) {
3670 return this->left->remove_pseudo() /
3671 this->right->remove_pseudo();
3672 }
3673 return this->shared_from_this();
3674 }
3675
3676//------------------------------------------------------------------------------
3682//------------------------------------------------------------------------------
3683 virtual shared_leaf<T, SAFE_MATH> to_vizgraph(std::stringstream &stream,
3684 jit::register_map &registers) {
3685 if (registers.find(this) == registers.end()) {
3686 const std::string name = jit::to_string('r', this);
3687 registers[this] = name;
3688 stream << " " << name
3689 << " [label = \"\\\\\", shape = oval, style = filled, fillcolor = blue, fontcolor = white];" << std::endl;
3690
3691 auto l = this->left->to_vizgraph(stream, registers);
3692 stream << " " << name << " -- " << registers[l.get()] << ";" << std::endl;
3693 auto r = this->right->to_vizgraph(stream, registers);
3694 stream << " " << name << " -- " << registers[r.get()] << ";" << std::endl;
3695 }
3696
3697 return this->shared_from_this();
3698 }
3699 };
3700
3701//------------------------------------------------------------------------------
3709//------------------------------------------------------------------------------
3710 template<jit::float_scalar T, bool SAFE_MATH=false>
3711 shared_leaf<T, SAFE_MATH> divide(shared_leaf<T, SAFE_MATH> l,
3712 shared_leaf<T, SAFE_MATH> r) {
3713 auto temp = std::make_shared<divide_node<T, SAFE_MATH>> (l, r)->reduce();
3714// Test for hash collisions.
3715 for (size_t i = temp->get_hash();
3716 i < std::numeric_limits<size_t>::max(); i++) {
3717 if (leaf_node<T, SAFE_MATH>::caches.nodes.find(i) ==
3718 leaf_node<T, SAFE_MATH>::caches.nodes.end()) {
3719 leaf_node<T, SAFE_MATH>::caches.nodes[i] = temp;
3720 return temp;
3721 } else if (temp->is_match(leaf_node<T, SAFE_MATH>::caches.nodes[i])) {
3722 return leaf_node<T, SAFE_MATH>::caches.nodes[i];
3723 }
3724 }
3725#if defined(__clang__) || defined(__GNUC__)
3726 __builtin_unreachable();
3727#else
3728 assert(false && "Should never reach.");
3729#endif
3730 }
3731
3732//------------------------------------------------------------------------------
3743//------------------------------------------------------------------------------
3744 template<jit::float_scalar T, bool SAFE_MATH=false>
3749
3750//------------------------------------------------------------------------------
3762//------------------------------------------------------------------------------
3763 template<jit::float_scalar T, jit::float_scalar L, bool SAFE_MATH=false>
3768
3769//------------------------------------------------------------------------------
3781//------------------------------------------------------------------------------
3782 template<jit::float_scalar T, jit::float_scalar R, bool SAFE_MATH=false>
3787
3789 template<jit::float_scalar T, bool SAFE_MATH=false>
3790 using shared_divide = std::shared_ptr<divide_node<T, SAFE_MATH>>;
3791
3792//------------------------------------------------------------------------------
3800//------------------------------------------------------------------------------
3801 template<jit::float_scalar T, bool SAFE_MATH=false>
3802 shared_divide<T, SAFE_MATH> divide_cast(shared_leaf<T, SAFE_MATH> x) {
3803 return std::dynamic_pointer_cast<divide_node<T, SAFE_MATH>> (x);
3804 }
3805
3806//******************************************************************************
3807// fused multiply add node.
3808//******************************************************************************
3809//------------------------------------------------------------------------------
3816//------------------------------------------------------------------------------
3817 template<jit::float_scalar T, bool SAFE_MATH=false>
3818 class fma_node final : public triple_node<T, SAFE_MATH> {
3819 private:
3820//------------------------------------------------------------------------------
3827//------------------------------------------------------------------------------
3828 shared_leaf<T, SAFE_MATH>
3829 reduce_nested_fma(shared_subtract<T, SAFE_MATH> sub) {
3830 auto temp = fma_cast(this->left);
3831 if (temp.get()) {
3832 if (is_constant_combinable(sub->get_right(), temp->get_left()) &&
3833 is_constant_combinable(sub->get_right(), temp->get_right()) &&
3834 is_constant_combinable(this->right, temp->get_right()) &&
3835 temp->get_middle()->is_match(sub->get_left())) {
3836 return fma(fma(temp->get_left(),
3837 sub->get_left(),
3838 temp->get_right() - temp->get_left()*sub->get_right()),
3839 sub->get_left(),
3840 this->right - temp->get_right()*sub->get_right());
3841 } else {
3842 if (temp->get_middle()->is_match(sub->get_left()) &&
3843 is_constant_combinable(sub->get_right(), this->right)) {
3844 auto temp2 = temp->reduce_nested_fma(sub);
3845 if (temp2.get()) {
3846 return fma(temp2,
3847 sub->get_left(),
3848 this->right - temp->get_right()*sub->get_right());
3849 }
3850 }
3851 }
3852 }
3853 return this->shared_from_this();
3854 }
3855
3856//------------------------------------------------------------------------------
3863//------------------------------------------------------------------------------
3864 static std::string to_string(leaf_node<T, SAFE_MATH> *l,
3865 leaf_node<T, SAFE_MATH> *m,
3866 leaf_node<T, SAFE_MATH> *r) {
3867 return "fma" + jit::format_to_string(reinterpret_cast<size_t> (l))
3868 + jit::format_to_string(reinterpret_cast<size_t> (m))
3869 + jit::format_to_string(reinterpret_cast<size_t> (r));
3870 }
3871
3872 public:
3873//------------------------------------------------------------------------------
3879//------------------------------------------------------------------------------
3880 fma_node(shared_leaf<T, SAFE_MATH> l,
3881 shared_leaf<T, SAFE_MATH> m,
3882 shared_leaf<T, SAFE_MATH> r) :
3883 triple_node<T, SAFE_MATH> (l, m, r, fma_node::to_string(l.get(),
3884 m.get(),
3885 r.get())) {}
3886
3887//------------------------------------------------------------------------------
3893//------------------------------------------------------------------------------
3894 virtual backend::buffer<T> evaluate() {
3895 backend::buffer<T> l_result = this->left->evaluate();
3896 backend::buffer<T> r_result = this->right->evaluate();
3897
3898// If all the elements on the left are zero, return the left side without
3899// reevaluating the right side.
3900 if (l_result.is_zero()) {
3901 return r_result;
3902 }
3903
3904 backend::buffer<T> m_result = this->middle->evaluate();
3905 return backend::fma(l_result, m_result, r_result);
3906 }
3907
3908//------------------------------------------------------------------------------
3912//------------------------------------------------------------------------------
3913 virtual shared_leaf<T, SAFE_MATH> reduce() {
3914 auto l = constant_cast(this->left);
3915 auto m = constant_cast(this->middle);
3916 auto r = constant_cast(this->right);
3917
3918 if ((l.get() && l->is(0)) ||
3919 (m.get() && m->is(0))) {
3920 return this->right;
3921 } else if (r.get() && r->is(0)) {
3922 return this->left*this->middle;
3923 } else if (l.get() && m.get() && r.get()) {
3924 return constant<T, SAFE_MATH> (this->evaluate());
3925 } else if (l.get() && m.get()) {
3926 return this->left*this->middle + this->right;
3927 } else if (l.get() && l->is(-1)) {
3928 return this->right - this->middle;
3929 } else if (m.get() && m->is(-1)) {
3930 return this->right - this->left;
3931 } else if (l.get() && l->is(1)) {
3932 return this->middle + this->right;
3933 } else if (m.get() && m->is(1)) {
3934 return this->left + this->right;
3935 }
3936
3937// Check if the left and middle are combinable. This will be constant merged in
3938// multiply reduction.
3939 if (is_constant_combinable(this->left, this->middle) ||
3940 is_variable_combinable(this->left, this->middle)) {
3941 return (this->left*this->middle) + this->right;
3942 }
3943
3944// fma(c2,c1,a) -> fma(c1,c2,a)
3945 if (is_constant_promotable(this->middle,
3946 this->left)) {
3947 return fma(this->middle, this->left, this->right);
3948 }
3949
3950// fma(a,b,a) -> a*(1 + b)
3951// fma(b,a,a) -> a*(1 + b)
3952 if (this->left->is_match(this->right)) {
3953 return this->left*(1.0 + this->middle);
3954 } else if (this->middle->is_match(this->right)) {
3955 return this->middle*(1.0 + this->left);
3956 }
3957
3958// fma(c1,c2 + a,c3) -> fma(c4,a,c5)
3959 auto ma = add_cast(this->middle);
3960 if (ma.get()) {
3961 if (is_constant_combinable(this->left, ma->get_left()) &&
3962 is_constant_combinable(this->left, this->right)) {
3963 return fma(this->left,
3964 ma->get_right(),
3965 fma(this->left, ma->get_left(), this->right));
3966 }
3967 }
3968
3969// fma(c1,c2 - a,c3) -> fma(-c1,a,c1*c2 + c3)
3970// fma(c1,a - c2,c3) -> fma(c1,a,c3 - c1*c2)
3971 auto ms = subtract_cast(this->middle);
3972 if (ms.get()) {
3973 if (is_constant_combinable(this->left, ms->get_left()) &&
3974 is_constant_combinable(this->left, this->right)) {
3975 return fma(-this->left, ms->get_right(),
3976 this->left*ms->get_left() + this->right);
3977 } else if (is_constant_combinable(this->left, ms->get_right()) &&
3978 is_constant_combinable(this->left, this->right)) {
3979 return fma(this->left, ms->get_left(),
3980 this->right - this->left*ms->get_right());
3981 }
3982
3983 auto temp = this->reduce_nested_fma(ms);
3984 if (temp.get() != this) {
3985 return temp;
3986 }
3987 }
3988
3989// Common factor reduction. If the left and right are both multiply nodes check
3990// for a common factor. So you can change a*b + (a*c) -> a*(b + c).
3991 auto lm = multiply_cast(this->left);
3992 auto mm = multiply_cast(this->middle);
3993 auto rm = multiply_cast(this->right);
3994 if (rm.get()) {
3995 if (rm->get_left()->is_match(this->left)) {
3996 return this->left*(this->middle + rm->get_right());
3997 } else if (rm->get_left()->is_match(this->middle)) {
3998 return this->middle*(this->left + rm->get_right());
3999 } else if (rm->get_right()->is_match(this->left)) {
4000 return this->left*(this->middle + rm->get_left());
4001 } else if (rm->get_right()->is_match(this->middle)) {
4002 return this->middle*(this->left + rm->get_left());
4003 }
4004
4005// Change case of
4006// fma(a,b,-c1*b) -> a*b - c1*b
4007 auto rmlc = constant_cast(rm->get_left());
4008 if (rmlc.get() && rmlc->evaluate().is_negative()) {
4009 return this->left*this->middle -
4010 (-1.0*rm->get_left())*rm->get_right();
4011 }
4012
4013// Change cases like
4014// fma(c1,a,c2*b) -> c1*fma(c3,b,a)
4015// fma(a,c1,c2*b) -> c1*fma(c3,b,a)
4016// fma(c1,a,b*c2) -> c1*fma(c3,b,a)
4017// fma(a,c1,b*c2) -> c1*fma(c3,b,a)
4018 if (is_constant_combinable(this->left,
4019 rm->get_left()) &&
4020 !this->left->has_constant_zero()) {
4021 auto temp = rm->get_left()/this->left;
4022 if (temp->is_normal()) {
4023 return this->left*fma(temp,
4024 rm->get_right(),
4025 this->middle);
4026 }
4027 }
4028 if (is_constant_combinable(this->middle,
4029 rm->get_left()) &&
4030 !this->middle->has_constant_zero()) {
4031 auto temp = rm->get_left()/this->middle;
4032 if (temp->is_normal()) {
4033 return this->middle*fma(temp,
4034 rm->get_right(),
4035 this->left);
4036 }
4037 }
4038 if (is_constant_combinable(this->left,
4039 rm->get_right()) &&
4040 !this->left->has_constant_zero()) {
4041 auto temp = rm->get_right()/this->left;
4042 if (temp->is_normal()) {
4043 return this->left*fma(temp,
4044 rm->get_left(),
4045 this->middle);
4046 }
4047 }
4048 if (is_constant_combinable(this->middle,
4049 rm->get_right()) &&
4050 !this->middle->has_constant_zero()) {
4051 auto temp = rm->get_right()/this->middle;
4052 if (temp->is_normal()) {
4053 return this->middle*fma(temp,
4054 rm->get_left(),
4055 this->left);
4056 }
4057 }
4058
4059// fma(a,b*c,b*d) -> b*fma(a,c,d)
4060// fma(a,c*b,b*d) -> b*fma(a,c,d)
4061// fma(a,b*c,d*b) -> b*fma(a,c,d)
4062// fma(a,c*b,d*b) -> b*fma(a,c,d)
4063 if (mm.get()) {
4064 if (mm->get_left()->is_match(rm->get_left())) {
4065 return mm->get_left()*fma(this->left,
4066 mm->get_right(),
4067 rm->get_right());
4068 } else if (mm->get_left()->is_match(rm->get_right())) {
4069 return mm->get_left()*fma(this->left,
4070 mm->get_right(),
4071 rm->get_left());
4072 } else if (mm->get_right()->is_match(rm->get_left())) {
4073 return mm->get_right()*fma(this->left,
4074 mm->get_left(),
4075 rm->get_right());
4076 } else if (mm->get_right()->is_match(rm->get_right())) {
4077 return mm->get_right()*fma(this->left,
4078 mm->get_left(),
4079 rm->get_left());
4080 }
4081 }
4082
4083// Convert fma(a*b,c,d*e) -> fma(d,e,a*b*c)
4084// Convert fma(a,b*c,d*e) -> fma(d,e,a*b*c)
4085 if ((lm.get() || mm.get()) &&
4086 (this->left->get_complexity() + this->middle->get_complexity() >
4087 this->right->get_complexity())) {
4088 return fma(rm->get_left(), rm->get_right(),
4089 this->left*this->middle);
4090 }
4091 }
4092
4093// Handle cases like.
4094// fma(c1*a,b,c2*d) -> c1*(a*b + c2/c1*d)
4095// fma(a*c1,b,c2*d) -> c1*(a*b + c2/c1*d)
4096// fma(c1*a,b,d*c2*d) -> c1*(a*b + c2/c1*d)
4097// fma(a*c1,b,d*c2*d) -> c1*(a*b + c2/c1*d)
4098 if (lm.get() && rm.get()) {
4099 if (is_constant_combinable(rm->get_left(),
4100 lm->get_left()) &&
4101 !lm->get_left()->has_constant_zero()) {
4102 auto temp = rm->get_left()/lm->get_left();
4103 if (temp->is_normal()){
4104 return lm->get_left()*fma(lm->get_right(),
4105 this->middle,
4106 temp*rm->get_right());
4107 }
4108 }
4109 if (is_constant_combinable(rm->get_left(),
4110 lm->get_right()) &&
4111 !lm->get_right()->has_constant_zero()) {
4112 auto temp = rm->get_left()/lm->get_right();
4113 if (temp->is_normal()){
4114 return lm->get_right()*fma(lm->get_left(),
4115 this->middle,
4116 temp*rm->get_right());
4117 }
4118 }
4119 if (is_constant_combinable(rm->get_right(),
4120 lm->get_left()) &&
4121 !lm->get_left()->has_constant_zero()) {
4122 auto temp = rm->get_right()/lm->get_left();
4123 if (temp->is_normal()) {
4124 return lm->get_left()*fma(lm->get_right(),
4125 this->middle,
4126 temp*rm->get_left());
4127 }
4128 }
4129 if (is_constant_combinable(rm->get_right(),
4130 lm->get_right()) &&
4131 !lm->get_right()->has_constant_zero()) {
4132 auto temp = rm->get_right()/lm->get_right();
4133 if (temp->is_normal()) {
4134 return lm->get_right()*fma(lm->get_left(),
4135 this->middle,
4136 temp*rm->get_left());
4137 }
4138 }
4139 }
4140
4141// Move constant multiplies to the left.
4142 if (lm.get()) {
4143// fma(c1*a,b,c) -> fma(c1,a*b,c)
4144 if (is_constant_promotable(lm->get_left(),
4145 lm->get_right())) {
4146 return fma(lm->get_left(),
4147 lm->get_right()*this->middle,
4148 this->right);
4149 }
4150 } else if (mm.get()) {
4151// fma(c1,c2*a,b) -> fma(c3,a,b)
4152// fma(c1,a*c2,b) -> fma(c3,a,b)
4153// fma(a,c1*b,c) -> fma(c1,a*b,c)
4154 if (is_constant_combinable(this->left,
4155 mm->get_left())) {
4156 auto temp = this->left*mm->get_left();
4157 if (temp->is_normal()) {
4158 return fma(temp,
4159 mm->get_right(),
4160 this->right);
4161 }
4162 }
4163 if (is_constant_combinable(this->left,
4164 mm->get_right())) {
4165 auto temp = this->left*mm->get_right();
4166 if (temp->is_normal()) {
4167 return fma(temp,
4168 mm->get_left(),
4169 this->right);
4170 }
4171 }
4172 if (is_constant_promotable(mm->get_left(),
4173 this->left)) {
4174 return fma(mm->get_left(),
4175 this->left*mm->get_right(),
4176 this->right);
4177 }
4178 }
4179
4180// fma(a,b*c,b) -> b*fma(a,c,1)
4181 if (mm.get()) {
4182 if (mm->get_left()->is_match(this->right)) {
4183 return mm->get_left()*fma(this->left,
4184 mm->get_right(),
4185 1.0);
4186 } else if (mm->get_right()->is_match(this->right)) {
4187 return mm->get_right()*fma(this->left,
4188 mm->get_left(),
4189 1.0);
4190 }
4191 }
4192
4193// fma(c1,a,c2/b) -> c1*(a + c3/b)
4194// fma(a,c1,c2/b) -> c1*(a + c3/b)
4195 auto rd = divide_cast(this->right);
4196 if (rd.get()) {
4197 if (is_constant_combinable(this->left,
4198 rd->get_left()) &&
4199 !this->left->has_constant_zero()) {
4200 auto temp = rd->get_left()/this->left;
4201 if (temp->is_normal()) {
4202 return this->left*(this->middle +
4203 temp/rd->get_right());
4204 }
4205 }
4206 if (is_constant_combinable(this->middle,
4207 rd->get_left()) &&
4208 !this->middle->has_constant_zero()) {
4209 auto temp = rd->get_left()/this->middle;
4210 if (temp->is_normal()) {
4211 return this->middle*(this->left +
4212 temp/rd->get_right());
4213 }
4214 }
4215 }
4216
4217// Reduce fma(a/b,b,c) -> a + c
4218// Reduce fma(a,b/a,c) -> b + c
4219 auto ld = divide_cast(this->left);
4220 if (ld.get() && ld->get_right()->is_match(this->middle)) {
4221 return ld->get_left() + this->right;
4222 }
4223 auto md = divide_cast(this->middle);
4224 if (md.get() && md->get_right()->is_match(this->left)) {
4225 return md->get_left() + this->right;
4226 }
4227
4228// Common denominator reductions.
4229 if (ld.get() && rd.get()) {
4230// fma(b/c,a,b,d) -> b(a/c + 1/d)
4231 if (ld->get_left()->is_match(rd->get_left())) {
4232 return ld->get_left()*(this->middle/ld->get_right() +
4233 1.0/rd->get_right());
4234 }
4235
4236// fma(a/(b*c),d,e/c) -> fma(a,d,e*b)/(b*c)
4237// fma(a/(c*b),d,e/c) -> fma(a,d,e*b)/(c*b)
4238// fma(a/c,d,e/(c*b)) -> fma(a*b,d,e)/(b*c)
4239// fma(a/c,d,e/(b*c)) -> fma(a*b,d,e)/(c*b)
4240 auto ldrm = multiply_cast(ld->get_right());
4241 auto rdrm = multiply_cast(rd->get_right());
4242
4243 if (ldrm.get()) {
4244 if (ldrm->get_right()->is_match(rd->get_right())) {
4245 return fma(ld->get_left(), this->middle,
4246 rd->get_left()*ldrm->get_left()) /
4247 ld->get_right();
4248 } else if (ldrm->get_left()->is_match(rd->get_right())) {
4249 return fma(ld->get_left(), this->middle,
4250 rd->get_left()*ldrm->get_right()) /
4251 ld->get_right();
4252 }
4253 } else if (rdrm.get()) {
4254 if (rdrm->get_right()->is_match(ld->get_right())) {
4255 return fma(ld->get_left()*rdrm->get_left(),
4256 this->middle, rd->get_left()) /
4257 rd->get_right();
4258 } else if (rdrm->get_left()->is_match(ld->get_right())) {
4259 return fma(ld->get_left()*rdrm->get_right(),
4260 this->middle, rd->get_left()) /
4261 rd->get_right();
4262 }
4263 }
4264 } else if (md.get() && rd.get()) {
4265// fma(a,d/(b*c),e/c) -> fma(a,d,e*b)/(b*c)
4266// fma(a,d/(c*b),e/c) -> fma(a,d,e*b)/(c*b)
4267// fma(a,d/c,e/(c*b)) -> fma(a,d*b,e)/(b*c)
4268// fma(a,d/c,e/(b*c)) -> fma(a,d*b,e)/(c*b)
4269 auto mdrm = multiply_cast(md->get_right());
4270 auto rdrm = multiply_cast(rd->get_right());
4271
4272 if (mdrm.get()) {
4273 if (mdrm->get_right()->is_match(rd->get_right())) {
4274 return fma(this->left, md->get_left(),
4275 rd->get_left()*mdrm->get_left()) /
4276 md->get_right();
4277 } else if (mdrm->get_left()->is_match(rd->get_right())) {
4278 return fma(this->left, md->get_left(),
4279 rd->get_left()*mdrm->get_right()) /
4280 md->get_right();
4281 }
4282 } else if (rdrm.get()) {
4283 if (rdrm->get_right()->is_match(md->get_right())) {
4284 return fma(this->left, md->get_left()*rdrm->get_left(),
4285 rd->get_left()) /
4286 rd->get_right();
4287 } else if (rdrm->get_left()->is_match(md->get_right())) {
4288 return fma(this->left, md->get_left()*rdrm->get_right(),
4289 rd->get_left()) /
4290 rd->get_right();
4291 }
4292 }
4293 }
4294
4295// Chained fma reductions.
4296 auto rfma = fma_cast(this->right);
4297 if (rfma.get()) {
4298// fma(a, b, fma(c, b, d)) -> fma(b, a + c, d)
4299// fma(b, a, fma(c, b, d)) -> fma(b, a + c, d)
4300// fma(a, b, fma(b, c, d)) -> fma(b, a + c, d)
4301// fma(b, a, fma(b, c, d)) -> fma(b, a + c, d)
4302 if (this->middle->is_match(rfma->get_middle())) {
4303 return fma(this->middle,
4304 this->left + rfma->get_left(),
4305 rfma->get_right());
4306 } else if (this->left->is_match(rfma->get_middle())) {
4307 return fma(this->left,
4308 this->middle + rfma->get_left(),
4309 rfma->get_right());
4310 } else if (this->middle->is_match(rfma->get_left())) {
4311 return fma(this->middle,
4312 this->left + rfma->get_middle(),
4313 rfma->get_right());
4314 } else if (this->left->is_match(rfma->get_left())) {
4315 return fma(this->left,
4316 this->middle + rfma->get_middle(),
4317 rfma->get_right());
4318 }
4319
4320 if (mm.get()) {
4321// fma(a, e*b, fma(c, b, d)) -> fma(b, fma(a, e, c), d)
4322// fma(a, b*e, fma(c, b, d)) -> fma(b, fma(a, e, c), d)
4323// fma(a, e*b, fma(b, c, d)) -> fma(b, fma(a, e, c), d)
4324// fma(a, b*e, fma(b, c, d)) -> fma(b, fma(a, e, c), d)
4325 if (mm->get_right()->is_match(rfma->get_middle())) {
4326 return fma(mm->get_right(),
4327 fma(this->left,
4328 mm->get_left(),
4329 rfma->get_left()),
4330 rfma->get_right());
4331 } else if (mm->get_left()->is_match(rfma->get_middle())) {
4332 return fma(mm->get_left(),
4333 fma(this->left,
4334 mm->get_right(),
4335 rfma->get_left()),
4336 rfma->get_right());
4337 } else if (mm->get_right()->is_match(rfma->get_left())) {
4338 return fma(mm->get_right(),
4339 fma(this->left,
4340 mm->get_left(),
4341 rfma->get_middle()),
4342 rfma->get_right());
4343 } else if (mm->get_left()->is_match(rfma->get_left())) {
4344 return fma(mm->get_left(),
4345 fma(this->left,
4346 mm->get_right(),
4347 rfma->get_middle()),
4348 rfma->get_right());
4349 }
4350 } else if (lm.get()) {
4351// fma(e*b, a, fma(c, b, d)) -> fma(b, fma(a, e, c), d)
4352// fma(b*e, a, fma(c, b, d)) -> fma(b, fma(a, e, c), d)
4353// fma(e*b, a, fma(b, c, d)) -> fma(b, fma(a, e, c), d)
4354// fma(e*d, a, fma(b, c, d)) -> fma(b, fma(a, e, c), d)
4355 if (lm->get_right()->is_match(rfma->get_middle())) {
4356 return fma(lm->get_right(),
4357 fma(this->middle,
4358 lm->get_left(),
4359 rfma->get_left()),
4360 rfma->get_right());
4361 } else if (lm->get_left()->is_match(rfma->get_middle())) {
4362 return fma(lm->get_left(),
4363 fma(this->middle,
4364 lm->get_right(),
4365 rfma->get_left()),
4366 rfma->get_right());
4367 } else if (lm->get_right()->is_match(rfma->get_left())) {
4368 return fma(lm->get_right(),
4369 fma(this->middle,
4370 lm->get_left(),
4371 rfma->get_middle()),
4372 rfma->get_right());
4373 } else if (lm->get_left()->is_match(rfma->get_left())) {
4374 return fma(lm->get_left(),
4375 fma(this->middle,
4376 lm->get_right(),
4377 rfma->get_middle()),
4378 rfma->get_right());
4379 }
4380 }
4381
4382 auto rfmamm = multiply_cast(rfma->get_middle());
4383 auto rfmalm = multiply_cast(rfma->get_left());
4384 if (rfmamm.get()) {
4385// fma(a, b, fma(c, e*b, d)) -> fma(b, fma(c, e, a), d)
4386// fma(b, a, fma(c, e*b, d)) -> fma(b, fma(c, e, a), d)
4387// fma(a, b, fma(c, b*e, d)) -> fma(b, fma(c, e, a), d)
4388// fma(b, a, fma(c, b*e, d)) -> fma(b, fma(c, e, a), d)
4389 if (rfmamm->get_right()->is_match(this->middle)) {
4390 return fma(this->middle,
4391 fma(rfma->get_left(),
4392 rfmamm->get_left(),
4393 this->left),
4394 rfma->get_right());
4395 } else if (rfmamm->get_right()->is_match(this->left)) {
4396 return fma(this->left,
4397 fma(rfma->get_left(),
4398 rfmamm->get_left(),
4399 this->middle),
4400 rfma->get_right());
4401 } else if (rfmamm->get_left()->is_match(this->middle)) {
4402 return fma(this->middle,
4403 fma(rfma->get_left(),
4404 rfmamm->get_right(),
4405 this->left),
4406 rfma->get_right());
4407 } else if (rfmamm->get_left()->is_match(this->left)) {
4408 return fma(this->left,
4409 fma(rfma->get_left(),
4410 rfmamm->get_right(),
4411 this->middle),
4412 rfma->get_right());
4413 }
4414 } else if (rfmalm.get()) {
4415// fma(a, b, fma(e*b, c, d)) -> fma(b, fma(c, e, a), d)
4416// fma(b, a, fma(e*b, c, d)) -> fma(b, fma(c, e, a), d)
4417// fma(a, b, fma(b*e, c, d)) -> fma(b, fma(c, e, a), d)
4418// fma(b, a, fma(b*e, c, d)) -> fma(b, fma(c, e, a), d)
4419 if (rfmalm->get_right()->is_match(this->middle)) {
4420 return fma(this->middle,
4421 fma(rfma->get_middle(),
4422 rfmalm->get_left(),
4423 this->left),
4424 rfma->get_right());
4425 } else if (rfmalm->get_right()->is_match(this->left)) {
4426 return fma(this->left,
4427 fma(rfma->get_middle(),
4428 rfmalm->get_left(),
4429 this->middle),
4430 rfma->get_right());
4431 } else if (rfmalm->get_left()->is_match(this->middle)) {
4432 return fma(this->middle,
4433 fma(rfma->get_middle(),
4434 rfmalm->get_right(),
4435 this->left),
4436 rfma->get_right());
4437 } else if (rfmalm->get_left()->is_match(this->left)) {
4438 return fma(this->left,
4439 fma(rfma->get_middle(),
4440 rfmalm->get_right(),
4441 this->middle),
4442 rfma->get_right());
4443 }
4444 }
4445
4446 if (mm.get() && rfmamm.get()) {
4447// fma(a, f*b, fma(c, e*b, d)) -> fma(b, fma(a, f, c*e), d)
4448// fma(a, b*f, fma(c, e*b, d)) -> fma(b, fma(a, f, c*e), d)
4449// fma(a, f*b, fma(c, b*e, d)) -> fma(b, fma(a, f, c*e), d)
4450// fma(a, b*f, fma(c, b*e, d)) -> fma(b, fma(a, f, c*e), d)
4451 if (mm->get_right()->is_match(rfmamm->get_right())) {
4452 return fma(mm->get_right(),
4453 fma(this->left,
4454 mm->get_left(),
4455 rfma->get_left()*rfmamm->get_left()),
4456 rfma->get_right());
4457 } else if (mm->get_left()->is_match(rfmamm->get_right())) {
4458 return fma(mm->get_left(),
4459 fma(this->left,
4460 mm->get_right(),
4461 rfma->get_left()*rfmamm->get_left()),
4462 rfma->get_right());
4463 } else if (mm->get_right()->is_match(rfmamm->get_left())) {
4464 return fma(mm->get_right(),
4465 fma(this->left,
4466 mm->get_left(),
4467 rfma->get_left()*rfmamm->get_right()),
4468 rfma->get_right());
4469 } else if (mm->get_left()->is_match(rfmamm->get_left())) {
4470 return fma(mm->get_left(),
4471 fma(this->left,
4472 mm->get_right(),
4473 rfma->get_left()*rfmamm->get_right()),
4474 rfma->get_right());
4475 }
4476 } else if (lm.get() && rfmamm.get()) {
4477// fma(f*b, a, fma(c, e*b, d)) -> fma(b, fma(a, f, c*e), d)
4478// fma(b*f, a, fma(c, e*b, d)) -> fma(b, fma(a, f, c*e), d)
4479// fma(f*b, a, fma(c, b*e, d)) -> fma(b, fma(a, f, c*e), d)
4480// fma(b*f, a, fma(c, b*e, d)) -> fma(b, fma(a, f, c*e), d)
4481 if (lm->get_right()->is_match(rfmamm->get_right())) {
4482 return fma(lm->get_right(),
4483 fma(this->middle,
4484 lm->get_left(),
4485 rfma->get_left()*rfmamm->get_left()),
4486 rfma->get_right());
4487 } else if (lm->get_left()->is_match(rfmamm->get_right())) {
4488 return fma(lm->get_left(),
4489 fma(this->middle,
4490 lm->get_right(),
4491 rfma->get_left()*rfmamm->get_left()),
4492 rfma->get_right());
4493 } else if (lm->get_right()->is_match(rfmamm->get_left())) {
4494 return fma(lm->get_right(),
4495 fma(this->middle,
4496 lm->get_left(),
4497 rfma->get_left()*rfmamm->get_right()),
4498 rfma->get_right());
4499 } else if (lm->get_left()->is_match(rfmamm->get_left())) {
4500 return fma(lm->get_left(),
4501 fma(this->middle,
4502 lm->get_right(),
4503 rfma->get_left()*rfmamm->get_right()),
4504 rfma->get_right());
4505 }
4506 } else if (mm.get() && rfmalm.get()) {
4507// fma(a, f*b, fma(e*b, c, d)) -> fma(b, fma(a, f, c*e), d)
4508// fma(a, b*f, fma(e*b, c, d)) -> fma(b, fma(a, f, c*e), d)
4509// fma(a, f*b, fma(b*e, c, d)) -> fma(b, fma(a, f, c*e), d)
4510// fma(a, b*f, fma(b*e, c, d)) -> fma(b, fma(a, f, c*e), d)
4511 if (mm->get_right()->is_match(rfmalm->get_right())) {
4512 return fma(mm->get_right(),
4513 fma(this->left,
4514 mm->get_left(),
4515 rfma->get_middle()*rfmalm->get_left()),
4516 rfma->get_right());
4517 } else if (mm->get_left()->is_match(rfmalm->get_right())) {
4518 return fma(mm->get_left(),
4519 fma(this->left,
4520 mm->get_right(),
4521 rfma->get_middle()*rfmalm->get_left()),
4522 rfma->get_right());
4523 } else if (mm->get_right()->is_match(rfmalm->get_left())) {
4524 return fma(mm->get_right(),
4525 fma(this->left,
4526 mm->get_left(),
4527 rfma->get_middle()*rfmalm->get_right()),
4528 rfma->get_right());
4529 } else if (mm->get_left()->is_match(rfmalm->get_left())) {
4530 return fma(mm->get_left(),
4531 fma(this->left,
4532 mm->get_right(),
4533 rfma->get_middle()*rfmalm->get_right()),
4534 rfma->get_right());
4535 }
4536 } else if (lm.get() && rfmalm.get()) {
4537// fma(f*b, a, fma(e*b, c, d)) -> fma(b, fma(a, f, c*e), d)
4538// fma(b*f, a, fma(e*b, c, d)) -> fma(b, fma(a, f, c*e), d)
4539// fma(f*b, a, fma(b*e, c, d)) -> fma(b, fma(a, f, c*e), d)
4540// fma(b*f, a, fma(b*e, c, d)) -> fma(b, fma(a, f, c*e), d)
4541 if (lm->get_right()->is_match(rfmalm->get_right())) {
4542 return fma(lm->get_right(),
4543 fma(this->middle,
4544 lm->get_left(),
4545 rfma->get_middle()*rfmalm->get_left()),
4546 rfma->get_right());
4547 } else if (lm->get_left()->is_match(rfmalm->get_right())) {
4548 return fma(lm->get_left(),
4549 fma(this->middle,
4550 lm->get_right(),
4551 rfma->get_middle()*rfmalm->get_left()),
4552 rfma->get_right());
4553 } else if (lm->get_right()->is_match(rfmalm->get_left())) {
4554 return fma(lm->get_right(),
4555 fma(this->middle,
4556 lm->get_left(),
4557 rfma->get_middle()*rfmalm->get_right()),
4558 rfma->get_right());
4559 } else if (lm->get_left()->is_match(rfmalm->get_left())) {
4560 return fma(lm->get_left(),
4561 fma(this->middle,
4562 lm->get_right(),
4563 rfma->get_middle()*rfmalm->get_right()),
4564 rfma->get_right());
4565 }
4566 }
4567
4568 if (is_variable_combinable(this->middle, rfma->get_middle())) {
4569 if (is_greater_exponent(this->middle, rfma->get_middle())) {
4570// fma(a,x^b,fma(c,x^d,e)) -> fma(x^d,fma(x^(d-b),a,c),e) if b > d
4571 return fma(rfma->get_middle(),
4572 fma(this->middle/rfma->get_middle(),
4573 this->left,
4574 rfma->get_left()),
4575 rfma->get_right());
4576 } else {
4577// fma(a,x^b,fma(c,x^d,e)) -> fma(x^b,fma(x^(d-b),c,a),e) if d > b
4578 return fma(this->middle,
4579 fma(rfma->get_middle()/this->middle,
4580 rfma->get_left(),
4581 this->left),
4582 rfma->get_right());
4583 }
4584 } else if (is_variable_combinable(this->left, rfma->get_middle())) {
4585 if (is_greater_exponent(this->left, rfma->get_middle())) {
4586// fma(x^b,a,fma(c,x^d,e)) -> fma(x^d,fma(x^(d-b),a,c),e) if b > d
4587 return fma(rfma->get_middle(),
4588 fma(this->left/rfma->get_middle(),
4589 this->middle,
4590 rfma->get_left()),
4591 rfma->get_right());
4592 } else {
4593// fma(x^b,a,fma(c,x^d,e)) -> fma(x^b,fma(x^(d-b),c,a),e) if d > b
4594 return fma(this->left,
4595 fma(rfma->get_middle()/this->left,
4596 rfma->get_left(),
4597 this->middle),
4598 rfma->get_right());
4599 }
4600 } else if (is_variable_combinable(this->middle, rfma->get_left())) {
4601 if (is_greater_exponent(this->middle, rfma->get_left())) {
4602// fma(a,x^b,fma(x^d,c,e)) -> fma(x^d,fma(x^(d-b),a,c),e) if b > d
4603 return fma(rfma->get_left(),
4604 fma(this->middle/rfma->get_left(),
4605 this->left,
4606 rfma->get_middle()),
4607 rfma->get_right());
4608 } else {
4609// fma(a,x^b,fma(x^d,c,e)) -> fma(x^b,fma(x^(d-b),c,a),e) if d > b
4610 return fma(this->middle,
4611 fma(rfma->get_left()/this->middle,
4612 rfma->get_middle(),
4613 this->left),
4614 rfma->get_right());
4615 }
4616 } else if (is_variable_combinable(this->left, rfma->get_left())) {
4617 if (is_greater_exponent(this->left, rfma->get_left())) {
4618// fma(x^b,a,fma(x^d,c,e)) -> fma(x^d,fma(x^(d-b),a,c),e) if b > d
4619 return fma(rfma->get_left(),
4620 fma(this->left/rfma->get_left(),
4621 this->middle,
4622 rfma->get_middle()),
4623 rfma->get_right());
4624 } else {
4625// fma(x^b,a,fma(x^d,c,e)) -> fma(x^b,fma(x^(d-b),c,a),e) if d > b
4626 return fma(this->left,
4627 fma(rfma->get_left()/this->left,
4628 rfma->get_middle(),
4629 this->middle),
4630 rfma->get_right());
4631 }
4632 }
4633
4634// fma(a,b,fma(a,b,c)) -> fma(2*a,b,c)
4635// fma(a,b,fma(b,a,c)) -> fma(2*a,b,c)
4636 if (this->left->is_match(rfma->get_left()) &&
4637 this->middle->is_match(rfma->get_middle())) {
4638 return fma(2.0*this->left, this->middle, rfma->get_right());
4639 } else if (this->left->is_match(rfma->get_middle()) &&
4640 this->middle->is_match(rfma->get_left())) {
4641 return fma(2.0*this->left, this->middle, rfma->get_right());
4642 }
4643
4644// fma(a,b/c,fma(e,f/c,g)) -> (a*b + e*f)/c + g
4645// fma(a,b/c,fma(e/c,f,g)) -> (a*b + e*f)/c + g
4646// fma(a/c,b,fma(e,f/c,g)) -> (a*b + e*f)/c + g
4647// fma(a/c,b,fma(e/c,f,g)) -> (a*b + e*f)/c + g
4648 auto fmald = divide_cast(rfma->get_left());
4649 auto fmamd = divide_cast(rfma->get_middle());
4650 if (ld.get()) {
4651 if (fmald.get() && ld->get_right()->is_match(fmald->get_right())) {
4652 return (ld->get_left()*this->middle +
4653 fmald->get_left()*rfma->get_middle())/ld->get_right() +
4654 rfma->get_right();
4655 } else if (fmamd.get() && ld->get_right()->is_match(fmamd->get_right())) {
4656 return (ld->get_left()*this->middle +
4657 fmamd->get_left()*rfma->get_left())/ld->get_right() +
4658 rfma->get_right();
4659 }
4660 } else if (md.get()) {
4661 if (fmald.get() && md->get_right()->is_match(fmald->get_right())) {
4662 return (md->get_left()*this->left +
4663 fmald->get_left()*rfma->get_middle())/md->get_right() +
4664 rfma->get_right();
4665 } else if (fmamd.get() && md->get_right()->is_match(fmamd->get_right())) {
4666 return (md->get_left()*this->left +
4667 fmamd->get_left()*rfma->get_left())/md->get_right() +
4668 rfma->get_right();
4669 }
4670 }
4671 }
4672
4673// Check to see if it is worth moving nodes out of a fma nodes. These should be
4674// restricted to variable like nodes. Only do this reduction if the complexity
4675// reduces.
4676 if (this->left->is_all_variables()) {
4677 auto rdl = this->right/this->left;
4678 if (rdl->get_complexity() < this->left->get_complexity() +
4679 this->right->get_complexity()) {
4680 return (this->middle + rdl)*this->left;
4681 }
4682 } else if (this->middle->is_all_variables()) {
4683 auto rdm = this->right/this->middle;
4684 auto rdmc = constant_cast(rdm->get_power_exponent());
4685 if ((rdm->get_complexity() < this->middle->get_complexity() +
4686 this->right->get_complexity()) &&
4687 !(rdmc.get() && rdmc->evaluate().is_negative())) {
4688 return (this->left + rdm)*this->middle;
4689 }
4690 }
4691
4692// Change negative exponents to divide so that can be factored out.
4693// fma(a,b^-c,d) = a/b^c + d
4694// fma(b^-c,a,d) = a/b^c + d
4695 auto lp = pow_cast(this->left);
4696 if (lp.get()) {
4697 auto exponent = constant_cast(lp->get_right());
4698 if (exponent.get() && exponent->evaluate().is_negative()) {
4699 return this->middle/pow(lp->get_left(), -lp->get_right()) +
4700 this->right;
4701 }
4702 }
4703 auto mp = pow_cast(this->middle);
4704 if (mp.get()) {
4705 auto exponent = constant_cast(mp->get_right());
4706 if (exponent.get() && exponent->evaluate().is_negative()) {
4707 return this->left/pow(mp->get_left(), -mp->get_right()) +
4708 this->right;
4709 }
4710
4711// fma(2,a^2,a) -> a*fma(2,a,1)
4712// Note this case is handled earlier. fma(2,a,a^2) -> a*fma(2,1,a)
4713 if (is_variable_combinable(this->middle,
4714 this->right)) {
4715 auto temp = this->right/this->middle;
4716 auto temp_exponent = constant_cast(temp->get_power_exponent());
4717 if (temp_exponent.get() && temp_exponent->evaluate().is_negative()) {
4718 return this->right*fma(this->left,
4719 this->middle/this->right,
4720 1.0);
4721 }
4722 }
4723 }
4724
4725// a^b*c^b + d -> (a*c)^b + d
4726 if (lp.get() && mp.get()) {
4727 if (lp->get_right()->is_match(mp->get_right())) {
4728 return pow(lp->get_left()*mp->get_left(),
4729 lp->get_right()) +
4730 this->right;
4731 }
4732 }
4733
4734// fma(2,(ab)^2,a^2b) -> a^2*fma(2, b^2, b)
4735 if (rm.get() && mp.get()) {
4736 auto mplm = multiply_cast(mp->get_left());
4737 if (mplm.get()) {
4738 if (is_variable_combinable(mplm->get_left(),
4739 rm->get_left())) {
4740 auto temp = pow(mplm->get_left(),
4741 mp->get_right());
4742 return temp*fma(this->left,
4743 this->middle/temp,
4744 this->right/temp);
4745 } else if (is_variable_combinable(mplm->get_right(),
4746 rm->get_left())) {
4747 auto temp = pow(mplm->get_right(),
4748 mp->get_right());
4749 return temp*fma(this->left,
4750 this->middle/temp,
4751 this->right/temp);
4752 }
4753 }
4754 }
4755// fma(2,(a*b)^2,fma(a^2,b,c)) -> fma(a^2,fma(2,b^2,b),c)
4756 if (rfma.get() && mp.get()) {
4757 auto mplm = multiply_cast(mp->get_left());
4758 if (mplm.get()) {
4759 if (is_variable_combinable(mplm->get_left(),
4760 rfma->get_left())) {
4761 auto temp = pow(mplm->get_left(),
4762 mp->get_right());
4763 return fma(temp,
4764 fma(this->left,
4765 this->middle/temp,
4766 rfma->get_middle()),
4767 rfma->get_right());
4768 } else if (is_variable_combinable(mplm->get_right(),
4769 rfma->get_left())) {
4770 auto temp = pow(mplm->get_right(),
4771 mp->get_right());
4772 return fma(temp,
4773 fma(this->left,
4774 this->middle/temp,
4775 rfma->get_middle()),
4776 rfma->get_right());
4777 }
4778
4779// fma(2,(a*b)^2,fma(3,a^2*b,c)) -> a^2*fma(2,b^2,fma(3,b,c))
4780 auto rfmamm = multiply_cast(rfma->get_middle());
4781 if (rfmamm.get()) {
4782 if (is_variable_combinable(mplm->get_left(),
4783 rfmamm->get_left())) {
4784 auto temp = pow(mplm->get_left(),
4785 mp->get_right());
4786 return temp*fma(this->left,
4787 this->middle/temp,
4788 fma(rfma->get_left(),
4789 rfma->get_middle()/temp,
4790 rfma->get_right()));
4791 }
4792 }
4793 }
4794 }
4795
4796// fma(a,b/c,b/d) -> b*(a/c + 1/d)
4797// fma(a,c/b,d/b) -> (a*c + d)/b
4798 if (md.get() && rd.get()) {
4799 if (md->get_left()->is_match(rd->get_left())) {
4800 return md->get_left()*(this->left/md->get_right() +
4801 1.0/rd->get_right());
4802 } else if (md->get_right()->is_match(rd->get_right())) {
4803 return (this->left*md->get_left() +
4804 rd->get_left())/md->get_right();
4805 }
4806 }
4807// fma(b/c,a,b/d) -> b*(a/c + 1/d)
4808// fma(c/b,a,d/b) -> (a*c + d)/b
4809 if (ld.get() && rd.get()) {
4810 if (ld->get_left()->is_match(rd->get_left())) {
4811 return ld->get_left()*(this->middle/ld->get_right() +
4812 1.0/rd->get_right());
4813 } else if (ld->get_right()->is_match(rd->get_right())) {
4814 return (this->middle*ld->get_left() +
4815 rd->get_left())/ld->get_right();
4816 }
4817 }
4818
4819// fma(a/b,c,(d/b)*e) -> fma(a,c,d*e)/b
4820// fma(a/b,c,e*(d/b)) -> fma(a,c,d*e)/b
4821 if (rm.get() && ld.get()) {
4822 auto rmld = divide_cast(rm->get_left());
4823 if (rmld.get() && ld->get_right()->is_match(rmld->get_right())) {
4824 return fma(ld->get_left(), this->middle, rmld->get_left()*rm->get_right())/ld->get_right();
4825 }
4826 auto rmrd = divide_cast(rm->get_right());
4827 if (rmrd.get() && ld->get_right()->is_match(rmrd->get_right())) {
4828 return fma(ld->get_left(), this->middle, rmrd->get_left()*rm->get_left())/ld->get_right();
4829 }
4830 }
4831// fma(a,c/b,(d/b)*e) -> fma(a,c,d*e)/b
4832// fma(a,c/b,e*(d/b)) -> fma(a,c,d*e)/b
4833 if (rm.get() && md.get()) {
4834 auto rmld = divide_cast(rm->get_left());
4835 if (rmld.get() && md->get_right()->is_match(rmld->get_right())) {
4836 return fma(this->left, md->get_left(), rmld->get_left()*rm->get_right())/md->get_right();
4837 }
4838 auto rmrd = divide_cast(rm->get_right());
4839 if (rmrd.get() && md->get_right()->is_match(rmrd->get_right())) {
4840 return fma(this->left, md->get_left(), rmrd->get_left()*rm->get_left())/md->get_right();
4841 }
4842 }
4843
4844// fma(a/b*c,d,e/b) -> fma(a*c,d,e)/b
4845// fma(a*c/b,d,e/b) -> fma(a*c,d,e)/b
4846 if (rd.get() && lm.get()) {
4847 auto lmld = divide_cast(lm->get_left());
4848 if (lmld.get() && rd->get_right()->is_match(lmld->get_right())) {
4849 return fma(lmld->get_left()*lm->get_right(), this->middle, rd->get_left())/rd->get_right();
4850 }
4851 auto lmrd = divide_cast(lm->get_right());
4852 if (lmrd.get() && rd->get_right()->is_match(lmrd->get_right())) {
4853 return fma(lmld->get_left()*lm->get_left(), this->middle, rd->get_left())/rd->get_right();
4854 }
4855 }
4856// fma(a,c/b*d,e/b) -> fma(a,c*d,e)/b
4857// fma(a,c*d/b,e/b) -> fma(a,c*d,e)/b
4858 if (rd.get() && mm.get()) {
4859 auto mmld = divide_cast(mm->get_left());
4860 if (mmld.get() && rd->get_right()->is_match(mmld->get_right())) {
4861 return fma(this->left, mmld->get_left()*mm->get_right(), rd->get_left())/rd->get_right();
4862 }
4863 auto mmrd = divide_cast(mm->get_right());
4864 if (mmrd.get() && rd->get_right()->is_match(mmrd->get_right())) {
4865 return fma(this->left, mmrd->get_left()*mm->get_left(), rd->get_left())/rd->get_right();
4866 }
4867 }
4868
4869// fma(a, b/c, ((f/c)*e)*d) -> fma(a, b, f*e*d)/c
4870// fma(a/c, b, ((f/c)*e)*d) -> fma(a, b, f*e*d)/c
4871// fma(a, b/c, (e*(f/c))*d) -> fma(a, b, f*e*d)/c
4872// fma(a/c, b, (e*(f/c))*d) -> fma(a, b, f*e*d)/c
4873// fma(a, b/c, d*((f/c)*e)) -> fma(a, b, f*e*d)/c
4874// fma(a/c, b, d*((f/c)*e)) -> fma(a, b, f*e*d)/c
4875// fma(a, b/c, d*(e*(f/c))) -> fma(a, b, f*e*d)/c
4876// fma(a/c, b, d*(e*(f/c))) -> fma(a, b, f*e*d)/c
4877 if (md.get() && rm.get()) {
4878 auto rmlm = multiply_cast(rm->get_left());
4879 if (rmlm.get()) {
4880 auto rmlmld = divide_cast(rmlm->get_left());
4881 if (rmlmld.get() && rmlmld->get_right()->is_match(md->get_right())) {
4882 return fma(this->left, md->get_left(),
4883 rmlmld->get_left()*rmlm->get_right()*rm->get_right())/md->get_right();
4884 }
4885 auto rmlmrd = divide_cast(rmlm->get_right());
4886 if (rmlmrd.get() && rmlmrd->get_right()->is_match(md->get_right())) {
4887 return fma(this->left, md->get_left(),
4888 rmlmrd->get_left()*rmlm->get_left()*rm->get_right())/md->get_right();
4889 }
4890 }
4891 auto rmrm = multiply_cast(rm->get_right());
4892 if (rmrm.get()) {
4893 auto rmrmld = divide_cast(rmrm->get_left());
4894 if (rmrmld.get() && rmrmld->get_right()->is_match(md->get_right())) {
4895 return fma(this->left, md->get_left(),
4896 rmrmld->get_left()*rmrm->get_right()*rm->get_left())/md->get_right();
4897 }
4898 auto rmrmrd = divide_cast(rmrm->get_right());
4899 if (rmrmrd.get() && rmrmrd->get_right()->is_match(md->get_right())) {
4900 return fma(this->left, md->get_left(),
4901 rmrmrd->get_left()*rmrm->get_left()*rm->get_left())/md->get_right();
4902 }
4903 }
4904 } else if (ld.get() && rm.get()) {
4905 auto rmlm = multiply_cast(rm->get_left());
4906 if (rmlm.get()) {
4907 auto rmlmld = divide_cast(rmlm->get_left());
4908 if (rmlmld.get() && rmlmld->get_right()->is_match(ld->get_right())) {
4909 return fma(ld->get_left(), this->middle,
4910 rmlmld->get_left()*rmlm->get_right()*rm->get_right())/ld->get_right();
4911 }
4912 auto rmlmrd = divide_cast(rmlm->get_right());
4913 if (rmlmrd.get() && rmlmrd->get_right()->is_match(ld->get_right())) {
4914 return fma(ld->get_left(), this->middle,
4915 rmlmrd->get_left()*rmlm->get_right()*rm->get_right())/ld->get_right();
4916 }
4917 }
4918 auto rmrm = multiply_cast(rm->get_right());
4919 if (rmrm.get()) {
4920 auto rmrmld = divide_cast(rmrm->get_left());
4921 if (rmrmld.get() && rmrmld->get_right()->is_match(ld->get_right())) {
4922 return fma(ld->get_left(), this->middle,
4923 rmrmld->get_left()*rmrm->get_right()*rm->get_left())/ld->get_right();
4924 }
4925 auto rmrmrd = divide_cast(rmrm->get_right());
4926 if (rmrmrd.get() && rmrmrd->get_right()->is_match(ld->get_right())) {
4927 return fma(ld->get_left(), this->middle,
4928 rmrmrd->get_left()*rmrm->get_left()*rm->get_left())/ld->get_right();
4929 }
4930 }
4931 }
4932
4933// fma(exp(a), exp(b), c) -> exp(a + b) + c
4934 auto le = exp_cast(this->left);
4935 auto me = exp_cast(this->middle);
4936 if (le.get() && me.get()) {
4937 return exp(le->get_arg() + me->get_arg()) + this->right;
4938 }
4939
4940// fma(exp(a), exp(b)*c, d) -> fma(exp(a)*exp(b), c, d)
4941// fma(exp(a), c*exp(b), d) -> fma(exp(a)*exp(b), c, d)
4942 if (mm.get() && le.get()) {
4943 auto mmle = exp_cast(mm->get_left());
4944 if (mmle.get()) {
4945 return fma(this->left*mm->get_left(),
4946 mm->get_right(),
4947 this->right);
4948 }
4949 auto mmre = exp_cast(mm->get_right());
4950 if (mmre.get()) {
4951 return fma(this->left*mm->get_right(),
4952 mm->get_left(),
4953 this->right);
4954 }
4955 }
4956// fma(exp(a)*c, exp(b), d) -> fma(exp(a)*exp(b), c, d)
4957// fma(c*exp(a), exp(b), d) -> fma(exp(a)*exp(b), c, d)
4958 if (lm.get() && me.get()) {
4959 auto lmle = exp_cast(lm->get_left());
4960 if (lmle.get()) {
4961 return fma(lm->get_left()*this->middle,
4962 lm->get_right(),
4963 this->right);
4964 }
4965 auto lmre = exp_cast(lm->get_right());
4966 if (lmre.get()) {
4967 return fma(lm->get_right()*this->middle,
4968 lm->get_left(),
4969 this->right);
4970 }
4971 }
4972
4973// fma(exp(a)*c, exp(b)*d, e) -> fma(exp(a)*exp(b), c*d, e)
4974// fma(exp(a)*c, d*exp(b), e) -> fma(exp(a)*exp(b), c*d, e)
4975// fma(c*exp(a), exp(b)*d, e) -> fma(exp(a)*exp(b), c*d, e)
4976// fma(c*exp(a), d*exp(b), e) -> fma(exp(a)*exp(b), c*d, e)
4977 if (lm.get() && mm.get()) {
4978 auto lmle = exp_cast(lm->get_left());
4979 if (lmle.get()) {
4980 auto mmle = exp_cast(mm->get_left());
4981 if (mmle.get()) {
4982 return fma(lm->get_left()*mm->get_left(),
4983 lm->get_right()*mm->get_right(),
4984 this->right);
4985 }
4986 auto mmre = exp_cast(mm->get_right());
4987 if (mmre.get()) {
4988 return fma(lm->get_left()*mm->get_right(),
4989 lm->get_right()*mm->get_left(),
4990 this->right);
4991 }
4992 }
4993 auto lmre = exp_cast(lm->get_right());
4994 if (lmre.get()) {
4995 auto mmle = exp_cast(mm->get_left());
4996 if (mmle.get()) {
4997 return fma(lm->get_right()*mm->get_left(),
4998 lm->get_left()*mm->get_right(),
4999 this->right);
5000 }
5001 auto mmre = exp_cast(mm->get_right());
5002 if (mmre.get()) {
5003 return fma(lm->get_right()*mm->get_right(),
5004 lm->get_left()*mm->get_left(),
5005 this->right);
5006 }
5007 }
5008 }
5009
5010// fma(exp(a)*c, exp(b)/d, e) -> fma(exp(a)*exp(b), c/d, e)
5011// fma(exp(a)*c, d/exp(b), e) -> fma(exp(a)/exp(b), c*d, e)
5012// fma(c*exp(a), exp(b)/d, e) -> fma(exp(a)*exp(b), c/d, e)
5013// fma(c*exp(a), d/exp(b), e) -> fma(exp(a)/exp(b), c*d, e)
5014 if (lm.get() && md.get()) {
5015 auto lmle = exp_cast(lm->get_left());
5016 if (lmle.get()) {
5017 auto mdle = exp_cast(md->get_left());
5018 if (mdle.get()) {
5019 return fma(lm->get_left()*md->get_left(),
5020 lm->get_right()/md->get_right(),
5021 this->right);
5022 }
5023 auto mdre = exp_cast(md->get_right());
5024 if (mdre.get()) {
5025 return fma(lm->get_left()/md->get_right(),
5026 lm->get_right()*md->get_left(),
5027 this->right);
5028 }
5029 }
5030 auto lmre = exp_cast(lm->get_right());
5031 if (lmre.get()) {
5032 auto mdle = exp_cast(md->get_left());
5033 if (mdle.get()) {
5034 return fma(lm->get_right()*md->get_left(),
5035 lm->get_left()/md->get_right(),
5036 this->right);
5037 }
5038 auto mdre = exp_cast(md->get_right());
5039 if (mdre.get()) {
5040 return fma(lm->get_right()/md->get_right(),
5041 lm->get_left()*md->get_left(),
5042 this->right);
5043 }
5044 }
5045 }
5046
5047// fma(exp(a)/c, exp(b)*d, e) -> fma(exp(a)*exp(b), d/c, e)
5048// fma(exp(a)/c, d*exp(b), e) -> fma(exp(a)*exp(b), d/c, e)
5049// fma(c/exp(a), exp(b)*d, e) -> fma(exp(b)/exp(a), c*d, e)
5050// fma(c/exp(a), d*exp(b), e) -> fma(exp(b)/exp(a), c*d, e)
5051 if (ld.get() && mm.get()) {
5052 auto ldle = exp_cast(ld->get_left());
5053 if (ldle.get()) {
5054 auto mmle = exp_cast(mm->get_left());
5055 if (mmle.get()) {
5056 return fma(ld->get_left()*mm->get_left(),
5057 mm->get_right()/ld->get_right(),
5058 this->right);
5059 }
5060 auto mmre = exp_cast(mm->get_right());
5061 if (mmre.get()) {
5062 return fma(ld->get_left()*mm->get_right(),
5063 mm->get_left()/ld->get_right(),
5064 this->right);
5065 }
5066 }
5067 auto ldre = exp_cast(ld->get_right());
5068 if (ldre.get()) {
5069 auto mmle = exp_cast(mm->get_left());
5070 if (mmle.get()) {
5071 return fma(mm->get_left()/ld->get_right(),
5072 ld->get_left()*mm->get_right(),
5073 this->right);
5074 }
5075 auto mmre = exp_cast(mm->get_right());
5076 if (mmre.get()) {
5077 return fma(mm->get_right()/ld->get_right(),
5078 ld->get_left()*mm->get_left(),
5079 this->right);
5080 }
5081 }
5082 }
5083
5084// fma(exp(a)/c, exp(b)/d, e) -> (exp(a)*exp(b))/(c*d) + e
5085// fma(exp(a)/c, d/exp(b), e) -> fma(exp(a)/exp(b), d/c, e)
5086// fma(c/exp(a), exp(b)/d, e) -> fma(exp(b)/exp(a), c/d, e)
5087// fma(c/exp(a), d/exp(b), e) -> (c*d)/(exp(a)*exp(b)) + e
5088 if (ld.get() && md.get()) {
5089 auto ldle = exp_cast(ld->get_left());
5090 if (ldle.get()) {
5091 auto mdle = exp_cast(md->get_left());
5092 if (mdle.get()) {
5093 return ((ld->get_left()*md->get_left()) /
5094 (ld->get_right()*md->get_right())) +
5095 this->right;
5096 }
5097 auto mdre = exp_cast(md->get_right());
5098 if (mdre.get()) {
5099 return fma(ld->get_left()/md->get_right(),
5100 md->get_left()/ld->get_right(),
5101 this->right);
5102 }
5103 }
5104 auto ldre = exp_cast(ld->get_right());
5105 if (ldre.get()) {
5106 auto mdle = exp_cast(md->get_left());
5107 if (mdle.get()) {
5108 return fma(md->get_left()/ld->get_right(),
5109 ld->get_left()/md->get_right(),
5110 this->right);
5111 }
5112 auto mdre = exp_cast(md->get_right());
5113 if (mdre.get()) {
5114 return ((ld->get_left()*md->get_left()) /
5115 (ld->get_right()*md->get_right())) +
5116 this->right;
5117 }
5118 }
5119 }
5120
5121// fma(sqrt(a),sqrt(b),c) -> sqrt(a*b) + c
5122 auto lsqr = sqrt_cast(this->left);
5123 auto msqr = sqrt_cast(this->middle);
5124 if (lsqr.get() && msqr.get()) {
5125 return sqrt(lsqr->get_arg()*msqr->get_arg()) + this->right;
5126 }
5127
5128 if constexpr (std::floating_point<T>) {
5129// fma(hypot(b,c),Sqrt(a),d) -> Sqrt((b^2 + c^2)*a) + d
5130// fma(Sqrt(a),hypot(b,c),d) -> Sqrt(a*(b^2 + c^2)) + d
5131 auto lhypot = hypot_cast(this->left);
5132 auto mhypot = hypot_cast(this->middle);
5133 if (lhypot.get() && msqr.get()) {
5134 return sqrt((pow(lhypot->get_left(), static_cast<T> (2)) +
5135 pow(lhypot->get_right(), static_cast<T> (2))) *
5136 msqr->get_arg()) + this->right;
5137 } else if (mhypot.get() && lsqr.get()) {
5138 return sqrt((pow(mhypot->get_left(), static_cast<T> (2)) +
5139 pow(mhypot->get_right(), static_cast<T> (2))) *
5140 lsqr->get_arg()) + this->right;
5141 }
5142
5143// fma(Sqrt(x^2),copysign(1,x),d) -> x
5144// fma(copysign(1,x),Sqrt(x^2),d) -> x
5145 auto lcsc = copysign_cast(this->left);
5146 auto mcsc = copysign_cast(this->middle);
5147 if (mcsc.get() && lsqr.get()) {
5148 auto mcsclc = constant_cast(mcsc->get_left());
5149 auto lsqrpc = pow_cast(lsqr->get_arg());
5150 if (lsqrpc.get() && mcsclc.get() && mcsclc->is(1) &&
5151 lsqrpc->get_left()->is_match(mcsc->get_right())) {
5152 return mcsc->get_right() + this->right;
5153 }
5154 } else if (lcsc.get() && msqr.get()) {
5155 auto lcsclc = constant_cast(lcsc->get_left());
5156 auto msqrpc = pow_cast(msqr->get_arg());
5157 if (msqrpc.get() && lcsclc.get() && lcsclc->is(1) &&
5158 msqrpc->get_left()->is_match(lcsc->get_right())) {
5159 return lcsc->get_right() + this->right;
5160 }
5161 }
5162 }
5163
5164 return this->shared_from_this();
5165 }
5166
5167//------------------------------------------------------------------------------
5174//------------------------------------------------------------------------------
5175 virtual shared_leaf<T, SAFE_MATH>
5176 df(shared_leaf<T, SAFE_MATH> x) {
5177 if (this->is_match(x)) {
5178 return one<T, SAFE_MATH> ();
5179 }
5180
5181 const size_t hash = reinterpret_cast<size_t> (x.get());
5182 if (this->df_cache.find(hash) == this->df_cache.end()) {
5183 auto temp_right = fma(this->left,
5184 this->middle->df(x),
5185 this->right->df(x));
5186
5187 this->df_cache[hash] = fma(this->left->df(x),
5188 this->middle,
5189 temp_right);
5190 }
5191 return this->df_cache[hash];
5192 }
5193
5194//------------------------------------------------------------------------------
5202//------------------------------------------------------------------------------
5203 virtual shared_leaf<T, SAFE_MATH>
5204 compile(std::ostringstream &stream,
5205 jit::register_map &registers,
5206 const jit::register_map &thread_mem,
5207 const jit::register_usage &usage) {
5208 if (registers.find(this) == registers.end()) {
5209 auto l = this->left->compile(stream, registers,
5210 thread_mem, usage);
5211 auto m = this->middle->compile(stream, registers,
5212 thread_mem, usage);
5213 auto r = this->right->compile(stream, registers,
5214 thread_mem, usage);
5215
5216 registers[this] = jit::to_string('r', this);
5217 stream << " const ";
5218 jit::add_type<T> (stream);
5219 stream << " " << registers[this] << " = ";
5220 if constexpr (SAFE_MATH) {
5221 stream << "(" << registers[l.get()] << " == ";
5222 if constexpr (jit::complex_scalar<T>) {
5223 jit::add_type<T> (stream);
5224 stream << "(0, 0)";
5225 } else {
5226 stream << "0";
5227 }
5228 stream << " || " << registers[m.get()] << " == ";
5229 if constexpr (jit::complex_scalar<T>) {
5230 jit::add_type<T> (stream);
5231 stream << "(0, 0)";
5232 } else {
5233 stream << "0";
5234 }
5235 stream << ") ? " << registers[r.get()] << " : ";
5236 }
5237 if constexpr (jit::complex_scalar<T>) {
5238 stream << registers[l.get()] << "*"
5239 << registers[m.get()] << " + "
5240 << registers[r.get()];
5241 } else {
5242 stream << "fma("
5243 << registers[l.get()] << ", "
5244 << registers[m.get()] << ", "
5245 << registers[r.get()] << ")";
5246 }
5247 this->endline(stream, usage);
5248 }
5249
5250 return this->shared_from_this();
5251 }
5252
5253//------------------------------------------------------------------------------
5258//------------------------------------------------------------------------------
5259 virtual bool is_match(shared_leaf<T, SAFE_MATH> x) {
5260 if (this == x.get()) {
5261 return true;
5262 }
5263
5264 auto x_cast = fma_cast(x);
5265 if (x_cast.get()) {
5266 return this->left->is_match(x_cast->get_left()) &&
5267 this->middle->is_match(x_cast->get_middle()) &&
5268 this->right->is_match(x_cast->get_right());
5269 }
5270
5271 return false;
5272 }
5273
5274//------------------------------------------------------------------------------
5276//------------------------------------------------------------------------------
5277 virtual void to_latex() const {
5278 std::cout << "\\left(";
5279 if (add_cast(this->left).get() ||
5280 subtract_cast(this->left).get()) {
5281 std::cout << "\\left(";
5282 this->left->to_latex();
5283 std::cout << "\\right)";
5284 } else {
5285 this->left->to_latex();
5286 }
5287 std::cout << " ";
5288 if (add_cast(this->middle).get() ||
5289 subtract_cast(this->middle).get()) {
5290 std::cout << "\\left(";
5291 this->middle->to_latex();
5292 std::cout << "\\right)";
5293 } else {
5294 this->middle->to_latex();
5295 }
5296 std::cout << "+";
5297 this->right->to_latex();
5298 std::cout << "\\right)";
5299 }
5300
5301//------------------------------------------------------------------------------
5305//------------------------------------------------------------------------------
5306 virtual shared_leaf<T, SAFE_MATH> remove_pseudo() {
5307 if (this->has_pseudo()) {
5308 return fma(this->left->remove_pseudo(),
5309 this->middle->remove_pseudo(),
5310 this->right->remove_pseudo());
5311 }
5312 return this->shared_from_this();
5313 }
5314
5315//------------------------------------------------------------------------------
5321//------------------------------------------------------------------------------
5322 virtual shared_leaf<T, SAFE_MATH> to_vizgraph(std::stringstream &stream,
5323 jit::register_map &registers) {
5324 if (registers.find(this) == registers.end()) {
5325 const std::string name = jit::to_string('r', this);
5326 registers[this] = name;
5327 stream << " " << name
5328 << " [label = \"fma\", shape = oval, style = filled, fillcolor = blue, fontcolor = white];" << std::endl;
5329
5330 auto l = this->left->to_vizgraph(stream, registers);
5331 stream << " " << name << " -- " << registers[l.get()] << ";" << std::endl;
5332 auto m = this->middle->to_vizgraph(stream, registers);
5333 stream << " " << name << " -- " << registers[m.get()] << ";" << std::endl;
5334 auto r = this->right->to_vizgraph(stream, registers);
5335 stream << " " << name << " -- " << registers[r.get()] << ";" << std::endl;
5336 }
5337
5338 return this->shared_from_this();
5339 }
5340 };
5341
5342//------------------------------------------------------------------------------
5351//------------------------------------------------------------------------------
5352 template<jit::float_scalar T, bool SAFE_MATH=false>
5353 shared_leaf<T, SAFE_MATH> fma(shared_leaf<T, SAFE_MATH> l,
5354 shared_leaf<T, SAFE_MATH> m,
5355 shared_leaf<T, SAFE_MATH> r) {
5356 auto temp = std::make_shared<fma_node<T, SAFE_MATH>> (l, m, r)->reduce();
5357// Test for hash collisions.
5358 for (size_t i = temp->get_hash();
5359 i < std::numeric_limits<size_t>::max(); i++) {
5360 if (leaf_node<T, SAFE_MATH>::caches.nodes.find(i) ==
5361 leaf_node<T, SAFE_MATH>::caches.nodes.end()) {
5362 leaf_node<T, SAFE_MATH>::caches.nodes[i] = temp;
5363 return temp;
5364 } else if (temp->is_match(leaf_node<T, SAFE_MATH>::caches.nodes[i])) {
5365 return leaf_node<T, SAFE_MATH>::caches.nodes[i];
5366 }
5367 }
5368#if defined(__clang__) || defined(__GNUC__)
5369 __builtin_unreachable();
5370#else
5371 assert(false && "Should never reach.");
5372#endif
5373 }
5374
5375//------------------------------------------------------------------------------
5388//------------------------------------------------------------------------------
5389 template<jit::float_scalar T, jit::float_scalar L, bool SAFE_MATH=false>
5395
5396//------------------------------------------------------------------------------
5409//------------------------------------------------------------------------------
5410 template<jit::float_scalar T, jit::float_scalar M, bool SAFE_MATH=false>
5416
5417//------------------------------------------------------------------------------
5430//------------------------------------------------------------------------------
5431 template<jit::float_scalar T, jit::float_scalar R, bool SAFE_MATH=false>
5437
5438//------------------------------------------------------------------------------
5452//------------------------------------------------------------------------------
5453 template<jit::float_scalar T, jit::float_scalar L, jit::float_scalar M, bool SAFE_MATH=false>
5454 shared_leaf<T, SAFE_MATH> fma(const L l,
5455 const M m,
5456 shared_leaf<T, SAFE_MATH> r) {
5457 return fma<T, SAFE_MATH> (constant<T, SAFE_MATH> (static_cast<T> (l)),
5458 constant<T, SAFE_MATH> (static_cast<T> (m)), r);
5459 }
5460
5461//------------------------------------------------------------------------------
5475//------------------------------------------------------------------------------
5476 template<jit::float_scalar T, jit::float_scalar M, jit::float_scalar R, bool SAFE_MATH=false>
5477 shared_leaf<T, SAFE_MATH> fma(shared_leaf<T, SAFE_MATH> l,
5478 const M m,
5479 const R r) {
5480 return fma<T, SAFE_MATH> (l, constant<T, SAFE_MATH> (static_cast<T> (m)),
5481 constant<T, SAFE_MATH> (static_cast<T> (r)));
5482 }
5483
5484//------------------------------------------------------------------------------
5498//------------------------------------------------------------------------------
5499 template<jit::float_scalar T, jit::float_scalar L, jit::float_scalar R, bool SAFE_MATH=false>
5500 shared_leaf<T, SAFE_MATH> fma(const L l,
5501 shared_leaf<T, SAFE_MATH> m,
5502 const R r) {
5503 return fma<T, SAFE_MATH> (constant<T, SAFE_MATH> (static_cast<T> (l)), m,
5504 constant<T, SAFE_MATH> (static_cast<T> (r)));
5505 }
5506
5508 template<jit::float_scalar T, bool SAFE_MATH=false>
5509 using shared_fma = std::shared_ptr<fma_node<T, SAFE_MATH>>;
5510
5511//------------------------------------------------------------------------------
5519//------------------------------------------------------------------------------
5520 template<jit::float_scalar T, bool SAFE_MATH=false>
5521 shared_fma<T, SAFE_MATH> fma_cast(shared_leaf<T, SAFE_MATH> x) {
5522 return std::dynamic_pointer_cast<fma_node<T, SAFE_MATH>> (x);
5523 }
5524
5525//******************************************************************************
5526// Modulo node.
5527//******************************************************************************
5528//------------------------------------------------------------------------------
5535//------------------------------------------------------------------------------
5536 template<std::floating_point T, bool SAFE_MATH=false>
5537 class modulo_node final : public no_derivative<T, SAFE_MATH, branch_node<T, SAFE_MATH>> {
5538 private:
5539//------------------------------------------------------------------------------
5545//------------------------------------------------------------------------------
5546 static std::string to_string(leaf_node<T, SAFE_MATH> *l,
5547 leaf_node<T, SAFE_MATH> *r) {
5548 return jit::format_to_string(reinterpret_cast<size_t> (l)) + "%" +
5549 jit::format_to_string(reinterpret_cast<size_t> (r));
5550 }
5551
5552 public:
5553//------------------------------------------------------------------------------
5558//------------------------------------------------------------------------------
5565
5566//------------------------------------------------------------------------------
5572//------------------------------------------------------------------------------
5573 virtual backend::buffer<T> evaluate() {
5574 backend::buffer<T> l_result = this->left->evaluate();
5575 backend::buffer<T> r_result = this->right->evaluate();
5576 return l_result % r_result;
5577 }
5578
5579//------------------------------------------------------------------------------
5583//------------------------------------------------------------------------------
5584 virtual shared_leaf<T, SAFE_MATH> reduce() {
5585// Constant reductions.
5586 auto l = constant_cast(this->left);
5587 auto r = constant_cast(this->right);
5588
5589 if (l.get() && r.get()) {
5590 return constant<T, SAFE_MATH> (this->evaluate());
5591 }
5592
5593 return this->shared_from_this();
5594 }
5595
5596//------------------------------------------------------------------------------
5604//------------------------------------------------------------------------------
5605 virtual shared_leaf<T, SAFE_MATH>
5606 compile(std::ostringstream &stream,
5607 jit::register_map &registers,
5608 const jit::register_map &thread_mem,
5609 const jit::register_usage &usage) {
5610 if (registers.find(this) == registers.end()) {
5611 auto l = this->left->compile(stream, registers,
5612 thread_mem, usage);
5613 auto r = this->right->compile(stream, registers,
5614 thread_mem, usage);
5615
5616 registers[this] = jit::to_string('r', this);
5617 stream << " const ";
5618 jit::add_type<T> (stream);
5619 stream << " " << registers[this] << " = fmod("
5620 << registers[l.get()] << ","
5621 << registers[r.get()] << ")";
5622 this->endline(stream, usage);
5623 }
5624
5625 return this->shared_from_this();
5626 }
5627
5628//------------------------------------------------------------------------------
5630//------------------------------------------------------------------------------
5631 virtual void to_latex() const {
5632 bool l_brackets = add_cast(this->left).get() ||
5633 subtract_cast(this->left).get();
5634 bool r_brackets = add_cast(this->right).get() ||
5635 subtract_cast(this->right).get();
5636 if (l_brackets) {
5637 std::cout << "\\left(";
5638 }
5639 this->left->to_latex();
5640 if (l_brackets) {
5641 std::cout << "\\right)";
5642 }
5643 std::cout << "\%";
5644 if (r_brackets) {
5645 std::cout << "\\left(";
5646 }
5647 this->right->to_latex();
5648 if (r_brackets) {
5649 std::cout << "\\right)";
5650 }
5651 }
5652
5653//------------------------------------------------------------------------------
5659//------------------------------------------------------------------------------
5660 virtual shared_leaf<T, SAFE_MATH> to_vizgraph(std::stringstream &stream,
5661 jit::register_map &registers) {
5662 if (registers.find(this) == registers.end()) {
5663 const std::string name = jit::to_string('r', this);
5664 registers[this] = name;
5665 stream << " " << name
5666 << " [label = \"%\", shape = oval, style = filled, fillcolor = blue, fontcolor = white];" << std::endl;
5667
5668 auto l = this->left->to_vizgraph(stream, registers);
5669 stream << " " << name << " -- " << registers[l.get()] << ";" << std::endl;
5670 auto r = this->right->to_vizgraph(stream, registers);
5671 stream << " " << name << " -- " << registers[r.get()] << ";" << std::endl;
5672 }
5673
5674 return this->shared_from_this();
5675 }
5676 };
5677
5678//------------------------------------------------------------------------------
5686//------------------------------------------------------------------------------
5687 template<std::floating_point T, bool SAFE_MATH=false>
5688 shared_leaf<T, SAFE_MATH> modulo(shared_leaf<T, SAFE_MATH> l,
5689 shared_leaf<T, SAFE_MATH> r) {
5690 auto temp = std::make_shared<modulo_node<T, SAFE_MATH>> (l, r)->reduce();
5691// Test for hash collisions.
5692 for (size_t i = temp->get_hash();
5693 i < std::numeric_limits<size_t>::max(); i++) {
5694 if (leaf_node<T, SAFE_MATH>::caches.nodes.find(i) ==
5695 leaf_node<T, SAFE_MATH>::caches.nodes.end()) {
5696 leaf_node<T, SAFE_MATH>::caches.nodes[i] = temp;
5697 return temp;
5698 } else if (temp->is_match(leaf_node<T, SAFE_MATH>::caches.nodes[i])) {
5699 return leaf_node<T, SAFE_MATH>::caches.nodes[i];
5700 }
5701 }
5702#if defined(__clang__) || defined(__GNUC__)
5703 __builtin_unreachable();
5704#else
5705 assert(false && "Should never reach.");
5706#endif
5707 }
5708
5709//------------------------------------------------------------------------------
5720//------------------------------------------------------------------------------
5721 template<std::floating_point T, bool SAFE_MATH=false>
5726
5727//------------------------------------------------------------------------------
5739//------------------------------------------------------------------------------
5740 template<std::floating_point T, std::floating_point L, bool SAFE_MATH=false>
5745
5746//------------------------------------------------------------------------------
5758//------------------------------------------------------------------------------
5759 template<std::floating_point T, std::floating_point R, bool SAFE_MATH=false>
5764
5766 template<std::floating_point T, bool SAFE_MATH=false>
5767 using shared_modulo = std::shared_ptr<modulo_node<T, SAFE_MATH>>;
5768
5769//------------------------------------------------------------------------------
5777//------------------------------------------------------------------------------
5778 template<std::floating_point T, bool SAFE_MATH=false>
5779 shared_modulo<T, SAFE_MATH> modulo_cast(shared_leaf<T, SAFE_MATH> x) {
5780 return std::dynamic_pointer_cast<modulo_node<T, SAFE_MATH>> (x);
5781 }
5782}
5783
5784#endif /* arithmetic_h */
backend::buffer
Class representing a generic buffer.
Definition backend.hpp:29
backend::buffer::subtract_row
void subtract_row(const buffer< T > &x)
Subtract row operation.
Definition backend.hpp:413
backend::buffer::multiply_row
void multiply_row(const buffer< T > &x)
Multiply row operation.
Definition backend.hpp:437
backend::buffer::add_col
void add_col(const buffer< T > &x)
Add col operation.
Definition backend.hpp:401
backend::buffer::divide_col
void divide_col(const buffer< T > &x)
Divide col operation.
Definition backend.hpp:473
backend::buffer::add_row
void add_row(const buffer< T > &x)
Add row operation.
Definition backend.hpp:356
backend::buffer::subtract_col
void subtract_col(const buffer< T > &x)
Subtract col operation.
Definition backend.hpp:425
backend::buffer::multiply_col
void multiply_col(const buffer< T > &x)
Multiply col operation.
Definition backend.hpp:449
backend::buffer::divide_row
void divide_row(const buffer< T > &x)
Divide row operation.
Definition backend.hpp:461
backend::buffer::is_zero
bool is_zero() const
Is every element zero.
Definition backend.hpp:126
graph::add_node
An addition node.
Definition arithmetic.hpp:132
graph::add_node::evaluate
virtual backend::buffer< T > evaluate()
Evaluate the results of addition.
Definition arithmetic.hpp:166
graph::add_node::to_vizgraph
virtual shared_leaf< T, SAFE_MATH > to_vizgraph(std::stringstream &stream, jit::register_map &registers)
Convert the node to vizgraph.
Definition arithmetic.hpp:733
graph::add_node::df
virtual shared_leaf< T, SAFE_MATH > df(shared_leaf< T, SAFE_MATH > x)
Transform node to derivative.
Definition arithmetic.hpp:621
graph::add_node::is_match
virtual bool is_match(shared_leaf< T, SAFE_MATH > x)
Query if the nodes match.
Definition arithmetic.hpp:671
graph::add_node::remove_pseudo
virtual shared_leaf< T, SAFE_MATH > remove_pseudo()
Remove pseudo variable nodes.
Definition arithmetic.hpp:718
graph::add_node::add_node
add_node(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > r)
Construct an addition node.
Definition arithmetic.hpp:154
graph::add_node::reduce
virtual shared_leaf< T, SAFE_MATH > reduce()
Reduce an addition node.
Definition arithmetic.hpp:177
graph::add_node::to_latex
virtual void to_latex() const
Convert the node to latex.
Definition arithmetic.hpp:691
graph::add_node::compile
virtual shared_leaf< T, SAFE_MATH > compile(std::ostringstream &stream, jit::register_map &registers, const jit::register_map &thread_mem, const jit::register_usage &usage)
Compile the node.
Definition arithmetic.hpp:643
graph::branch_node
Class representing a branch node.
Definition node.hpp:1337
graph::branch_node::right
shared_leaf< T, SAFE_MATH > right
Right branch of the tree.
Definition node.hpp:1342
graph::branch_node::left
shared_leaf< T, SAFE_MATH > left
Left branch of the tree.
Definition node.hpp:1340
graph::divide_node
A division node.
Definition arithmetic.hpp:2857
graph::divide_node::reduce
virtual shared_leaf< T, SAFE_MATH > reduce()
Reduce an division node.
Definition arithmetic.hpp:2910
graph::divide_node::remove_pseudo
virtual shared_leaf< T, SAFE_MATH > remove_pseudo()
Remove pseudo variable nodes.
Definition arithmetic.hpp:3668
graph::divide_node::divide_node
divide_node(shared_leaf< T, SAFE_MATH > n, shared_leaf< T, SAFE_MATH > d)
Construct an addition node.
Definition arithmetic.hpp:2879
graph::divide_node::is_match
virtual bool is_match(shared_leaf< T, SAFE_MATH > x)
Query if the nodes match.
Definition arithmetic.hpp:3638
graph::divide_node::evaluate
virtual backend::buffer< T > evaluate()
Evaluate the results of division.
Definition arithmetic.hpp:2891
graph::divide_node::to_vizgraph
virtual shared_leaf< T, SAFE_MATH > to_vizgraph(std::stringstream &stream, jit::register_map &registers)
Convert the node to vizgraph.
Definition arithmetic.hpp:3683
graph::divide_node::to_latex
virtual void to_latex() const
Convert the node to latex.
Definition arithmetic.hpp:3655
graph::divide_node::compile
virtual shared_leaf< T, SAFE_MATH > compile(std::ostringstream &stream, jit::register_map &registers, const jit::register_map &thread_mem, const jit::register_usage &usage)
Compile the node.
Definition arithmetic.hpp:3594
graph::divide_node::df
virtual shared_leaf< T, SAFE_MATH > df(shared_leaf< T, SAFE_MATH > x)
Transform node to derivative.
Definition arithmetic.hpp:3571
graph::fma_node
A fused multiply add node.
Definition arithmetic.hpp:3818
graph::fma_node::remove_pseudo
virtual shared_leaf< T, SAFE_MATH > remove_pseudo()
Remove pseudo variable nodes.
Definition arithmetic.hpp:5306
graph::fma_node::df
virtual shared_leaf< T, SAFE_MATH > df(shared_leaf< T, SAFE_MATH > x)
Transform node to derivative.
Definition arithmetic.hpp:5176
graph::fma_node::compile
virtual shared_leaf< T, SAFE_MATH > compile(std::ostringstream &stream, jit::register_map &registers, const jit::register_map &thread_mem, const jit::register_usage &usage)
Compile the node.
Definition arithmetic.hpp:5204
graph::fma_node::evaluate
virtual backend::buffer< T > evaluate()
Evaluate the results of fused multiply add.
Definition arithmetic.hpp:3894
graph::fma_node::to_vizgraph
virtual shared_leaf< T, SAFE_MATH > to_vizgraph(std::stringstream &stream, jit::register_map &registers)
Convert the node to vizgraph.
Definition arithmetic.hpp:5322
graph::fma_node::fma_node
fma_node(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > m, shared_leaf< T, SAFE_MATH > r)
Construct a fused multiply add node.
Definition arithmetic.hpp:3880
graph::fma_node::to_latex
virtual void to_latex() const
Convert the node to latex.
Definition arithmetic.hpp:5277
graph::fma_node::reduce
virtual shared_leaf< T, SAFE_MATH > reduce()
Reduce a fused multiply add node.
Definition arithmetic.hpp:3913
graph::fma_node::is_match
virtual bool is_match(shared_leaf< T, SAFE_MATH > x)
Query if the nodes match.
Definition arithmetic.hpp:5259
graph::leaf_node
Class representing a node leaf.
Definition node.hpp:371
graph::leaf_node::compile
virtual std::shared_ptr< leaf_node< T, SAFE_MATH > > compile(std::ostringstream &stream, jit::register_map &registers, const jit::register_map &thread_mem, const jit::register_usage &usage)=0
Compile the node.
graph::leaf_node::endline
virtual void endline(std::ostringstream &stream, const jit::register_usage &usage, const char end=';') const
End a line in the kernel source.
Definition node.hpp:653
graph::leaf_node::to_vizgraph
virtual std::shared_ptr< leaf_node< T, SAFE_MATH > > to_vizgraph(std::stringstream &stream, jit::register_map &registers)=0
Convert the node to vizgraph.
graph::leaf_node::evaluate
virtual backend::buffer< T > evaluate()=0
Evaluate method.
graph::leaf_node::has_pseudo
virtual bool has_pseudo() const
Query if the node contains pseudo variables.
Definition node.hpp:633
graph::leaf_node::hash
const size_t hash
Hash for node.
Definition node.hpp:374
graph::leaf_node::df_cache
std::unordered_map< size_t, std::shared_ptr< leaf_node< T, SAFE_MATH > > > df_cache
Cache derivative terms.
Definition node.hpp:378
graph::modulo_node
An Modulo node.
Definition arithmetic.hpp:5537
graph::modulo_node::modulo_node
modulo_node(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > r)
Construct an modulo node.
Definition arithmetic.hpp:5559
graph::modulo_node::to_latex
virtual void to_latex() const
Convert the node to latex.
Definition arithmetic.hpp:5631
graph::modulo_node::to_vizgraph
virtual shared_leaf< T, SAFE_MATH > to_vizgraph(std::stringstream &stream, jit::register_map &registers)
Convert the node to vizgraph.
Definition arithmetic.hpp:5660
graph::modulo_node::compile
virtual shared_leaf< T, SAFE_MATH > compile(std::ostringstream &stream, jit::register_map &registers, const jit::register_map &thread_mem, const jit::register_usage &usage)
Compile the node.
Definition arithmetic.hpp:5606
graph::modulo_node::evaluate
virtual backend::buffer< T > evaluate()
Evaluate the results of modulo.
Definition arithmetic.hpp:5573
graph::modulo_node::reduce
virtual shared_leaf< T, SAFE_MATH > reduce()
Reduce an modulo node.
Definition arithmetic.hpp:5584
graph::multiply_node
A multiplication node.
Definition arithmetic.hpp:1724
graph::multiply_node::evaluate
virtual backend::buffer< T > evaluate()
Evaluate the results of multiplication.
Definition arithmetic.hpp:1863
graph::multiply_node::compile
virtual shared_leaf< T, SAFE_MATH > compile(std::ostringstream &stream, jit::register_map &registers, const jit::register_map &thread_mem, const jit::register_usage &usage)
Compile the node.
Definition arithmetic.hpp:2608
graph::multiply_node::df
virtual shared_leaf< T, SAFE_MATH > df(shared_leaf< T, SAFE_MATH > x)
Transform node to derivative.
Definition arithmetic.hpp:2585
graph::multiply_node::remove_pseudo
virtual shared_leaf< T, SAFE_MATH > remove_pseudo()
Remove pseudo variable nodes.
Definition arithmetic.hpp:2709
graph::multiply_node::reduce
virtual shared_leaf< T, SAFE_MATH > reduce()
Reduce an multiplication node.
Definition arithmetic.hpp:1888
graph::multiply_node::multiply_node
multiply_node(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > r)
Construct a multiplication node.
Definition arithmetic.hpp:1852
graph::multiply_node::is_match
virtual bool is_match(shared_leaf< T, SAFE_MATH > x)
Query if the nodes match.
Definition arithmetic.hpp:2660
graph::multiply_node::to_vizgraph
virtual shared_leaf< T, SAFE_MATH > to_vizgraph(std::stringstream &stream, jit::register_map &registers)
Convert the node to vizgraph.
Definition arithmetic.hpp:2724
graph::multiply_node::to_latex
virtual void to_latex() const
Convert the node to latex.
Definition arithmetic.hpp:2682
graph::no_derivative
Nodes without derivatives.
Definition node.hpp:1713
graph::subtract_node
A subtraction node.
Definition arithmetic.hpp:871
graph::subtract_node::compile
virtual shared_leaf< T, SAFE_MATH > compile(std::ostringstream &stream, jit::register_map &registers, const jit::register_map &thread_mem, const jit::register_usage &usage)
Compile the node.
Definition arithmetic.hpp:1483
graph::subtract_node::subtract_node
subtract_node(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > r)
Construct a subtraction node.
Definition arithmetic.hpp:893
graph::subtract_node::df
virtual shared_leaf< T, SAFE_MATH > df(shared_leaf< T, SAFE_MATH > x)
Transform node to derivative.
Definition arithmetic.hpp:1461
graph::subtract_node::remove_pseudo
virtual shared_leaf< T, SAFE_MATH > remove_pseudo()
Remove pseudo variable nodes.
Definition arithmetic.hpp:1555
graph::subtract_node::is_match
virtual bool is_match(shared_leaf< T, SAFE_MATH > x)
Query if the nodes match.
Definition arithmetic.hpp:1511
graph::subtract_node::reduce
virtual shared_leaf< T, SAFE_MATH > reduce()
Reduce an subtraction node.
Definition arithmetic.hpp:916
graph::subtract_node::to_vizgraph
virtual shared_leaf< T, SAFE_MATH > to_vizgraph(std::stringstream &stream, jit::register_map &registers)
Convert the node to vizgraph.
Definition arithmetic.hpp:1570
graph::subtract_node::evaluate
virtual backend::buffer< T > evaluate()
Evaluate the results of subtraction.
Definition arithmetic.hpp:905
graph::subtract_node::to_latex
virtual void to_latex() const
Convert the node to latex.
Definition arithmetic.hpp:1528
graph::triple_node
Class representing a triple branch node.
Definition node.hpp:1467
graph::triple_node::middle
shared_leaf< T, SAFE_MATH > middle
Middle branch of the tree.
Definition node.hpp:1470
jit::complex_scalar
Complex scalar concept.
Definition register.hpp:25
assert
subroutine assert(test, message)
Assert check.
Definition f_binding_test.f90:38
backend::fma
buffer< T > fma(buffer< T > &a, buffer< T > &b, buffer< T > &c)
Fused multiply add operation.
Definition backend.hpp:1219
graph
Name space for graph nodes.
Definition arithmetic.hpp:13
graph::piecewise_2D_cast
shared_piecewise_2D< T, SAFE_MATH > piecewise_2D_cast(shared_leaf< T, SAFE_MATH > x)
Cast to a piecewise 2D node.
Definition piecewise.hpp:1508
graph::operator*
shared_leaf< T, SAFE_MATH > operator*(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > r)
Build multiply node from two leaves.
Definition arithmetic.hpp:2786
graph::shared_add
std::shared_ptr< add_node< T, SAFE_MATH > > shared_add
Convenience type alias for shared add nodes.
Definition arithmetic.hpp:843
graph::pow_cast
shared_pow< T, SAFE_MATH > pow_cast(shared_leaf< T, SAFE_MATH > x)
Cast to a power node.
Definition math.hpp:1448
graph::subtract_cast
shared_subtract< T, SAFE_MATH > subtract_cast(shared_leaf< T, SAFE_MATH > x)
Cast to a subtract node.
Definition arithmetic.hpp:1710
graph::modulo
shared_leaf< T, SAFE_MATH > modulo(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > r)
Build modulo node from two leaves.
Definition arithmetic.hpp:5688
graph::zero
constexpr shared_leaf< T, SAFE_MATH > zero()
Forward declare for zero.
Definition node.hpp:1154
graph::operator/
shared_leaf< T, SAFE_MATH > operator/(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > r)
Build divide node from two leaves.
Definition arithmetic.hpp:3745
graph::sqrt_cast
shared_sqrt< T, SAFE_MATH > sqrt_cast(shared_leaf< T, SAFE_MATH > x)
Cast to a sqrt node.
Definition math.hpp:331
graph::multiply
shared_leaf< T, SAFE_MATH > multiply(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > r)
Build multiply node from two leaves.
Definition arithmetic.hpp:2752
graph::add_cast
shared_add< T, SAFE_MATH > add_cast(shared_leaf< T, SAFE_MATH > x)
Cast to a add node.
Definition arithmetic.hpp:855
graph::is_greater_exponent
bool is_greater_exponent(shared_leaf< T, SAFE_MATH > a, shared_leaf< T, SAFE_MATH > b)
Check if the exponent is greater than the other.
Definition arithmetic.hpp:111
graph::operator+
shared_leaf< T, SAFE_MATH > operator+(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > r)
Build add node from two leaves.
Definition arithmetic.hpp:798
graph::shared_divide
std::shared_ptr< divide_node< T, SAFE_MATH > > shared_divide
Convenience type alias for shared divide nodes.
Definition arithmetic.hpp:3790
graph::pow
shared_leaf< T, SAFE_MATH > pow(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > r)
Build power node.
Definition math.hpp:1384
graph::sin_cast
shared_sine< T, SAFE_MATH > sin_cast(shared_leaf< T, SAFE_MATH > x)
Cast to a sine node.
Definition trigonometry.hpp:258
graph::divide_cast
shared_divide< T, SAFE_MATH > divide_cast(shared_leaf< T, SAFE_MATH > x)
Cast to a divide node.
Definition arithmetic.hpp:3802
graph::piecewise_1D_cast
shared_piecewise_1D< T, SAFE_MATH > piecewise_1D_cast(shared_leaf< T, SAFE_MATH > x)
Cast to a piecewise 1D node.
Definition piecewise.hpp:849
graph::operator-
shared_leaf< T, SAFE_MATH > operator-(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > r)
Build subtract node from two leaves.
Definition arithmetic.hpp:1634
graph::divide
shared_leaf< T, SAFE_MATH > divide(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > r)
Build divide node from two leaves.
Definition arithmetic.hpp:3711
graph::multiply_cast
shared_multiply< T, SAFE_MATH > multiply_cast(shared_leaf< T, SAFE_MATH > x)
Cast to a multiply node.
Definition arithmetic.hpp:2843
graph::shared_modulo
std::shared_ptr< modulo_node< T, SAFE_MATH > > shared_modulo
Convenience type alias for shared modulo nodes.
Definition arithmetic.hpp:5767
graph::modulo_cast
shared_modulo< T, SAFE_MATH > modulo_cast(shared_leaf< T, SAFE_MATH > x)
Cast to a modulo node.
Definition arithmetic.hpp:5779
graph::exp
shared_leaf< T, SAFE_MATH > exp(shared_leaf< T, SAFE_MATH > x)
Define exp convenience function.
Definition math.hpp:558
graph::exp_cast
shared_exp< T, SAFE_MATH > exp_cast(shared_leaf< T, SAFE_MATH > x)
Cast to a exp node.
Definition math.hpp:592
graph::is_constant_promotable
bool is_constant_promotable(shared_leaf< T, SAFE_MATH > a, shared_leaf< T, SAFE_MATH > b)
Check if the constants are promotable.
Definition arithmetic.hpp:53
graph::fma
shared_leaf< T, SAFE_MATH > fma(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > m, shared_leaf< T, SAFE_MATH > r)
Build fused multiply add node.
Definition arithmetic.hpp:5353
graph::constant_cast
shared_constant< T, SAFE_MATH > constant_cast(shared_leaf< T, SAFE_MATH > x)
Cast to a constant node.
Definition node.hpp:1202
graph::is_variable_combinable
bool is_variable_combinable(shared_leaf< T, SAFE_MATH > a, shared_leaf< T, SAFE_MATH > b)
Check if the variable is combinable.
Definition arithmetic.hpp:75
graph::i
constexpr T i
Convenience type for imaginary constant.
Definition node.hpp:1186
graph::subtract
shared_leaf< T, SAFE_MATH > subtract(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > r)
Build subtract node from two leaves.
Definition arithmetic.hpp:1599
graph::shared_fma
std::shared_ptr< fma_node< T, SAFE_MATH > > shared_fma
Convenience type alias for shared add nodes.
Definition arithmetic.hpp:5509
graph::sqrt
shared_leaf< T, SAFE_MATH > sqrt(shared_leaf< T, SAFE_MATH > x)
Define sqrt convenience function.
Definition math.hpp:297
graph::shared_leaf
std::shared_ptr< leaf_node< T, SAFE_MATH > > shared_leaf
Convenience type alias for shared leaf nodes.
Definition node.hpp:691
graph::hypot_cast
shared_hypot< T, SAFE_MATH > hypot_cast(shared_leaf< T, SAFE_MATH > x)
Cast to a hypot node.
Definition math.hpp:2287
graph::operator%
shared_leaf< T, SAFE_MATH > operator%(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > r)
Build modulo node from two leaves.
Definition arithmetic.hpp:5722
graph::is_variable_promotable
bool is_variable_promotable(shared_leaf< T, SAFE_MATH > a, shared_leaf< T, SAFE_MATH > b)
Check if the variable is variable is promotable.
Definition arithmetic.hpp:91
graph::shared_multiply
std::shared_ptr< multiply_node< T, SAFE_MATH > > shared_multiply
Convenience type alias for shared multiply nodes.
Definition arithmetic.hpp:2831
graph::add
shared_leaf< T, SAFE_MATH > add(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > r)
Build add node from two leaves.
Definition arithmetic.hpp:764
graph::fma_cast
shared_fma< T, SAFE_MATH > fma_cast(shared_leaf< T, SAFE_MATH > x)
Cast to a fma node.
Definition arithmetic.hpp:5521
graph::copysign_cast
shared_copysign< T, SAFE_MATH > copysign_cast(shared_leaf< T, SAFE_MATH > x)
Cast to a copysign node.
Definition math.hpp:2597
graph::cos_cast
shared_cosine< T, SAFE_MATH > cos_cast(shared_leaf< T, SAFE_MATH > x)
Cast to a cosine node.
Definition trigonometry.hpp:512
graph::is_constant_combinable
bool is_constant_combinable(shared_leaf< T, SAFE_MATH > a, shared_leaf< T, SAFE_MATH > b)
Check if nodes are constant combinable.
Definition arithmetic.hpp:25
graph::shared_subtract
std::shared_ptr< subtract_node< T, SAFE_MATH > > shared_subtract
Convenience type alias for shared subtract nodes.
Definition arithmetic.hpp:1698
jit::format_to_string
std::string format_to_string(const T value)
Convert a value to a string while avoiding locale.
Definition register.hpp:213
jit::register_usage
std::unordered_map< void *, size_t > register_usage
Type alias for counting register usage.
Definition register.hpp:261
jit::register_map
std::unordered_map< void *, std::string > register_map
Type alias for mapping node pointers to register names.
Definition register.hpp:259
jit::to_string
std::string to_string(const char prefix, const NODE *pointer)
Convert a graph::leaf_node pointer to a string.
Definition register.hpp:247
node.hpp
Base nodes of graph computation framework.
piecewise_1D
void piecewise_1D()
Tests for 1D piecewise nodes.
Definition piecewise_test.cpp:82
piecewise_2D
void piecewise_2D()
Tests for 2D piecewise nodes.
Definition piecewise_test.cpp:321