DataTree.cc 15.9 KB
Newer Older
sebastien's avatar
sebastien committed
1
/*
2
 * Copyright (C) 2003-2011 Dynare Team
sebastien's avatar
sebastien committed
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 *
 * This file is part of Dynare.
 *
 * Dynare is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Dynare is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 */

20
#include <cstdlib>
21
#include <cassert>
22
#include <iostream>
sebastien's avatar
sebastien committed
23

24
25
#include "DataTree.hh"

26
27
28
DataTree::DataTree(SymbolTable &symbol_table_arg,
                   NumericalConstants &num_constants_arg,
                   ExternalFunctionsTable &external_functions_table_arg) :
29
  symbol_table(symbol_table_arg),
sebastien's avatar
sebastien committed
30
  num_constants(num_constants_arg),
31
  external_functions_table(external_functions_table_arg),
32
  node_counter(0)
sebastien's avatar
sebastien committed
33
{
34
35
36
  Zero = AddNonNegativeConstant("0");
  One = AddNonNegativeConstant("1");
  Two = AddNonNegativeConstant("2");
sebastien's avatar
sebastien committed
37
38

  MinusOne = AddUMinus(One);
39

40
41
  NaN = AddNonNegativeConstant("NaN");
  Infinity = AddNonNegativeConstant("Inf");
42
  MinusInfinity = AddUMinus(Infinity);
sebastien's avatar
sebastien committed
43

44
  Pi = AddNonNegativeConstant("3.141592653589793");
sebastien's avatar
sebastien committed
45
46
}

47
48
DataTree::~DataTree()
{
49
  for (node_list_t::iterator it = node_list.begin(); it != node_list.end(); it++)
sebastien's avatar
sebastien committed
50
51
52
    delete *it;
}

53
expr_t
54
DataTree::AddNonNegativeConstant(const string &value) throw (NumericalConstants::InvalidFloatingPointNumberException)
sebastien's avatar
sebastien committed
55
{
56
  int id = num_constants.AddNonNegativeConstant(value);
sebastien's avatar
sebastien committed
57

58
  num_const_node_map_t::iterator it = num_const_node_map.find(id);
sebastien's avatar
sebastien committed
59
60
61
62
63
64
  if (it != num_const_node_map.end())
    return it->second;
  else
    return new NumConstNode(*this, id);
}

sebastien's avatar
sebastien committed
65
66
VariableNode *
DataTree::AddVariableInternal(int symb_id, int lag)
sebastien's avatar
sebastien committed
67
{
68
  variable_node_map_t::iterator it = variable_node_map.find(make_pair(symb_id, lag));
sebastien's avatar
sebastien committed
69
70
71
  if (it != variable_node_map.end())
    return it->second;
  else
sebastien's avatar
sebastien committed
72
    return new VariableNode(*this, symb_id, lag);
sebastien's avatar
sebastien committed
73
}
sebastien's avatar
sebastien committed
74
75
76

VariableNode *
DataTree::AddVariable(int symb_id, int lag)
sebastien's avatar
sebastien committed
77
{
78
  assert(lag == 0);
sebastien's avatar
sebastien committed
79
  return AddVariableInternal(symb_id, lag);
sebastien's avatar
sebastien committed
80
81
}

82
83
expr_t
DataTree::AddPlus(expr_t iArg1, expr_t iArg2)
sebastien's avatar
sebastien committed
84
85
86
87
88
{
  if (iArg1 != Zero && iArg2 != Zero)
    {
      // Simplify x+(-y) in x-y
      UnaryOpNode *uarg2 = dynamic_cast<UnaryOpNode *>(iArg2);
sebastien's avatar
sebastien committed
89
90
      if (uarg2 != NULL && uarg2->get_op_code() == oUminus)
        return AddMinus(iArg1, uarg2->get_arg());
sebastien's avatar
sebastien committed
91
92
93
94
95

      // To treat commutativity of "+"
      // Nodes iArg1 and iArg2 are sorted by index
      if (iArg1->idx > iArg2->idx)
        {
96
          expr_t tmp = iArg1;
sebastien's avatar
sebastien committed
97
98
99
100
101
102
103
104
105
106
107
108
109
          iArg1 = iArg2;
          iArg2 = tmp;
        }
      return AddBinaryOp(iArg1, oPlus, iArg2);
    }
  else if (iArg1 != Zero)
    return iArg1;
  else if (iArg2 != Zero)
    return iArg2;
  else
    return Zero;
}

110
111
expr_t
DataTree::AddMinus(expr_t iArg1, expr_t iArg2)
sebastien's avatar
sebastien committed
112
{
113
  if (iArg2 == Zero)
sebastien's avatar
sebastien committed
114
    return iArg1;
115
116

  if (iArg1 == Zero)
sebastien's avatar
sebastien committed
117
    return AddUMinus(iArg2);
118
119

  if (iArg1 == iArg2)
sebastien's avatar
sebastien committed
120
    return Zero;
121
122

  return AddBinaryOp(iArg1, oMinus, iArg2);
sebastien's avatar
sebastien committed
123
124
}

125
126
expr_t
DataTree::AddUMinus(expr_t iArg1)
sebastien's avatar
sebastien committed
127
128
129
130
131
{
  if (iArg1 != Zero)
    {
      // Simplify -(-x) in x
      UnaryOpNode *uarg = dynamic_cast<UnaryOpNode *>(iArg1);
sebastien's avatar
sebastien committed
132
133
      if (uarg != NULL && uarg->get_op_code() == oUminus)
        return uarg->get_arg();
sebastien's avatar
sebastien committed
134
135
136
137
138
139
140

      return AddUnaryOp(oUminus, iArg1);
    }
  else
    return Zero;
}

141
142
expr_t
DataTree::AddTimes(expr_t iArg1, expr_t iArg2)
sebastien's avatar
sebastien committed
143
144
145
146
147
148
149
150
151
152
153
{
  if (iArg1 == MinusOne)
    return AddUMinus(iArg2);
  else if (iArg2 == MinusOne)
    return AddUMinus(iArg1);
  else if (iArg1 != Zero && iArg1 != One && iArg2 != Zero && iArg2 != One)
    {
      // To treat commutativity of "*"
      // Nodes iArg1 and iArg2 are sorted by index
      if (iArg1->idx > iArg2->idx)
        {
154
          expr_t tmp = iArg1;
sebastien's avatar
sebastien committed
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
          iArg1 = iArg2;
          iArg2 = tmp;
        }
      return AddBinaryOp(iArg1, oTimes, iArg2);
    }
  else if (iArg1 != Zero && iArg1 != One && iArg2 == One)
    return iArg1;
  else if (iArg2 != Zero && iArg2 != One && iArg1 == One)
    return iArg2;
  else if (iArg2 == One && iArg1 == One)
    return One;
  else
    return Zero;
}

170
171
expr_t
DataTree::AddDivide(expr_t iArg1, expr_t iArg2)
sebastien's avatar
sebastien committed
172
{
173
  if (iArg2 == One)
sebastien's avatar
sebastien committed
174
    return iArg1;
175
176
177

  // This test should be before the next two, otherwise 0/0 won't be rejected
  if (iArg2 == Zero)
178
    {
179
      cerr << "ERROR: Division by zero!" << endl;
180
      exit(EXIT_FAILURE);
181
    }
182
183
184
185
186
187
188
189

  if (iArg1 == Zero)
    return Zero;

  if (iArg1 == iArg2)
    return One;

  return AddBinaryOp(iArg1, oDivide, iArg2);
190
}
sebastien's avatar
sebastien committed
191

192
193
expr_t
DataTree::AddLess(expr_t iArg1, expr_t iArg2)
ferhat's avatar
ferhat committed
194
195
196
197
{
  return AddBinaryOp(iArg1, oLess, iArg2);
}

198
199
expr_t
DataTree::AddGreater(expr_t iArg1, expr_t iArg2)
ferhat's avatar
ferhat committed
200
201
202
203
{
  return AddBinaryOp(iArg1, oGreater, iArg2);
}

204
205
expr_t
DataTree::AddLessEqual(expr_t iArg1, expr_t iArg2)
ferhat's avatar
ferhat committed
206
207
208
209
{
  return AddBinaryOp(iArg1, oLessEqual, iArg2);
}

210
211
expr_t
DataTree::AddGreaterEqual(expr_t iArg1, expr_t iArg2)
ferhat's avatar
ferhat committed
212
213
214
215
{
  return AddBinaryOp(iArg1, oGreaterEqual, iArg2);
}

216
217
expr_t
DataTree::AddEqualEqual(expr_t iArg1, expr_t iArg2)
ferhat's avatar
ferhat committed
218
219
220
221
{
  return AddBinaryOp(iArg1, oEqualEqual, iArg2);
}

222
223
expr_t
DataTree::AddDifferent(expr_t iArg1, expr_t iArg2)
ferhat's avatar
ferhat committed
224
225
226
227
{
  return AddBinaryOp(iArg1, oDifferent, iArg2);
}

228
229
expr_t
DataTree::AddPower(expr_t iArg1, expr_t iArg2)
sebastien's avatar
sebastien committed
230
{
231
  if (iArg1 != Zero && iArg2 != Zero && iArg1 != One && iArg2 != One)
sebastien's avatar
sebastien committed
232
    return AddBinaryOp(iArg1, oPower, iArg2);
233
234
  else if (iArg1 == One)
    return One;
sebastien's avatar
sebastien committed
235
236
237
238
239
240
241
242
  else if (iArg2 == One)
    return iArg1;
  else if (iArg2 == Zero)
    return One;
  else
    return Zero;
}

243
244
245
246
247
248
249
expr_t
DataTree::AddPowerDeriv(expr_t iArg1, expr_t iArg2, int powerDerivOrder)
{
  assert(powerDerivOrder > 0);
  return AddBinaryOp(iArg1, oPowerDeriv, iArg2, powerDerivOrder);
}

250
251
expr_t
DataTree::AddExp(expr_t iArg1)
sebastien's avatar
sebastien committed
252
253
254
255
256
257
258
{
  if (iArg1 != Zero)
    return AddUnaryOp(oExp, iArg1);
  else
    return One;
}

259
260
expr_t
DataTree::AddLog(expr_t iArg1)
sebastien's avatar
sebastien committed
261
262
263
264
265
266
267
{
  if (iArg1 != Zero && iArg1 != One)
    return AddUnaryOp(oLog, iArg1);
  else if (iArg1 == One)
    return Zero;
  else
    {
268
      cerr << "ERROR: log(0) not defined!" << endl;
269
      exit(EXIT_FAILURE);
sebastien's avatar
sebastien committed
270
271
272
    }
}

273
274
expr_t
DataTree::AddLog10(expr_t iArg1)
sebastien's avatar
sebastien committed
275
276
277
278
279
280
281
{
  if (iArg1 != Zero && iArg1 != One)
    return AddUnaryOp(oLog10, iArg1);
  else if (iArg1 == One)
    return Zero;
  else
    {
282
      cerr << "ERROR: log10(0) not defined!" << endl;
283
      exit(EXIT_FAILURE);
sebastien's avatar
sebastien committed
284
285
286
    }
}

287
288
expr_t
DataTree::AddCos(expr_t iArg1)
sebastien's avatar
sebastien committed
289
290
291
292
293
294
295
{
  if (iArg1 != Zero)
    return AddUnaryOp(oCos, iArg1);
  else
    return One;
}

296
297
expr_t
DataTree::AddSin(expr_t iArg1)
sebastien's avatar
sebastien committed
298
299
300
301
302
303
304
{
  if (iArg1 != Zero)
    return AddUnaryOp(oSin, iArg1);
  else
    return Zero;
}

305
306
expr_t
DataTree::AddTan(expr_t iArg1)
sebastien's avatar
sebastien committed
307
308
309
310
311
312
313
{
  if (iArg1 != Zero)
    return AddUnaryOp(oTan, iArg1);
  else
    return Zero;
}

314
315
expr_t
DataTree::AddAcos(expr_t iArg1)
sebastien's avatar
sebastien committed
316
317
318
319
320
321
322
{
  if (iArg1 != One)
    return AddUnaryOp(oAcos, iArg1);
  else
    return Zero;
}

323
324
expr_t
DataTree::AddAsin(expr_t iArg1)
sebastien's avatar
sebastien committed
325
326
327
328
329
330
331
{
  if (iArg1 != Zero)
    return AddUnaryOp(oAsin, iArg1);
  else
    return Zero;
}

332
333
expr_t
DataTree::AddAtan(expr_t iArg1)
sebastien's avatar
sebastien committed
334
335
336
337
338
339
340
{
  if (iArg1 != Zero)
    return AddUnaryOp(oAtan, iArg1);
  else
    return Zero;
}

341
342
expr_t
DataTree::AddCosh(expr_t iArg1)
sebastien's avatar
sebastien committed
343
344
345
346
347
348
349
{
  if (iArg1 != Zero)
    return AddUnaryOp(oCosh, iArg1);
  else
    return One;
}

350
351
expr_t
DataTree::AddSinh(expr_t iArg1)
sebastien's avatar
sebastien committed
352
353
354
355
356
357
358
{
  if (iArg1 != Zero)
    return AddUnaryOp(oSinh, iArg1);
  else
    return Zero;
}

359
360
expr_t
DataTree::AddTanh(expr_t iArg1)
sebastien's avatar
sebastien committed
361
362
363
364
365
366
367
{
  if (iArg1 != Zero)
    return AddUnaryOp(oTanh, iArg1);
  else
    return Zero;
}

368
369
expr_t
DataTree::AddAcosh(expr_t iArg1)
sebastien's avatar
sebastien committed
370
371
372
373
374
375
376
{
  if (iArg1 != One)
    return AddUnaryOp(oAcosh, iArg1);
  else
    return Zero;
}

377
378
expr_t
DataTree::AddAsinh(expr_t iArg1)
sebastien's avatar
sebastien committed
379
380
381
382
383
384
385
{
  if (iArg1 != Zero)
    return AddUnaryOp(oAsinh, iArg1);
  else
    return Zero;
}

386
387
expr_t
DataTree::AddAtanh(expr_t iArg1)
sebastien's avatar
sebastien committed
388
389
390
391
392
393
394
{
  if (iArg1 != Zero)
    return AddUnaryOp(oAtanh, iArg1);
  else
    return Zero;
}

395
396
expr_t
DataTree::AddSqrt(expr_t iArg1)
sebastien's avatar
sebastien committed
397
398
399
400
401
402
403
{
  if (iArg1 != Zero)
    return AddUnaryOp(oSqrt, iArg1);
  else
    return Zero;
}

404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
expr_t
DataTree::AddAbs(expr_t iArg1)
{
  if (iArg1 == Zero)
    return Zero;
  if (iArg1 == One)
    return One;
  else
    return AddUnaryOp(oAbs, iArg1);
}

expr_t
DataTree::AddSign(expr_t iArg1)
{
  if (iArg1 == Zero)
    return Zero;
  if (iArg1 == One)
    return One;
  else
    return AddUnaryOp(oSign, iArg1);
}

426
427
expr_t
DataTree::AddErf(expr_t iArg1)
428
429
430
431
432
433
434
{
  if (iArg1 != Zero)
    return AddUnaryOp(oErf, iArg1);
  else
    return Zero;
}

435
436
expr_t
DataTree::AddMax(expr_t iArg1, expr_t iArg2)
437
{
438
  return AddBinaryOp(iArg1, oMax, iArg2);
439
440
}

441
442
expr_t
DataTree::AddMin(expr_t iArg1, expr_t iArg2)
443
{
444
  return AddBinaryOp(iArg1, oMin, iArg2);
445
446
}

447
448
expr_t
DataTree::AddNormcdf(expr_t iArg1, expr_t iArg2, expr_t iArg3)
449
450
{
  return AddTrinaryOp(iArg1, oNormcdf, iArg2, iArg3);
451
452
}

453
454
expr_t
DataTree::AddNormpdf(expr_t iArg1, expr_t iArg2, expr_t iArg3)
455
456
{
  return AddTrinaryOp(iArg1, oNormpdf, iArg2, iArg3);
457
458
}

459
460
expr_t
DataTree::AddSteadyState(expr_t iArg1)
461
462
463
464
{
  return AddUnaryOp(oSteadyState, iArg1);
}

465
466
467
468
469
470
471
472
473
474
475
476
expr_t
DataTree::AddSteadyStateParamDeriv(expr_t iArg1, int param_symb_id)
{
  return AddUnaryOp(oSteadyStateParamDeriv, iArg1, 0, "", param_symb_id);
}

expr_t
DataTree::AddSteadyStateParam2ndDeriv(expr_t iArg1, int param1_symb_id, int param2_symb_id)
{
  return AddUnaryOp(oSteadyStateParam2ndDeriv, iArg1, 0, "", param1_symb_id, param2_symb_id);
}

477
478
expr_t
DataTree::AddExpectation(int iArg1, expr_t iArg2)
479
480
481
482
{
  return AddUnaryOp(oExpectation, iArg2, iArg1);
}

483
484
expr_t
DataTree::AddExpectation(string *iArg1, expr_t iArg2)
485
486
487
488
{
  return AddUnaryOp(oExpectation, iArg2, 0, *iArg1);
}

489
490
expr_t
DataTree::AddEqual(expr_t iArg1, expr_t iArg2)
sebastien's avatar
sebastien committed
491
492
493
494
495
{
  return AddBinaryOp(iArg1, oEqual, iArg2);
}

void
496
DataTree::AddLocalVariable(int symb_id, expr_t value) throw (LocalVariableException)
sebastien's avatar
sebastien committed
497
{
498
  assert(symbol_table.getType(symb_id) == eModelLocalVariable);
sebastien's avatar
sebastien committed
499

sebastien's avatar
sebastien committed
500
  // Throw an exception if symbol already declared
501
  map<int, expr_t>::iterator it = local_variables_table.find(symb_id);
sebastien's avatar
sebastien committed
502
  if (it != local_variables_table.end())
503
    throw LocalVariableException(symbol_table.getName(symb_id));
sebastien's avatar
sebastien committed
504

505
  local_variables_table[symb_id] = value;
sebastien's avatar
sebastien committed
506
}
sebastien's avatar
sebastien committed
507

508
509
expr_t
DataTree::AddExternalFunction(int symb_id, const vector<expr_t> &arguments)
510
511
{
  assert(symbol_table.getType(symb_id) == eExternalFunction);
512

513
  external_function_node_map_t::iterator it = external_function_node_map.find(make_pair(arguments, symb_id));
514
515
516
  if (it != external_function_node_map.end())
    return it->second;

517
518
  return new ExternalFunctionNode(*this, symb_id, arguments);
}
sebastien's avatar
sebastien committed
519

520
521
expr_t
DataTree::AddFirstDerivExternalFunctionNode(int top_level_symb_id, const vector<expr_t> &arguments, int input_index)
522
523
{
  assert(symbol_table.getType(top_level_symb_id) == eExternalFunction);
524

525
526
527
  first_deriv_external_function_node_map_t::iterator it
    = first_deriv_external_function_node_map.find(make_pair(make_pair(arguments, input_index),
                                                            top_level_symb_id));
528
529
530
  if (it != first_deriv_external_function_node_map.end())
    return it->second;

531
532
  return new FirstDerivExternalFunctionNode(*this, top_level_symb_id, arguments, input_index);
}
533

534
535
expr_t
DataTree::AddSecondDerivExternalFunctionNode(int top_level_symb_id, const vector<expr_t> &arguments, int input_index1, int input_index2)
536
537
{
  assert(symbol_table.getType(top_level_symb_id) == eExternalFunction);
538

539
540
541
542
  second_deriv_external_function_node_map_t::iterator it
    = second_deriv_external_function_node_map.find(make_pair(make_pair(arguments,
                                                                       make_pair(input_index1, input_index2)),
                                                             top_level_symb_id));
543
544
545
  if (it != second_deriv_external_function_node_map.end())
    return it->second;

546
  return new SecondDerivExternalFunctionNode(*this, top_level_symb_id, arguments, input_index1, input_index2);
sebastien's avatar
sebastien committed
547
}
sebastien's avatar
sebastien committed
548

549
550
551
bool
DataTree::isSymbolUsed(int symb_id) const
{
552
  for (variable_node_map_t::const_iterator it = variable_node_map.begin();
553
       it != variable_node_map.end(); it++)
554
555
556
557
558
559
560
561
562
    if (it->first.first == symb_id)
      return true;

  if (local_variables_table.find(symb_id) != local_variables_table.end())
    return true;

  return false;
}

563
564
565
566
567
568
int
DataTree::getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException)
{
  throw UnknownDerivIDException();
}

569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
SymbolType
DataTree::getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException)
{
  throw UnknownDerivIDException();
}

int
DataTree::getLagByDerivID(int deriv_id) const throw (UnknownDerivIDException)
{
  throw UnknownDerivIDException();
}

int
DataTree::getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException)
{
  throw UnknownDerivIDException();
}

void
DataTree::addAllParamDerivId(set<int> &deriv_id_set)
{
}

592
593
594
595
596
int
DataTree::getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDException)
{
  throw UnknownDerivIDException();
}
sebastien's avatar
sebastien committed
597
598

bool
Sébastien Villemot's avatar
Sébastien Villemot committed
599
DataTree::isUnaryOpUsed(UnaryOpcode opcode) const
sebastien's avatar
sebastien committed
600
{
601
  for (unary_op_node_map_t::const_iterator it = unary_op_node_map.begin();
602
       it != unary_op_node_map.end(); it++)
603
    if (it->first.first.second == opcode)
Sébastien Villemot's avatar
Sébastien Villemot committed
604
605
606
607
608
609
610
611
      return true;

  return false;
}

bool
DataTree::isBinaryOpUsed(BinaryOpcode opcode) const
{
612
  for (binary_op_node_map_t::const_iterator it = binary_op_node_map.begin();
Sébastien Villemot's avatar
Sébastien Villemot committed
613
614
615
616
617
618
619
620
621
622
       it != binary_op_node_map.end(); it++)
    if (it->first.second == opcode)
      return true;

  return false;
}

bool
DataTree::isTrinaryOpUsed(TrinaryOpcode opcode) const
{
623
  for (trinary_op_node_map_t::const_iterator it = trinary_op_node_map.begin();
Sébastien Villemot's avatar
Sébastien Villemot committed
624
625
       it != trinary_op_node_map.end(); it++)
    if (it->first.second == opcode)
sebastien's avatar
sebastien committed
626
627
628
629
      return true;

  return false;
}
630
631
632
633

bool
DataTree::isExternalFunctionUsed(int symb_id) const
{
634
  for (external_function_node_map_t::const_iterator it = external_function_node_map.begin();
635
636
637
638
639
640
641
642
643
644
       it != external_function_node_map.end(); it++)
    if (it->first.second == symb_id)
      return true;

  return false;
}

bool
DataTree::isFirstDerivExternalFunctionUsed(int symb_id) const
{
645
  for (first_deriv_external_function_node_map_t::const_iterator it = first_deriv_external_function_node_map.begin();
646
647
648
649
650
651
652
653
654
655
       it != first_deriv_external_function_node_map.end(); it++)
    if (it->first.second == symb_id)
      return true;

  return false;
}

bool
DataTree::isSecondDerivExternalFunctionUsed(int symb_id) const
{
656
  for (second_deriv_external_function_node_map_t::const_iterator it = second_deriv_external_function_node_map.begin();
657
658
659
660
661
662
       it != second_deriv_external_function_node_map.end(); it++)
    if (it->first.second == symb_id)
      return true;

  return false;
}
Sébastien Villemot's avatar
Sébastien Villemot committed
663
664
665
666
667
668
669
670
671
672
673

int
DataTree::minLagForSymbol(int symb_id) const
{
  int r = 0;
  for (variable_node_map_t::const_iterator it = variable_node_map.begin();
       it != variable_node_map.end(); ++it)
    if (it->first.first == symb_id && it->first.second < r)
      r = it->first.second;
  return r;
}
674
675
676
677
678
679
680
681
682
683
684

void
DataTree::writePowerDerivCHeader(ostream &output) const
{
  if (isBinaryOpUsed(oPowerDeriv))
    output << "double getPowerDeriv(double, double, int);" << endl;
}

void
DataTree::writePowerDeriv(ostream &output, bool use_dll) const
{
685
  if (use_dll && isBinaryOpUsed(oPowerDeriv))
686
687
688
689
690
    output << "/*" << endl
           << " * The k-th derivative of x^p" << endl
           << " */" << endl
           << "double getPowerDeriv(double x, double p, int k)" << endl
           << "{" << endl
691
692
693
           << "#ifdef _MSC_VER" << endl
           << "# define nearbyint(x) (fabs((x)-floor(x)) < fabs((x)-ceil(x)) ? floor(x) : ceil(x))" << endl
           << "#endif" << endl
694
695
696
697
698
699
700
701
702
703
704
705
           << "  if ( fabs(x) < " << NEAR_ZERO << " && p > 0 && k >= p && fabs(p-nearbyint(p)) < " << NEAR_ZERO << " )" << endl
           << "    return 0.0;" << endl
           << "  else" << endl
           << "    {" << endl
           << "      int i = 0;" << endl
           << "      double dxp = pow(x, p-k);" << endl
           << "      for (; i<k; i++)" << endl
           << "        dxp *= p--;" << endl
           << "      return dxp;" << endl
           << "    }" << endl
           << "}" << endl;
}