[BACK]Return to math.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / builtin

Diff for /OpenXM_contrib2/asir2000/builtin/math.c between version 1.12 and 1.13

version 1.12, 2015/08/14 13:51:54 version 1.13, 2018/03/29 01:32:50
Line 45 
Line 45 
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,   * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.   * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *   *
  * $OpenXM: OpenXM_contrib2/asir2000/builtin/math.c,v 1.11 2015/08/06 10:01:52 fujimoto Exp $   * $OpenXM: OpenXM_contrib2/asir2000/builtin/math.c,v 1.12 2015/08/14 13:51:54 fujimoto Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include <math.h>  #include <math.h>
Line 58  void Pdsqrt(),Pdsin(),Pdcos(),Pdtan(),Pdasin(),Pdacos(
Line 58  void Pdsqrt(),Pdsin(),Pdcos(),Pdtan(),Pdasin(),Pdacos(
 void Pabs(),Pdfloor(),Pdceil(),Pdrint(),Pdisnan();  void Pabs(),Pdfloor(),Pdceil(),Pdrint(),Pdisnan();
   
 struct ftab math_tab[] = {  struct ftab math_tab[] = {
         {"dsqrt",Pdsqrt,1},    {"dsqrt",Pdsqrt,1},
         {"dabs",Pabs,1},    {"dabs",Pabs,1},
         {"dsin",Pdsin,1},    {"dsin",Pdsin,1},
         {"dcos",Pdcos,1},    {"dcos",Pdcos,1},
         {"dtan",Pdtan,1},    {"dtan",Pdtan,1},
         {"dlog",Pdlog,1},    {"dlog",Pdlog,1},
         {"dexp",Pdexp,1},    {"dexp",Pdexp,1},
         {"dasin",Pdasin,1},    {"dasin",Pdasin,1},
         {"dacos",Pdacos,1},    {"dacos",Pdacos,1},
         {"datan",Pdatan,1},    {"datan",Pdatan,1},
         {"floor",Pdfloor,1},    {"floor",Pdfloor,1},
         {"dfloor",Pdfloor,1},    {"dfloor",Pdfloor,1},
         {"ceil",Pdceil,1},    {"ceil",Pdceil,1},
         {"dceil",Pdceil,1},    {"dceil",Pdceil,1},
         {"rint",Pdrint,1},    {"rint",Pdrint,1},
         {"drint",Pdrint,1},    {"drint",Pdrint,1},
         {"disnan",Pdisnan,1},    {"disnan",Pdisnan,1},
         {0,0,0},    {0,0,0},
 };  };
   
 void get_ri(Num z,double *r,double *i)  void get_ri(Num z,double *r,double *i)
 {  {
         if ( !z ) {    if ( !z ) {
                 *r = 0; *i = 0; return;      *r = 0; *i = 0; return;
         }    }
         if ( OID(z) != O_N )    if ( OID(z) != O_N )
                 error("get_ri : invalid argument");      error("get_ri : invalid argument");
         switch ( NID(z) ) {    switch ( NID(z) ) {
                 case N_Q: case N_R: case N_B:      case N_Q: case N_R: case N_B:
                         *r = ToReal(z); *i = 0;        *r = ToReal(z); *i = 0;
                         break;        break;
                 case N_C:      case N_C:
                         *r = ToReal(((C)z)->r);        *r = ToReal(((C)z)->r);
                         *i = ToReal(((C)z)->i);        *i = ToReal(((C)z)->i);
                         break;        break;
                 default:      default:
                         error("get_ri : invalid argument");        error("get_ri : invalid argument");
                         break;        break;
         }    }
 }  }
   
 void Pabs(arg,rp)  void Pabs(arg,rp)
 NODE arg;  NODE arg;
 Real *rp;  Real *rp;
 {  {
         double s,r,i;    double s,r,i;
   
         if ( !ARG0(arg) ) {    if ( !ARG0(arg) ) {
                 *rp = 0; return;      *rp = 0; return;
         }    }
         get_ri((Num)ARG0(arg),&r,&i);    get_ri((Num)ARG0(arg),&r,&i);
         if ( i == 0 )    if ( i == 0 )
                 s = fabs(r);      s = fabs(r);
         else if ( r == 0 )    else if ( r == 0 )
                 s = fabs(i);      s = fabs(i);
         else    else
                 s = sqrt(r*r+i*i);      s = sqrt(r*r+i*i);
         MKReal(s,*rp);    MKReal(s,*rp);
 }  }
   
 void Pdsqrt(arg,rp)  void Pdsqrt(arg,rp)
 NODE arg;  NODE arg;
 Num *rp;  Num *rp;
 {  {
         double s,r,i,a;    double s,r,i,a;
         C z;    C z;
         Real real;    Real real;
   
         if ( !ARG0(arg) ) {    if ( !ARG0(arg) ) {
                 *rp = 0; return;      *rp = 0; return;
         }    }
         get_ri((Num)ARG0(arg),&r,&i);    get_ri((Num)ARG0(arg),&r,&i);
         if ( i == 0 )    if ( i == 0 )
                 if ( r > 0 ) {      if ( r > 0 ) {
                         s = sqrt(r);        s = sqrt(r);
                         MKReal(s,real);        MKReal(s,real);
                         *rp = (Num)real;        *rp = (Num)real;
                 } else {      } else {
                         NEWC(z);        NEWC(z);
                         z->r = 0;        z->r = 0;
                         s = sqrt(-r); MKReal(s,real); z->i = (Num)real;        s = sqrt(-r); MKReal(s,real); z->i = (Num)real;
                         *rp = (Num)z;        *rp = (Num)z;
                 }      }
         else {    else {
                 a = sqrt(r*r+i*i);      a = sqrt(r*r+i*i);
                 NEWC(z);      NEWC(z);
                 s = sqrt((r+a)/2); MKReal(s,real); z->r = (Num)real;      s = sqrt((r+a)/2); MKReal(s,real); z->r = (Num)real;
                 s = i>0?sqrt((-r+a)/2):-sqrt((-r+a)/2);      s = i>0?sqrt((-r+a)/2):-sqrt((-r+a)/2);
                 MKReal(s,real); z->i = (Num)real;      MKReal(s,real); z->i = (Num)real;
                 *rp = (Num)z;      *rp = (Num)z;
         }    }
 }  }
   
 void Pdsin(arg,rp)  void Pdsin(arg,rp)
 NODE arg;  NODE arg;
 Real *rp;  Real *rp;
 {  {
         double s;    double s;
   
         s = sin(ToReal(ARG0(arg)));    s = sin(ToReal(ARG0(arg)));
         MKReal(s,*rp);    MKReal(s,*rp);
 }  }
   
 void Pdcos(arg,rp)  void Pdcos(arg,rp)
 NODE arg;  NODE arg;
 Real *rp;  Real *rp;
 {  {
         double s;    double s;
   
         s = cos(ToReal(ARG0(arg)));    s = cos(ToReal(ARG0(arg)));
         MKReal(s,*rp);    MKReal(s,*rp);
 }  }
   
 void Pdtan(arg,rp)  void Pdtan(arg,rp)
 NODE arg;  NODE arg;
 Real *rp;  Real *rp;
 {  {
         double s;    double s;
   
         s = tan(ToReal(ARG0(arg)));    s = tan(ToReal(ARG0(arg)));
         MKReal(s,*rp);    MKReal(s,*rp);
 }  }
   
 void Pdasin(arg,rp)  void Pdasin(arg,rp)
 NODE arg;  NODE arg;
 Real *rp;  Real *rp;
 {  {
         double s;    double s;
   
         s = asin(ToReal(ARG0(arg)));    s = asin(ToReal(ARG0(arg)));
         MKReal(s,*rp);    MKReal(s,*rp);
 }  }
   
 void Pdacos(arg,rp)  void Pdacos(arg,rp)
 NODE arg;  NODE arg;
 Real *rp;  Real *rp;
 {  {
         double s;    double s;
   
         s = acos(ToReal(ARG0(arg)));    s = acos(ToReal(ARG0(arg)));
         MKReal(s,*rp);    MKReal(s,*rp);
 }  }
   
 void Pdatan(arg,rp)  void Pdatan(arg,rp)
 NODE arg;  NODE arg;
 Real *rp;  Real *rp;
 {  {
         double s;    double s;
   
         s = atan(ToReal(ARG0(arg)));    s = atan(ToReal(ARG0(arg)));
         MKReal(s,*rp);    MKReal(s,*rp);
 }  }
   
 void Pdlog(arg,rp)  void Pdlog(arg,rp)
 NODE arg;  NODE arg;
 Real *rp;  Real *rp;
 {  {
         double s;    double s;
   
         s = log(ToReal(ARG0(arg)));    s = log(ToReal(ARG0(arg)));
         MKReal(s,*rp);    MKReal(s,*rp);
 }  }
   
 void Pdexp(arg,rp)  void Pdexp(arg,rp)
 NODE arg;  NODE arg;
 Real *rp;  Real *rp;
 {  {
         double s;    double s;
   
         s = exp(ToReal(ARG0(arg)));    s = exp(ToReal(ARG0(arg)));
         MKReal(s,*rp);    MKReal(s,*rp);
 }  }
   
 void Pdfloor(arg,rp)  void Pdfloor(arg,rp)
 NODE arg;  NODE arg;
 Q *rp;  Q *rp;
 {  {
         L a;    L a;
         unsigned int au,al;    unsigned int au,al;
         int sgn;    int sgn;
         Q q;    Q q;
         double d;    double d;
   
         if ( !ARG0(arg) ) {    if ( !ARG0(arg) ) {
                 *rp = 0;      *rp = 0;
                 return;      return;
         }    }
         d = floor(ToReal(ARG0(arg)));    d = floor(ToReal(ARG0(arg)));
         if ( d < -9.223372036854775808e18 || d >= 9.223372036854775808e18 )    if ( d < -9.223372036854775808e18 || d >= 9.223372036854775808e18 )
                 error("dfloor : OverFlow");      error("dfloor : OverFlow");
         if ( !d ) {    if ( !d ) {
                 *rp = 0;      *rp = 0;
                 return;      return;
         }    }
         a = (L)d;    a = (L)d;
         if ( a < 0 ) {    if ( a < 0 ) {
                 sgn = -1;      sgn = -1;
                 a = -a;      a = -a;
         } else    } else
                 sgn = 1;      sgn = 1;
 #if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__MINGW32__) || defined(__x86_64)  #if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__MINGW32__) || defined(__x86_64)
         au = ((unsigned int *)&a)[1];    au = ((unsigned int *)&a)[1];
         al = ((unsigned int *)&a)[0];    al = ((unsigned int *)&a)[0];
 #else  #else
         al = ((unsigned int *)&a)[1];    al = ((unsigned int *)&a)[1];
         au = ((unsigned int *)&a)[0];    au = ((unsigned int *)&a)[0];
 #endif  #endif
         if ( au ) {    if ( au ) {
                 NEWQ(q); SGN(q) = sgn; NM(q)=NALLOC(2); DN(q)=0;      NEWQ(q); SGN(q) = sgn; NM(q)=NALLOC(2); DN(q)=0;
                 PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;      PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;
         } else {    } else {
                 UTOQ(al,q); SGN(q) = sgn;      UTOQ(al,q); SGN(q) = sgn;
         }    }
         *rp = q;    *rp = q;
 }  }
   
 void Pdceil(arg,rp)  void Pdceil(arg,rp)
 NODE arg;  NODE arg;
 Q *rp;  Q *rp;
 {  {
         L a;    L a;
         unsigned int au,al;    unsigned int au,al;
         int sgn;    int sgn;
         Q q;    Q q;
         double d;    double d;
   
         if ( !ARG0(arg) ) {    if ( !ARG0(arg) ) {
                 *rp = 0;      *rp = 0;
                 return;      return;
         }    }
         d = ceil(ToReal(ARG0(arg)));    d = ceil(ToReal(ARG0(arg)));
         if ( d < -9.223372036854775808e18 || d >= 9.223372036854775808e18 )    if ( d < -9.223372036854775808e18 || d >= 9.223372036854775808e18 )
                 error("dceil : OverFlow");      error("dceil : OverFlow");
         if ( !d ) {    if ( !d ) {
                 *rp = 0;      *rp = 0;
                 return;      return;
         }    }
         a = (L)d;    a = (L)d;
         if ( a < 0 ) {    if ( a < 0 ) {
                 sgn = -1;      sgn = -1;
                 a = -a;      a = -a;
         } else    } else
                 sgn = 1;      sgn = 1;
 #if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__MINGW32__) || defined(__x86_64)  #if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__MINGW32__) || defined(__x86_64)
         au = ((unsigned int *)&a)[1];    au = ((unsigned int *)&a)[1];
         al = ((unsigned int *)&a)[0];    al = ((unsigned int *)&a)[0];
 #else  #else
         al = ((unsigned int *)&a)[1];    al = ((unsigned int *)&a)[1];
         au = ((unsigned int *)&a)[0];    au = ((unsigned int *)&a)[0];
 #endif  #endif
         if ( au ) {    if ( au ) {
                 NEWQ(q); SGN(q) = sgn; NM(q)=NALLOC(2); DN(q)=0;      NEWQ(q); SGN(q) = sgn; NM(q)=NALLOC(2); DN(q)=0;
                 PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;      PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;
         } else {    } else {
                 UTOQ(al,q); SGN(q) = sgn;      UTOQ(al,q); SGN(q) = sgn;
         }    }
         *rp = q;    *rp = q;
 }  }
   
 void Pdrint(arg,rp)  void Pdrint(arg,rp)
 NODE arg;  NODE arg;
 Q *rp;  Q *rp;
 {  {
         L a;    L a;
         unsigned int au,al;    unsigned int au,al;
         int sgn;    int sgn;
         Q q;    Q q;
         double d;    double d;
   
         if ( !ARG0(arg) ) {    if ( !ARG0(arg) ) {
                 *rp = 0;      *rp = 0;
                 return;      return;
         }    }
 #if defined(VISUAL) || defined(__MINGW32__)  #if defined(VISUAL) || defined(__MINGW32__)
         d = ToReal(ARG0(arg));    d = ToReal(ARG0(arg));
         if ( d > 0 )    if ( d > 0 )
                 d = floor(d+0.5);      d = floor(d+0.5);
         else    else
                 d = ceil(d-0.5);      d = ceil(d-0.5);
 #else  #else
         d = rint(ToReal(ARG0(arg)));    d = rint(ToReal(ARG0(arg)));
 #endif  #endif
         if ( d < -9.223372036854775808e18 || d >= 9.223372036854775808e18 )    if ( d < -9.223372036854775808e18 || d >= 9.223372036854775808e18 )
                 error("drint : OverFlow");      error("drint : OverFlow");
         a = (L)d;    a = (L)d;
         if ( a < 0 ) {    if ( a < 0 ) {
                 sgn = -1;      sgn = -1;
                 a = -a;      a = -a;
         } else    } else
                 sgn = 1;      sgn = 1;
 #if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__MINGW32__) || defined(__x86_64)  #if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__MINGW32__) || defined(__x86_64)
         au = ((unsigned int *)&a)[1];    au = ((unsigned int *)&a)[1];
         al = ((unsigned int *)&a)[0];    al = ((unsigned int *)&a)[0];
 #else  #else
         al = ((unsigned int *)&a)[1];    al = ((unsigned int *)&a)[1];
         au = ((unsigned int *)&a)[0];    au = ((unsigned int *)&a)[0];
 #endif  #endif
         if ( au ) {    if ( au ) {
                 NEWQ(q); SGN(q) = sgn; NM(q)=NALLOC(2); DN(q)=0;      NEWQ(q); SGN(q) = sgn; NM(q)=NALLOC(2); DN(q)=0;
                 PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;      PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;
         } else if ( al ) {    } else if ( al ) {
                 UTOQ(al,q); SGN(q) = sgn;      UTOQ(al,q); SGN(q) = sgn;
         } else    } else
                 q = 0;      q = 0;
         *rp = q;    *rp = q;
 }  }
   
 void Pdisnan(NODE arg,Q *rp)  void Pdisnan(NODE arg,Q *rp)
 {  {
         Real r;    Real r;
         double d;    double d;
 #if defined(VISUAL) || defined(__MINGW32__)  #if defined(VISUAL) || defined(__MINGW32__)
     int c;      int c;
 #endif  #endif
   
         r = (Real)ARG0(arg);    r = (Real)ARG0(arg);
         if ( !r || !NUM(r) || !REAL(r) ) {    if ( !r || !NUM(r) || !REAL(r) ) {
                 *rp = 0;      *rp = 0;
                 return;      return;
         }    }
         d = ToReal(r);    d = ToReal(r);
 #if defined(VISUAL) || defined(__MINGW32__)  #if defined(VISUAL) || defined(__MINGW32__)
         c = _fpclass(d);    c = _fpclass(d);
         if ( c == _FPCLASS_SNAN || c == _FPCLASS_QNAN ) *rp = ONE;    if ( c == _FPCLASS_SNAN || c == _FPCLASS_QNAN ) *rp = ONE;
         else if ( c == _FPCLASS_PINF || c == _FPCLASS_NINF ) STOQ(2,*rp);    else if ( c == _FPCLASS_PINF || c == _FPCLASS_NINF ) STOQ(2,*rp);
 #else  #else
         if ( isnan(d) ) *rp = ONE;    if ( isnan(d) ) *rp = ONE;
         else if ( isinf(d) ) STOQ(2,*rp);    else if ( isinf(d) ) STOQ(2,*rp);
 #endif  #endif
         else *rp = 0;    else *rp = 0;
 }  }
Legend:
Removed from v.1.12  
changed lines
  Added in v.1.13
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>