|
Open64 (mfef90, whirl2f, and IR tools)
TAG: version-openad; SVN changeset: 916
|
|
Open64 (mfef90, whirl2f, and IR tools)
TAG: version-openad; SVN changeset: 916
|
00001 /* 00002 00003 Copyright (C) 2000, 2001 Silicon Graphics, Inc. All Rights Reserved. 00004 00005 This program is free software; you can redistribute it and/or modify it 00006 under the terms of version 2 of the GNU General Public License as 00007 published by the Free Software Foundation. 00008 00009 This program is distributed in the hope that it would be useful, but 00010 WITHOUT ANY WARRANTY; without even the implied warranty of 00011 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 00012 00013 Further, this software is distributed without any warranty that it is 00014 free of the rightful claim of any third person regarding infringement 00015 or the like. Any license provided herein, whether implied or 00016 otherwise, applies only to this software file. Patent licenses, if 00017 any, provided herein do not apply to combinations of this program with 00018 other software, or any other product whatsoever. 00019 00020 You should have received a copy of the GNU General Public License along 00021 with this program; if not, write the Free Software Foundation, Inc., 59 00022 Temple Place - Suite 330, Boston MA 02111-1307, USA. 00023 00024 Contact information: Silicon Graphics, Inc., 1600 Amphitheatre Pky, 00025 Mountain View, CA 94043, or: 00026 00027 http://www.sgi.com 00028 00029 For further information regarding this notice, see: 00030 00031 http://oss.sgi.com/projects/GenInfo/NoticeExplan 00032 00033 */ 00034 00035 00036 #include <stdio.h> 00037 #if defined _CRAY && !defined _CRAYMPP 00038 00039 /* Native library routines are never used on Cray systems */ 00040 00041 #else /* ! _CRAY */ 00042 00043 #include <signal.h> 00044 #include <setjmp.h> 00045 #include <errno.h> 00046 00047 #include "arith.internal.h" 00048 #include "int64.h" 00049 00050 int ar_rounding_modes = 0xf; /* All rounding modes allowed */ 00051 int ar_underflow_modes = 1<<AR_UNDERFLOW_TO_DENORM; 00052 00053 /* 00054 * Solaris vendor compiler will evaluate false on all the following 00055 * three condtions. Let GNU compiler do so, too. In case of Solaris, 00056 * Use native math library for intrinsic function evaluation 00057 */ 00058 00059 /*#if defined(_CRAYMPP) || defined(__mips) 00060 "ar_power" is not correct here use other one 00061 for __mips 00062 */ 00063 #if defined(_CRAYMPP) 00064 00065 /* Call native, F90-compiled routines to evaluate all functions */ 00066 00067 #if !defined(__mips) 00068 # define NULL 0 00069 #endif 00070 #define IEEE_FLOAT_32 (UNROUNDED_TYPE(AR_Float_IEEE_NR_32)) 00071 #define IEEE_FLOAT_64 (UNROUNDED_TYPE(AR_Float_IEEE_NR_64)) 00072 #define IEEE_FLOAT_128 (UNROUNDED_TYPE(AR_Float_IEEE_NR_128)) 00073 #define IEEE_COMPLEX_32 (UNROUNDED_TYPE(AR_Complex_IEEE_NR_32)) 00074 #define IEEE_COMPLEX_64 (UNROUNDED_TYPE(AR_Complex_IEEE_NR_64)) 00075 #define IEEE_COMPLEX_128 (UNROUNDED_TYPE(AR_Complex_IEEE_NR_128)) 00076 00077 int status; 00078 00079 static int ar_native1(void (function)(), ar_data *result, const AR_TYPE *resulttype, const ar_data *opnd); 00080 00081 static int ar_native2(void (function)(), ar_data *result, const AR_TYPE *resulttype, const ar_data *opnd1, const ar_data *opnd2); 00082 00083 /* Fortran character index */ 00084 int 00085 ar_index (ar_data *result, const AR_TYPE *resulttype, 00086 const char *str1, long len1, const char *str2, long len2, long backward) 00087 { 00088 long index; 00089 long back = backward; 00090 int status; 00091 00092 extern long _F90_INDEX(const char *st1, const char *st2, long *back, 00093 int len1, int len2); 00094 00095 index = _F90_INDEX(str1, str2, &back, len1, len2); 00096 result->ar_i64.part1 = 0; 00097 result->ar_i64.part2 = 0; 00098 result->ar_i64.part3 = index>>16; 00099 result->ar_i64.part4 = index & 0xffff; 00100 00101 status = AR_STAT_OK; 00102 switch (*resulttype) { 00103 case AR_Int_8_S: 00104 if (( INT8_SIGN(result) && !IS_INT8_UPPER_ONES(result)) || 00105 (!INT8_SIGN(result) && !IS_INT8_UPPER_ZERO(result))) { 00106 return AR_STAT_OVERFLOW; 00107 } 00108 break; 00109 00110 case AR_Int_16_S: 00111 if (( INT16_SIGN(result) && !IS_INT16_UPPER_ONES(result)) || 00112 (!INT16_SIGN(result) && !IS_INT16_UPPER_ZERO(result))) { 00113 return AR_STAT_OVERFLOW; 00114 } 00115 break; 00116 00117 case AR_Int_32_S: 00118 if (( INT32_SIGN(result) && !IS_INT32_UPPER_ONES(result)) || 00119 (!INT32_SIGN(result) && !IS_INT32_UPPER_ZERO(result))) { 00120 return AR_STAT_OVERFLOW; 00121 } 00122 break; 00123 00124 case AR_Int_64_S: 00125 break; 00126 00127 default: 00128 return AR_STAT_INVALID_TYPE; 00129 } 00130 WORD_SWAP(result->ar_i64); 00131 return AR_status((AR_DATA*)result, resulttype); 00132 } 00133 00134 00135 /* Fortran character scan */ 00136 int 00137 ar_scan (ar_data *result, const AR_TYPE *resulttype, 00138 const char *str1, long len1, const char *str2, long len2, long backward) 00139 { 00140 long index; 00141 long back = backward; 00142 int status; 00143 00144 extern long _F90_SCAN(const char *st1, const char *st2, long *back, 00145 int len1, int len2); 00146 00147 00148 index = _F90_SCAN(str1, str2, &back, len1, len2); 00149 result->ar_i64.part1 = 0; 00150 result->ar_i64.part2 = 0; 00151 result->ar_i64.part3 = index>>16; 00152 result->ar_i64.part4 = index & 0xffff; 00153 00154 status = AR_STAT_OK; 00155 switch (*resulttype) { 00156 case AR_Int_8_S: 00157 if (( INT8_SIGN(result) && !IS_INT8_UPPER_ONES(result)) || 00158 (!INT8_SIGN(result) && !IS_INT8_UPPER_ZERO(result))) { 00159 return AR_STAT_OVERFLOW; 00160 } 00161 break; 00162 00163 case AR_Int_16_S: 00164 if (( INT16_SIGN(result) && !IS_INT16_UPPER_ONES(result)) || 00165 (!INT16_SIGN(result) && !IS_INT16_UPPER_ZERO(result))) { 00166 return AR_STAT_OVERFLOW; 00167 } 00168 break; 00169 00170 case AR_Int_32_S: 00171 if (( INT32_SIGN(result) && !IS_INT32_UPPER_ONES(result)) || 00172 (!INT32_SIGN(result) && !IS_INT32_UPPER_ZERO(result))) { 00173 return AR_STAT_OVERFLOW; 00174 } 00175 break; 00176 00177 case AR_Int_64_S: 00178 break; 00179 00180 default: 00181 return AR_STAT_INVALID_TYPE; 00182 } 00183 00184 WORD_SWAP(result->ar_i64); 00185 return AR_status((AR_DATA*)result, resulttype); 00186 } 00187 00188 00189 /* Fortran character verify */ 00190 int 00191 ar_verify (ar_data *result, const AR_TYPE *resulttype, 00192 const char *str1, long len1, const char *str2, long len2, long backward) 00193 { 00194 long index; 00195 long back = backward; 00196 00197 extern long _F90_VERIFY(const char *st1, const char *st2, long *back, 00198 int len1, int len2); 00199 00200 index = _F90_VERIFY(str1, str2, &back, len1, len2); 00201 result->ar_i64.part1 = 0; 00202 result->ar_i64.part2 = 0; 00203 result->ar_i64.part3 = index>>16; 00204 result->ar_i64.part4 = index & 0xffff; 00205 00206 status = AR_STAT_OK; 00207 switch (*resulttype) { 00208 case AR_Int_8_S: 00209 if (( INT8_SIGN(result) && !IS_INT8_UPPER_ONES(result)) || 00210 (!INT8_SIGN(result) && !IS_INT8_UPPER_ZERO(result))) { 00211 return AR_STAT_OVERFLOW; 00212 } 00213 break; 00214 00215 case AR_Int_16_S: 00216 if (( INT16_SIGN(result) && !IS_INT16_UPPER_ONES(result)) || 00217 (!INT16_SIGN(result) && !IS_INT16_UPPER_ZERO(result))) { 00218 return AR_STAT_OVERFLOW; 00219 } 00220 break; 00221 00222 case AR_Int_32_S: 00223 if (( INT32_SIGN(result) && !IS_INT32_UPPER_ONES(result)) || 00224 (!INT32_SIGN(result) && !IS_INT32_UPPER_ZERO(result))) { 00225 return AR_STAT_OVERFLOW; 00226 } 00227 break; 00228 00229 case AR_Int_64_S: 00230 break; 00231 00232 default: 00233 return AR_STAT_INVALID_TYPE; 00234 } 00235 00236 WORD_SWAP(result->ar_i64); 00237 return AR_status((AR_DATA*)result, resulttype); 00238 } 00239 00240 00241 /* Fortran-90 reshape */ 00242 int 00243 ar_reshape (void *result, const void *source, const void *shape, 00244 const void *pad, const void *order) 00245 { 00246 extern void _RESHAPE(); 00247 00248 if(*(char**)source == NULL || *(char**)shape == NULL) 00249 return AR_STAT_UNDEFINED; 00250 00251 _RESHAPE(result, source, shape, pad, order); 00252 00253 return AR_STAT_OK; 00254 } 00255 00256 00257 /* Fortran-90 transfer */ 00258 int 00259 ar_transfer (void *result, const void *source, const void *mold, long *length) 00260 { 00261 long size; 00262 00263 extern void _TRANSFER(); 00264 00265 if(*(char**)source == NULL || *(char**)mold == NULL) 00266 return AR_STAT_UNDEFINED; 00267 00268 if(length != NULL) { 00269 size = *length; 00270 _TRANSFER(result, source, mold, &size); 00271 } 00272 else 00273 _TRANSFER(result, source, mold, (int*)NULL); 00274 00275 return AR_STAT_OK; 00276 } 00277 00278 00279 /* Fortran-90 modulo */ 00280 int 00281 ar_modulo (ar_data *result, const AR_TYPE *resulttype, 00282 const ar_data *opnd1, const AR_TYPE *opnd1type, 00283 const ar_data *opnd2, const AR_TYPE *opnd2type) 00284 { 00285 ar_data result_32; 00286 ar_data opnd1_32; 00287 ar_data opnd2_32; 00288 AR_TYPE int32type = AR_Int_32_S; 00289 int status; 00290 00291 #if !defined _CRAYMPP 00292 #define ARMOD armod_ 00293 #define ARMODD armodd_ 00294 #define ARMODXD armodxd_ 00295 #define ARMODI armodi_ 00296 #define ARMODJ armodj_ 00297 #endif 00298 extern void ARMOD (ar_data *res, 00299 const ar_data *arg1, const ar_data *arg2); 00300 extern void ARMODD (ar_data *res, 00301 const ar_data *arg1, const ar_data *arg2); 00302 extern void ARMODXD(ar_data *res, 00303 const ar_data *arg1, const ar_data *arg2); 00304 extern void ARMODI (ar_data *res, 00305 const ar_data *arg1, const ar_data *arg2); 00306 extern void ARMODJ (ar_data *res, 00307 const ar_data *arg1, const ar_data *arg2); 00308 00309 switch (UNROUNDED_TYPE(*resulttype)) { 00310 00311 case IEEE_FLOAT_32: 00312 return ar_native2(ARMOD, result, resulttype, opnd1, opnd2); 00313 00314 case IEEE_FLOAT_64: 00315 return ar_native2(ARMODD, result, resulttype, opnd1, opnd2); 00316 00317 case IEEE_FLOAT_128: 00318 return ar_native2(ARMODXD, result, resulttype, opnd1, opnd2); 00319 00320 default: 00321 switch (*resulttype) { 00322 00323 case AR_Int_8_S: 00324 case AR_Int_16_S: 00325 ZERO_INT32_ALL(&opnd1_32); 00326 ZERO_INT32_ALL(&opnd2_32); 00327 ZERO_INT32_ALL(&result_32); 00328 (void) AR_convert((AR_DATA*) &opnd1_32, &int32type, 00329 (AR_DATA*)opnd1, resulttype); 00330 (void) AR_convert((AR_DATA*) &opnd2_32, &int32type, 00331 (AR_DATA*)opnd2, resulttype); 00332 status = ar_native2(ARMODI, 00333 &result_32, &int32type, 00334 &opnd1_32, &opnd2_32); 00335 if (IS_ERROR_STATUS(status)) 00336 return status; 00337 return AR_convert((AR_DATA*) result, resulttype, 00338 (AR_DATA*) &result_32, &int32type); 00339 00340 case AR_Int_32_S: 00341 return ar_native2(ARMODI, 00342 result, resulttype, opnd1, opnd2); 00343 00344 case AR_Int_64_S: 00345 return ar_native2(ARMODJ, 00346 result, resulttype, opnd1, opnd2); 00347 00348 default: 00349 return AR_STAT_INVALID_TYPE; 00350 } 00351 00352 } 00353 } 00354 00355 00356 /* Fortran-90 selected_real_kind */ 00357 int 00358 ar_selected_real_kind (ar_data *result, const AR_TYPE *resulttype, 00359 const ar_data *opnd1, const AR_TYPE *opnd1type, 00360 const ar_data *opnd2, const AR_TYPE *opnd2type) 00361 { 00362 int status; 00363 AR_TYPE int32type = AR_Int_32_S; 00364 00365 #if !defined _CRAYMPP 00366 #define ARSELRK arselrk_ 00367 #endif 00368 extern void ARSELRK(ar_data *res, const ar_data *arg1, const ar_data *arg2); 00369 00370 ZERO_INT32_UPPER(result); 00371 status = ar_native2(ARSELRK, result, &int32type, opnd1, opnd2); 00372 00373 switch (*resulttype) { 00374 00375 case AR_Int_8_S: 00376 case AR_Int_16_S: 00377 if (!IS_ERROR_STATUS(status)) 00378 status = AR_convert((AR_DATA*) result, resulttype, 00379 (AR_DATA*) result, &int32type); 00380 break; 00381 00382 case AR_Int_32_S: 00383 break; 00384 00385 case AR_Int_64_S: 00386 if(result->ar_i64.part3 & 0x8000) 00387 result->ar_i64.part1 = result->ar_i64.part2 = 0xFFFF; 00388 break; 00389 00390 default: 00391 return AR_STAT_INVALID_TYPE; 00392 } 00393 00394 return status; 00395 } 00396 00397 00398 /* Square root */ 00399 int 00400 ar_sqrt (ar_data *result, const AR_TYPE *resulttype, 00401 const ar_data *opnd, const AR_TYPE *opndtype) 00402 { 00403 #if !defined _CRAYMPP 00404 #define ARSQRT arsqrt_ 00405 #define ARDSQRT ardsqrt_ 00406 #define ARXDSQRT arxdsqrt_ 00407 #define ARCSQRT arcsqrt_ 00408 #define ARCDSQRT arcdsqrt_ 00409 #define ARCXDSQRT arcxdsqrt_ 00410 #endif 00411 extern void ARSQRT (ar_data *res, const ar_data *arg); 00412 extern void ARDSQRT (ar_data *res, const ar_data *arg); 00413 extern void ARXDSQRT (ar_data *res, const ar_data *arg); 00414 extern void ARCSQRT (ar_data *res, const ar_data *arg); 00415 extern void ARCDSQRT (ar_data *res, const ar_data *arg); 00416 extern void ARCXDSQRT(ar_data *res, const ar_data *arg); 00417 00418 switch (UNROUNDED_TYPE(*resulttype)) { 00419 00420 case IEEE_FLOAT_32: 00421 return ar_native1(ARSQRT, result, resulttype, opnd); 00422 00423 case IEEE_FLOAT_64: 00424 return ar_native1(ARDSQRT, result, resulttype, opnd); 00425 00426 case IEEE_FLOAT_128: 00427 return ar_native1(ARXDSQRT, result, resulttype, opnd); 00428 00429 case IEEE_COMPLEX_32: 00430 return ar_native1(ARCSQRT, result, resulttype, opnd); 00431 00432 case IEEE_COMPLEX_64: 00433 return ar_native1(ARCDSQRT, result, resulttype, opnd); 00434 00435 case IEEE_COMPLEX_128: 00436 return ar_native1(ARCXDSQRT, result, resulttype, opnd); 00437 00438 default: 00439 return AR_STAT_INVALID_TYPE; 00440 } 00441 } 00442 00443 00444 /* Natural (base "e") logarithm */ 00445 int 00446 ar_log (ar_data *result, const AR_TYPE *resulttype, 00447 const ar_data *opnd, const AR_TYPE *opndtype) 00448 { 00449 #if !defined _CRAYMPP 00450 #define ARLOG arlog_ 00451 #define ARDLOG ardlog_ 00452 #define ARXDLOG arxdlog_ 00453 #define ARCLOG arclog_ 00454 #define ARCDLOG arcdlog_ 00455 #define ARCXDLOG arcxdlog_ 00456 #endif 00457 extern void ARLOG (ar_data *res, const ar_data *arg); 00458 extern void ARDLOG (ar_data *res, const ar_data *arg); 00459 extern void ARXDLOG (ar_data *res, const ar_data *arg); 00460 extern void ARCLOG (ar_data *res, const ar_data *arg); 00461 extern void ARCDLOG (ar_data *res, const ar_data *arg); 00462 extern void ARCXDLOG(ar_data *res, const ar_data *arg); 00463 00464 switch (UNROUNDED_TYPE(*resulttype)) { 00465 00466 case IEEE_FLOAT_32: 00467 return ar_native1(ARLOG, result, resulttype, opnd); 00468 00469 case IEEE_FLOAT_64: 00470 return ar_native1(ARDLOG, result, resulttype, opnd); 00471 00472 case IEEE_FLOAT_128: 00473 return ar_native1(ARXDLOG, result, resulttype, opnd); 00474 00475 case IEEE_COMPLEX_32: 00476 return ar_native1(ARCLOG, result, resulttype, opnd); 00477 00478 case IEEE_COMPLEX_64: 00479 return ar_native1(ARCDLOG, result, resulttype, opnd); 00480 00481 case IEEE_COMPLEX_128: 00482 return ar_native1(ARCXDLOG, result, resulttype, opnd); 00483 00484 default: 00485 return AR_STAT_INVALID_TYPE; 00486 } 00487 } 00488 00489 00490 /* Exponential ("e" ** x) function */ 00491 int 00492 ar_exp (ar_data *result, const AR_TYPE *resulttype, 00493 const ar_data *opnd, const AR_TYPE *opndtype) 00494 { 00495 #if !defined _CRAYMPP 00496 #define AREXP arexp_ 00497 #define ARDEXP ardexp_ 00498 #define ARXDEXP arxdexp_ 00499 #define ARCEXP arcexp_ 00500 #define ARCDEXP arcdexp_ 00501 #define ARCXDEXP arcxdexp_ 00502 #endif 00503 extern void AREXP (ar_data *res, const ar_data *arg); 00504 extern void ARDEXP (ar_data *res, const ar_data *arg); 00505 extern void ARXDEXP (ar_data *res, const ar_data *arg); 00506 extern void ARCEXP (ar_data *res, const ar_data *arg); 00507 extern void ARCDEXP (ar_data *res, const ar_data *arg); 00508 extern void ARCXDEXP(ar_data *res, const ar_data *arg); 00509 00510 switch (UNROUNDED_TYPE(*resulttype)) { 00511 00512 case IEEE_FLOAT_32: 00513 return ar_native1(AREXP, result, resulttype, opnd); 00514 00515 case IEEE_FLOAT_64: 00516 return ar_native1(ARDEXP, result, resulttype, opnd); 00517 00518 case IEEE_FLOAT_128: 00519 return ar_native1(ARXDEXP, result, resulttype, opnd); 00520 00521 case IEEE_COMPLEX_32: 00522 return ar_native1(ARCEXP, result, resulttype, opnd); 00523 00524 case IEEE_COMPLEX_64: 00525 return ar_native1(ARCDEXP, result, resulttype, opnd); 00526 00527 case IEEE_COMPLEX_128: 00528 return ar_native1(ARCXDEXP, result, resulttype, opnd); 00529 00530 default: 00531 return AR_STAT_INVALID_TYPE; 00532 } 00533 } 00534 00535 00536 /* Complex absolute value */ 00537 int 00538 ar_cabs (ar_data *result, const AR_TYPE *resulttype, 00539 const ar_data *opnd, const AR_TYPE *opndtype) 00540 { 00541 #if !defined _CRAYMPP 00542 #define ARCABS arcabs_ 00543 #define ARCDABS arcdabs_ 00544 #define ARCXDABS arcxdabs_ 00545 #endif 00546 extern void ARCABS (ar_data *res, const ar_data *arg); 00547 extern void ARCDABS (ar_data *res, const ar_data *arg); 00548 extern void ARCXDABS(ar_data *res, const ar_data *arg); 00549 00550 if (UNROUNDED_TYPE(*resulttype) == IEEE_FLOAT_32 && 00551 UNROUNDED_TYPE(*opndtype) == IEEE_COMPLEX_32) 00552 return ar_native1(ARCABS, result, resulttype, opnd); 00553 00554 if (UNROUNDED_TYPE(*resulttype) == IEEE_FLOAT_64 && 00555 UNROUNDED_TYPE(*opndtype) == IEEE_COMPLEX_64) 00556 return ar_native1(ARCDABS, result, resulttype, opnd); 00557 00558 if (UNROUNDED_TYPE(*resulttype) == IEEE_FLOAT_128 && 00559 UNROUNDED_TYPE(*opndtype) == IEEE_COMPLEX_128) 00560 return ar_native1(ARCXDABS, result, resulttype, opnd); 00561 00562 return AR_STAT_INVALID_TYPE; 00563 } 00564 00565 /* Exponentiation */ 00566 int 00567 ar_power(ar_data *result, const AR_TYPE *resulttype, 00568 const ar_data *base, const AR_TYPE *basetype, 00569 const ar_data *power, const AR_TYPE *powertype) 00570 { 00571 int status; 00572 ar_data tbase, tpow; 00573 AR_TYPE btype, ptype; 00574 ar_data base_32; 00575 ar_data tbase_32; 00576 ar_data tpow_32; 00577 ar_data result_32; 00578 AR_TYPE int32type = AR_Int_32_S; 00579 00580 #if !defined _CRAYMPP 00581 #define ARPOWGG arpowgg_ 00582 #define ARPOWHH arpowhh_ 00583 #define ARPOWII arpowii_ 00584 #define ARPOWJJ arpowjj_ 00585 #define ARPOWRG arpowrg_ 00586 #define ARPOWRH arpowrh_ 00587 #define ARPOWRI arpowri_ 00588 #define ARPOWRJ arpowrj_ 00589 #define ARPOWDG arpowdg_ 00590 #define ARPOWDH arpowdh_ 00591 #define ARPOWDI arpowdi_ 00592 #define ARPOWDJ arpowdj_ 00593 #define ARPOWXDG arpowxdg_ 00594 #define ARPOWXDH arpowxdh_ 00595 #define ARPOWXDI arpowxdi_ 00596 #define ARPOWXDJ arpowxdj_ 00597 #define ARPOWCG arpowcg_ 00598 #define ARPOWCH arpowch_ 00599 #define ARPOWCI arpowci_ 00600 #define ARPOWCJ arpowcj_ 00601 #define ARPOWGR arpowgr_ 00602 #define ARPOWHR arpowhr_ 00603 #define ARPOWIR arpowir_ 00604 #define ARPOWJR arpowjr_ 00605 #define ARPOWRR arpowrr_ 00606 #define ARPOWDR arpowdr_ 00607 #define ARPOWXDR arpowxdr_ 00608 #define ARPOWCR arpowcr_ 00609 #define ARPOWDD arpowdd_ 00610 #define ARPOWXDXD arpowxdxd_ 00611 #define ARPOWCC arpowcc_ 00612 #define ARPOWCDG arpowcdg_ 00613 #define ARPOWCDH arpowcdh_ 00614 #define ARPOWCDI arpowcdi_ 00615 #define ARPOWCDJ arpowcdj_ 00616 #define ARPOWCXDG arpowcxdg_ 00617 #define ARPOWCXDH arpowcxdh_ 00618 #define ARPOWCXDI arpowcxdi_ 00619 #define ARPOWCXDJ arpowcxdj_ 00620 #define ARPOWCDCD arpowcdcd_ 00621 #define ARPOWCXDCXD arpowcxdcxd_ 00622 #endif 00623 00624 extern void ARPOWGG (ar_data *res, const ar_data *base, 00625 const ar_data *power); 00626 extern void ARPOWHH (ar_data *res, const ar_data *base, 00627 const ar_data *power); 00628 extern void ARPOWII (ar_data *res, const ar_data *base, 00629 const ar_data *power); 00630 extern void ARPOWJJ (ar_data *res, const ar_data *base, 00631 const ar_data *power); 00632 extern void ARPOWRG (ar_data *res, const ar_data *base, 00633 const ar_data *power); 00634 extern void ARPOWRH (ar_data *res, const ar_data *base, 00635 const ar_data *power); 00636 extern void ARPOWRI (ar_data *res, const ar_data *base, 00637 const ar_data *power); 00638 extern void ARPOWRJ (ar_data *res, const ar_data *base, 00639 const ar_data *power); 00640 extern void ARPOWDG (ar_data *res, const ar_data *base, 00641 const ar_data *power); 00642 extern void ARPOWDH (ar_data *res, const ar_data *base, 00643 const ar_data *power); 00644 extern void ARPOWDI (ar_data *res, const ar_data *base, 00645 const ar_data *power); 00646 extern void ARPOWDJ (ar_data *res, const ar_data *base, 00647 const ar_data *power); 00648 extern void ARPOWXDG (ar_data *res, const ar_data *base, 00649 const ar_data *power); 00650 extern void ARPOWXDH (ar_data *res, const ar_data *base, 00651 const ar_data *power); 00652 extern void ARPOWXDI (ar_data *res, const ar_data *base, 00653 const ar_data *power); 00654 extern void ARPOWXDJ (ar_data *res, const ar_data *base, 00655 const ar_data *power); 00656 extern void ARPOWCG (ar_data *res, const ar_data *base, 00657 const ar_data *power); 00658 extern void ARPOWCH (ar_data *res, const ar_data *base, 00659 const ar_data *power); 00660 extern void ARPOWCI (ar_data *res, const ar_data *base, 00661 const ar_data *power); 00662 extern void ARPOWCJ (ar_data *res, const ar_data *base, 00663 const ar_data *power); 00664 extern void ARPOWGR (ar_data *res, const ar_data *base, 00665 const ar_data *power); 00666 extern void ARPOWHR (ar_data *res, const ar_data *base, 00667 const ar_data *power); 00668 extern void ARPOWIR (ar_data *res, const ar_data *base, 00669 const ar_data *power); 00670 extern void ARPOWJR (ar_data *res, const ar_data *base, 00671 const ar_data *power); 00672 extern void ARPOWRR (ar_data *res, const ar_data *base, 00673 const ar_data *power); 00674 extern void ARPOWDR (ar_data *res, const ar_data *base, 00675 const ar_data *power); 00676 extern void ARPOWXDR (ar_data *res, const ar_data *base, 00677 const ar_data *power); 00678 extern void ARPOWCR (ar_data *res, const ar_data *base, 00679 const ar_data *power); 00680 extern void ARPOWDD (ar_data *res, const ar_data *base, 00681 const ar_data *power); 00682 extern void ARPOWXDXD (ar_data *res, const ar_data *base, 00683 const ar_data *power); 00684 extern void ARPOWCC (ar_data *res, const ar_data *base, 00685 const ar_data *power); 00686 extern void ARPOWCDG (ar_data *res, const ar_data *base, 00687 const ar_data *power); 00688 extern void ARPOWCDH (ar_data *res, const ar_data *base, 00689 const ar_data *power); 00690 extern void ARPOWCDI (ar_data *res, const ar_data *base, 00691 const ar_data *power); 00692 extern void ARPOWCDJ (ar_data *res, const ar_data *base, 00693 const ar_data *power); 00694 extern void ARPOWCXDG (ar_data *res, const ar_data *base, 00695 const ar_data *power); 00696 extern void ARPOWCXDH (ar_data *res, const ar_data *base, 00697 const ar_data *power); 00698 extern void ARPOWCXDI (ar_data *res, const ar_data *base, 00699 const ar_data *power); 00700 extern void ARPOWCXDJ (ar_data *res, const ar_data *base, 00701 const ar_data *power); 00702 extern void ARPOWCDCD (ar_data *res, const ar_data *base, 00703 const ar_data *power); 00704 extern void ARPOWCXDCXD (ar_data *res, const ar_data *base, 00705 const ar_data *power); 00706 00707 /* Prepare for power function evalution by converting 00708 * base and power operand types to values supported by 00709 * the native power functions. 00710 */ 00711 00712 if(AR_CLASS(*basetype) == AR_CLASS_INT) { 00713 00714 if(UNROUNDED_TYPE(*powertype) == IEEE_FLOAT_32) { 00715 if(*resulttype == *powertype) { 00716 switch (*basetype) { 00717 case AR_Int_8_S: 00718 case AR_Int_16_S: 00719 ZERO_INT32_ALL(&base_32); 00720 (void) AR_convert((AR_DATA*) &base_32, 00721 &int32type, 00722 (AR_DATA*)base, 00723 basetype); 00724 return ar_native2(ARPOWIR, 00725 result, resulttype, 00726 &base_32, power); 00727 case AR_Int_32_S: 00728 return ar_native2(ARPOWIR, 00729 result, resulttype, 00730 base, power); 00731 case AR_Int_64_S: 00732 return ar_native2(ARPOWJR, 00733 result, resulttype, 00734 base, power); 00735 default: 00736 return AR_STAT_INVALID_TYPE; 00737 } 00738 } 00739 else 00740 return AR_STAT_INVALID_TYPE; 00741 } 00742 00743 btype = ptype = *powertype; 00744 } 00745 else if(AR_CLASS(*powertype) == AR_CLASS_INT || 00746 (AR_FLOAT_SIZE(*powertype) == AR_FLOAT_32 && 00747 AR_FLOAT_IS_COMPLEX(*powertype) != AR_FLOAT_COMPLEX)) { 00748 00749 /* base**I or base**R power functions */ 00750 00751 btype = *basetype; 00752 ptype = *powertype; 00753 } 00754 00755 /* Otherwise, process arg types to simulate power function with 00756 * base type == power type using the greatest precision and/or 00757 * generality required. 00758 */ 00759 00760 else { 00761 00762 /* Convert to base type == power type using the greatest 00763 * precision and/or generality required. 00764 */ 00765 00766 if(AR_FLOAT_SIZE(*basetype) >= AR_FLOAT_SIZE(*powertype)) 00767 btype = (AR_TYPE) 00768 (*basetype | AR_FLOAT_IS_COMPLEX(*powertype)); 00769 else 00770 btype = (AR_TYPE) 00771 (*powertype | AR_FLOAT_IS_COMPLEX(*basetype)); 00772 00773 ptype = btype; 00774 } 00775 00776 /* 00777 * Verify that the resulttype matches the simulated function's 00778 * return type given by the expanded base type. 00779 */ 00780 00781 if(*resulttype != btype) 00782 return AR_STAT_INVALID_TYPE; 00783 00784 /* 00785 * Setup the operands to the power function converting to the 00786 * correct (expanded) type if necessary. 00787 */ 00788 00789 status = AR_STAT_OK; 00790 if(*basetype != btype) 00791 status = AR_convert((AR_DATA*)&tbase, &btype, 00792 (AR_DATA*)base, basetype); 00793 else 00794 tbase = *base; 00795 00796 if(*powertype != ptype) 00797 status = AR_convert((AR_DATA*)&tpow, &ptype, 00798 (AR_DATA*)power, powertype); 00799 else 00800 tpow = *power; 00801 00802 if (IS_ERROR_STATUS(status)) 00803 return status; 00804 00805 /* 00806 * Call the correct native power function determined by 00807 * the (expanded) base and power types. 00808 */ 00809 00810 switch (UNROUNDED_TYPE(btype)) { 00811 00812 case IEEE_FLOAT_32: 00813 if(AR_CLASS(ptype) == AR_CLASS_INT) 00814 switch (ptype) { 00815 case AR_Int_8_S: 00816 case AR_Int_16_S: 00817 ZERO_INT32_ALL(&tpow_32); 00818 (void) AR_convert((AR_DATA*) &tpow_32, 00819 &int32type, 00820 (AR_DATA*) &tpow, 00821 &ptype); 00822 return ar_native2(ARPOWRI, 00823 result, resulttype, 00824 &tbase, &tpow_32); 00825 case AR_Int_32_S: 00826 return ar_native2(ARPOWRI, 00827 result, resulttype, 00828 &tbase, &tpow); 00829 case AR_Int_64_S: 00830 return ar_native2(ARPOWRJ, 00831 result, resulttype, 00832 &tbase, &tpow); 00833 default: 00834 return AR_STAT_INVALID_TYPE; 00835 } 00836 else 00837 return ar_native2(ARPOWRR, 00838 result, resulttype, 00839 &tbase, &tpow); 00840 00841 case IEEE_FLOAT_64: 00842 if(AR_CLASS(ptype) == AR_CLASS_INT) 00843 switch (ptype) { 00844 case AR_Int_8_S: 00845 case AR_Int_16_S: 00846 ZERO_INT32_ALL(&tpow_32); 00847 (void) AR_convert((AR_DATA*) &tpow_32, 00848 &int32type, 00849 (AR_DATA*) &tpow, 00850 &ptype); 00851 return ar_native2(ARPOWDI, 00852 result, resulttype, 00853 &tbase, &tpow_32); 00854 case AR_Int_32_S: 00855 return ar_native2(ARPOWDI, 00856 result, resulttype, 00857 &tbase, &tpow); 00858 case AR_Int_64_S: 00859 return ar_native2(ARPOWDJ, 00860 result, resulttype, 00861 &tbase, &tpow); 00862 default: 00863 return AR_STAT_INVALID_TYPE; 00864 } 00865 else if(AR_FLOAT_SIZE(ptype) == AR_FLOAT_32) 00866 return ar_native2(ARPOWDR, 00867 result, resulttype, 00868 &tbase, &tpow); 00869 else 00870 return ar_native2(ARPOWDD, 00871 result, resulttype, 00872 &tbase, &tpow); 00873 00874 case IEEE_FLOAT_128: 00875 if(AR_CLASS(ptype) == AR_CLASS_INT) 00876 switch (ptype) { 00877 case AR_Int_8_S: 00878 case AR_Int_16_S: 00879 ZERO_INT32_ALL(&tpow_32); 00880 (void) AR_convert((AR_DATA*) &tpow_32, 00881 &int32type, 00882 (AR_DATA*) &tpow, 00883 &ptype); 00884 return ar_native2(ARPOWXDI, 00885 result, resulttype, 00886 &tbase, &tpow_32); 00887 case AR_Int_32_S: 00888 return ar_native2(ARPOWXDI, 00889 result, resulttype, 00890 &tbase, &tpow); 00891 case AR_Int_64_S: 00892 return ar_native2(ARPOWXDJ, 00893 result, resulttype, 00894 &tbase, &tpow); 00895 default: 00896 return AR_STAT_INVALID_TYPE; 00897 } 00898 else if(AR_FLOAT_SIZE(ptype) == AR_FLOAT_32) 00899 return ar_native2(ARPOWXDR, 00900 result, resulttype, 00901 &tbase, &tpow); 00902 else 00903 return ar_native2(ARPOWXDXD, 00904 result, resulttype, 00905 &tbase, &tpow); 00906 00907 case IEEE_COMPLEX_32: 00908 if(AR_CLASS(ptype) == AR_CLASS_INT) 00909 switch (ptype) { 00910 case AR_Int_8_S: 00911 case AR_Int_16_S: 00912 ZERO_INT32_ALL(&tpow_32); 00913 (void) AR_convert((AR_DATA*) &tpow_32, 00914 &int32type, 00915 (AR_DATA*) &tpow, 00916 &ptype); 00917 return ar_native2(ARPOWCI, 00918 result, resulttype, 00919 &tbase, &tpow_32); 00920 case AR_Int_32_S: 00921 return ar_native2(ARPOWCI, 00922 result, resulttype, 00923 &tbase, &tpow); 00924 case AR_Int_64_S: 00925 return ar_native2(ARPOWCJ, 00926 result, resulttype, 00927 &tbase, &tpow); 00928 default: 00929 return AR_STAT_INVALID_TYPE; 00930 } 00931 else if(AR_FLOAT_IS_COMPLEX(ptype) != AR_FLOAT_COMPLEX) 00932 return ar_native2(ARPOWCR, 00933 result, resulttype, 00934 &tbase, &tpow); 00935 else 00936 return ar_native2(ARPOWCC, 00937 result, resulttype, 00938 &tbase, &tpow); 00939 00940 case IEEE_COMPLEX_64: 00941 if(AR_CLASS(ptype) == AR_CLASS_INT) 00942 switch (ptype) { 00943 case AR_Int_8_S: 00944 case AR_Int_16_S: 00945 ZERO_INT32_ALL(&tpow_32); 00946 (void) AR_convert((AR_DATA*) &tpow_32, 00947 &int32type, 00948 (AR_DATA*) &tpow, 00949 &ptype); 00950 return ar_native2(ARPOWCDI, 00951 result, resulttype, 00952 &tbase, &tpow_32); 00953 case AR_Int_32_S: 00954 return ar_native2(ARPOWCDI, 00955 result, resulttype, 00956 &tbase, &tpow); 00957 case AR_Int_64_S: 00958 return ar_native2(ARPOWCDJ, 00959 result, resulttype, 00960 &tbase, &tpow); 00961 default: 00962 return AR_STAT_INVALID_TYPE; 00963 } 00964 else 00965 return ar_native2(ARPOWCDCD, 00966 result, resulttype, 00967 &tbase, &tpow); 00968 00969 case IEEE_COMPLEX_128: 00970 if(AR_CLASS(ptype) == AR_CLASS_INT) 00971 switch (ptype) { 00972 case AR_Int_8_S: 00973 case AR_Int_16_S: 00974 ZERO_INT32_ALL(&tpow_32); 00975 (void) AR_convert((AR_DATA*) &tpow_32, 00976 &int32type, 00977 (AR_DATA*) &tpow, 00978 &ptype); 00979 return ar_native2(ARPOWCXDI, 00980 result, resulttype, 00981 &tbase, &tpow_32); 00982 case AR_Int_32_S: 00983 return ar_native2(ARPOWCXDI, 00984 result, resulttype, 00985 &tbase, &tpow); 00986 case AR_Int_64_S: 00987 return ar_native2(ARPOWCXDJ, 00988 result, resulttype, 00989 &tbase, &tpow); 00990 default: 00991 return AR_STAT_INVALID_TYPE; 00992 } 00993 else 00994 return ar_native2(ARPOWCXDCXD, 00995 result, resulttype, 00996 &tbase, &tpow); 00997 00998 default: 00999 switch (btype) { 01000 case AR_Int_8_S: 01001 case AR_Int_16_S: 01002 ZERO_INT32_ALL(&tbase_32); 01003 ZERO_INT32_ALL(&tpow_32); 01004 ZERO_INT32_ALL(&result_32); 01005 (void) AR_convert((AR_DATA*) &tbase_32, 01006 &int32type, 01007 (AR_DATA*) &tbase, 01008 &btype); 01009 (void) AR_convert((AR_DATA*) &tpow_32, 01010 &int32type, 01011 (AR_DATA*) &tpow, 01012 &ptype); 01013 status = ar_native2(ARPOWII, 01014 &result_32, &int32type, 01015 &tbase_32, &tpow_32); 01016 if (IS_ERROR_STATUS(status)) 01017 return status; 01018 return AR_convert((AR_DATA*) result, resulttype, 01019 (AR_DATA*) &result_32, &int32type); 01020 case AR_Int_32_S: 01021 ZERO_INT32_UPPER(result); 01022 return ar_native2(ARPOWII, 01023 result, resulttype, 01024 &tbase, &tpow); 01025 case AR_Int_64_S: 01026 return ar_native2(ARPOWJJ, 01027 result, resulttype, 01028 &tbase, &tpow); 01029 default: 01030 return AR_STAT_INVALID_TYPE; 01031 } 01032 } 01033 } 01034 01035 01036 /* Native complex division */ 01037 int 01038 ar_divide_complex (ar_data *result, const AR_TYPE *resulttype, 01039 const ar_data *opnd1, const AR_TYPE *opnd1type, 01040 const ar_data *opnd2, const AR_TYPE *opnd2type) 01041 { 01042 #if !defined _CRAYMPP 01043 #define ARCDIV arcdiv_ 01044 #define ARCDDIV arcddiv_ 01045 #define ARCXDDIV arcxddiv_ 01046 #endif 01047 extern void ARCDIV (ar_data *res, const ar_data *num, 01048 const ar_data *den); 01049 extern void ARCDDIV(ar_data *res, const ar_data *num, 01050 const ar_data *den); 01051 extern void ARCXDDIV(ar_data *res, const ar_data *num, 01052 const ar_data *den); 01053 01054 switch(UNROUNDED_TYPE(*resulttype)) { 01055 01056 case IEEE_COMPLEX_32: 01057 return ar_native2(ARCDIV, result, resulttype, opnd1, opnd2); 01058 01059 case IEEE_COMPLEX_64: 01060 return ar_native2(ARCDDIV, result, resulttype, opnd1, opnd2); 01061 01062 case IEEE_COMPLEX_128: 01063 return ar_native2(ARCXDDIV, result, resulttype, opnd1, opnd2); 01064 } 01065 01066 return AR_STAT_INVALID_TYPE; 01067 } 01068 01069 static int calling_math_lib = 0; 01070 static jmp_buf math_lib_jmpbuf; 01071 01072 /* Catch floating-point exceptions */ 01073 static volatile int fp_error = 0; 01074 static void fptrap (int sig) { 01075 fp_error = 1; 01076 signal (SIGFPE, fptrap); 01077 } 01078 01079 /* Catch math library detected errors */ 01080 int _lerror() { 01081 errno = ERANGE; 01082 return 0; 01083 } 01084 01085 /* These macros are used to protect ourselves when calling math libraries. */ 01086 01087 #define ARMOR_ON \ 01088 oldfpe = signal (SIGFPE, fptrap); \ 01089 fp_error = 0; \ 01090 errno = 0; \ 01091 calling_math_lib = 1; \ 01092 if (setjmp (math_lib_jmpbuf)) \ 01093 errno = 1; \ 01094 else 01095 01096 #define ARMOR_OFF(status, result, resulttype) \ 01097 calling_math_lib = 0; \ 01098 signal (SIGFPE, oldfpe); \ 01099 if (fp_error | errno) { \ 01100 if (fp_error > 0) \ 01101 status = AR_STAT_OVERFLOW; \ 01102 else if (fp_error == 0) \ 01103 status = AR_STAT_UNDEFINED; \ 01104 } else \ 01105 status |= AR_status((AR_DATA*)result, resulttype); 01106 01107 static int 01108 ar_native1(void (function)(), ar_data *result, const AR_TYPE *resulttype, const ar_data *opnd) 01109 { 01110 int status = AR_STAT_OK; 01111 void (*oldfpe)(); 01112 01113 ARMOR_ON 01114 (function) (result, opnd); 01115 ARMOR_OFF (status, result, resulttype); 01116 01117 return status; 01118 01119 } 01120 01121 static int 01122 ar_native2(void (function)(), ar_data *result, const AR_TYPE *resulttype, const ar_data *opnd1, const ar_data *opnd2) 01123 { 01124 int status = AR_STAT_OK; 01125 void (*oldfpe)(); 01126 01127 ARMOR_ON 01128 (function) (result, opnd1, opnd2); 01129 ARMOR_OFF (status, result, resulttype); 01130 01131 return status; 01132 01133 } 01134 01135 #else /* !(defined(_CRAYMPP) || defined(__mips)) */ 01136 01137 /* Use native math library for intrinsic function evaluation */ 01138 01139 #include <math.h> 01140 #include <signal.h> 01141 #include <errno.h> 01142 #include <setjmp.h> 01143 #include <stdio.h> 01144 01145 #undef double_t 01146 #undef complex_t 01147 #undef NAT_ROUND 01148 #undef NAT_COMPLEX 01149 01150 #define double_t double 01151 01152 #if defined(__sparc__) || defined(__mips) 01153 01154 #define NAT_ROUND AR_Float_IEEE_NR_32 01155 #define NAT_COMPLEX AR_Complex_IEEE_NR_32 01156 01157 #else /* IEEE 754 native arithmetic assumed */ 01158 01159 #define NAT_ROUND AR_Float_IEEE_NR_64 01160 #define NAT_COMPLEX AR_Complex_IEEE_NR_64 01161 01162 #if __svr4__ 01163 #include <complex.h> 01164 #define complex_t double complex 01165 #endif 01166 01167 #endif 01168 01169 static AR_TYPE native_round_single = NAT_ROUND; 01170 static AR_TYPE native_complex = NAT_COMPLEX; 01171 static AR_TYPE native_double = AR_Float_IEEE_NR_64; 01172 static AR_TYPE native_long_double = AR_Float_IEEE_NR_64; 01173 static AR_TYPE native_dbl_complex = AR_Complex_IEEE_NR_64; 01174 01175 01176 /* Fortran character index */ 01177 int 01178 ar_index (ar_data *result, const AR_TYPE *resulttype, 01179 const char *str1, const char *str2, const ar_data *backward) 01180 { 01181 return AR_STAT_INVALID_TYPE; /* Not available */ 01182 } 01183 01184 01185 /* Fortran character scan */ 01186 int 01187 ar_scan (ar_data *result, const AR_TYPE *resulttype, 01188 const char *str1, const char *str2, const ar_data *backward) 01189 { 01190 return AR_STAT_INVALID_TYPE; /* Not available */ 01191 } 01192 01193 01194 /* Fortran character verify */ 01195 int 01196 ar_verify (ar_data *result, const AR_TYPE *resulttype, 01197 const char *str1, const char *str2, const ar_data *backward) 01198 { 01199 return AR_STAT_INVALID_TYPE; /* Not available */ 01200 } 01201 01202 01203 /* Fortran-90 reshape */ 01204 int 01205 ar_reshape (void *result, const void *source, const void *shape, 01206 const void *pad, const void *order) 01207 { 01208 return AR_STAT_INVALID_TYPE; /* Not available */ 01209 } 01210 01211 01212 /* Fortran-90 transfer */ 01213 int 01214 ar_transfer (void *result, const void *source, const void *mold, long *length) 01215 { 01216 return AR_STAT_INVALID_TYPE; /* Not available */ 01217 } 01218 01219 01220 /* Fortran-90 modulo */ 01221 int 01222 ar_modulo (ar_data *result, const AR_TYPE *resulttype, 01223 const ar_data *opnd1, const AR_TYPE *opnd1type, 01224 const ar_data *opnd2, const AR_TYPE *opnd2type) 01225 { 01226 return AR_STAT_INVALID_TYPE; /* Not available */ 01227 } 01228 01229 /* Fortran-90 selected_real_kind */ 01230 int 01231 ar_selected_real_kind (ar_data *result, const AR_TYPE *resulttype, 01232 const ar_data *opnd1, const AR_TYPE *opnd1type, 01233 const ar_data *opnd2, const AR_TYPE *opnd2type) 01234 { 01235 return AR_STAT_INVALID_TYPE; /* Not available */ 01236 } 01237 01238 /* Square root */ 01239 int 01240 ar_sqrt (ar_data *result, const AR_TYPE *resulttype, 01241 const ar_data *opnd, const AR_TYPE *opndtype) 01242 { 01243 int status; 01244 ar_data nat_opnd, nat_result; 01245 01246 /* Check whether operation exists in native library */ 01247 01248 if (*resulttype == native_double) { 01249 *(double_t *) result = (sqrt) (*(double_t *) opnd); 01250 return AR_status((AR_DATA*)result, resulttype); 01251 } 01252 01253 #ifdef complex_t 01254 if (*resulttype == native_complex) { 01255 *(complex_t *) result = (csqrt) (*(complex_t *) opnd); 01256 return AR_status((AR_DATA*)result, resulttype); 01257 } 01258 #endif 01259 01260 /* The operation is not a native library operation. 01261 * Convert, approximate with native arithmetic, and convert back. 01262 */ 01263 01264 if (AR_FLOAT_IS_COMPLEX (*resulttype) == AR_FLOAT_COMPLEX) { 01265 #ifdef complex_t 01266 status = AR_convert ((AR_DATA*)&nat_opnd, &native_complex, 01267 (AR_DATA*)opnd, opndtype); 01268 status |= ar_sqrt (&nat_result, &native_complex, 01269 &nat_opnd, &native_complex); 01270 status &= ~(AR_NEGATIVE | AR_ZERO); 01271 status |= AR_convert ((AR_DATA*)result, resulttype, 01272 (AR_DATA*)&nat_result, &native_complex); 01273 #else 01274 status = AR_STAT_INVALID_TYPE; 01275 #endif 01276 } else { 01277 status = AR_convert ((AR_DATA*)&nat_opnd, &native_long_double, 01278 (AR_DATA*)opnd, opndtype); 01279 status |= ar_sqrt (&nat_result, &native_long_double, 01280 &nat_opnd, &native_long_double); 01281 status &= ~(AR_NEGATIVE | AR_ZERO); 01282 status |= AR_convert ((AR_DATA*)result, resulttype, 01283 (AR_DATA*)&nat_result, &native_long_double); 01284 } 01285 01286 return status; 01287 } 01288 01289 01290 /* Natural (base "e") logarithm */ 01291 int 01292 ar_log (ar_data *result, const AR_TYPE *resulttype, 01293 const ar_data *opnd, const AR_TYPE *opndtype) 01294 { 01295 int status; 01296 ar_data nat_opnd, nat_result; 01297 01298 /* Check whether operation exists in native library */ 01299 01300 if (*resulttype == native_double) { 01301 *(double_t *) result = (log) (*(double_t *) opnd); 01302 return AR_status((AR_DATA*)result, resulttype); 01303 } 01304 #ifdef complex_t 01305 if (*resulttype == native_complex) { 01306 *(complex_t *) result = (clog) (*(complex_t *) opnd); 01307 return AR_status((AR_DATA*)result, resulttype); 01308 } 01309 #endif 01310 01311 /* The operation is not a native library operation. 01312 * Convert, approximate with native arithmetic, and convert back. 01313 */ 01314 01315 if (AR_FLOAT_IS_COMPLEX (*resulttype) == AR_FLOAT_COMPLEX) { 01316 #ifdef complex_t 01317 status = AR_convert ((AR_DATA*)&nat_opnd, &native_complex, 01318 (AR_DATA*)opnd, opndtype); 01319 status |= ar_log (&nat_result, &native_complex, 01320 &nat_opnd, &native_complex); 01321 status &= ~(AR_NEGATIVE | AR_ZERO); 01322 status |= AR_convert ((AR_DATA*)result, resulttype, 01323 (AR_DATA*)&nat_result, &native_complex); 01324 #else 01325 status = AR_STAT_INVALID_TYPE; 01326 #endif 01327 } else { 01328 status = AR_convert ((AR_DATA*)&nat_opnd, &native_long_double, 01329 (AR_DATA*)opnd, opndtype); 01330 status |= ar_log (&nat_result, &native_long_double, 01331 &nat_opnd, &native_long_double); 01332 status &= ~(AR_NEGATIVE | AR_ZERO); 01333 status |= AR_convert ((AR_DATA*)result, resulttype, 01334 (AR_DATA*)&nat_result, &native_long_double); 01335 } 01336 01337 return status; 01338 } 01339 01340 01341 /* Exponential ("e" ** x) function */ 01342 int 01343 ar_exp (ar_data *result, const AR_TYPE *resulttype, 01344 const ar_data *opnd, const AR_TYPE *opndtype) 01345 { 01346 int status; 01347 ar_data nat_opnd, nat_result; 01348 01349 if (*resulttype == native_double) { 01350 *(double_t *) result = (exp) (*(double_t *) opnd); 01351 return AR_status((AR_DATA*)result, resulttype); 01352 } 01353 #ifdef complex_t 01354 if (*resulttype == native_complex) { 01355 *(complex_t *) result = (cexp) (*(complex_t *) opnd); 01356 return AR_status((AR_DATA*)result, resulttype); 01357 } 01358 #endif 01359 01360 /* The operation is not a native library operation. 01361 * Convert, approximate with native arithmetic, and convert back. 01362 */ 01363 01364 if (AR_FLOAT_IS_COMPLEX (*resulttype) == AR_FLOAT_COMPLEX) { 01365 #ifdef complex_t 01366 status = AR_convert ((AR_DATA*)&nat_opnd, &native_complex, 01367 (AR_DATA*)opnd, opndtype); 01368 status |= ar_exp (&nat_result, &native_complex, 01369 &nat_opnd, &native_complex); 01370 status &= ~(AR_NEGATIVE | AR_ZERO); 01371 status |= AR_convert ((AR_DATA*)result, resulttype, 01372 (AR_DATA*)&nat_result, &native_complex); 01373 #else 01374 status = AR_STAT_INVALID_TYPE; 01375 #endif 01376 } else { 01377 status = AR_convert ((AR_DATA*)&nat_opnd, &native_long_double, 01378 (AR_DATA*)opnd, opndtype); 01379 status |= ar_exp (&nat_result, &native_long_double, 01380 &nat_opnd, &native_long_double); 01381 status &= ~(AR_NEGATIVE | AR_ZERO); 01382 status |= AR_convert ((AR_DATA*)result, resulttype, 01383 (AR_DATA*)&nat_result, &native_long_double); 01384 } 01385 01386 return status; 01387 } 01388 01389 01390 /* Complex absolute value */ 01391 int 01392 ar_cabs (ar_data *result, const AR_TYPE *resulttype, 01393 const ar_data *opnd, const AR_TYPE *opndtype) 01394 { 01395 int status; 01396 ar_data re, im, resq, imsq, sum, nat_opnd, nat_result; 01397 AR_TYPE reimtype; 01398 01399 #if defined(complex_t) 01400 if (*resulttype == native_double && 01401 *opndtype == native_complex) { 01402 *(double_t *) result = (cabs) (*(complex_t *) opnd); 01403 return AR_status((AR_DATA*)result, resulttype); 01404 } 01405 01406 /* The operation is not a native library operation. 01407 * Convert, execute, and convert back. 01408 */ 01409 01410 status = AR_convert ((AR_DATA*)&nat_opnd, &native_complex, 01411 (AR_DATA*)opnd, opndtype); 01412 status |= ar_cabs (&nat_result, &native_round_single, 01413 &nat_opnd, &native_complex); 01414 status &= ~(AR_NEGATIVE | AR_ZERO); 01415 status |= AR_convert ((AR_DATA*)result, resulttype, 01416 (AR_DATA*)&nat_result, &native_round_single); 01417 #else 01418 01419 /* If no native complex operations, compute sqrt (real**2 * imag**2). */ 01420 01421 status = ar_decompose_complex (&re, &im, &reimtype, 01422 opnd, opndtype); 01423 status |= AR_multiply ((AR_DATA*)&resq, &reimtype, 01424 (AR_DATA*)&re, &reimtype, 01425 (AR_DATA*)&re, &reimtype); 01426 status |= AR_multiply ((AR_DATA*)&imsq, &reimtype, 01427 (AR_DATA*)&im, &reimtype, 01428 (AR_DATA*)&im, &reimtype); 01429 status |= AR_add ((AR_DATA*)&sum, &reimtype, 01430 (AR_DATA*)&resq, &reimtype, 01431 (AR_DATA*)&imsq, &reimtype); 01432 status &= ~(AR_STAT_ZERO | AR_STAT_NEGATIVE); 01433 status |= AR_sqrt ((AR_DATA*)result, resulttype, (AR_DATA*)&sum, &reimtype); 01434 01435 #endif 01436 01437 return status; 01438 } 01439 01440 01441 /* Exponentiation */ 01442 int 01443 ar_power(ar_data *result, const AR_TYPE *resulttype, 01444 const ar_data *base, const AR_TYPE *basetype, 01445 const ar_data *power, const AR_TYPE *powertype) 01446 { 01447 int i; 01448 int status; 01449 int status2; 01450 ar_data tbase, tpow, tpow2, temp, one; 01451 AR_TYPE inttype = AR_Int_64_S; 01452 static int loopchk = 0; 01453 01454 #ifdef complex_t 01455 if (*basetype == native_complex || *powertype == native_complex) { 01456 if (*resulttype != native_complex) 01457 return AR_STAT_INVALID_TYPE; 01458 if (*basetype != native_complex) { 01459 status = AR_convert ((AR_DATA*)&tbase, &native_complex, 01460 (AR_DATA*)base, basetype); 01461 if (IS_ERROR_STATUS(status)) 01462 return AR_STAT_INVALID_TYPE; 01463 } else 01464 tbase = *base; 01465 if (*powertype != native_complex) { 01466 status = AR_convert ((AR_DATA*)&tpow, &native_complex, 01467 (AR_DATA*)power, powertype); 01468 if (IS_ERROR_STATUS(status)) 01469 return AR_STAT_INVALID_TYPE; 01470 } else 01471 tpow = *power; 01472 *(complex_t *) result = (cpow) (*(complex_t *) &tbase, 01473 *(complex_t *) &tpow); 01474 return AR_status((AR_DATA*)result, resulttype); 01475 } 01476 #endif 01477 01478 if (*basetype == native_double || 01479 *powertype == native_double) { 01480 if (*resulttype != native_double) 01481 return AR_STAT_INVALID_TYPE; 01482 if (*basetype != native_double) { 01483 status = AR_convert ((AR_DATA*)&tbase, &native_double, 01484 (AR_DATA*)base, basetype); 01485 if (IS_ERROR_STATUS(status)) 01486 return AR_STAT_INVALID_TYPE; 01487 } else 01488 tbase = *base; 01489 if (*powertype != native_double) { 01490 status = AR_convert ((AR_DATA*)&tpow, &native_double, 01491 (AR_DATA*)power, powertype); 01492 if (IS_ERROR_STATUS(status)) 01493 return AR_STAT_INVALID_TYPE; 01494 } else 01495 tpow = *power; 01496 *(double_t *) result = (pow) (*(double_t *) &tbase, 01497 *(double_t *) &tpow); 01498 return AR_status((AR_DATA*)result, resulttype); 01499 } 01500 01501 /* The operation is not a native library operation. 01502 * Convert, approximate with native arithmetic, and convert back. 01503 */ 01504 01505 /* Use native exponentiation routines if appropriate. */ 01506 01507 /* 01508 * Native routines will not be used directly. 01509 */ 01510 01511 /* Get "1" of the proper type */ 01512 if (AR_one ((AR_DATA*)&one, basetype) & AR_STAT_INVALID_TYPE) 01513 return AR_STAT_INVALID_TYPE; 01514 01515 status = AR_status ((AR_DATA*)base, basetype); 01516 01517 /* 0 ** 0 and 0 ** (neg) are invalid; 0 ** (pos) == 0. */ 01518 if (status & AR_STAT_ZERO) { 01519 status = AR_status ((AR_DATA*)power, powertype); 01520 if (status & (AR_STAT_ZERO | AR_STAT_NEGATIVE)) 01521 return AR_STAT_UNDEFINED; 01522 return AR_convert ((AR_DATA*)result, resulttype, (AR_DATA*)base, basetype); 01523 } 01524 01525 /* 1 ** (anything) == 1 */ 01526 if (AR_compare ((AR_DATA*)base, basetype, (AR_DATA*)&one, basetype) == AR_Compare_EQ) 01527 return AR_convert ((AR_DATA*)result, resulttype, (AR_DATA*)&one, basetype); 01528 01529 /* Exponentiation by integer powers */ 01530 if (AR_CLASS (*powertype) == AR_CLASS_INT) { 01531 01532 if (*resulttype != *basetype) 01533 return AR_STAT_INVALID_TYPE; 01534 01535 /* (integer) ** 1 == (integer) */ 01536 if (AR_CLASS (*basetype) == AR_CLASS_INT && 01537 AR_compare ((AR_DATA*)power, powertype, (AR_DATA*)&one, powertype) == 01538 AR_Compare_EQ) 01539 return AR_convert ((AR_DATA*)result, resulttype, (AR_DATA*)base, basetype); 01540 01541 /* (-1) ** (even) == 1, (-1) ** (odd) == -1 */ 01542 AR_negate ((AR_DATA*)&tbase, basetype, (AR_DATA*)&one, basetype); 01543 if (AR_compare ((AR_DATA*)base, basetype, (AR_DATA*)&tbase, basetype) == 01544 AR_Compare_EQ) { 01545 if (power->ar_i64.part4 & 1) 01546 return AR_convert ((AR_DATA*)result, resulttype, 01547 (AR_DATA*)base, basetype); 01548 return AR_convert ((AR_DATA*)result, resulttype, 01549 (AR_DATA*)&one, basetype); 01550 } 01551 01552 /* Compute negative integer powers by computing a power 01553 * of a reciprocal. 01554 */ 01555 if (INT_SIGN(*powertype, power)) { 01556 01557 /* (int other than 0, 1, -1) ** (negative) == 0 */ 01558 if (AR_CLASS (*basetype) == AR_CLASS_INT) 01559 return AR_convert ((AR_DATA*)result, resulttype, 01560 (AR_DATA*)&AR_const_zero, &inttype); 01561 01562 status = AR_divide ((AR_DATA*)&tbase, basetype, 01563 (AR_DATA*)&one, basetype, 01564 (AR_DATA*)base, basetype); 01565 status &= AR_ERROR_STATUS; 01566 if (status) { 01567 AR_convert ((AR_DATA*)result, resulttype, 01568 &AR_const_zero, &inttype); 01569 return status; 01570 } 01571 AR_negate ((AR_DATA*)&tpow, powertype, (AR_DATA*)power, powertype); 01572 01573 } else { 01574 tbase = *base; 01575 tpow = *power; 01576 } 01577 01578 /* Perform exponentiation by repeated multiplication */ 01579 status = AR_convert ((AR_DATA*)result, resulttype, (AR_DATA*)&one, basetype); 01580 status2 = AR_STAT_OK; 01581 01582 /* 01583 * Solaris porting 01584 * The following block contains a "undefined symbol" - opnd1type 01585 * it seems opnd1type should be basetype, since in this function (power) 01586 * there are only base/power these two operands. 01587 */ 01588 switch (AR_INT_SIZE (*basetype)) { 01589 case AR_INT_SIZE_8: 01590 AR_convert ((AR_DATA*) &tpow2, &inttype, 01591 (AR_DATA*) &tpow, basetype); 01592 break; 01593 01594 case AR_INT_SIZE_16: 01595 AR_convert ((AR_DATA*) &tpow2, &inttype, 01596 (AR_DATA*) &tpow, basetype); 01597 break; 01598 01599 case AR_INT_SIZE_32: 01600 AR_convert ((AR_DATA*) &tpow2, &inttype, 01601 (AR_DATA*) &tpow, basetype); 01602 break; 01603 01604 case AR_INT_SIZE_46: 01605 case AR_INT_SIZE_64: 01606 tpow2 = tpow; 01607 break; 01608 01609 default: 01610 return (AR_STAT_INVALID_TYPE); 01611 } 01612 01613 while (!IS_INT64_ZERO(&tpow2)) { 01614 if (IS_ERROR_STATUS(status2)) { 01615 AR_convert ((AR_DATA*)result, resulttype, 01616 &AR_const_zero, &inttype); 01617 return status2 & AR_ERROR_STATUS; 01618 } 01619 if (tpow2.ar_i64.part4 & 1) { 01620 status = AR_multiply ((AR_DATA*)result, 01621 resulttype, 01622 (AR_DATA*)result, 01623 resulttype, 01624 (AR_DATA*)&tbase, 01625 basetype); 01626 if (IS_ERROR_STATUS(status)) { 01627 AR_convert ((AR_DATA*)result, 01628 resulttype, 01629 &AR_const_zero, 01630 &inttype); 01631 return status & AR_ERROR_STATUS; 01632 } 01633 } 01634 status2 = AR_multiply ((AR_DATA*)&tbase, basetype, 01635 (AR_DATA*)&tbase, basetype, 01636 (AR_DATA*)&tbase, basetype); 01637 SHRIGHT64 (tpow2.ar_i64); 01638 } 01639 01640 return status; 01641 } 01642 01643 /* 01644 * Exponentiation by non-native floating powers 01645 */ 01646 01647 if(loopchk > 0) 01648 return AR_STAT_INVALID_TYPE; 01649 01650 if (AR_CLASS (*basetype) == AR_CLASS_FLOAT && 01651 AR_FLOAT_IS_COMPLEX (*basetype) == AR_FLOAT_COMPLEX || 01652 AR_CLASS (*powertype) == AR_CLASS_FLOAT && 01653 AR_FLOAT_IS_COMPLEX (*powertype) == AR_FLOAT_COMPLEX) { 01654 01655 #ifndef complex_t 01656 return AR_STAT_INVALID_TYPE; 01657 #else 01658 /* Compute (native complex) ** (native complex) */ 01659 status = AR_convert ((AR_DATA*)&tbase, &native_complex, (AR_DATA*)base, basetype); 01660 status |= AR_convert ((AR_DATA*)&tpow, &native_complex, (AR_DATA*)power, powertype); 01661 if (!IS_ERROR_STATUS(status)) { 01662 loopchk++; 01663 status = ar_power (&temp, &native_complex, 01664 &tbase, &native_complex, 01665 &tpow, &native_complex); 01666 loopchk--; 01667 if (!IS_ERROR_STATUS(status)) 01668 status = AR_convert ((AR_DATA*)result, resulttype, 01669 (AR_DATA*)&temp, &native_complex); 01670 } 01671 return status; 01672 #endif 01673 } 01674 01675 /* Compute as (native double prec) ** (native double prec) */ 01676 status = AR_convert ((AR_DATA*)&tbase, &native_long_double, (AR_DATA*)base, basetype); 01677 status |= AR_convert ((AR_DATA*)&tpow, &native_long_double, (AR_DATA*)power, powertype); 01678 if (!IS_ERROR_STATUS(status)) { 01679 loopchk++; 01680 status = ar_power (&temp, &native_long_double, 01681 &tbase, &native_long_double, 01682 &tpow, &native_long_double); 01683 loopchk--; 01684 if (!IS_ERROR_STATUS(status)) 01685 status = AR_convert ((AR_DATA*)result, resulttype, 01686 (AR_DATA*)&temp, &native_long_double); 01687 } 01688 return status; 01689 } 01690 01691 01692 /* Native complex division */ 01693 int 01694 ar_divide_complex (ar_data *result, const AR_TYPE *resulttype, 01695 const ar_data *opnd1, const AR_TYPE *opnd1type, 01696 const ar_data *opnd2, const AR_TYPE *opnd2type) 01697 { 01698 01699 /* Assume types all match (see logic in AR_divide) */ 01700 01701 AR_DATA a, b, c, d, ac, bd, bc, ad, cc, dd, acbd, bcad, ccdd, re, im; 01702 AR_TYPE reimtype1, reimtype2, temptype; 01703 int status, restat, imstat; 01704 01705 status = ar_decompose_complex ((ar_data*)&a, (ar_data*)&b, &reimtype1, opnd1, opnd1type); 01706 status |= ar_decompose_complex ((ar_data*)&c, (ar_data*)&d, &reimtype2, opnd2, opnd2type); 01707 01708 /* PDGCS requests that a different sequence be used when the 01709 * imaginary part of the denominator is zero. A meeting of 01710 * managers on 11/30/93 decided in favor of this expediency. 01711 * Note that we do NOT apply special-case processing when 01712 * the real part of the denominator is zero. 01713 */ 01714 01715 imstat = AR_status (&d, &reimtype2); 01716 if (imstat & AR_STAT_ZERO) { 01717 01718 /* zero imaginary part, use short sequence */ 01719 restat = AR_divide (&re, &reimtype1, 01720 &a, &reimtype1, &c, &reimtype2); 01721 imstat = AR_divide (&im, &reimtype1, 01722 &b, &reimtype1, &c, &reimtype2); 01723 01724 } else { 01725 01726 /* 01727 * general sequence: 01728 * 01729 * a + bi (a + bi)(c - di) (ac + bd) (bc - ad)i 01730 * ------ = ---------------- = --------- + ---------- 01731 * c + di (c + di)(c - di) c*c + d*d c*c + d*d 01732 */ 01733 01734 status |= AR_multiply (&ac, &reimtype1, 01735 &a, &reimtype1, &c, &reimtype2); 01736 status |= AR_multiply (&bd, &reimtype1, &b, 01737 &reimtype1, &d, &reimtype2); 01738 status |= AR_multiply (&bc, &reimtype1, 01739 &b, &reimtype1, &c, &reimtype2); 01740 status |= AR_multiply (&ad, &reimtype1, 01741 &a, &reimtype1, &d, &reimtype2); 01742 status |= AR_multiply (&cc, &reimtype2, 01743 &c, &reimtype2, &c, &reimtype2); 01744 status |= AR_multiply (&dd, &reimtype2, 01745 &d, &reimtype2, &d, &reimtype2); 01746 status |= AR_add (&acbd, &reimtype1, 01747 &ac, &reimtype1, &bd, &reimtype1); 01748 status |= AR_subtract (&bcad, &reimtype1, 01749 &bc, &reimtype1, &ad, &reimtype1); 01750 status |= AR_add (&ccdd, &reimtype1, 01751 &cc, &reimtype1, &dd, &reimtype1); 01752 01753 restat = AR_divide (&re, &reimtype1, 01754 &acbd, &reimtype1, &ccdd, &reimtype1); 01755 imstat = AR_divide (&im, &reimtype1, 01756 &bcad, &reimtype1, &ccdd, &reimtype1); 01757 } 01758 01759 status |= ar_compose_complex (result, &temptype, (ar_data*)&re, (ar_data*)&im, &reimtype1); 01760 status |= restat | imstat; 01761 status &= ~(AR_STAT_ZERO | AR_STAT_NEGATIVE); 01762 status |= restat & imstat & AR_STAT_ZERO; 01763 return status; 01764 } 01765 01766 #endif /* _Solaris || defined(_CRAYMPP) || defined(__mips) */ 01767 01768 /* string -> floating point */ 01769 int 01770 ar_convert_str_to_float (ar_data *result, const AR_TYPE *resulttype, 01771 const char *str) 01772 { 01773 return ar_cvt_str_to_float (result, resulttype, str); 01774 } 01775 01776 #endif /* _CRAY */ 01777 01778 01779 static char USMID [] = "\n%Z%%M% %I% %G% %U%\n";
1.7.6.1 
