LORENE
cmp_arithm.C
1 /*
2  * Arithmetical operations for class Cmp
3  *
4  */
5 
6 /*
7  * Copyright (c) 1999-2000 Jean-Alain Marck
8  * Copyright (c) 1999-2001 Eric Gourgoulhon
9  * Copyright (c) 1999-2001 Philippe Grandclement
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 char cmp_arithm_C[] = "$Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_arithm.C,v 1.5 2014/10/13 08:52:47 j_novak Exp $" ;
31 
32 /*
33  * $Id: cmp_arithm.C,v 1.5 2014/10/13 08:52:47 j_novak Exp $
34  * $Log: cmp_arithm.C,v $
35  * Revision 1.5 2014/10/13 08:52:47 j_novak
36  * Lorene classes and functions now belong to the namespace Lorene.
37  *
38  * Revision 1.4 2014/10/06 15:13:03 j_novak
39  * Modified #include directives to use c++ syntax.
40  *
41  * Revision 1.3 2003/06/20 13:59:59 f_limousin
42  * L'assert pour le mapping est realise a partir du mapping lui meme et non a partir du pointeur sur le mapping.
43  *
44  * Revision 1.2 2002/10/16 14:36:34 j_novak
45  * Reorganization of #include instructions of standard C++, in order to
46  * use experimental version 3 of gcc.
47  *
48  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
49  * LORENE
50  *
51  * Revision 2.11 2001/05/26 15:07:47 eric
52  * Ajout de operator% : multiplication de deux Cmp avec desaliasage
53  *
54  * Revision 2.10 2000/02/15 13:39:57 phil
55  * correction -= et +=
56  *
57  * Revision 2.9 2000/01/28 16:34:57 eric
58  * Utilisation des nouvelles fonctions Cmp::check_dzpuis et Cmp::dz_nonzero
59  * dans les tests sur dzpuis.
60  *
61  * Revision 2.8 1999/12/10 16:32:52 eric
62  * Dans l'arithmetique membre (+=, -=, *=), on n'appelle desormais
63  * del_deriv() que tout a la fin.
64  *
65  * Revision 2.7 1999/11/26 14:23:38 eric
66  * Ajout du membre dzpuis.
67  *
68  * Revision 2.6 1999/11/12 17:08:35 eric
69  * Ajout de la partie manquante de l'arithmetique.
70  *
71  * Revision 2.5 1999/10/28 13:23:43 phil
72  * Deverouillage des ETATNONDEF
73  *
74  * Revision 2.4 1999/10/27 08:45:21 eric
75  * ntroduction du membre Valeur va.
76  *
77  * Revision 2.3 1999/10/22 08:14:58 eric
78  * La fonction annule() est rebaptisee annule_hard().
79  *
80  * Revision 2.2 1999/09/14 17:19:44 phil
81  * ajout de Cmp operator* (double, const Cmp&)
82  *
83  * Revision 2.1 1999/03/03 11:13:56 hyc
84  * *** empty log message ***
85  *
86  *
87  * $Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_arithm.C,v 1.5 2014/10/13 08:52:47 j_novak Exp $
88  *
89  */
90 
91 // headers C
92 #include <cassert>
93 #include <cstdlib>
94 
95 // headers Lorene
96 #include "cmp.h"
97 #include "type_parite.h"
98 
99  //********************//
100  // OPERATEURS UNAIRES //
101  //********************//
102 
103 namespace Lorene {
104 Cmp operator+(const Cmp & ci) {
105  return ci ;
106 }
107 
108 Cmp operator-(const Cmp & ci) {
109 
110  // Cas particulier
111  if ((ci.get_etat() == ETATZERO) || (ci.get_etat() == ETATNONDEF)) {
112  return ci ;
113  }
114 
115  // Cas general
116  assert(ci.get_etat() == ETATQCQ) ; // sinon...
117  Cmp r(ci.get_mp()) ; // Cmp resultat
118  r.set_etat_qcq() ;
119  r.va = - ci.va ;
120  r.set_dzpuis( ci.get_dzpuis() ) ;
121 
122  // Termine
123  return r ;
124 }
125 
126  //**********//
127  // ADDITION //
128  //**********//
129 // Cmp + Cmp
130 // ---------
131 Cmp operator+(const Cmp & c1, const Cmp & c2) {
132 
133  if (c1.get_etat() == ETATNONDEF)
134  return c1 ;
135  if (c2.get_etat() == ETATNONDEF)
136  return c2 ;
137  assert(c1.get_mp() == c2.get_mp()) ;
138 
139  // Cas particuliers
140  if (c1.get_etat() == ETATZERO) {
141  return c2 ;
142  }
143  if (c2.get_etat() == ETATZERO) {
144  return c1 ;
145  }
146  assert(c1.get_etat() == ETATQCQ) ; // sinon...
147  assert(c2.get_etat() == ETATQCQ) ; // sinon...
148 
149  // Cas general
150 
151  if ( c1.dz_nonzero() && c2.dz_nonzero() ) {
152  if ( c1.get_dzpuis() != c2.get_dzpuis() ) {
153  cout << "Operation Cmp + Cmp forbidden in the external " << endl;
154  cout << " compactified domain ! " << endl ;
155  abort() ;
156  }
157  }
158 
159  Cmp r(c1) ; // Le resultat
160  r.va += c2.va ;
161 
162  if (c1.dz_nonzero()) {
163  r.set_dzpuis( c1.get_dzpuis() ) ;
164  }
165  else{
166  r.set_dzpuis( c2.get_dzpuis() ) ;
167  }
168 
169  // Termine
170  return r ;
171 }
172 
173 // Cmp + double
174 // ------------
175 Cmp operator+(const Cmp& t1, double x)
176 {
177  // Protections
178  assert(t1.get_etat() != ETATNONDEF) ;
179 
180  // Cas particuliers
181  if (x == double(0)) {
182  return t1 ;
183  }
184 
185  assert( t1.check_dzpuis(0) ) ;
186 
187  Cmp resu(t1) ;
188 
189  if (t1.get_etat() == ETATZERO) {
190  resu = x ;
191  }
192  else{
193  assert(resu.get_etat() == ETATQCQ) ; // sinon ...
194  resu.va = resu.va + x ;
195  }
196 
197  resu.set_dzpuis(0) ;
198 
199  return resu ;
200 }
201 
202 // double + Cmp
203 // ------------
204 Cmp operator+(double x, const Cmp& t1)
205 {
206  return t1 + x ;
207 }
208 
209 // Cmp + int
210 // ---------
211 Cmp operator+(const Cmp& t1, int m)
212 {
213  return t1 + double(m) ;
214 }
215 
216 // int + Cmp
217 // ---------
218 Cmp operator+(int m, const Cmp& t1)
219 {
220  return t1 + double(m) ;
221 }
222 
223 
224 
225 
226 
227  //**************//
228  // SOUSTRACTION //
229  //**************//
230 
231 // Cmp - Cmp
232 // ---------
233 Cmp operator-(const Cmp & c1, const Cmp & c2) {
234 
235  if (c1.get_etat() == ETATNONDEF)
236  return c1 ;
237  if (c2.get_etat() == ETATNONDEF)
238  return c2 ;
239 
240  assert(c1.get_mp() == c2.get_mp()) ;
241 
242  // Cas particuliers
243  if (c1.get_etat() == ETATZERO) {
244  return -c2 ;
245  }
246  if (c2.get_etat() == ETATZERO) {
247  return c1 ;
248  }
249  assert(c1.get_etat() == ETATQCQ) ; // sinon...
250  assert(c2.get_etat() == ETATQCQ) ; // sinon...
251 
252  // Cas general
253  if ( c1.dz_nonzero() && c2.dz_nonzero() ) {
254  if ( c1.get_dzpuis() != c2.get_dzpuis() ) {
255  cout << "Operation Cmp - Cmp forbidden in the external " << endl;
256  cout << " compactified domain ! " << endl ;
257  abort() ;
258  }
259  }
260 
261  Cmp r(c1) ; // Le resultat
262  r.va -= c2.va ;
263 
264  if (c1.dz_nonzero()) {
265  r.set_dzpuis( c1.get_dzpuis() ) ;
266  }
267  else{
268  r.set_dzpuis( c2.get_dzpuis() ) ;
269  }
270 
271  // Termine
272  return r ;
273 }
274 
275 // Cmp - double
276 // ------------
277 Cmp operator-(const Cmp& t1, double x)
278 {
279  // Protections
280  assert(t1.get_etat() != ETATNONDEF) ;
281 
282  // Cas particuliers
283  if (x == double(0)) {
284  return t1 ;
285  }
286 
287  assert( t1.check_dzpuis(0) ) ;
288 
289  Cmp resu(t1) ;
290 
291  if (t1.get_etat() == ETATZERO) {
292  resu = - x ;
293  }
294  else{
295  assert(resu.get_etat() == ETATQCQ) ; // sinon ...
296  resu.va = resu.va - x ;
297  }
298 
299  resu.set_dzpuis(0) ;
300 
301  return resu ;
302 }
303 
304 // double - Cmp
305 // ------------
306 Cmp operator-(double x, const Cmp& t1)
307 {
308  return - (t1 - x) ;
309 }
310 
311 // Cmp - int
312 // ---------
313 Cmp operator-(const Cmp& t1, int m)
314 {
315  return t1 - double(m) ;
316 }
317 
318 // int - Cmp
319 // ---------
320 Cmp operator-(int m, const Cmp& t1)
321 {
322  return double(m) - t1 ;
323 }
324 
325 
326 
327 
328 
329 
330  //****************//
331  // MULTIPLICATION //
332  //****************//
333 
334 // Cmp * Cmp
335 // ---------
336 Cmp operator*(const Cmp& c1, const Cmp& c2) {
337 
338 
339  // Cas particuliers
340  if ((c1.get_etat() == ETATZERO) || (c1.get_etat() == ETATNONDEF)){
341  return c1 ;
342  }
343  if ((c2.get_etat() == ETATZERO)|| (c2.get_etat() == ETATNONDEF)) {
344  return c2 ;
345  }
346  assert(c1.get_etat() == ETATQCQ) ; // sinon...
347  assert(c2.get_etat() == ETATQCQ) ; // sinon...
348 
349  // Protection
350  assert(*(c1.get_mp()) == *(c2.get_mp())) ;
351 
352  // Cas general
353  Cmp r(c1) ; // Le resultat
354  r.va *= c2.va ;
355 
356  r.set_dzpuis( c1.get_dzpuis() + c2.get_dzpuis() ) ;
357 
358  // Termine
359  return r ;
360 }
361 
362 // Cmp % Cmp (multiplication with desaliasing)
363 // -------------------------------------------
364 Cmp operator%(const Cmp& c1, const Cmp& c2) {
365 
366 
367  // Cas particuliers
368  if ((c1.get_etat() == ETATZERO) || (c1.get_etat() == ETATNONDEF)){
369  return c1 ;
370  }
371  if ((c2.get_etat() == ETATZERO)|| (c2.get_etat() == ETATNONDEF)) {
372  return c2 ;
373  }
374  assert(c1.get_etat() == ETATQCQ) ; // sinon...
375  assert(c2.get_etat() == ETATQCQ) ; // sinon...
376 
377  // Protection
378  assert(c1.get_mp() == c2.get_mp()) ;
379 
380  // Cas general
381  Cmp r( c1.get_mp() ) ; // Le resultat
382  r.set_etat_qcq() ;
383  r.va = c1.va % c2.va ;
384 
385  r.set_dzpuis( c1.get_dzpuis() + c2.get_dzpuis() ) ;
386 
387  // Termine
388  return r ;
389 }
390 
391 
392 
393 
394 // double * Cmp
395 // ------------
396 Cmp operator*(double a, const Cmp& c1) {
397 
398  // Cas particuliers
399  if ((c1.get_etat() == ETATZERO) || (c1.get_etat() == ETATNONDEF)) {
400  return c1 ;
401  }
402 
403  assert(c1.get_etat() == ETATQCQ) ; // sinon...
404 
405  // Cas general
406  Cmp r(c1.get_mp()) ;
407  r.set_dzpuis( c1.get_dzpuis() ) ;
408 
409  if ( a == double(0) ) {
410  r.set_etat_zero() ;
411  }
412  else {
413  r.set_etat_qcq() ;
414  r.va = a * c1.va ;
415  }
416 
417 
418  // Termine
419  return r ;
420 }
421 
422 
423 // Cmp * double
424 // ------------
425 Cmp operator*(const Cmp& t1, double x)
426 {
427  return x * t1 ;
428 }
429 
430 // Cmp * int
431 // ---------
432 Cmp operator*(const Cmp& t1, int m)
433 {
434  return t1 * double(m) ;
435 }
436 
437 // int * Cmp
438 // ---------
439 Cmp operator*(int m, const Cmp& t1)
440 {
441  return double(m) * t1 ;
442 }
443 
444 
445 
446 
447 
448 
449 
450  //**********//
451  // DIVISION //
452  //**********//
453 
454 
455 // Cmp / Cmp
456 // ---------
457 Cmp operator/(const Cmp& c1, const Cmp& c2) {
458 
459  // Protections
460  assert(c1.get_etat() != ETATNONDEF) ;
461  assert(c2.get_etat() != ETATNONDEF) ;
462  assert(c1.get_mp() == c2.get_mp()) ;
463 
464  // Cas particuliers
465  if (c2.get_etat() == ETATZERO) {
466  cout << "Division by 0 in Cmp / Cmp !" << endl ;
467  abort() ;
468  }
469  if (c1.get_etat() == ETATZERO) {
470  return c1 ;
471  }
472 
473  // Cas general
474 
475  assert(c1.get_etat() == ETATQCQ) ; // sinon...
476  assert(c2.get_etat() == ETATQCQ) ; // sinon...
477 
478  Cmp r(c1.get_mp()) ; // Le resultat
479 
480  r.set_etat_qcq() ;
481  r.va = c1.va / c2.va ;
482 
483  r.set_dzpuis( c1.get_dzpuis() - c2.get_dzpuis() ) ;
484 
485  // Termine
486  return r ;
487 }
488 
489 // Cmp / double
490 // -------------
491 Cmp operator/(const Cmp& c1, double x) {
492 
493  if (c1.get_etat() == ETATNONDEF)
494  return c1 ;
495 
496  // Cas particuliers
497  if ( x == double(0) ) {
498  cout << "Division by 0 in Cmp / double !" << endl ;
499  abort() ;
500  }
501  if (c1.get_etat() == ETATZERO) {
502  return c1 ;
503  }
504 
505  assert(c1.get_etat() == ETATQCQ) ; // sinon...
506 
507  Cmp r(c1.get_mp()) ; // Le resultat
508 
509  r.set_etat_qcq() ;
510  r.va = c1.va / x ;
511 
512  r.set_dzpuis( c1.get_dzpuis() ) ;
513 
514  // Termine
515  return r ;
516 }
517 
518 
519 // double / Cmp
520 // ------------
521 Cmp operator/(double x, const Cmp& c2) {
522 
523  if (c2.get_etat() == ETATNONDEF)
524  return c2 ;
525 
526  if (c2.get_etat() == ETATZERO) {
527  cout << "Division by 0 in Cmp / Cmp !" << endl ;
528  abort() ;
529  }
530 
531 
532  assert(c2.get_etat() == ETATQCQ) ; // sinon...
533 
534  Cmp r(c2.get_mp()) ; // Le resultat
535  r.set_dzpuis( - c2.get_dzpuis() ) ;
536 
537  if ( x == double(0) ) {
538  r.set_etat_zero() ;
539  }
540  else {
541  r.set_etat_qcq() ;
542  r.va = x / c2.va ;
543  }
544 
545  // Termine
546  return r ;
547 }
548 
549 
550 // Cmp / int
551 // ---------
552 Cmp operator/(const Cmp& c1, int m) {
553 
554  return c1 / double(m) ;
555 
556 }
557 
558 
559 // int / Cmp
560 // ---------
561 Cmp operator/(int m, const Cmp& c2) {
562 
563  return double(m) / c2 ;
564 
565 }
566 
567  //*******************//
568  // operateurs +=,... //
569  //*******************//
570 
571 //---------
572 // += Cmp
573 //---------
574 
575 void Cmp::operator+=(const Cmp & ci) {
576 
577  // Protection
578  assert(mp == ci.get_mp()) ; // meme mapping
579  if (etat == ETATNONDEF)
580  return ;
581 
582  // Cas particulier
583  if (ci.get_etat() == ETATZERO) {
584  return ;
585  }
586 
587  if (ci.get_etat() == ETATNONDEF) {
588  set_etat_nondef() ;
589  return ;
590  }
591 
592  // Cas general
593 
594 
595  if ( dz_nonzero() && ci.dz_nonzero() ) {
596  if ( dzpuis != ci.dzpuis ) {
597  cout << "Operation += Cmp forbidden in the external " << endl;
598  cout << " compactified domain ! " << endl ;
599  abort() ;
600  }
601  }
602 
603  if (etat == ETATZERO) {
604  (*this) = ci ;
605  }
606  else {
607  va += ci.va ;
608 
609  if( ci.dz_nonzero() ) {
610  set_dzpuis(ci.dzpuis) ;
611  }
612  }
613  // Menage (a ne faire qu'a la fin seulement)
614  del_deriv() ;
615 
616 
617 }
618 
619 //---------
620 // -= Cmp
621 //---------
622 
623 void Cmp::operator-=(const Cmp & ci) {
624 
625  // Protection
626  assert(mp == ci.get_mp()) ; // meme mapping
627  if (etat == ETATNONDEF)
628  return ;
629 
630  // Cas particulier
631  if (ci.get_etat() == ETATZERO) {
632  return ;
633  }
634 
635  if (ci.get_etat() == ETATNONDEF) {
636  set_etat_nondef() ;
637  return ;
638  }
639 
640  // Cas general
641  if ( dz_nonzero() && ci.dz_nonzero() ) {
642  if ( dzpuis != ci.dzpuis ) {
643  cout << "Operation -= Cmp forbidden in the external " << endl;
644  cout << " compactified domain ! " << endl ;
645  abort() ;
646  }
647  }
648 
649 
650  if (etat == ETATZERO) {
651  (*this) = -ci ;
652  }
653  else {
654  va -= ci.va ;
655 
656  if( ci.dz_nonzero() ) {
657  set_dzpuis(ci.dzpuis) ;
658  }
659  }
660  // Menage (a ne faire qu'a la fin seulement)
661  del_deriv() ;
662 }
663 
664 //---------
665 // *= Cmp
666 //---------
667 
668 void Cmp::operator*=(const Cmp & ci) {
669 
670  // Protection
671  assert(mp == ci.get_mp()) ; // meme mapping
672  if (etat == ETATNONDEF)
673  return ;
674 
675  // Cas particulier
676  if (ci.get_etat() == ETATZERO) {
677  set_etat_zero() ;
678  return ;
679  }
680 
681  if (etat == ETATZERO) {
682  return ;
683  }
684 
685  if (ci.get_etat() == ETATNONDEF) {
686  set_etat_nondef() ;
687  return ;
688  }
689 
690  // Cas general
691 
692  assert(etat == ETATQCQ) ; // sinon....
693 
694  va *= ci.va ;
695 
696  dzpuis += ci.dzpuis ;
697 
698  // Menage (a ne faire qu'a la fin seulement)
699  del_deriv() ;
700 
701 }
702 }
Lorene::Cmp::set_etat_nondef
void set_etat_nondef()
Sets the logical state to ETATNONDEF (undefined).
Definition: cmp.C:297
Lorene::Cmp::check_dzpuis
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition: cmp.C:715
Lorene::Cmp::del_deriv
void del_deriv()
Logical destructor of the derivatives.
Definition: cmp.C:265
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::operator*
Base_val operator*(const Base_val &, const Base_val &)
This operator is used when calling multiplication or division of Valeur .
Definition: base_val_mult.C:102
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp
Definition: cmp.h:464
Lorene::Cmp::get_etat
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
Lorene::Cmp::etat
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition: cmp.h:454
Lorene::Cmp::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: cmp.C:289
Lorene::operator/
Cmp operator/(const Cmp &, const Cmp &)
Cmp / Cmp.
Definition: cmp_arithm.C:457
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Lorene::Cmp::operator+=
void operator+=(const Cmp &)
+= Cmp
Definition: cmp_arithm.C:575
Lorene::operator-
Cmp operator-(const Cmp &)
- Cmp
Definition: cmp_arithm.C:108
Lorene::Cmp::dzpuis
int dzpuis
Power of r by which the quantity represented by this must be divided in the external compactified z...
Definition: cmp.h:461
Lorene::Cmp::operator-=
void operator-=(const Cmp &)
-= Cmp
Definition: cmp_arithm.C:623
Lorene::Cmp::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition: cmp.h:903
Lorene::operator%
Cmp operator%(const Cmp &, const Cmp &)
Cmp * Cmp with desaliasing.
Definition: cmp_arithm.C:364
Lorene::Cmp::operator*=
void operator*=(const Cmp &)
*= Cmp
Definition: cmp_arithm.C:668
Lorene::operator+
Cmp operator+(const Cmp &)
Definition: cmp_arithm.C:104
Lorene::Cmp::mp
const Map * mp
Reference mapping.
Definition: cmp.h:451
Lorene::Cmp::dz_nonzero
bool dz_nonzero() const
Returns true if the last domain is compactified and *this is not zero in this domain.
Definition: cmp.C:660
Lorene::Cmp::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: cmp.C:304
Lorene::Cmp::get_mp
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
Lorene::Cmp::set_dzpuis
void set_dzpuis(int)
Set a value to dzpuis.
Definition: cmp.C:654