LORENE
tbl_math.C
1 /*
2  * Mathematical functions for class Tbl
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 tbl_math_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tbl/tbl_math.C,v 1.4 2014/10/13 08:53:41 j_novak Exp $" ;
31 
32 /*
33  * $Id: tbl_math.C,v 1.4 2014/10/13 08:53:41 j_novak Exp $
34  * $Log: tbl_math.C,v $
35  * Revision 1.4 2014/10/13 08:53:41 j_novak
36  * Lorene classes and functions now belong to the namespace Lorene.
37  *
38  * Revision 1.3 2014/10/06 15:13:18 j_novak
39  * Modified #include directives to use c++ syntax.
40  *
41  * Revision 1.2 2012/01/17 10:38:48 j_penner
42  * added a Heaviside function
43  *
44  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
45  * LORENE
46  *
47  * Revision 2.4 1999/12/02 17:47:48 phil
48  * ajout de racine_cubique
49  *
50  * Revision 2.3 1999/10/29 15:05:32 eric
51  * Ajout de fonctions mathematiques (abs, norme, etc...).
52  *
53  * Revision 2.2 1999/03/01 15:10:19 eric
54  * *** empty log message ***
55  *
56  * Revision 2.1 1999/02/24 15:26:13 hyc
57  * *** empty log message ***
58  *
59  *
60  * $Header: /cvsroot/Lorene/C++/Source/Tbl/tbl_math.C,v 1.4 2014/10/13 08:53:41 j_novak Exp $
61  *
62  */
63 
64 // Headers C
65 // ---------
66 #include <cmath>
67 #include <cstdlib>
68 
69 // Headers Lorene
70 // --------------
71 #include "tbl.h"
72 
73  //-------//
74  // Sinus //
75  //-------//
76 
77 namespace Lorene {
78 Tbl sin(const Tbl& ti)
79 {
80  // Protection
81  assert(ti.get_etat() != ETATNONDEF) ;
82 
83  // Cas ETATZERO
84  if (ti.get_etat() == ETATZERO) {
85  return ti ;
86  }
87 
88  // Cas general
89  assert(ti.get_etat() == ETATQCQ) ; // sinon...
90  Tbl to(ti.dim) ; // Tbl resultat
91  to.set_etat_qcq() ;
92  int taille = ti.get_taille() ;
93  for (int i=0 ; i<taille ; i++) {
94  to.t[i] = sin(ti.t[i]) ;
95  }
96  return to ;
97 }
98 
99  //---------//
100  // Cosinus //
101  //---------//
102 
103 Tbl cos(const Tbl& ti)
104 {
105  // Protection
106  assert(ti.get_etat() != ETATNONDEF) ;
107 
108  Tbl to(ti.dim) ; // Tbl resultat
109  to.set_etat_qcq() ;
110 
111  // Cas ETATZERO
112  if (ti.get_etat() == ETATZERO) {
113  int taille = ti.get_taille() ;
114  for (int i=0 ; i<taille ; i++) {
115  to.t[i] = 1 ;
116  }
117  return to ;
118  }
119 
120  // Cas general
121  assert(ti.get_etat() == ETATQCQ) ; // sinon...
122  int taille = ti.get_taille() ;
123  for (int i=0 ; i<taille ; i++) {
124  to.t[i] = cos(ti.t[i]) ;
125  }
126  return to ;
127 }
128 
129  //----------//
130  // Tangente //
131  //----------//
132 
133 Tbl tan(const Tbl& ti)
134 {
135  // Protection
136  assert(ti.get_etat() != ETATNONDEF) ;
137 
138  // Cas ETATZERO
139  if (ti.get_etat() == ETATZERO) {
140  return ti ;
141  }
142 
143  // Cas general
144  assert(ti.get_etat() == ETATQCQ) ; // sinon...
145  Tbl to(ti.dim) ; // Tbl resultat
146  to.set_etat_qcq() ;
147  int taille = ti.get_taille() ;
148  for (int i=0 ; i<taille ; i++) {
149  to.t[i] = tan(ti.t[i]) ;
150  }
151  return to ;
152 }
153 
154  //----------//
155  // ArcSinus //
156  //----------//
157 
158 Tbl asin(const Tbl& ti)
159 {
160  // Protection
161  assert(ti.get_etat() != ETATNONDEF) ;
162 
163  // Cas ETATZERO
164  if (ti.get_etat() == ETATZERO) {
165  return ti ;
166  }
167 
168  // Cas general
169  assert(ti.get_etat() == ETATQCQ) ; // sinon...
170  Tbl to(ti.dim) ; // Tbl resultat
171  to.set_etat_qcq() ;
172  int taille = ti.get_taille() ;
173  for (int i=0 ; i<taille ; i++) {
174  to.t[i] = asin(ti.t[i]) ;
175  }
176  return to ;
177 }
178 
179  //------------//
180  // ArcCosinus //
181  //------------//
182 
183 Tbl acos(const Tbl& ti)
184 {
185  // Protection
186  assert(ti.get_etat() != ETATNONDEF) ;
187 
188  Tbl to(ti.dim) ; // Tbl resultat
189  to.set_etat_qcq() ;
190 
191  // Cas ETATZERO
192  if (ti.get_etat() == ETATZERO) {
193  int taille = ti.get_taille() ;
194  for (int i=0 ; i<taille ; i++) {
195  to.t[i] = M_PI * .5 ;
196  }
197  return to ;
198  }
199 
200  // Cas general
201  assert(ti.get_etat() == ETATQCQ) ; // sinon...
202  int taille = ti.get_taille() ;
203  for (int i=0 ; i<taille ; i++) {
204  to.t[i] = acos(ti.t[i]) ;
205  }
206  return to ;
207 }
208 
209  //-------------//
210  // ArcTangente //
211  //-------------//
212 
213 Tbl atan(const Tbl& ti)
214 {
215  // Protection
216  assert(ti.get_etat() != ETATNONDEF) ;
217 
218  // Cas ETATZERO
219  if (ti.get_etat() == ETATZERO) {
220  return ti ;
221  }
222 
223  // Cas general
224  assert(ti.get_etat() == ETATQCQ) ; // sinon...
225  Tbl to(ti.dim) ; // Tbl resultat
226  to.set_etat_qcq() ;
227  int taille = ti.get_taille() ;
228  for (int i=0 ; i<taille ; i++) {
229  to.t[i] = atan(ti.t[i]) ;
230  }
231  return to ;
232 }
233 
234  //------//
235  // Sqrt //
236  //------//
237 
238 Tbl sqrt(const Tbl& ti)
239 {
240  // Protection
241  assert(ti.get_etat() != ETATNONDEF) ;
242 
243  // Cas ETATZERO
244  if (ti.get_etat() == ETATZERO) {
245  return ti ;
246  }
247 
248  // Cas general
249  assert(ti.get_etat() == ETATQCQ) ; // sinon...
250  Tbl to(ti.dim) ; // Tbl resultat
251  to.set_etat_qcq() ;
252  int taille = ti.get_taille() ;
253  for (int i=0 ; i<taille ; i++) {
254  to.t[i] = sqrt(ti.t[i]) ;
255  }
256  return to ;
257 }
258 
259  //---------------//
260  // Exponantielle //
261  //---------------//
262 
263 Tbl exp(const Tbl& ti)
264 {
265  // Protection
266  assert(ti.get_etat() != ETATNONDEF) ;
267 
268  Tbl to(ti.dim) ; // Tbl resultat
269  to.set_etat_qcq() ;
270 
271  // Cas ETATZERO
272  if (ti.get_etat() == ETATZERO) {
273  int taille = ti.get_taille() ;
274  for (int i=0 ; i<taille ; i++) {
275  to.t[i] = 1 ;
276  }
277  return to ;
278  }
279 
280  // Cas general
281  assert(ti.get_etat() == ETATQCQ) ; // sinon...
282  int taille = ti.get_taille() ;
283  for (int i=0 ; i<taille ; i++) {
284  to.t[i] = exp(ti.t[i]) ;
285  }
286  return to ;
287 }
288 
289  //--------------------//
290  // Heaviside Function //
291  //--------------------//
292 
293 Tbl Heaviside(const Tbl& ti)
294 {
295  // Protection
296  assert(ti.get_etat() != ETATNONDEF) ;
297 
298  Tbl to(ti.dim) ; // Tbl resultat
299  to.set_etat_qcq() ;
300 
301  // Cas ETATZERO
302  if (ti.get_etat() == ETATZERO) {
303  int taille = ti.get_taille() ;
304  for (int i=0 ; i<taille ; i++) {
305  to.t[i] = 0 ;
306  }
307  return to ;
308  }
309 
310  // Cas general
311  assert(ti.get_etat() == ETATQCQ) ; // Otherwise
312  int taille = ti.get_taille() ;
313  for (int i=0 ; i<taille ; i++) {
314  if(ti.t[i] >= 0)
315  to.t[i] = 1 ;
316  else
317  to.t[i] = 0 ;
318  }
319  return to ;
320 }
321 
322  //-------------//
323  // Log naturel //
324  //-------------//
325 
326 Tbl log(const Tbl& ti)
327 {
328  // Protection
329  assert(ti.get_etat() != ETATNONDEF) ;
330 
331  // Cas ETATZERO
332  if (ti.get_etat() == ETATZERO) {
333  cout << "Tbl log: log(ETATZERO) !" << endl ;
334  abort () ;
335  }
336 
337  // Cas general
338  assert(ti.get_etat() == ETATQCQ) ; // sinon...
339  Tbl to(ti.dim) ; // Tbl resultat
340  to.set_etat_qcq() ;
341  int taille = ti.get_taille() ;
342  for (int i=0 ; i<taille ; i++) {
343  to.t[i] = log(ti.t[i]) ;
344  }
345  return to ;
346 }
347 
348  //-------------//
349  // Log decimal //
350  //-------------//
351 
352 Tbl log10(const Tbl& ti)
353 {
354  // Protection
355  assert(ti.get_etat() != ETATNONDEF) ;
356 
357  // Cas ETATZERO
358  if (ti.get_etat() == ETATZERO) {
359  cout << "Tbl log10: log10(ETATZERO) !" << endl ;
360  abort () ;
361  }
362 
363  // Cas general
364  assert(ti.get_etat() == ETATQCQ) ; // sinon...
365  Tbl to(ti.dim) ; // Tbl resultat
366  to.set_etat_qcq() ;
367  int taille = ti.get_taille() ;
368  for (int i=0 ; i<taille ; i++) {
369  to.t[i] = log10(ti.t[i]) ;
370  }
371  return to ;
372 }
373 
374  //--------------//
375  // Power entier //
376  //--------------//
377 
378 Tbl pow(const Tbl& ti, int n)
379 {
380  // Protection
381  assert(ti.get_etat() != ETATNONDEF) ;
382 
383  // Cas ETATZERO
384  if (ti.get_etat() == ETATZERO) {
385  if (n > 0) {
386  return ti ;
387  }
388  else {
389  cout << "Tbl pow: ETATZERO^n avec n<=0 ! "<< endl ;
390  abort () ;
391  }
392  }
393 
394  // Cas general
395  assert(ti.get_etat() == ETATQCQ) ; // sinon...
396  Tbl to(ti.dim) ; // Tbl resultat
397  to.set_etat_qcq() ;
398  double x = n ;
399  int taille = ti.get_taille() ;
400  for (int i=0 ; i<taille ; i++) {
401  to.t[i] = pow(ti.t[i], x) ;
402  }
403  return to ;
404 }
405 
406  //--------------//
407  // Power double //
408  //--------------//
409 
410 Tbl pow(const Tbl& ti, double x)
411 {
412  // Protection
413  assert(ti.get_etat() != ETATNONDEF) ;
414 
415  // Cas ETATZERO
416  if (ti.get_etat() == ETATZERO) {
417  if (x > 0) {
418  return ti ;
419  }
420  else {
421  cout << "Tbl pow: ETATZERO^x avec x<=0 !" << endl ;
422  abort () ;
423  }
424  }
425 
426  // Cas general
427  assert(ti.get_etat() == ETATQCQ) ; // sinon...
428  Tbl to(ti.dim) ; // Tbl resultat
429  to.set_etat_qcq() ;
430  int taille = ti.get_taille() ;
431  for (int i=0 ; i<taille ; i++) {
432  to.t[i] = pow(ti.t[i], x) ;
433  }
434  return to ;
435 }
436 
437  //----------------//
438  // Absolute value //
439  //----------------//
440 
441 Tbl abs(const Tbl& ti)
442 {
443  // Protection
444  assert(ti.get_etat() != ETATNONDEF) ;
445 
446  // Cas ETATZERO
447  if (ti.get_etat() == ETATZERO) {
448  return ti ;
449  }
450 
451  // Cas general
452  assert(ti.get_etat() == ETATQCQ) ; // sinon...
453 
454  Tbl to(ti.dim) ; // Tbl resultat
455  to.set_etat_qcq() ;
456 
457  const double* xi = ti.t ;
458  double* xo = to.t ;
459  int taille = ti.get_taille() ;
460 
461  for (int i=0 ; i<taille ; i++) {
462  xo[i] = fabs( xi[i] ) ;
463  }
464 
465  return to ;
466 }
467  //----------------//
468  // Cubic //
469  //----------------//
470 
472 {
473  // Protection
474  assert(ti.get_etat() != ETATNONDEF) ;
475 
476  // Cas ETATZERO
477  if (ti.get_etat() == ETATZERO) {
478  return ti ;
479  }
480 
481  // Cas general
482  assert(ti.get_etat() == ETATQCQ) ; // sinon...
483 
484  Tbl absolute(abs(ti)) ;
485  Tbl res (pow(absolute, 1./3.)) ;
486 
487  for (int i=0 ; i<ti.get_taille() ; i++)
488  if (ti.t[i] < 0)
489  res.t[i] *= -1 ;
490 
491  return res ;
492 }
493 
494  //-------------------------------//
495  // max //
496  //-------------------------------//
497 
498 double max(const Tbl& ti) {
499 
500  // Protection
501  assert(ti.get_etat() != ETATNONDEF) ;
502 
503  // Cas particulier
504  if (ti.get_etat() == ETATZERO) {
505  return double(0) ;
506  }
507 
508  // Cas general
509  assert(ti.get_etat() == ETATQCQ) ; // sinon....
510 
511  const double* x = ti.t ;
512  double resu = x[0] ;
513  for (int i=1; i<ti.get_taille(); i++) {
514  if ( x[i] > resu ) resu = x[i] ;
515  }
516 
517  return resu ;
518 }
519 
520  //-------------------------------//
521  // min //
522  //-------------------------------//
523 
524 double min(const Tbl& ti) {
525 
526  // Protection
527  assert(ti.get_etat() != ETATNONDEF) ;
528 
529  // Cas particulier
530  if (ti.get_etat() == ETATZERO) {
531  return double(0) ;
532  }
533 
534  // Cas general
535  assert(ti.get_etat() == ETATQCQ) ; // sinon....
536 
537  const double* x = ti.t ;
538  double resu = x[0] ;
539  for (int i=1; i<ti.get_taille(); i++) {
540  if ( x[i] < resu ) resu = x[i] ;
541  }
542 
543  return resu ;
544 }
545 
546  //-------------------------------//
547  // norme //
548  //-------------------------------//
549 
550 double norme(const Tbl& ti) {
551 
552  // Protection
553  assert(ti.get_etat() != ETATNONDEF) ;
554 
555  double resu = 0 ;
556 
557  if (ti.get_etat() != ETATZERO) { // on n'effectue la somme que si necessaire
558 
559  assert(ti.get_etat() == ETATQCQ) ; // sinon....
560  const double* x = ti.t ;
561  for (int i=0; i<ti.get_taille(); i++) {
562  resu += fabs( x[i] ) ;
563  }
564 
565  }
566 
567  return resu ;
568 }
569 
570  //-------------------------------//
571  // diffrel //
572  //-------------------------------//
573 
574 double diffrel(const Tbl& t1, const Tbl& t2) {
575 
576  // Protections
577  assert(t1.get_etat() != ETATNONDEF) ;
578  assert(t2.get_etat() != ETATNONDEF) ;
579 
580  double norm2 = norme(t2) ;
581  double normdiff = norme(t1-t2) ;
582  double resu ;
583  if ( norm2 == double(0) ) {
584  resu = normdiff ;
585  }
586  else {
587  resu = normdiff / norm2 ;
588  }
589 
590  return resu ;
591 
592 }
593 
594  //-------------------------------//
595  // diffrelmax //
596  //-------------------------------//
597 
598 double diffrelmax(const Tbl& t1, const Tbl& t2) {
599 
600  // Protections
601  assert(t1.get_etat() != ETATNONDEF) ;
602  assert(t2.get_etat() != ETATNONDEF) ;
603 
604  double max2 = max(abs(t2)) ;
605  double maxdiff = max(abs(t1-t2)) ;
606  double resu ;
607  if ( max2 == double(0) ) {
608  resu = maxdiff ;
609  }
610  else {
611  resu = maxdiff / max2 ;
612  }
613 
614  return resu ;
615 
616 }
617 }
Lorene::racine_cubique
Cmp racine_cubique(const Cmp &)
Cube root.
Definition: cmp_math.C:245
Lorene::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
Lorene::Heaviside
Mtbl Heaviside(const Mtbl &)
Heaviside function.
Definition: mtbl_math.C:317
Lorene::exp
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::diffrelmax
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:539
Lorene::norme
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:481
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
Lorene::abs
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:410
Lorene::atan
Cmp atan(const Cmp &)
Arctangent.
Definition: cmp_math.C:195
Lorene::Tbl::dim
Dim_tbl dim
Number of dimensions, size,...
Definition: tbl.h:172
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Lorene::Tbl::t
double * t
The array of double.
Definition: tbl.h:173
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Lorene::asin
Cmp asin(const Cmp &)
Arcsine.
Definition: cmp_math.C:144
Lorene::min
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition: cmp_math.C:458
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::log10
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:322
Lorene::acos
Cmp acos(const Cmp &)
Arccosine.
Definition: cmp_math.C:169
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:94
Lorene::Tbl::get_taille
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:397
Lorene::tan
Cmp tan(const Cmp &)
Tangent.
Definition: cmp_math.C:120
Lorene::Tbl::get_etat
int get_etat() const
Gives the logical state.
Definition: tbl.h:394
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:69