Open64: osprey/arith/native.c Source File

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