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