1 /*************************************************************************** 2 math_fun.cpp - mathematical GDL library function 3 ------------------- 4 begin : July 22 2002 5 copyright : (C) 2002 by Marc Schellens 6 email : m_schellens@users.sf.net 7 ***************************************************************************/ 8 9 /*************************************************************************** 10 * * 11 * This program is free software; you can redistribute it and/or modify * 12 * it under the terms of the GNU General Public License as published by * 13 * the Free Software Foundation; either version 2 of the License, or * 14 * (at your option) any later version. * 15 * * 16 ***************************************************************************/ 17 18 19 #include "includefirst.hpp" 20 21 #include <memory> 22 #include <complex> 23 #include <cmath> 24 25 #include <gsl/gsl_sf.h> 26 #include <gsl/gsl_sf_laguerre.h> 27 #include <gsl/gsl_randist.h> 28 #include <gsl/gsl_cdf.h> 29 #include <gsl/gsl_linalg.h> 30 31 #include "objects.hpp" 32 #include "datatypes.hpp" 33 #include "envt.hpp" 34 #include "math_utl.hpp" 35 #include "math_fun.hpp" 36 37 //#define GDL_DEBUG 38 #undef GDL_DEBUG 39 40 #ifdef _MSC_VER 41 #define round(f) floor(f+0.5) 42 #endif 43 44 namespace lib { 45 46 using namespace std; 47 48 template< typename srcT, typename destT> TransposeFromToGSL(srcT * src,destT * dest,SizeT srcStride1,SizeT nEl)49 void TransposeFromToGSL( srcT* src, destT* dest, SizeT srcStride1, SizeT nEl) 50 { 51 for( SizeT d = 0, ix = 0, srcDim0 = 0; d<nEl; ++d) 52 { 53 dest[ d] = src[ ix]; 54 ix += srcStride1; 55 if( ix >= nEl) 56 ix = ++srcDim0; 57 } 58 } 59 60 template< typename srcT, typename destT> FromToGSL(srcT * src,destT * dest,SizeT nEl)61 void FromToGSL( srcT* src, destT* dest, SizeT nEl) 62 { 63 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 64 { 65 #pragma omp for 66 for( SizeT d = 0; d<nEl; ++d) 67 { 68 dest[ d] = src[ d]; 69 } 70 } 71 } 72 svdc(EnvT * e)73 void svdc( EnvT* e) 74 { 75 e->NParam( 4); 76 77 static int doubleKWIx = e->KeywordIx( "DOUBLE"); 78 bool doubleKW = e->KeywordSet( doubleKWIx); 79 80 BaseGDL* A = e->GetParDefined( 0); 81 doubleKW = doubleKW || (A->Type() == GDL_DOUBLE) || (A->Type() == GDL_COMPLEXDBL); 82 83 if( doubleKW) 84 { 85 A = e->GetParAs< DDoubleGDL>( 0); 86 } 87 else 88 { 89 A = e->GetParAs< DFloatGDL>( 0); 90 } 91 if( A->Rank() != 2) 92 e->Throw( "Argument must be a 2-D matrix: "+e->GetParString(0)); 93 94 e->AssureGlobalPar( 1); // W 95 e->AssureGlobalPar( 2); // U 96 e->AssureGlobalPar( 3); // V 97 98 static int columnKWIx = e->KeywordIx( "COLUMN"); 99 bool columnKW = e->KeywordSet( columnKWIx); 100 static int itmaxKWIx = e->KeywordIx( "ITMAX"); 101 DLong itMax = 30; 102 e->AssureLongScalarKWIfPresent( itmaxKWIx, itMax); 103 104 DLong n; 105 DLong m; 106 if( columnKW) 107 { 108 n = A->Dim( 1); 109 m = A->Dim( 0); 110 } 111 else 112 { 113 n = A->Dim( 0); 114 m = A->Dim( 1); 115 } 116 if( m < n) 117 e->Throw( "SVD of NxM matrix with N>M is not implemented yet."); 118 119 DLong nEl = A->N_Elements(); 120 121 if( doubleKW) 122 { 123 DDoubleGDL* AA = static_cast<DDoubleGDL*>( A); 124 125 gsl_matrix *aGSL = gsl_matrix_alloc( m, n); 126 GDLGuard<gsl_matrix> g1( aGSL, gsl_matrix_free); 127 if( !columnKW) 128 memcpy(aGSL->data, &(*AA)[0], nEl*sizeof( DDouble)); 129 else 130 TransposeFromToGSL< DDouble, double>( &(*AA)[0], aGSL->data, AA->Dim( 0), nEl); 131 132 gsl_matrix *vGSL = gsl_matrix_alloc( n, n); 133 GDLGuard<gsl_matrix> g2( vGSL, gsl_matrix_free); 134 gsl_vector *wGSL = gsl_vector_alloc( n); 135 GDLGuard<gsl_vector> g3( wGSL, gsl_vector_free); 136 137 gsl_vector *work = gsl_vector_alloc( n); 138 GDLGuard<gsl_vector> g4( work, gsl_vector_free); 139 gsl_linalg_SV_decomp( aGSL, vGSL, wGSL, work); 140 // gsl_vector_free( work); 141 142 // aGSL -> uGSL 143 gsl_matrix *uGSL = aGSL; // why? 144 145 // U 146 DDoubleGDL* U = new DDoubleGDL( AA->Dim(), BaseGDL::NOZERO); 147 if( !columnKW) 148 memcpy( &(*U)[0], uGSL->data, nEl*sizeof( DDouble)); 149 else 150 TransposeFromToGSL< double, DDouble>( uGSL->data, &(*U)[0], U->Dim( 1), nEl); 151 // gsl_matrix_free( uGSL); 152 e->SetPar( 2, U); 153 154 // V 155 DDoubleGDL* V = new DDoubleGDL( dimension( n, n), BaseGDL::NOZERO); 156 if( !columnKW) 157 memcpy( &(*V)[0], vGSL->data, n*n*sizeof( DDouble)); 158 else 159 TransposeFromToGSL< double, DDouble>( vGSL->data, &(*V)[0], n, n*n); 160 // gsl_matrix_free( vGSL); 161 e->SetPar( 3, V); 162 163 // W 164 DDoubleGDL* W = new DDoubleGDL( dimension( n), BaseGDL::NOZERO); 165 memcpy( &(*W)[0], wGSL->data, n*sizeof( DDouble)); 166 // gsl_vector_free( wGSL); 167 e->SetPar( 1, W); 168 } 169 else // float 170 { 171 DFloatGDL* AA = static_cast<DFloatGDL*>( A); 172 173 gsl_matrix *aGSL = gsl_matrix_alloc( m, n); 174 GDLGuard<gsl_matrix> g1( aGSL, gsl_matrix_free); 175 if( !columnKW) 176 FromToGSL< DFloat, double>( &(*AA)[0], aGSL->data, nEl); 177 else 178 TransposeFromToGSL< DFloat, double>( &(*AA)[0], aGSL->data, AA->Dim( 0), nEl); 179 180 gsl_matrix *vGSL = gsl_matrix_alloc( n, n); 181 GDLGuard<gsl_matrix> g2( vGSL, gsl_matrix_free); 182 gsl_vector *wGSL = gsl_vector_alloc( n); 183 GDLGuard<gsl_vector> g3( wGSL, gsl_vector_free); 184 185 gsl_vector *work = gsl_vector_alloc( n); 186 GDLGuard<gsl_vector> g4( work, gsl_vector_free); 187 gsl_linalg_SV_decomp( aGSL, vGSL, wGSL, work); 188 // gsl_vector_free( work); 189 190 // aGSL -> uGSL 191 gsl_matrix *uGSL = aGSL; // why? 192 193 // U 194 DFloatGDL* U = new DFloatGDL( AA->Dim(), BaseGDL::NOZERO); 195 if( !columnKW) 196 FromToGSL< double, DFloat>( uGSL->data, &(*U)[0], nEl); 197 else 198 TransposeFromToGSL< double, DFloat>( uGSL->data, &(*U)[0], U->Dim( 1), nEl); 199 // gsl_matrix_free( uGSL); 200 e->SetPar( 2, U); 201 202 // V 203 DFloatGDL* V = new DFloatGDL( dimension( n, n), BaseGDL::NOZERO); 204 if( !columnKW) 205 FromToGSL< double, DFloat>( vGSL->data, &(*V)[0], n*n); 206 else 207 TransposeFromToGSL< double, DFloat>( vGSL->data, &(*V)[0], n, n*n); 208 // gsl_matrix_free( vGSL); 209 e->SetPar( 3, V); 210 211 // W 212 DFloatGDL* W = new DFloatGDL( dimension( n), BaseGDL::NOZERO); 213 FromToGSL< double, DFloat>( wGSL->data, &(*W)[0], n); 214 // gsl_vector_free( wGSL); 215 e->SetPar( 1, W); 216 } 217 } 218 219 ////These functions return LONG or LONG64. 220 //// Note: IDL adds the supplementary optization that they are not called at all if they are applied on integer types, as, e.g., round(integer)==integer. 221 //// TBD: This trick must be performed at interpreter level. 222 /// O3 compiled code speed is OK wrt IDL in all cases. 223 224 template< typename T> //ONLY USED FOR FLOATS AND DOUBLE 225 BaseGDL* round_fun_template(BaseGDL * p0,bool isKWSetL64)226 round_fun_template (BaseGDL* p0, bool isKWSetL64) 227 { 228 T* p0C = static_cast<T*> (p0); 229 SizeT nEl = p0->N_Elements (); 230 // L64 keyword support 231 if (isKWSetL64) 232 { 233 DLong64GDL* res = new DLong64GDL (p0C->Dim (), BaseGDL::NOZERO); 234 if (nEl == 1) 235 { 236 (*res)[ 0] = round ((*p0C)[ 0]); 237 return res; 238 } 239 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 240 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = round ((*p0C)[ i]); 241 return res; 242 } 243 else 244 { 245 DLongGDL* res = new DLongGDL (p0C->Dim (), BaseGDL::NOZERO); 246 if (nEl == 1) 247 { 248 (*res)[ 0] = round ((*p0C)[ 0]); 249 return res; 250 } 251 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 252 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = round ((*p0C)[ i]); 253 254 return res; 255 } 256 } 257 258 BaseGDL* round_fun(EnvT * e)259 round_fun (EnvT* e) { 260 e->NParam (1); //, "ROUND"); 261 BaseGDL* p0 = e->GetParDefined (0); //, "ROUND"); 262 263 SizeT nEl = p0->N_Elements (); 264 if (nEl == 0) e->Throw ("Variable is undefined: " + e->GetParString (0)); 265 if (!(NumericType (p0->Type ()))) e->Throw (p0->TypeStr () + " expression: not allowed in this context: " + e->GetParString (0)); 266 267 //L64 means it: output IS ALWAYS L64. 268 if (e->KeywordSet (0)) 269 { //L64 270 if (p0->Type () == GDL_COMPLEX) 271 { 272 DComplexGDL* p0C = static_cast<DComplexGDL*> (p0); 273 SizeT nEl = p0->N_Elements (); 274 DLong64GDL* res = new DLong64GDL (p0C->Dim (), BaseGDL::NOZERO); 275 if (nEl == 1) 276 { 277 (*res)[ 0] = round ((*p0C)[ 0].real ()); 278 return res; 279 } 280 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 281 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = round ((*p0C)[ i].real ()); 282 return res; 283 } 284 else if (p0->Type () == GDL_COMPLEXDBL) 285 { 286 DComplexDblGDL* p0C = static_cast<DComplexDblGDL*> (p0); 287 SizeT nEl = p0->N_Elements (); 288 DLong64GDL* res = new DLong64GDL (p0C->Dim (), BaseGDL::NOZERO); 289 if (nEl == 1) 290 { 291 (*res)[ 0] = round ((*p0C)[ 0].real ()); 292 return res; 293 } 294 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 295 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = round ((*p0C)[ i].real ()); 296 return res; 297 } 298 else if (p0->Type () == GDL_DOUBLE) return round_fun_template< DDoubleGDL>(p0, true); 299 else if (p0->Type () == GDL_FLOAT) return round_fun_template< DFloatGDL>(p0, true); 300 else if (p0->Type () == GDL_LONG64) return p0->Dup (); 301 else return p0->Convert2 (GDL_LONG64, BaseGDL::COPY); 302 } 303 else 304 { 305 if (p0->Type () == GDL_COMPLEX) 306 { 307 DComplexGDL* p0C = static_cast<DComplexGDL*> (p0); 308 SizeT nEl = p0->N_Elements (); 309 DLongGDL* res = new DLongGDL (p0C->Dim (), BaseGDL::NOZERO); 310 if (nEl == 1) 311 { 312 (*res)[ 0] = round ((*p0C)[ 0].real ()); 313 return res; 314 } 315 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 316 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = round ((*p0C)[ i].real ()); 317 return res; 318 } 319 else if (p0->Type () == GDL_COMPLEXDBL) 320 { 321 DComplexDblGDL* p0C = static_cast<DComplexDblGDL*> (p0); 322 SizeT nEl = p0->N_Elements (); 323 DLongGDL* res = new DLongGDL (p0C->Dim (), BaseGDL::NOZERO); 324 if (nEl == 1) 325 { 326 (*res)[ 0] = round ((*p0C)[ 0].real ()); 327 return res; 328 } 329 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 330 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = round ((*p0C)[ i].real ()); 331 return res; 332 } 333 else if (p0->Type () == GDL_DOUBLE) return round_fun_template< DDoubleGDL>(p0, false); 334 else if (p0->Type () == GDL_FLOAT) return round_fun_template< DFloatGDL>(p0, false); 335 else return p0->Dup (); 336 } 337 } 338 339 template< typename T> ceil_fun_template(BaseGDL * p0,bool isKWSetL64)340 BaseGDL* ceil_fun_template(BaseGDL* p0, bool isKWSetL64) 341 { 342 T* p0C = static_cast<T*> (p0); 343 SizeT nEl = p0->N_Elements (); 344 // L64 keyword support 345 if (isKWSetL64) 346 { 347 DLong64GDL* res = new DLong64GDL (p0C->Dim (), BaseGDL::NOZERO); 348 if (nEl == 1) 349 { 350 (*res)[ 0] = ceil ((*p0C)[ 0]); 351 return res; 352 } 353 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 354 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = ceil ((*p0C)[ i]); 355 return res; 356 } 357 else 358 { 359 DLongGDL* res = new DLongGDL (p0C->Dim (), BaseGDL::NOZERO); 360 if (nEl == 1) 361 { 362 (*res)[ 0] = ceil ((*p0C)[ 0]); 363 return res; 364 } 365 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 366 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = ceil ((*p0C)[ i]); 367 368 return res; 369 } 370 } 371 ceil_fun(EnvT * e)372 BaseGDL* ceil_fun(EnvT* e) 373 { 374 e->NParam (1); //, "CEIL"); 375 BaseGDL* p0 = e->GetParDefined (0); //, "CEIL"); 376 377 SizeT nEl = p0->N_Elements (); 378 if (nEl == 0) e->Throw ("Variable is undefined: " + e->GetParString (0)); 379 if (!(NumericType (p0->Type ()))) e->Throw (p0->TypeStr () + " expression: not allowed in this context: " + e->GetParString (0)); 380 381 //L64 means it: output IS ALWAYS L64. 382 if (e->KeywordSet (0)) 383 { //L64 384 if (p0->Type () == GDL_COMPLEX) 385 { 386 DComplexGDL* p0C = static_cast<DComplexGDL*> (p0); 387 SizeT nEl = p0->N_Elements (); 388 DLong64GDL* res = new DLong64GDL (p0C->Dim (), BaseGDL::NOZERO); 389 if (nEl == 1) 390 { 391 (*res)[ 0] = ceil ((*p0C)[ 0].real ()); 392 return res; 393 } 394 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 395 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = ceil ((*p0C)[ i].real ()); 396 return res; 397 } 398 else if (p0->Type () == GDL_COMPLEXDBL) 399 { 400 DComplexDblGDL* p0C = static_cast<DComplexDblGDL*> (p0); 401 SizeT nEl = p0->N_Elements (); 402 DLong64GDL* res = new DLong64GDL (p0C->Dim (), BaseGDL::NOZERO); 403 if (nEl == 1) 404 { 405 (*res)[ 0] = ceil ((*p0C)[ 0].real ()); 406 return res; 407 } 408 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 409 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = ceil ((*p0C)[ i].real ()); 410 return res; 411 } 412 else if (p0->Type () == GDL_DOUBLE) return ceil_fun_template< DDoubleGDL>(p0, true); 413 else if (p0->Type () == GDL_FLOAT) return ceil_fun_template< DFloatGDL>(p0, true); 414 else if (p0->Type () == GDL_LONG64) return p0->Dup (); 415 else return p0->Convert2 (GDL_LONG64, BaseGDL::COPY); 416 } 417 else 418 { 419 if (p0->Type () == GDL_COMPLEX) 420 { 421 DComplexGDL* p0C = static_cast<DComplexGDL*> (p0); 422 SizeT nEl = p0->N_Elements (); 423 DLongGDL* res = new DLongGDL (p0C->Dim (), BaseGDL::NOZERO); 424 if (nEl == 1) 425 { 426 (*res)[ 0] = ceil ((*p0C)[ 0].real ()); 427 return res; 428 } 429 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 430 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = ceil ((*p0C)[ i].real ()); 431 return res; 432 } 433 else if (p0->Type () == GDL_COMPLEXDBL) 434 { 435 DComplexDblGDL* p0C = static_cast<DComplexDblGDL*> (p0); 436 SizeT nEl = p0->N_Elements (); 437 DLongGDL* res = new DLongGDL (p0C->Dim (), BaseGDL::NOZERO); 438 if (nEl == 1) 439 { 440 (*res)[ 0] = ceil ((*p0C)[ 0].real ()); 441 return res; 442 } 443 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 444 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = ceil ((*p0C)[ i].real ()); 445 return res; 446 } 447 else if (p0->Type () == GDL_DOUBLE) return ceil_fun_template< DDoubleGDL>(p0, false); 448 else if (p0->Type () == GDL_FLOAT) return ceil_fun_template< DFloatGDL>(p0, false); 449 else return p0->Dup (); 450 } 451 } 452 453 template< typename T> floor_fun_template(BaseGDL * p0,bool isKWSetL64)454 BaseGDL* floor_fun_template(BaseGDL* p0, bool isKWSetL64) 455 { 456 T* p0C = static_cast<T*> (p0); 457 SizeT nEl = p0->N_Elements (); 458 // L64 keyword support 459 if (isKWSetL64) 460 { 461 DLong64GDL* res = new DLong64GDL (p0C->Dim (), BaseGDL::NOZERO); 462 if (nEl == 1) 463 { 464 (*res)[ 0] = floor ((*p0C)[ 0]); 465 return res; 466 } 467 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 468 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = floor ((*p0C)[ i]); 469 return res; 470 } 471 else 472 { 473 DLongGDL* res = new DLongGDL (p0C->Dim (), BaseGDL::NOZERO); 474 if (nEl == 1) 475 { 476 (*res)[ 0] = floor ((*p0C)[ 0]); 477 return res; 478 } 479 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 480 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = floor ((*p0C)[ i]); 481 482 return res; 483 } 484 } 485 floor_fun(EnvT * e)486 BaseGDL* floor_fun(EnvT* e) 487 { 488 e->NParam (1); //, "floor"); 489 BaseGDL* p0 = e->GetParDefined (0); //, "floor"); 490 491 SizeT nEl = p0->N_Elements (); 492 if (nEl == 0) e->Throw ("Variable is undefined: " + e->GetParString (0)); 493 if (!(NumericType (p0->Type ()))) e->Throw (p0->TypeStr () + " expression: not allowed in this context: " + e->GetParString (0)); 494 495 //L64 means it: output IS ALWAYS L64. 496 if (e->KeywordSet (0)) 497 { //L64 498 if (p0->Type () == GDL_COMPLEX) 499 { 500 DComplexGDL* p0C = static_cast<DComplexGDL*> (p0); 501 SizeT nEl = p0->N_Elements (); 502 DLong64GDL* res = new DLong64GDL (p0C->Dim (), BaseGDL::NOZERO); 503 if (nEl == 1) 504 { 505 (*res)[ 0] = floor ((*p0C)[ 0].real ()); 506 return res; 507 } 508 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 509 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = floor ((*p0C)[ i].real ()); 510 return res; 511 } 512 else if (p0->Type () == GDL_COMPLEXDBL) 513 { 514 DComplexDblGDL* p0C = static_cast<DComplexDblGDL*> (p0); 515 SizeT nEl = p0->N_Elements (); 516 DLong64GDL* res = new DLong64GDL (p0C->Dim (), BaseGDL::NOZERO); 517 if (nEl == 1) 518 { 519 (*res)[ 0] = floor ((*p0C)[ 0].real ()); 520 return res; 521 } 522 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 523 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = floor ((*p0C)[ i].real ()); 524 return res; 525 } 526 else if (p0->Type () == GDL_DOUBLE) return floor_fun_template< DDoubleGDL>(p0, true); 527 else if (p0->Type () == GDL_FLOAT) return floor_fun_template< DFloatGDL>(p0, true); 528 else if (p0->Type () == GDL_LONG64) return p0->Dup (); 529 else return p0->Convert2 (GDL_LONG64, BaseGDL::COPY); 530 } 531 else 532 { 533 if (p0->Type () == GDL_COMPLEX) 534 { 535 DComplexGDL* p0C = static_cast<DComplexGDL*> (p0); 536 SizeT nEl = p0->N_Elements (); 537 DLongGDL* res = new DLongGDL (p0C->Dim (), BaseGDL::NOZERO); 538 if (nEl == 1) 539 { 540 (*res)[ 0] = floor ((*p0C)[ 0].real ()); 541 return res; 542 } 543 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 544 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = floor ((*p0C)[ i].real ()); 545 return res; 546 } 547 else if (p0->Type () == GDL_COMPLEXDBL) 548 { 549 DComplexDblGDL* p0C = static_cast<DComplexDblGDL*> (p0); 550 SizeT nEl = p0->N_Elements (); 551 DLongGDL* res = new DLongGDL (p0C->Dim (), BaseGDL::NOZERO); 552 if (nEl == 1) 553 { 554 (*res)[ 0] = floor ((*p0C)[ 0].real ()); 555 return res; 556 } 557 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 558 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = floor ((*p0C)[ i].real ()); 559 return res; 560 } 561 else if (p0->Type () == GDL_DOUBLE) return floor_fun_template< DDoubleGDL>(p0, false); 562 else if (p0->Type () == GDL_FLOAT) return floor_fun_template< DFloatGDL>(p0, false); 563 else return p0->Dup (); 564 } 565 } 566 567 // GDL Direct functions (no new environment created because the function has no keywords and only one parameter-->no overheads) 568 //RETURNS: float (all 32 bits + strings), double, complex, double complex outputs. 569 570 //SQRT 571 template< typename T> sqrt_fun_template(BaseGDL * p0)572 BaseGDL* sqrt_fun_template(BaseGDL* p0) { 573 T* p0C = static_cast<T*> (p0); 574 T* res = new T(p0C->Dim(), BaseGDL::NOZERO); 575 SizeT nEl = p0->N_Elements(); 576 if (nEl == 1) { 577 (*res)[0] = sqrt((*p0C)[0]); 578 return res; 579 } 580 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 581 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = sqrt((*p0C)[ i]); 582 return res; 583 } 584 585 template< typename T> sqrt_fun_template_grab(BaseGDL * p0)586 BaseGDL* sqrt_fun_template_grab(BaseGDL* p0) { 587 T* p0C = static_cast<T*> (p0); 588 SizeT nEl = p0->N_Elements(); 589 if (nEl == 1) { 590 (*p0C)[0] = sqrt((*p0C)[0]); 591 return p0; 592 } 593 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 594 for (SizeT i = 0; i < nEl; ++i) (*p0C)[ i] = sqrt((*p0C)[ i]); 595 return p0; 596 } 597 sqrt_fun(BaseGDL * p0,bool isReference)598 BaseGDL* sqrt_fun(BaseGDL* p0, bool isReference) { 599 SizeT nEl = p0->N_Elements(); 600 DType p0Type = p0->Type(); 601 if (p0Type == GDL_COMPLEX) 602 if (isReference) return sqrt_fun_template< DComplexGDL>(p0); 603 else return sqrt_fun_template_grab< DComplexGDL>(p0); 604 else if (p0Type == GDL_COMPLEXDBL) 605 if (isReference) return sqrt_fun_template< DComplexDblGDL>(p0); 606 else return sqrt_fun_template_grab< DComplexDblGDL>(p0); 607 else if (p0Type == GDL_DOUBLE) 608 if (isReference) return sqrt_fun_template< DDoubleGDL>(p0); 609 else return sqrt_fun_template_grab< DDoubleGDL>(p0); 610 else if (p0Type == GDL_FLOAT) 611 if (isReference) return sqrt_fun_template< DFloatGDL>(p0); 612 else return sqrt_fun_template_grab< DFloatGDL>(p0); 613 else { 614 DFloatGDL* res = static_cast<DFloatGDL*> (p0->Convert2(GDL_FLOAT, BaseGDL::COPY)); 615 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 616 for (SizeT i = 0; i < nEl; ++i)(*res)[ i] = sqrt((*res)[ i]); 617 return res; 618 } 619 } 620 621 //SIN 622 template< typename T> sin_fun_template(BaseGDL * p0)623 BaseGDL* sin_fun_template (BaseGDL* p0) 624 { 625 T* p0C = static_cast<T*> (p0); 626 T* res = new T (p0C->Dim (), BaseGDL::NOZERO); 627 SizeT nEl = p0->N_Elements (); 628 if (nEl == 1) 629 { 630 (*res)[0] = sin ((*p0C)[0]); 631 return res; 632 } 633 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 634 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = sin ((*p0C)[ i]); 635 return res; 636 } 637 638 template< typename T> sin_fun_template_grab(BaseGDL * p0)639 BaseGDL* sin_fun_template_grab (BaseGDL* p0) 640 { 641 T* p0C = static_cast<T*> (p0); 642 SizeT nEl = p0->N_Elements (); 643 if (nEl == 1) 644 { 645 (*p0C)[0] = sin ((*p0C)[0]); 646 return p0; 647 } 648 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 649 for (SizeT i = 0; i < nEl; ++i) (*p0C)[ i] = sin ((*p0C)[ i]); 650 return p0; 651 } 652 sin_fun(BaseGDL * p0,bool isReference)653 BaseGDL* sin_fun (BaseGDL* p0, bool isReference) 654 { 655 SizeT nEl = p0->N_Elements (); 656 DType p0Type = p0->Type (); 657 if (p0Type == GDL_COMPLEX) 658 if (isReference) return sin_fun_template< DComplexGDL>(p0); else return sin_fun_template_grab< DComplexGDL>(p0); 659 else if (p0Type == GDL_COMPLEXDBL) 660 if (isReference) return sin_fun_template< DComplexDblGDL>(p0); else return sin_fun_template_grab< DComplexDblGDL>(p0); 661 else if (p0Type == GDL_DOUBLE) 662 if (isReference) return sin_fun_template< DDoubleGDL>(p0); else return sin_fun_template_grab< DDoubleGDL>(p0); 663 else if (p0Type == GDL_FLOAT) 664 if (isReference) return sin_fun_template< DFloatGDL>(p0); else return sin_fun_template_grab< DFloatGDL>(p0); 665 else 666 { 667 DFloatGDL* res = static_cast<DFloatGDL*> (p0->Convert2 (GDL_FLOAT, BaseGDL::COPY)); 668 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 669 for (SizeT i = 0; i < nEl; ++i)(*res)[ i] = sin ((*res)[ i]); 670 return res; 671 } 672 } 673 674 //COS 675 template< typename T> cos_fun_template(BaseGDL * p0)676 BaseGDL* cos_fun_template(BaseGDL* p0) { 677 T* p0C = static_cast<T*> (p0); 678 T* res = new T(p0C->Dim(), BaseGDL::NOZERO); 679 SizeT nEl = p0->N_Elements(); 680 if (nEl == 1) { 681 (*res)[0] = cos((*p0C)[0]); 682 return res; 683 } 684 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 685 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = cos((*p0C)[ i]); 686 return res; 687 } 688 689 template< typename T> cos_fun_template_grab(BaseGDL * p0)690 BaseGDL* cos_fun_template_grab(BaseGDL* p0) { 691 T* p0C = static_cast<T*> (p0); 692 SizeT nEl = p0->N_Elements(); 693 if (nEl == 1) { 694 (*p0C)[0] = cos((*p0C)[0]); 695 return p0; 696 } 697 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 698 for (SizeT i = 0; i < nEl; ++i) (*p0C)[ i] = cos((*p0C)[ i]); 699 return p0; 700 } 701 cos_fun(BaseGDL * p0,bool isReference)702 BaseGDL* cos_fun(BaseGDL* p0, bool isReference) { 703 SizeT nEl = p0->N_Elements(); 704 DType p0Type = p0->Type(); 705 if (p0Type == GDL_COMPLEX) 706 if (isReference) return cos_fun_template< DComplexGDL>(p0); 707 else return cos_fun_template_grab< DComplexGDL>(p0); 708 else if (p0Type == GDL_COMPLEXDBL) 709 if (isReference) return cos_fun_template< DComplexDblGDL>(p0); 710 else return cos_fun_template_grab< DComplexDblGDL>(p0); 711 else if (p0Type == GDL_DOUBLE) 712 if (isReference) return cos_fun_template< DDoubleGDL>(p0); 713 else return cos_fun_template_grab< DDoubleGDL>(p0); 714 else if (p0Type == GDL_FLOAT) 715 if (isReference) return cos_fun_template< DFloatGDL>(p0); 716 else return cos_fun_template_grab< DFloatGDL>(p0); 717 else { 718 DFloatGDL* res = static_cast<DFloatGDL*> (p0->Convert2(GDL_FLOAT, BaseGDL::COPY)); 719 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 720 for (SizeT i = 0; i < nEl; ++i)(*res)[ i] = cos((*res)[ i]); 721 return res; 722 } 723 } 724 725 //TAN 726 template< typename T> tan_fun_template(BaseGDL * p0)727 BaseGDL* tan_fun_template(BaseGDL* p0) { 728 T* p0C = static_cast<T*> (p0); 729 T* res = new T(p0C->Dim(), BaseGDL::NOZERO); 730 SizeT nEl = p0->N_Elements(); 731 if (nEl == 1) { 732 (*res)[0] = tan((*p0C)[0]); 733 return res; 734 } 735 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 736 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = tan((*p0C)[ i]); 737 return res; 738 } 739 740 template< typename T> tan_fun_template_grab(BaseGDL * p0)741 BaseGDL* tan_fun_template_grab(BaseGDL* p0) { 742 T* p0C = static_cast<T*> (p0); 743 SizeT nEl = p0->N_Elements(); 744 if (nEl == 1) { 745 (*p0C)[0] = tan((*p0C)[0]); 746 return p0; 747 } 748 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 749 for (SizeT i = 0; i < nEl; ++i) (*p0C)[ i] = tan((*p0C)[ i]); 750 return p0; 751 } 752 tan_fun(BaseGDL * p0,bool isReference)753 BaseGDL* tan_fun(BaseGDL* p0, bool isReference) { 754 SizeT nEl = p0->N_Elements(); 755 DType p0Type = p0->Type(); 756 if (p0Type == GDL_COMPLEX) 757 if (isReference) return tan_fun_template< DComplexGDL>(p0); 758 else return tan_fun_template_grab< DComplexGDL>(p0); 759 else if (p0Type == GDL_COMPLEXDBL) 760 if (isReference) return tan_fun_template< DComplexDblGDL>(p0); 761 else return tan_fun_template_grab< DComplexDblGDL>(p0); 762 else if (p0Type == GDL_DOUBLE) 763 if (isReference) return tan_fun_template< DDoubleGDL>(p0); 764 else return tan_fun_template_grab< DDoubleGDL>(p0); 765 else if (p0Type == GDL_FLOAT) 766 if (isReference) return tan_fun_template< DFloatGDL>(p0); 767 else return tan_fun_template_grab< DFloatGDL>(p0); 768 else { 769 DFloatGDL* res = static_cast<DFloatGDL*> (p0->Convert2(GDL_FLOAT, BaseGDL::COPY)); 770 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 771 for (SizeT i = 0; i < nEl; ++i)(*res)[ i] = tan((*res)[ i]); 772 return res; 773 } 774 } 775 776 //SINH 777 template< typename T> sinh_fun_template(BaseGDL * p0)778 BaseGDL* sinh_fun_template(BaseGDL* p0) { 779 T* p0C = static_cast<T*> (p0); 780 T* res = new T(p0C->Dim(), BaseGDL::NOZERO); 781 SizeT nEl = p0->N_Elements(); 782 if (nEl == 1) { 783 (*res)[0] = sinh((*p0C)[0]); 784 return res; 785 } 786 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 787 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = sinh((*p0C)[ i]); 788 return res; 789 } 790 791 template< typename T> sinh_fun_template_grab(BaseGDL * p0)792 BaseGDL* sinh_fun_template_grab(BaseGDL* p0) { 793 T* p0C = static_cast<T*> (p0); 794 SizeT nEl = p0->N_Elements(); 795 if (nEl == 1) { 796 (*p0C)[0] = sinh((*p0C)[0]); 797 return p0; 798 } 799 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 800 for (SizeT i = 0; i < nEl; ++i) (*p0C)[ i] = sinh((*p0C)[ i]); 801 return p0; 802 } 803 sinh_fun(BaseGDL * p0,bool isReference)804 BaseGDL* sinh_fun(BaseGDL* p0, bool isReference) { 805 SizeT nEl = p0->N_Elements(); 806 DType p0Type = p0->Type(); 807 if (p0Type == GDL_COMPLEX) 808 if (isReference) return sinh_fun_template< DComplexGDL>(p0); 809 else return sinh_fun_template_grab< DComplexGDL>(p0); 810 else if (p0Type == GDL_COMPLEXDBL) 811 if (isReference) return sinh_fun_template< DComplexDblGDL>(p0); 812 else return sinh_fun_template_grab< DComplexDblGDL>(p0); 813 else if (p0Type == GDL_DOUBLE) 814 if (isReference) return sinh_fun_template< DDoubleGDL>(p0); 815 else return sinh_fun_template_grab< DDoubleGDL>(p0); 816 else if (p0Type == GDL_FLOAT) 817 if (isReference) return sinh_fun_template< DFloatGDL>(p0); 818 else return sinh_fun_template_grab< DFloatGDL>(p0); 819 else { 820 DFloatGDL* res = static_cast<DFloatGDL*> (p0->Convert2(GDL_FLOAT, BaseGDL::COPY)); 821 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 822 for (SizeT i = 0; i < nEl; ++i)(*res)[ i] = sinh((*res)[ i]); 823 return res; 824 } 825 } 826 827 //COSH 828 template< typename T> cosh_fun_template(BaseGDL * p0)829 BaseGDL* cosh_fun_template(BaseGDL* p0) { 830 T* p0C = static_cast<T*> (p0); 831 T* res = new T(p0C->Dim(), BaseGDL::NOZERO); 832 SizeT nEl = p0->N_Elements(); 833 if (nEl == 1) { 834 (*res)[0] = cosh((*p0C)[0]); 835 return res; 836 } 837 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 838 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = cosh((*p0C)[ i]); 839 return res; 840 } 841 842 template< typename T> cosh_fun_template_grab(BaseGDL * p0)843 BaseGDL* cosh_fun_template_grab(BaseGDL* p0) { 844 T* p0C = static_cast<T*> (p0); 845 SizeT nEl = p0->N_Elements(); 846 if (nEl == 1) { 847 (*p0C)[0] = cosh((*p0C)[0]); 848 return p0; 849 } 850 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 851 for (SizeT i = 0; i < nEl; ++i) (*p0C)[ i] = cosh((*p0C)[ i]); 852 return p0; 853 } 854 cosh_fun(BaseGDL * p0,bool isReference)855 BaseGDL* cosh_fun(BaseGDL* p0, bool isReference) { 856 SizeT nEl = p0->N_Elements(); 857 DType p0Type = p0->Type(); 858 if (p0Type == GDL_COMPLEX) 859 if (isReference) return cosh_fun_template< DComplexGDL>(p0); 860 else return cosh_fun_template_grab< DComplexGDL>(p0); 861 else if (p0Type == GDL_COMPLEXDBL) 862 if (isReference) return cosh_fun_template< DComplexDblGDL>(p0); 863 else return cosh_fun_template_grab< DComplexDblGDL>(p0); 864 else if (p0Type == GDL_DOUBLE) 865 if (isReference) return cosh_fun_template< DDoubleGDL>(p0); 866 else return cosh_fun_template_grab< DDoubleGDL>(p0); 867 else if (p0Type == GDL_FLOAT) 868 if (isReference) return cosh_fun_template< DFloatGDL>(p0); 869 else return cosh_fun_template_grab< DFloatGDL>(p0); 870 else { 871 DFloatGDL* res = static_cast<DFloatGDL*> (p0->Convert2(GDL_FLOAT, BaseGDL::COPY)); 872 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 873 for (SizeT i = 0; i < nEl; ++i)(*res)[ i] = cosh((*res)[ i]); 874 return res; 875 } 876 } 877 878 //TANH 879 template< typename T> tanh_fun_template(BaseGDL * p0)880 BaseGDL* tanh_fun_template(BaseGDL* p0) { 881 T* p0C = static_cast<T*> (p0); 882 T* res = new T(p0C->Dim(), BaseGDL::NOZERO); 883 SizeT nEl = p0->N_Elements(); 884 if (nEl == 1) { 885 (*res)[0] = tanh((*p0C)[0]); 886 return res; 887 } 888 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 889 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = tanh((*p0C)[ i]); 890 return res; 891 } 892 893 template< typename T> tanh_fun_template_grab(BaseGDL * p0)894 BaseGDL* tanh_fun_template_grab(BaseGDL* p0) { 895 T* p0C = static_cast<T*> (p0); 896 SizeT nEl = p0->N_Elements(); 897 if (nEl == 1) { 898 (*p0C)[0] = tanh((*p0C)[0]); 899 return p0; 900 } 901 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 902 for (SizeT i = 0; i < nEl; ++i) (*p0C)[ i] = tanh((*p0C)[ i]); 903 return p0; 904 } 905 tanh_fun(BaseGDL * p0,bool isReference)906 BaseGDL* tanh_fun(BaseGDL* p0, bool isReference) { 907 SizeT nEl = p0->N_Elements(); 908 DType p0Type = p0->Type(); 909 if (p0Type == GDL_COMPLEX) 910 if (isReference) return tanh_fun_template< DComplexGDL>(p0); 911 else return tanh_fun_template_grab< DComplexGDL>(p0); 912 else if (p0Type == GDL_COMPLEXDBL) 913 if (isReference) return tanh_fun_template< DComplexDblGDL>(p0); 914 else return tanh_fun_template_grab< DComplexDblGDL>(p0); 915 else if (p0Type == GDL_DOUBLE) 916 if (isReference) return tanh_fun_template< DDoubleGDL>(p0); 917 else return tanh_fun_template_grab< DDoubleGDL>(p0); 918 else if (p0Type == GDL_FLOAT) 919 if (isReference) return tanh_fun_template< DFloatGDL>(p0); 920 else return tanh_fun_template_grab< DFloatGDL>(p0); 921 else { 922 DFloatGDL* res = static_cast<DFloatGDL*> (p0->Convert2(GDL_FLOAT, BaseGDL::COPY)); 923 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 924 for (SizeT i = 0; i < nEl; ++i)(*res)[ i] = tanh((*res)[ i]); 925 return res; 926 } 927 } 928 929 template< typename T> asin_fun_template(BaseGDL * p0)930 BaseGDL* asin_fun_template(BaseGDL* p0) { 931 T* p0C = static_cast<T*> (p0); 932 T* res = new T(p0C->Dim(), BaseGDL::NOZERO); 933 SizeT nEl = p0->N_Elements(); 934 if (nEl == 1) { 935 (*res)[0] = asin((*p0C)[0]); 936 return res; 937 } 938 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 939 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = asin((*p0C)[ i]); 940 return res; 941 } 942 943 template< typename T> asin_fun_template_grab(BaseGDL * p0)944 BaseGDL* asin_fun_template_grab(BaseGDL* p0) { 945 T* p0C = static_cast<T*> (p0); 946 SizeT nEl = p0->N_Elements(); 947 if (nEl == 1) { 948 (*p0C)[0] = asin((*p0C)[0]); 949 return p0; 950 } 951 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 952 for (SizeT i = 0; i < nEl; ++i) (*p0C)[ i] = asin((*p0C)[ i]); 953 return p0; 954 } 955 asin_fun(BaseGDL * p0,bool isReference)956 BaseGDL* asin_fun (BaseGDL* p0, bool isReference) 957 { 958 SizeT nEl = p0->N_Elements (); 959 DType p0Type = p0->Type (); 960 if (p0Type == GDL_COMPLEX) 961 if (isReference) return asin_fun_template< DComplexGDL>(p0); else return asin_fun_template_grab< DComplexGDL>(p0); 962 else if (p0Type == GDL_COMPLEXDBL) 963 if (isReference) return asin_fun_template< DComplexDblGDL>(p0); else return asin_fun_template_grab< DComplexDblGDL>(p0); 964 else if (p0Type == GDL_DOUBLE) 965 if (isReference) return asin_fun_template< DDoubleGDL>(p0); else return asin_fun_template_grab< DDoubleGDL>(p0); 966 else if (p0Type == GDL_FLOAT) 967 if (isReference) return asin_fun_template< DFloatGDL>(p0); else return asin_fun_template_grab< DFloatGDL>(p0); 968 else 969 { 970 DFloatGDL* res = static_cast<DFloatGDL*> (p0->Convert2 (GDL_FLOAT, BaseGDL::COPY)); 971 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 972 for (SizeT i = 0; i < nEl; ++i)(*res)[ i] = asin ((*res)[ i]); 973 return res; 974 } 975 } 976 977 template< typename T> acos_fun_template(BaseGDL * p0)978 BaseGDL* acos_fun_template(BaseGDL* p0) { 979 T* p0C = static_cast<T*> (p0); 980 T* res = new T(p0C->Dim(), BaseGDL::NOZERO); 981 SizeT nEl = p0->N_Elements(); 982 if (nEl == 1) { 983 (*res)[0] = acos((*p0C)[0]); 984 return res; 985 } 986 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 987 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = acos((*p0C)[ i]); 988 return res; 989 } 990 991 template< typename T> acos_fun_template_grab(BaseGDL * p0)992 BaseGDL* acos_fun_template_grab(BaseGDL* p0) { 993 T* p0C = static_cast<T*> (p0); 994 SizeT nEl = p0->N_Elements(); 995 if (nEl == 1) { 996 (*p0C)[0] = acos((*p0C)[0]); 997 return p0; 998 } 999 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1000 for (SizeT i = 0; i < nEl; ++i) (*p0C)[ i] = acos((*p0C)[ i]); 1001 return p0; 1002 } 1003 acos_fun(BaseGDL * p0,bool isReference)1004 BaseGDL* acos_fun (BaseGDL* p0, bool isReference) 1005 { 1006 SizeT nEl = p0->N_Elements (); 1007 DType p0Type = p0->Type (); 1008 if (p0Type == GDL_COMPLEX) 1009 if (isReference) return acos_fun_template< DComplexGDL>(p0); else return acos_fun_template_grab< DComplexGDL>(p0); 1010 else if (p0Type == GDL_COMPLEXDBL) 1011 if (isReference) return acos_fun_template< DComplexDblGDL>(p0); else return acos_fun_template_grab< DComplexDblGDL>(p0); 1012 else if (p0Type == GDL_DOUBLE) 1013 if (isReference) return acos_fun_template< DDoubleGDL>(p0); else return acos_fun_template_grab< DDoubleGDL>(p0); 1014 else if (p0Type == GDL_FLOAT) 1015 if (isReference) return acos_fun_template< DFloatGDL>(p0); else return acos_fun_template_grab< DFloatGDL>(p0); 1016 else 1017 { 1018 DFloatGDL* res = static_cast<DFloatGDL*> (p0->Convert2 (GDL_FLOAT, BaseGDL::COPY)); 1019 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1020 for (SizeT i = 0; i < nEl; ++i)(*res)[ i] = acos ((*res)[ i]); 1021 return res; 1022 } 1023 } 1024 1025 //EXP 1026 template< typename T> exp_fun_template(BaseGDL * p0)1027 BaseGDL* exp_fun_template(BaseGDL* p0) { 1028 T* p0C = static_cast<T*> (p0); 1029 T* res = new T(p0C->Dim(), BaseGDL::NOZERO); 1030 SizeT nEl = p0->N_Elements(); 1031 if (nEl == 1) { 1032 (*res)[0] = exp((*p0C)[0]); 1033 return res; 1034 } 1035 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1036 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = exp((*p0C)[ i]); 1037 return res; 1038 } 1039 1040 template< typename T> exp_fun_template_grab(BaseGDL * p0)1041 BaseGDL* exp_fun_template_grab(BaseGDL* p0) { 1042 T* p0C = static_cast<T*> (p0); 1043 SizeT nEl = p0->N_Elements(); 1044 if (nEl == 1) { 1045 (*p0C)[0] = exp((*p0C)[0]); 1046 return p0; 1047 } 1048 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1049 for (SizeT i = 0; i < nEl; ++i) (*p0C)[ i] = exp((*p0C)[ i]); 1050 return p0; 1051 } 1052 exp_fun(BaseGDL * p0,bool isReference)1053 BaseGDL* exp_fun(BaseGDL* p0, bool isReference) { 1054 SizeT nEl = p0->N_Elements(); 1055 DType p0Type = p0->Type(); 1056 if (p0Type == GDL_COMPLEX) 1057 if (isReference) return exp_fun_template< DComplexGDL>(p0); 1058 else return exp_fun_template_grab< DComplexGDL>(p0); 1059 else if (p0Type == GDL_COMPLEXDBL) 1060 if (isReference) return exp_fun_template< DComplexDblGDL>(p0); 1061 else return exp_fun_template_grab< DComplexDblGDL>(p0); 1062 else if (p0Type == GDL_DOUBLE) 1063 if (isReference) return exp_fun_template< DDoubleGDL>(p0); 1064 else return exp_fun_template_grab< DDoubleGDL>(p0); 1065 else if (p0Type == GDL_FLOAT) 1066 if (isReference) return exp_fun_template< DFloatGDL>(p0); 1067 else return exp_fun_template_grab< DFloatGDL>(p0); 1068 else { 1069 DFloatGDL* res = static_cast<DFloatGDL*> (p0->Convert2(GDL_FLOAT, BaseGDL::COPY)); 1070 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1071 for (SizeT i = 0; i < nEl; ++i)(*res)[ i] = exp((*res)[ i]); 1072 return res; 1073 } 1074 } 1075 1076 1077 // BaseGDL* alog_fun( EnvT* e) alog_fun(BaseGDL * p0,bool isReference)1078 BaseGDL* alog_fun( BaseGDL* p0, bool isReference) 1079 { 1080 if( !isReference) //e->StealLocalPar( 0)) 1081 { 1082 return p0->LogThis(); 1083 } 1084 return p0->Log(); 1085 } 1086 1087 // BaseGDL* alog10_fun( EnvT* e) alog10_fun(BaseGDL * p0,bool isReference)1088 BaseGDL* alog10_fun( BaseGDL* p0, bool isReference) 1089 { 1090 assert( p0 != NULL); 1091 assert( p0->N_Elements() > 0); 1092 if( !isReference) //e->StealLocalPar( 0)) 1093 { 1094 return p0->Log10This(); 1095 } 1096 return p0->Log10(); 1097 } 1098 1099 //Following produce Float or doubles for complex. 1100 1101 template< typename T> abs_fun_template(BaseGDL * p0)1102 BaseGDL* abs_fun_template (BaseGDL* p0) 1103 { 1104 T* p0C = static_cast<T*> (p0); 1105 T* res = new T (p0C->Dim (), BaseGDL::NOZERO); 1106 SizeT nEl = p0->N_Elements (); 1107 if (nEl == 1) 1108 { 1109 (*res)[ 0] = abs ((*p0C)[ 0]); 1110 return res; 1111 } 1112 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1113 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = abs ((*p0C)[ i]); 1114 return res; 1115 } 1116 abs_fun(BaseGDL * p0,bool isReference)1117 BaseGDL* abs_fun (BaseGDL* p0, bool isReference) 1118 { 1119 if (p0->Type () == GDL_COMPLEX) 1120 { 1121 DComplexGDL* p0C = static_cast<DComplexGDL*> (p0); 1122 DFloatGDL* res = new DFloatGDL (p0C->Dim (), BaseGDL::NOZERO); 1123 SizeT nEl = p0->N_Elements (); 1124 if (nEl == 1) 1125 { 1126 (*res)[ 0] = abs ((*p0C)[ 0]); 1127 return res; 1128 } 1129 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1130 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = abs ((*p0C)[ i]); //sqrt(Creal*Creal + Cimag*Cimag); 1131 return res; 1132 } 1133 else if (p0->Type () == GDL_COMPLEXDBL) 1134 { 1135 DComplexDblGDL* p0C = static_cast<DComplexDblGDL*> (p0); 1136 DDoubleGDL* res = new DDoubleGDL (p0C->Dim (), BaseGDL::NOZERO); 1137 SizeT nEl = p0->N_Elements (); 1138 if (nEl == 1) 1139 { 1140 (*res)[ 0] = abs ((*p0C)[ 0]); 1141 return res; 1142 } 1143 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1144 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = abs ((*p0C)[ i]); //sqrt(Creal*Creal + Cimag*Cimag); 1145 return res; 1146 } 1147 else if (p0->Type () == GDL_DOUBLE) 1148 return abs_fun_template< DDoubleGDL>(p0); 1149 else if (p0->Type () == GDL_FLOAT) 1150 return abs_fun_template< DFloatGDL>(p0); 1151 else if (p0->Type () == GDL_LONG64) 1152 return abs_fun_template< DLong64GDL>(p0); 1153 else if (p0->Type () == GDL_LONG) 1154 return abs_fun_template< DLongGDL>(p0); 1155 else if (p0->Type () == GDL_INT) 1156 return abs_fun_template< DIntGDL>(p0); 1157 else if (isReference) 1158 { 1159 if (p0->Type () == GDL_ULONG64) 1160 return p0->Dup (); 1161 else if (p0->Type () == GDL_ULONG) 1162 return p0->Dup (); 1163 else if (p0->Type () == GDL_UINT) 1164 return p0->Dup (); 1165 else if (p0->Type () == GDL_BYTE) 1166 return p0->Dup (); 1167 } 1168 else 1169 { 1170 if (p0->Type () == GDL_ULONG64) 1171 return p0; 1172 else if (p0->Type () == GDL_ULONG) 1173 return p0; 1174 else if (p0->Type () == GDL_UINT) 1175 return p0; 1176 else if (p0->Type () == GDL_BYTE) 1177 return p0; 1178 } 1179 DFloatGDL* res = static_cast<DFloatGDL*> 1180 (p0->Convert2 (GDL_FLOAT, BaseGDL::COPY)); 1181 SizeT nEl = p0->N_Elements (); 1182 if (nEl == 1) 1183 { 1184 (*res)[ 0] = abs ((*res)[ 0]); 1185 return res; 1186 } 1187 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1188 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = abs ((*res)[ i]); 1189 return res; 1190 } 1191 1192 //COMPLEX SPECIALS conj_fun(BaseGDL * p0,bool isReference)1193 BaseGDL* conj_fun (BaseGDL* p0, bool isReference)//( EnvT* e) 1194 { 1195 SizeT nEl = p0->N_Elements (); 1196 1197 if (p0->Type () == GDL_COMPLEX) 1198 { 1199 DComplexGDL* res; 1200 if (isReference) res = static_cast<DComplexGDL*> (p0)->NewResult (); else res=static_cast<DComplexGDL*>(p0); 1201 DComplexGDL* p0C = static_cast<DComplexGDL*> (p0); 1202 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1203 for (SizeT i = 0; i < nEl; ++i) (*res)[i] = std::conj ((*p0C)[i]); 1204 return res; 1205 } 1206 if (p0->Type () == GDL_COMPLEXDBL) 1207 { 1208 DComplexDblGDL* res; 1209 if (isReference) res = static_cast<DComplexDblGDL*> (p0)->NewResult (); else res=res=static_cast<DComplexDblGDL*>(p0); 1210 DComplexDblGDL* p0C = static_cast<DComplexDblGDL*> (p0); 1211 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1212 for (SizeT i = 0; i < nEl; ++i) (*res)[i] = std::conj ((*p0C)[i]); 1213 return res; 1214 } 1215 if (p0->Type () == GDL_DOUBLE || 1216 p0->Type () == GDL_LONG64 || 1217 p0->Type () == GDL_ULONG64) 1218 { 1219 DComplexDblGDL* res = static_cast<DComplexDblGDL*> (p0->Convert2 (GDL_COMPLEXDBL, BaseGDL::COPY)); 1220 return res; 1221 } 1222 1223 // all other types 1224 return static_cast<DComplexGDL*> (p0->Convert2 (GDL_COMPLEX, BaseGDL::COPY)); 1225 } 1226 //returns Double or floats imaginary_fun(BaseGDL * p0,bool isReference)1227 BaseGDL* imaginary_fun (BaseGDL* p0, bool isReference)//( EnvT* e) 1228 { 1229 SizeT nEl = p0->N_Elements (); 1230 // complex types, return imaginary part 1231 if (p0->Type () == GDL_COMPLEX) 1232 { 1233 DComplexGDL* c0 = static_cast<DComplexGDL*> (p0); 1234 DFloatGDL* res = new DFloatGDL (c0->Dim (), BaseGDL::NOZERO); 1235 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1236 for (SizeT i = 0; i < nEl; ++i) (*res)[i] = imag ((*c0)[i]); 1237 return res; 1238 } 1239 if (p0->Type () == GDL_COMPLEXDBL) 1240 { 1241 DComplexDblGDL* c0 = static_cast<DComplexDblGDL*> (p0); 1242 DDoubleGDL* res = new DDoubleGDL (c0->Dim (), BaseGDL::NOZERO); 1243 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1244 for (SizeT i = 0; i < nEl; ++i) (*res)[i] = imag ((*c0)[i]); 1245 return res; 1246 } 1247 1248 // forbidden types 1249 DType t = p0->Type (); 1250 if (t == GDL_STRING) 1251 throw GDLException ("String expression not allowed in this context."); 1252 if (t == GDL_STRUCT) 1253 throw GDLException ("Struct expression not allowed in this context."); 1254 if (t == GDL_PTR) 1255 throw GDLException ("Pointer expression not allowed in this context."); 1256 if (t == GDL_OBJ) 1257 throw GDLException ("Object reference not allowed in this context."); 1258 1259 // all other types (return array of zeros) 1260 return new DFloatGDL (p0->Dim (), BaseGDL::ZERO); // ZERO 1261 } 1262 real_part_fun(BaseGDL * p0,bool isReference)1263 BaseGDL* real_part_fun (BaseGDL* p0, bool isReference) 1264 { 1265 SizeT nEl = p0->N_Elements (); 1266 // complex types, return real part 1267 if (p0->Type () == GDL_COMPLEX) 1268 { 1269 DComplexGDL* c0 = static_cast<DComplexGDL*> (p0); 1270 DFloatGDL* res = new DFloatGDL (c0->Dim (), BaseGDL::NOZERO); 1271 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1272 for (SizeT i = 0; i < nEl; ++i) (*res)[i] = real ((*c0)[i]); 1273 return res; 1274 } 1275 if (p0->Type () == GDL_COMPLEXDBL) 1276 { 1277 DComplexDblGDL* c0 = static_cast<DComplexDblGDL*> (p0); 1278 DDoubleGDL* res = new DDoubleGDL (c0->Dim (), BaseGDL::NOZERO); 1279 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1280 for (SizeT i = 0; i < nEl; ++i) (*res)[i] = real ((*c0)[i]); 1281 return res; 1282 } 1283 1284 DType t = p0->Type (); 1285 // avoid forbidden types 1286 if (t == GDL_STRUCT) 1287 throw GDLException ("Struct expression not allowed in this context."); 1288 if (t == GDL_PTR) 1289 throw GDLException ("Pointer expression not allowed in this context."); 1290 if (t == GDL_OBJ) 1291 throw GDLException ("Object reference not allowed in this context."); 1292 // Doubles to double, copy 1293 if (t == GDL_DOUBLE){ 1294 if (isReference) return p0->Dup (); else return p0; 1295 } 1296 // Floats to float, copy 1297 if (t == GDL_FLOAT){ 1298 if (isReference) return p0->Dup (); else return p0; 1299 } 1300 // all other types to float 1301 return static_cast<DFloatGDL*> (p0->Convert2 (GDL_FLOAT, BaseGDL::COPY)); 1302 } 1303 1304 1305 //atan() is different as complex is different. Note that atan(complex(0.888,0)) does not give exactly 0 for its imaginary part, as IDL does. 1306 1307 // atan() for complex 1308 // .. is now in C++11 1309 1310 // template< typename C> 1311 // inline C atanC(const C& c) 1312 // { 1313 // // double x = c.real(); 1314 // // double x2 = x * x; 1315 // // double y = c.imag(); 1316 // // return C(0.5 * atan2(2.0*x, 1.0 - x2 - y*y), 0.25 * log( (x2 + (y+1)*(y+1)) / (x2 + (y-1)*(y-1)) )); 1317 // const C i(0.0, 1.0); 1318 // const C one(1.0, 0.0); 1319 // return log((one + i * c) / (one - i * c)) / (C(2.0, 0.0) * i); 1320 // } 1321 1322 template< typename C> atanC(const C & c1,const C & c2)1323 inline C atanC(const C& c1, const C& c2) 1324 { 1325 const C i(0.0, 1.0); 1326 //const C one(1.0,0.0); 1327 // return -i * log((c2 + i * c1) / (sqrt(pow(c2, 2) + pow(c1, 2)))); 1328 return -i * log((c2 + i * c1) / sqrt((c2 * c2) + (c1 * c1))); 1329 } 1330 atan_fun(EnvT * e)1331 BaseGDL* atan_fun(EnvT* e) 1332 { 1333 //lots of different guards defined here to insure they do their work only at exit of atan_fun. 1334 Guard< DDoubleGDL> guardp0D; 1335 Guard< DDoubleGDL> guardp1D; 1336 Guard< DFloatGDL> guardp0F; 1337 Guard< DFloatGDL> guardp1F; 1338 Guard< DComplexDblGDL> guardp0DC; 1339 Guard< DComplexDblGDL> guardp1DC; 1340 Guard< DComplexGDL> guardp0C; 1341 Guard< DComplexGDL> guardp1C; 1342 1343 SizeT nParam = e->NParam(1); //, "ATAN"); 1344 1345 BaseGDL* p0 = e->GetParDefined(0); //, "ATAN"); 1346 1347 SizeT nEl = p0->N_Elements(); 1348 if (nEl == 0) 1349 e->Throw( 1350 "Variable is undefined: " + e->GetParString(0)); 1351 1352 if (nParam == 2) { 1353 BaseGDL* p1 = e->GetPar(1); 1354 if (p1 == NULL) 1355 e->Throw( 1356 "Variable is undefined: " + e->GetParString(1)); 1357 SizeT nEl1 = p1->N_Elements(); 1358 if (nEl1 == 0) 1359 e->Throw( 1360 "Variable is undefined: " + e->GetParString(1)); 1361 1362 DType t = (DTypeOrder[ p0->Type()] > DTypeOrder[ p1->Type()]) ? p0->Type() : p1->Type(); 1363 1364 // WRT. the previous version, written to insure that zero-dimension values are taken as of size as long as the other argument, 1365 // we hopefully keep the same behaviour but permit the speedup of parallelism by using local loop variables. 1366 dimension dim; 1367 int cas = 0; 1368 SizeT nElMin; 1369 if ((p0->Rank() == 0 && p1->Rank() == 0) || (p0->Rank() != 0 && p1->Rank() != 0)) { 1370 cas = 0; 1371 dim = (nEl1 > nEl) ? p0->Dim() : p1->Dim(); 1372 nElMin = (nEl1 > nEl) ? nEl : nEl1; 1373 } else if (p0->Rank() != 0 && p1->Rank() == 0) { 1374 cas = 1; 1375 dim = p0->Dim(); 1376 nElMin = nEl; 1377 } else if (p0->Rank() == 0 && p1->Rank() != 0) { 1378 cas = 2; 1379 dim = p1->Dim(); 1380 nElMin = nEl1; 1381 } 1382 1383 if (t == GDL_DOUBLE) { 1384 DDoubleGDL* p0D; 1385 DDoubleGDL* p1D; 1386 if (p0->Type() == GDL_DOUBLE) p0D = static_cast<DDoubleGDL*> (p0); 1387 else { 1388 p0D = static_cast<DDoubleGDL*> (p0->Convert2(GDL_DOUBLE, BaseGDL::COPY)); 1389 guardp0D.Reset(p0D); 1390 } 1391 if (p1->Type() == GDL_DOUBLE) p1D = static_cast<DDoubleGDL*> (p1); 1392 else { 1393 p1D = static_cast<DDoubleGDL*> (p1->Convert2(GDL_DOUBLE, BaseGDL::COPY)); 1394 guardp1D.Reset(p1D); 1395 } 1396 DDoubleGDL* res = new DDoubleGDL(dim, BaseGDL::NOZERO); 1397 if (nElMin == 1) { 1398 (*res)[ 0] = atan2((*p0D)[0], (*p1D)[0]); 1399 return res; 1400 } 1401 switch (cas) { 1402 case 0: 1403 #pragma omp parallel for if (nElMin >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nElMin)) 1404 for (SizeT i = 0; i < nElMin; ++i) (*res)[i] = atan2((*p0D)[i], (*p1D)[i]); 1405 return res; 1406 case 1: 1407 #pragma omp parallel for if (nElMin >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nElMin)) 1408 for (SizeT i = 0; i < nElMin; ++i) (*res)[i] = atan2((*p0D)[i], (*p1D)[0]); 1409 return res; 1410 case 2: 1411 #pragma omp parallel for if (nElMin >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nElMin)) 1412 for (SizeT i = 0; i < nElMin; ++i) (*res)[i] = atan2((*p0D)[0], (*p1D)[i]); 1413 return res; 1414 } 1415 } else if (t == GDL_FLOAT) { 1416 DFloatGDL* p0F; 1417 DFloatGDL* p1F; 1418 if (p0->Type() == GDL_FLOAT) p0F = static_cast<DFloatGDL*> (p0); 1419 else { 1420 p0F = static_cast<DFloatGDL*> (p0->Convert2(GDL_FLOAT, BaseGDL::COPY)); 1421 guardp0F.Reset(p0F); 1422 } 1423 if (p1->Type() == GDL_FLOAT) p1F = static_cast<DFloatGDL*> (p1); 1424 else { 1425 p1F = static_cast<DFloatGDL*> (p1->Convert2(GDL_FLOAT, BaseGDL::COPY)); 1426 guardp1F.Reset(p1F); 1427 } 1428 DFloatGDL* res = new DFloatGDL(dim, BaseGDL::NOZERO); 1429 if (nElMin == 1) { 1430 (*res)[ 0] = atan2((*p0F)[0], (*p1F)[0]); 1431 return res; 1432 } 1433 switch (cas) { 1434 case 0: 1435 #pragma omp parallel for if (nElMin >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nElMin)) 1436 for (SizeT i = 0; i < nElMin; ++i) (*res)[i] = atan2((*p0F)[i], (*p1F)[i]); 1437 return res; 1438 case 1: 1439 #pragma omp parallel for if (nElMin >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nElMin)) 1440 for (SizeT i = 0; i < nElMin; ++i) (*res)[i] = atan2((*p0F)[i], (*p1F)[0]); 1441 return res; 1442 case 2: 1443 #pragma omp parallel for if (nElMin >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nElMin)) 1444 for (SizeT i = 0; i < nElMin; ++i) (*res)[i] = atan2((*p0F)[0], (*p1F)[i]); 1445 return res; 1446 } 1447 } else if (t == GDL_COMPLEX) { 1448 DComplexGDL* p0C; 1449 DComplexGDL* p1C; 1450 if (p0->Type() == GDL_COMPLEX) p0C = static_cast<DComplexGDL*> (p0); 1451 else { 1452 p0C = static_cast<DComplexGDL*> (p0->Convert2(GDL_COMPLEX, BaseGDL::COPY)); 1453 guardp0C.reset(p0C); 1454 } 1455 if (p1->Type() == GDL_COMPLEX) p1C = static_cast<DComplexGDL*> (p1); 1456 else { 1457 p1C = static_cast<DComplexGDL*> (p1->Convert2(GDL_COMPLEX, BaseGDL::COPY)); 1458 guardp1C.Reset(p1C); 1459 } 1460 DComplexGDL* res = new DComplexGDL(dim, BaseGDL::NOZERO); 1461 if (nElMin == 1) { 1462 (*res)[ 0] = atanC((*p0C)[0], (*p1C)[0]); 1463 return res; 1464 } 1465 switch (cas) { 1466 case 0: 1467 #pragma omp parallel for if (nElMin >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nElMin)) 1468 for (SizeT i = 0; i < nElMin; ++i) (*res)[i] = atanC((*p0C)[i], (*p1C)[i]); 1469 return res; 1470 case 1: 1471 #pragma omp parallel for if (nElMin >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nElMin)) 1472 for (SizeT i = 0; i < nElMin; ++i) (*res)[i] = atanC((*p0C)[i], (*p1C)[0]); 1473 return res; 1474 case 2: 1475 #pragma omp parallel for if (nElMin >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nElMin)) 1476 for (SizeT i = 0; i < nElMin; ++i) (*res)[i] = atanC((*p0C)[0], (*p1C)[i]); 1477 return res; 1478 } 1479 } else if (t == GDL_COMPLEXDBL) { 1480 DComplexDblGDL* p0DC; 1481 DComplexDblGDL* p1DC; 1482 if (p0->Type() == GDL_COMPLEXDBL) p0DC = static_cast<DComplexDblGDL*> (p0); 1483 else { 1484 p0DC = static_cast<DComplexDblGDL*> (p0->Convert2(GDL_COMPLEXDBL, BaseGDL::COPY)); 1485 guardp0DC.Reset(p0DC); 1486 } 1487 if (p1->Type() == GDL_COMPLEXDBL) p1DC = static_cast<DComplexDblGDL*> (p1); 1488 else { 1489 p1DC = static_cast<DComplexDblGDL*> (p1->Convert2(GDL_COMPLEXDBL, BaseGDL::COPY)); 1490 guardp1DC.Reset(p1DC); 1491 } 1492 DComplexDblGDL* res = new DComplexDblGDL(dim, BaseGDL::NOZERO); 1493 if (nElMin == 1) { 1494 (*res)[ 0] = atanC((*p0DC)[0], (*p1DC)[0]); 1495 return res; 1496 } 1497 switch (cas) { 1498 case 0: 1499 #pragma omp parallel for if (nElMin >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nElMin)) 1500 for (SizeT i = 0; i < nElMin; ++i) (*res)[i] = atanC((*p0DC)[i], (*p1DC)[i]); 1501 return res; 1502 case 1: 1503 #pragma omp parallel for if (nElMin >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nElMin)) 1504 for (SizeT i = 0; i < nElMin; ++i) (*res)[i] = atanC((*p0DC)[i], (*p1DC)[0]); 1505 return res; 1506 case 2: 1507 #pragma omp parallel for if (nElMin >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nElMin)) 1508 for (SizeT i = 0; i < nElMin; ++i) (*res)[i] = atanC((*p0DC)[0], (*p1DC)[i]); 1509 return res; 1510 } 1511 } else { 1512 DDoubleGDL* p0D = static_cast<DDoubleGDL*> (p0->Convert2(GDL_DOUBLE, BaseGDL::COPY)); 1513 guardp0D.Reset(p0D); 1514 DDoubleGDL* p1D = static_cast<DDoubleGDL*> (p1->Convert2(GDL_DOUBLE, BaseGDL::COPY)); 1515 guardp1D.Reset(p1D); 1516 DFloatGDL* res = new DFloatGDL(dim, BaseGDL::NOZERO); 1517 if (nElMin == 1) { 1518 (*res)[ 0] = atan2((*p0D)[0], (*p1D)[0]); 1519 return res; 1520 } 1521 switch (cas) { 1522 case 0: 1523 #pragma omp parallel for if (nElMin >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nElMin)) 1524 for (SizeT i = 0; i < nElMin; ++i) (*res)[i] = atan2((*p0D)[i], (*p1D)[i]); 1525 return res; 1526 case 1: 1527 #pragma omp parallel for if (nElMin >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nElMin)) 1528 for (SizeT i = 0; i < nElMin; ++i) (*res)[i] = atan2((*p0D)[i], (*p1D)[0]); 1529 return res; 1530 case 2: 1531 #pragma omp parallel for if (nElMin >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nElMin)) 1532 for (SizeT i = 0; i < nElMin; ++i) (*res)[i] = atan2((*p0D)[0], (*p1D)[i]); 1533 return res; 1534 } 1535 } 1536 } else { 1537 static int phaseIx = e->KeywordIx("PHASE"); 1538 if (e->KeywordSet(phaseIx) && (p0->Type() == GDL_COMPLEX || p0->Type() == GDL_COMPLEXDBL)) { //special for phase 1539 if (p0->Type() == GDL_COMPLEX) { 1540 DComplexGDL* p0C = static_cast<DComplexGDL*> (p0); 1541 DFloatGDL* res = new DFloatGDL(p0C->Dim(), BaseGDL::NOZERO); 1542 if (nEl == 1) { 1543 (*res)[ 0] = atan2(((*p0C)[ 0]).imag(), ((*p0C)[ 0]).real()); 1544 return res; 1545 } 1546 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1547 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = atan2(((*p0C)[i]).imag(), ((*p0C)[i]).real()); 1548 return res; 1549 } else if (p0->Type() == GDL_COMPLEXDBL) { 1550 DComplexDblGDL* p0C = static_cast<DComplexDblGDL*> (p0); 1551 DDoubleGDL* res = new DDoubleGDL(p0C->Dim(), BaseGDL::NOZERO); 1552 if (nEl == 1) { 1553 (*res)[ 0] = atan2(((*p0C)[ 0]).imag(), ((*p0C)[ 0]).real()); 1554 return res; 1555 } 1556 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1557 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = atan2(((*p0C)[i]).imag(), ((*p0C)[i]).real()); 1558 return res; 1559 } 1560 } else { 1561 if (p0->Type() == GDL_DOUBLE) { 1562 DDoubleGDL* p0D = static_cast<DDoubleGDL*> (p0); 1563 DDoubleGDL* res = new DDoubleGDL(p0D->Dim(), BaseGDL::NOZERO); 1564 if (nEl == 1) { 1565 (*res)[ 0] = atan((*p0D)[ 0]); 1566 return res; 1567 } 1568 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1569 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = atan((*p0D)[i]); 1570 return res; 1571 } else if (p0->Type() == GDL_FLOAT) { 1572 DFloatGDL* p0F = static_cast<DFloatGDL*> (p0); 1573 DFloatGDL* res = new DFloatGDL(p0F->Dim(), BaseGDL::NOZERO); 1574 if (nEl == 1) { 1575 (*res)[ 0] = atan((*p0F)[ 0]); 1576 return res; 1577 } 1578 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1579 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = atan((*p0F)[i]); 1580 return res; 1581 } else if (p0->Type() == GDL_COMPLEX) { 1582 DComplexGDL* p0C = static_cast<DComplexGDL*> (p0); 1583 DComplexGDL* res = new DComplexGDL(p0C->Dim(), BaseGDL::NOZERO); 1584 if (nEl == 1) { 1585 (*res)[ 0] = std::atan((*p0C)[ 0]); 1586 return res; 1587 } 1588 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1589 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = std::atan((*p0C)[ i]); 1590 return res; 1591 } else if (p0->Type() == GDL_COMPLEXDBL) { 1592 DComplexDblGDL* p0C = static_cast<DComplexDblGDL*> (p0); 1593 DComplexDblGDL* res = new DComplexDblGDL(p0C->Dim(), BaseGDL::NOZERO); 1594 if (nEl == 1) { 1595 (*res)[ 0] = std::atan((*p0C)[ 0]); 1596 return res; 1597 } 1598 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1599 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = std::atan((*p0C)[ i]); 1600 return res; 1601 } else { 1602 DFloatGDL* res = static_cast<DFloatGDL*> (p0->Convert2(GDL_FLOAT, BaseGDL::COPY)); //use same allocation 1603 if (nEl == 1) { 1604 (*res)[ 0] = atan((*res)[ 0]); 1605 return res; 1606 } 1607 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) 1608 for (SizeT i = 0; i < nEl; ++i) (*res)[ i] = atan((*res)[i]); 1609 return res; 1610 } 1611 } 1612 } 1613 assert(false); 1614 return NULL; //neverreached end of non-void function 1615 } 1616 1617 // by medericboquien@users.sourceforge.net gauss_pdf(EnvT * e)1618 BaseGDL* gauss_pdf(EnvT* e) 1619 { 1620 // SizeT nParam = e->NParam(1); 1621 DDoubleGDL* v = static_cast<DDoubleGDL*>(e->GetParDefined(0)-> 1622 Convert2(GDL_DOUBLE,BaseGDL::COPY)); 1623 SizeT nv = v->N_Elements(); 1624 1625 for (int count = 0;count < nv;++count) 1626 (*v)[count] = gsl_cdf_ugaussian_P((*v)[count]); 1627 1628 if (e->GetParDefined(0)->Type() == GDL_DOUBLE) 1629 return v; 1630 else 1631 return v->Convert2(GDL_FLOAT,BaseGDL::CONVERT); 1632 return new DByteGDL(0); 1633 } 1634 1635 // by medericboquien@users.sourceforge.net gauss_cvf(EnvT * e)1636 BaseGDL* gauss_cvf(EnvT* e) 1637 { 1638 // SizeT nParam = e->NParam(1); 1639 DDoubleGDL* p = static_cast<DDoubleGDL*>(e->GetParDefined(0)-> 1640 Convert2(GDL_DOUBLE,BaseGDL::COPY)); 1641 1642 if (p->N_Elements() != 1) 1643 e->Throw("Parameter must be scalar or one element array: "+ 1644 e->GetParString(0)); 1645 if ((*p)[0] < 0. || (*p)[0] > 1.) 1646 e->Throw("Parameter must be in [0,1]: "+e->GetParString(0)); 1647 1648 (*p)[0] = gsl_cdf_ugaussian_Qinv((*p)[0]); 1649 1650 if (e->GetParDefined(0)->Type() == GDL_DOUBLE) 1651 return p; 1652 else 1653 return p->Convert2(GDL_FLOAT,BaseGDL::CONVERT); 1654 return new DByteGDL(0); 1655 } 1656 1657 // by medericboquien@users.sourceforge.net t_pdf(EnvT * e)1658 BaseGDL* t_pdf(EnvT* e) 1659 { 1660 // SizeT nParam = e->NParam(2); 1661 DDoubleGDL* v = e->GetParAs<DDoubleGDL>(0); 1662 DDoubleGDL* df = e->GetParAs<DDoubleGDL>(1); 1663 DDoubleGDL* res; 1664 1665 SizeT nv = v->N_Elements(); 1666 SizeT ndf = df->N_Elements(); 1667 1668 for (int i=0;i<ndf;++i) 1669 if ((*df)[i] <= 0.) 1670 e->Throw("Degrees of freedom must be positive."); 1671 1672 if (nv == 1 && ndf == 1) { 1673 res = new DDoubleGDL(dimension(1), BaseGDL::NOZERO); 1674 (*res)[0] = gsl_cdf_tdist_P((*v)[0],(*df)[0]); 1675 } else if (nv > 1 && ndf == 1) { 1676 res = new DDoubleGDL(dimension(nv), BaseGDL::NOZERO); 1677 for (SizeT count = 0; count < nv; ++count) 1678 (*res)[count] = gsl_cdf_tdist_P((*v)[count],(*df)[0]); 1679 } else if (nv == 1 && ndf > 1) { 1680 res = new DDoubleGDL(dimension(ndf), BaseGDL::NOZERO); 1681 for (SizeT count = 0; count < ndf; ++count) 1682 (*res)[count] = gsl_cdf_tdist_P((*v)[0],(*df)[count]); 1683 } else { 1684 SizeT nreturn = nv>ndf?ndf:nv; 1685 res = new DDoubleGDL(dimension(nreturn), BaseGDL::NOZERO); 1686 for (SizeT count = 0; count < nreturn; ++count) 1687 (*res)[count] = gsl_cdf_tdist_P((*v)[count],(*df)[count]); 1688 } 1689 1690 if (e->GetParDefined(0)->Type() != GDL_DOUBLE && e->GetParDefined(0)->Type() != GDL_DOUBLE) 1691 return res->Convert2(GDL_FLOAT,BaseGDL::CONVERT); 1692 else 1693 return res; 1694 return new DByteGDL(0); 1695 } 1696 1697 // by medericboquien@users.sourceforge.net 1698 laguerre(EnvT * e)1699 BaseGDL* laguerre(EnvT* e) 1700 { 1701 SizeT nParam = e->NParam(2); 1702 1703 DDoubleGDL* xvals = e->GetParAs<DDoubleGDL>(0); 1704 if (e->GetParDefined(0)->Type() == GDL_COMPLEX || e->GetParDefined(0)->Type() == GDL_COMPLEXDBL) 1705 e->Throw("Complex Laguerre not implemented: "); 1706 1707 DIntGDL* nval = e->GetParAs<DIntGDL>(1); 1708 if (nval->N_Elements() != 1) 1709 e->Throw("N and K must be scalars."); 1710 if ((*nval)[0] < 0) 1711 e->Throw("Argument N must be greater than or equal to zero."); 1712 1713 DDoubleGDL* kval; 1714 Guard<DDoubleGDL> kval_guard; 1715 if (nParam > 2) { 1716 kval = e->GetParAs<DDoubleGDL>(2); 1717 if (kval->N_Elements() != 1) 1718 e->Throw("N and K must be scalars."); 1719 if ((*kval)[0] < 0.) 1720 e->Throw("Argument K must be greater than or equal to zero."); 1721 } else { 1722 kval = new DDoubleGDL(0); 1723 kval_guard.Reset(kval); 1724 } 1725 1726 DDoubleGDL* res = new DDoubleGDL(xvals->Dim(), BaseGDL::NOZERO); 1727 DDouble k = (*kval)[0]; 1728 DInt n = (*nval)[0]; 1729 SizeT nEx = xvals->N_Elements(); 1730 1731 //is this really useful? parallelizing Laguerre? Do not forget openmp uses time. 1732 //#pragma omp parallel for if (nEx >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEx)) 1733 for (SizeT count = 0; count < nEx; ++count) (*res)[count] = gsl_sf_laguerre_n(n, k, (*xvals)[count]); 1734 1735 static DInt doubleKWIx = e->KeywordIx("DOUBLE"); 1736 static DInt coefKWIx = e->KeywordIx("COEFFICIENTS"); 1737 1738 if (e->KeywordPresent(coefKWIx)) { 1739 double gamma_kn1 = gsl_sf_gamma(k + n + 1.); 1740 DDoubleGDL* coefKW = new DDoubleGDL(dimension(n + 1), BaseGDL::NOZERO); 1741 1742 //GD: parallelizing this complicated loop cannot be done just like that. 1743 //I further doubt Laguerre is going to be used on large arrays. 1744 //#pragma omp parallel for if (n >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= n)) 1745 for (SizeT count = 0; count <= n; ++count) { 1746 double dcount = static_cast<double> (count); 1747 (*coefKW)[count] = ((count & 0x0001) ? -1.0 : 1.0) * gamma_kn1 / 1748 (gsl_sf_gamma(n - dcount + 1.) * gsl_sf_gamma(k + dcount + 1.) * 1749 gsl_sf_gamma(dcount + 1.)); 1750 } 1751 1752 if (e->GetParDefined(0)->Type() != GDL_DOUBLE && !e->KeywordSet(doubleKWIx)) 1753 coefKW = static_cast<DDoubleGDL*> (coefKW-> 1754 Convert2(GDL_FLOAT, BaseGDL::CONVERT)); 1755 e->SetKW(coefKWIx, coefKW); 1756 } 1757 1758 //convert things back 1759 if (e->GetParDefined(0)->Type() != GDL_DOUBLE && !e->KeywordSet(doubleKWIx)) 1760 return res->Convert2(GDL_FLOAT, BaseGDL::CONVERT); 1761 else 1762 return res; 1763 1764 } 1765 1766 // SA: based on equations 5-5 & 5-6 from Snyder (1987) USGS report no 1395 (page 31) 1767 // available for download at: http://pubs.er.usgs.gov/djvu/PP/pp_1395.djvu ll_arc_distance_helper(T c,T Az,T phi1,T l0,T & phi,T & l,bool degrees)1768 template <typename T> inline void ll_arc_distance_helper( 1769 T c, T Az, T phi1, T l0, T& phi, T& l, bool degrees) 1770 { 1771 // temporary variables 1772 T pi = 4 * atan((T)1.), 1773 dtor = degrees ? pi / 180. : 1, 1774 sin_c = sin(c), 1775 cos_c = cos(c), 1776 cos_Az = cos(Az * dtor), 1777 sin_phi1 = sin(phi1 * dtor), 1778 cos_phi1 = cos(phi1 * dtor); 1779 // computing the results 1780 phi = asin(sin_phi1 * cos_c + cos_phi1 * sin_c * cos_Az) / dtor; 1781 l = l0 * dtor + atan2( 1782 sin_c * sin(Az * dtor), (cos_phi1 * cos_c - sin_phi1 * sin_c * cos_Az) 1783 ); 1784 // placing the result in (-pi, pi) 1785 while (l < -pi) l += 2 * pi; 1786 while (l > pi) l -= 2 * pi; 1787 // converting to degrees if needed 1788 l /= dtor; 1789 } ll_arc_distance(EnvT * e)1790 BaseGDL* ll_arc_distance(EnvT* e) 1791 { 1792 // sanity check (for number of parameters) 1793 // SizeT nParam = e->NParam(); 1794 1795 // 1-st argument : longitude/latitude values pair (in radians unless DEGREE kw. present) 1796 BaseGDL* p0 = e->GetNumericParDefined(0); 1797 1798 // 2-nd argument : arc distance (in radians regardless of DEGREE kw. presence) 1799 BaseGDL* p1 = e->GetNumericParDefined(1); 1800 if (p1->N_Elements() != 1) 1801 e->Throw("second argument is expected to be a scalar or 1-element array"); 1802 1803 // 3-rd argument : azimuth (in radians unless DEGREE kw. present) 1804 BaseGDL* p2 = e->GetNumericParDefined(2); 1805 if (p2->N_Elements() != 1) 1806 e->Throw("third argument is expected to be a scalar or 1-element array"); 1807 1808 // chosing a type for the return value 1809 bool args_complexdbl = 1810 (p0->Type() == GDL_COMPLEXDBL || p1->Type() == GDL_COMPLEXDBL || p2->Type() == GDL_COMPLEXDBL); 1811 bool args_complex = args_complexdbl ? false : 1812 (p0->Type() == GDL_COMPLEX || p1->Type() == GDL_COMPLEX || p2->Type() == GDL_COMPLEX); 1813 DType type = ( 1814 p0->Type() == GDL_DOUBLE || p1->Type() == GDL_DOUBLE || p2->Type() == GDL_DOUBLE || args_complexdbl 1815 ) ? GDL_DOUBLE : GDL_FLOAT; 1816 1817 // converting datatypes if neccesarry 1818 if (p0->Type() != type) p0 = p0->Convert2(type, BaseGDL::COPY); 1819 if (p1->Type() != type) p1 = p1->Convert2(type, BaseGDL::COPY); 1820 if (p2->Type() != type) p2 = p2->Convert2(type, BaseGDL::COPY); 1821 1822 // calculating (by calling a helper template function for float/double versions) 1823 BaseGDL* rt = p0->New(dimension(2, BaseGDL::NOZERO)); 1824 if (type == GDL_FLOAT) 1825 { 1826 ll_arc_distance_helper( 1827 (*static_cast<DFloatGDL*>(p1))[0], 1828 (*static_cast<DFloatGDL*>(p2))[0], 1829 (*static_cast<DFloatGDL*>(p0))[1], 1830 (*static_cast<DFloatGDL*>(p0))[0], 1831 (*static_cast<DFloatGDL*>(rt))[1], 1832 (*static_cast<DFloatGDL*>(rt))[0], 1833 e->KeywordSet(0) //DEGREES (sole option) 1834 ); 1835 } 1836 else 1837 { 1838 ll_arc_distance_helper( 1839 (*static_cast<DDoubleGDL*>(p1))[0], 1840 (*static_cast<DDoubleGDL*>(p2))[0], 1841 (*static_cast<DDoubleGDL*>(p0))[1], 1842 (*static_cast<DDoubleGDL*>(p0))[0], 1843 (*static_cast<DDoubleGDL*>(rt))[1], 1844 (*static_cast<DDoubleGDL*>(rt))[0], 1845 e->KeywordSet(0) 1846 ); 1847 } 1848 1849 // handling complex/dcomplex conversion 1850 return rt->Convert2( 1851 args_complexdbl ? GDL_COMPLEXDBL : args_complex ? GDL_COMPLEX : type, 1852 BaseGDL::CONVERT 1853 ); 1854 } 1855 crossp(EnvT * e)1856 BaseGDL* crossp(EnvT* e) 1857 { 1858 BaseGDL* p0 = e->GetNumericParDefined(0); 1859 BaseGDL* p1 = e->GetNumericParDefined(1); 1860 if (p0->N_Elements() != 3 || p1->N_Elements() != 3) 1861 e->Throw("Both arguments must have 3 elements"); 1862 1863 BaseGDL *a, *b, *c; 1864 1865 a = (DTypeOrder[p0->Type()] >= DTypeOrder[p1->Type()] ? p0 : p1)->New(dimension(3), BaseGDL::ZERO); 1866 // a = 0 1867 // .--mem: new a (with the type and shape of the result) 1868 b = p0->CShift(-1)->Convert2(a->Type(), BaseGDL::CONVERT); 1869 // | .--mem: new b 1870 a->Add(b); // | | a = shift(p0, -1) 1871 delete b; // | `--mem: del b 1872 b = p1->CShift(-2)->Convert2(a->Type(), BaseGDL::CONVERT); 1873 // | .--mem: new b 1874 a->Mult(b); // | | a = shift(p0, -1) * shift(p1, -2) 1875 b->Sub(b); // | | b = 0 1876 c = p0->CShift(1)->Convert2(a->Type(), BaseGDL::CONVERT); 1877 // | | .--mem: new c 1878 b->Sub(c); // | | | b = - shift(p0, 1) 1879 delete c; // | | `--mem: del c 1880 c = p1->CShift(2)->Convert2(a->Type(), BaseGDL::CONVERT); 1881 // | | .--mem: new c 1882 b->Mult(c); // | | | b = - shift(p0, 1) * shift(p1, 2) 1883 delete c; // | | `--mem: del c 1884 a->Add(b); // | | a = shift(p0, -1) * shift(p1, -2) - shift(p0, 1) * shift(p1, 2) 1885 delete b; // | `--mem: del b 1886 return a; // `---> 1887 } 1888 1889 1890 // SA: adapted from the GPL-licensed GNU plotutils (plotutils-2.5/ode/specfun.c) 1891 // ----------------------------------------------------------------------------- 1892 template <typename T> inverf(T p)1893 T inverf (T p) /* Inverse Error Function */ 1894 { 1895 /* 1896 * Source: This routine was derived (using f2c) from the Fortran 1897 * subroutine MERFI found in ACM Algorithm 602, obtained from netlib. 1898 * 1899 * MDNRIS code is copyright 1978 by IMSL, Inc. Since MERFI has been 1900 * submitted to netlib, it may be used with the restrictions that it may 1901 * only be used for noncommercial purposes, and that IMSL be acknowledged 1902 * as the copyright-holder of the code. 1903 */ 1904 1905 /* Initialized data */ 1906 static T a1 = -.5751703, a2 = -1.896513, a3 = -.05496261, 1907 b0 = -.113773, b1 = -3.293474, b2 = -2.374996, b3 = -1.187515, 1908 c0 = -.1146666, c1 = -.1314774, c2 = -.2368201, c3 = .05073975, 1909 d0 = -44.27977, d1 = 21.98546, d2 = -7.586103, 1910 e0 = -.05668422, e1 = .3937021, e2 = -.3166501, e3 = .06208963, 1911 f0 = -6.266786, f1 = 4.666263, f2 = -2.962883, 1912 g0 = 1.851159e-4, g1 = -.002028152, g2 = -.1498384, g3 = .01078639, 1913 h0 = .09952975, h1 = .5211733, h2 = -.06888301; 1914 1915 /* Local variables */ 1916 static T a, b, f, w, x, y, z, sigma, z2, sd, wi, sn; 1917 1918 x = p; 1919 1920 /* determine sign of x */ 1921 sigma = (x > 0 ? 1.0 : -1.0); 1922 1923 /* Note: -1.0 < x < 1.0 */ 1924 1925 z = abs(x); 1926 1927 /* z between 0.0 and 0.85, approx. f by a 1928 rational function in z */ 1929 1930 if (z <= 0.85) 1931 { 1932 z2 = z * z; 1933 f = z + z * (b0 + a1 * z2 / (b1 + z2 + a2 / (b2 + z2 + a3 / (b3 + z2)))); 1934 } 1935 else /* z greater than 0.85 */ 1936 { 1937 a = 1.0 - z; 1938 b = z; 1939 1940 /* reduced argument is in (0.85,1.0), obtain the transformed variable */ 1941 1942 w = sqrt(-(T)log(a + a * b)); 1943 1944 if (w >= 4.0) 1945 /* w greater than 4.0, approx. f by a rational function in 1.0 / w */ 1946 { 1947 wi = 1.0 / w; 1948 sn = ((g3 * wi + g2) * wi + g1) * wi; 1949 sd = ((wi + h2) * wi + h1) * wi + h0; 1950 f = w + w * (g0 + sn / sd); 1951 } 1952 else if (w < 4.0 && w > 2.5) 1953 /* w between 2.5 and 4.0, approx. f by a rational function in w */ 1954 { 1955 sn = ((e3 * w + e2) * w + e1) * w; 1956 sd = ((w + f2) * w + f1) * w + f0; 1957 f = w + w * (e0 + sn / sd); 1958 1959 /* w between 1.13222 and 2.5, approx. f by 1960 a rational function in w */ 1961 } 1962 else if (w <= 2.5 && w > 1.13222) 1963 { 1964 sn = ((c3 * w + c2) * w + c1) * w; 1965 sd = ((w + d2) * w + d1) * w + d0; 1966 f = w + w * (c0 + sn / sd); 1967 } 1968 } 1969 y = sigma * f; 1970 1971 return y; 1972 } 1973 // ----------------------------------------------------------------------------- 1974 gdl_erfinv_fun(EnvT * e)1975 BaseGDL* gdl_erfinv_fun(EnvT* e) 1976 { 1977 BaseGDL* p0 = e->GetNumericParDefined(0); 1978 SizeT n = p0->N_Elements(); 1979 static int doubleIx = e->KeywordIx("DOUBLE"); 1980 if (e->KeywordSet(doubleIx) || p0->Type() == GDL_DOUBLE) 1981 { 1982 DDoubleGDL *ret = new DDoubleGDL(dimension(n)), *p0d = e->GetParAs<DDoubleGDL>(0); 1983 while (n != 0) --n, (*ret)[n] = inverf((*p0d)[n]); 1984 return ret; 1985 } 1986 else 1987 { 1988 DFloatGDL *ret = new DFloatGDL(dimension(n)), *p0f = e->GetParAs<DFloatGDL>(0); 1989 while (n != 0) --n, (*ret)[n] = inverf((*p0f)[n]); 1990 return ret; 1991 } 1992 } 1993 1994 } // namespace 1995