Annotation of OpenXM_contrib2/asir2000/builtin/math.c, Revision 1.13
1.2 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.3 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.13 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/math.c,v 1.12 2015/08/14 13:51:54 fujimoto Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include <math.h>
52: #include "parse.h"
1.12 fujimoto 53: #if defined(VISUAL) || defined(__MINGW32__)
1.10 ohara 54: #include <float.h>
55: #endif
1.1 noro 56:
57: void Pdsqrt(),Pdsin(),Pdcos(),Pdtan(),Pdasin(),Pdacos(),Pdatan(),Pdlog(),Pdexp();
1.9 noro 58: void Pabs(),Pdfloor(),Pdceil(),Pdrint(),Pdisnan();
1.1 noro 59:
60: struct ftab math_tab[] = {
1.13 ! noro 61: {"dsqrt",Pdsqrt,1},
! 62: {"dabs",Pabs,1},
! 63: {"dsin",Pdsin,1},
! 64: {"dcos",Pdcos,1},
! 65: {"dtan",Pdtan,1},
! 66: {"dlog",Pdlog,1},
! 67: {"dexp",Pdexp,1},
! 68: {"dasin",Pdasin,1},
! 69: {"dacos",Pdacos,1},
! 70: {"datan",Pdatan,1},
! 71: {"floor",Pdfloor,1},
! 72: {"dfloor",Pdfloor,1},
! 73: {"ceil",Pdceil,1},
! 74: {"dceil",Pdceil,1},
! 75: {"rint",Pdrint,1},
! 76: {"drint",Pdrint,1},
! 77: {"disnan",Pdisnan,1},
! 78: {0,0,0},
1.1 noro 79: };
80:
1.8 noro 81: void get_ri(Num z,double *r,double *i)
82: {
1.13 ! noro 83: if ( !z ) {
! 84: *r = 0; *i = 0; return;
! 85: }
! 86: if ( OID(z) != O_N )
! 87: error("get_ri : invalid argument");
! 88: switch ( NID(z) ) {
! 89: case N_Q: case N_R: case N_B:
! 90: *r = ToReal(z); *i = 0;
! 91: break;
! 92: case N_C:
! 93: *r = ToReal(((C)z)->r);
! 94: *i = ToReal(((C)z)->i);
! 95: break;
! 96: default:
! 97: error("get_ri : invalid argument");
! 98: break;
! 99: }
1.8 noro 100: }
101:
1.1 noro 102: void Pabs(arg,rp)
103: NODE arg;
104: Real *rp;
105: {
1.13 ! noro 106: double s,r,i;
1.1 noro 107:
1.13 ! noro 108: if ( !ARG0(arg) ) {
! 109: *rp = 0; return;
! 110: }
! 111: get_ri((Num)ARG0(arg),&r,&i);
! 112: if ( i == 0 )
! 113: s = fabs(r);
! 114: else if ( r == 0 )
! 115: s = fabs(i);
! 116: else
! 117: s = sqrt(r*r+i*i);
! 118: MKReal(s,*rp);
1.1 noro 119: }
120:
121: void Pdsqrt(arg,rp)
122: NODE arg;
1.8 noro 123: Num *rp;
1.1 noro 124: {
1.13 ! noro 125: double s,r,i,a;
! 126: C z;
! 127: Real real;
! 128:
! 129: if ( !ARG0(arg) ) {
! 130: *rp = 0; return;
! 131: }
! 132: get_ri((Num)ARG0(arg),&r,&i);
! 133: if ( i == 0 )
! 134: if ( r > 0 ) {
! 135: s = sqrt(r);
! 136: MKReal(s,real);
! 137: *rp = (Num)real;
! 138: } else {
! 139: NEWC(z);
! 140: z->r = 0;
! 141: s = sqrt(-r); MKReal(s,real); z->i = (Num)real;
! 142: *rp = (Num)z;
! 143: }
! 144: else {
! 145: a = sqrt(r*r+i*i);
! 146: NEWC(z);
! 147: s = sqrt((r+a)/2); MKReal(s,real); z->r = (Num)real;
! 148: s = i>0?sqrt((-r+a)/2):-sqrt((-r+a)/2);
! 149: MKReal(s,real); z->i = (Num)real;
! 150: *rp = (Num)z;
! 151: }
1.1 noro 152: }
153:
154: void Pdsin(arg,rp)
155: NODE arg;
156: Real *rp;
157: {
1.13 ! noro 158: double s;
1.1 noro 159:
1.13 ! noro 160: s = sin(ToReal(ARG0(arg)));
! 161: MKReal(s,*rp);
1.1 noro 162: }
163:
164: void Pdcos(arg,rp)
165: NODE arg;
166: Real *rp;
167: {
1.13 ! noro 168: double s;
1.1 noro 169:
1.13 ! noro 170: s = cos(ToReal(ARG0(arg)));
! 171: MKReal(s,*rp);
1.1 noro 172: }
173:
174: void Pdtan(arg,rp)
175: NODE arg;
176: Real *rp;
177: {
1.13 ! noro 178: double s;
1.1 noro 179:
1.13 ! noro 180: s = tan(ToReal(ARG0(arg)));
! 181: MKReal(s,*rp);
1.1 noro 182: }
183:
184: void Pdasin(arg,rp)
185: NODE arg;
186: Real *rp;
187: {
1.13 ! noro 188: double s;
1.1 noro 189:
1.13 ! noro 190: s = asin(ToReal(ARG0(arg)));
! 191: MKReal(s,*rp);
1.1 noro 192: }
193:
194: void Pdacos(arg,rp)
195: NODE arg;
196: Real *rp;
197: {
1.13 ! noro 198: double s;
1.1 noro 199:
1.13 ! noro 200: s = acos(ToReal(ARG0(arg)));
! 201: MKReal(s,*rp);
1.1 noro 202: }
203:
204: void Pdatan(arg,rp)
205: NODE arg;
206: Real *rp;
207: {
1.13 ! noro 208: double s;
1.1 noro 209:
1.13 ! noro 210: s = atan(ToReal(ARG0(arg)));
! 211: MKReal(s,*rp);
1.1 noro 212: }
213:
214: void Pdlog(arg,rp)
215: NODE arg;
216: Real *rp;
217: {
1.13 ! noro 218: double s;
1.1 noro 219:
1.13 ! noro 220: s = log(ToReal(ARG0(arg)));
! 221: MKReal(s,*rp);
1.1 noro 222: }
223:
224: void Pdexp(arg,rp)
225: NODE arg;
226: Real *rp;
227: {
1.13 ! noro 228: double s;
1.1 noro 229:
1.13 ! noro 230: s = exp(ToReal(ARG0(arg)));
! 231: MKReal(s,*rp);
1.1 noro 232: }
233:
234: void Pdfloor(arg,rp)
235: NODE arg;
236: Q *rp;
237: {
1.13 ! noro 238: L a;
! 239: unsigned int au,al;
! 240: int sgn;
! 241: Q q;
! 242: double d;
! 243:
! 244: if ( !ARG0(arg) ) {
! 245: *rp = 0;
! 246: return;
! 247: }
! 248: d = floor(ToReal(ARG0(arg)));
! 249: if ( d < -9.223372036854775808e18 || d >= 9.223372036854775808e18 )
! 250: error("dfloor : OverFlow");
! 251: if ( !d ) {
! 252: *rp = 0;
! 253: return;
! 254: }
! 255: a = (L)d;
! 256: if ( a < 0 ) {
! 257: sgn = -1;
! 258: a = -a;
! 259: } else
! 260: sgn = 1;
1.12 fujimoto 261: #if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__MINGW32__) || defined(__x86_64)
1.13 ! noro 262: au = ((unsigned int *)&a)[1];
! 263: al = ((unsigned int *)&a)[0];
1.1 noro 264: #else
1.13 ! noro 265: al = ((unsigned int *)&a)[1];
! 266: au = ((unsigned int *)&a)[0];
1.1 noro 267: #endif
1.13 ! noro 268: if ( au ) {
! 269: NEWQ(q); SGN(q) = sgn; NM(q)=NALLOC(2); DN(q)=0;
! 270: PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;
! 271: } else {
! 272: UTOQ(al,q); SGN(q) = sgn;
! 273: }
! 274: *rp = q;
1.1 noro 275: }
276:
277: void Pdceil(arg,rp)
278: NODE arg;
279: Q *rp;
280: {
1.13 ! noro 281: L a;
! 282: unsigned int au,al;
! 283: int sgn;
! 284: Q q;
! 285: double d;
! 286:
! 287: if ( !ARG0(arg) ) {
! 288: *rp = 0;
! 289: return;
! 290: }
! 291: d = ceil(ToReal(ARG0(arg)));
! 292: if ( d < -9.223372036854775808e18 || d >= 9.223372036854775808e18 )
! 293: error("dceil : OverFlow");
! 294: if ( !d ) {
! 295: *rp = 0;
! 296: return;
! 297: }
! 298: a = (L)d;
! 299: if ( a < 0 ) {
! 300: sgn = -1;
! 301: a = -a;
! 302: } else
! 303: sgn = 1;
1.12 fujimoto 304: #if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__MINGW32__) || defined(__x86_64)
1.13 ! noro 305: au = ((unsigned int *)&a)[1];
! 306: al = ((unsigned int *)&a)[0];
1.1 noro 307: #else
1.13 ! noro 308: al = ((unsigned int *)&a)[1];
! 309: au = ((unsigned int *)&a)[0];
1.1 noro 310: #endif
1.13 ! noro 311: if ( au ) {
! 312: NEWQ(q); SGN(q) = sgn; NM(q)=NALLOC(2); DN(q)=0;
! 313: PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;
! 314: } else {
! 315: UTOQ(al,q); SGN(q) = sgn;
! 316: }
! 317: *rp = q;
1.1 noro 318: }
319:
320: void Pdrint(arg,rp)
321: NODE arg;
322: Q *rp;
323: {
1.13 ! noro 324: L a;
! 325: unsigned int au,al;
! 326: int sgn;
! 327: Q q;
! 328: double d;
! 329:
! 330: if ( !ARG0(arg) ) {
! 331: *rp = 0;
! 332: return;
! 333: }
1.12 fujimoto 334: #if defined(VISUAL) || defined(__MINGW32__)
1.13 ! noro 335: d = ToReal(ARG0(arg));
! 336: if ( d > 0 )
! 337: d = floor(d+0.5);
! 338: else
! 339: d = ceil(d-0.5);
1.1 noro 340: #else
1.13 ! noro 341: d = rint(ToReal(ARG0(arg)));
1.1 noro 342: #endif
1.13 ! noro 343: if ( d < -9.223372036854775808e18 || d >= 9.223372036854775808e18 )
! 344: error("drint : OverFlow");
! 345: a = (L)d;
! 346: if ( a < 0 ) {
! 347: sgn = -1;
! 348: a = -a;
! 349: } else
! 350: sgn = 1;
1.12 fujimoto 351: #if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__MINGW32__) || defined(__x86_64)
1.13 ! noro 352: au = ((unsigned int *)&a)[1];
! 353: al = ((unsigned int *)&a)[0];
1.1 noro 354: #else
1.13 ! noro 355: al = ((unsigned int *)&a)[1];
! 356: au = ((unsigned int *)&a)[0];
1.1 noro 357: #endif
1.13 ! noro 358: if ( au ) {
! 359: NEWQ(q); SGN(q) = sgn; NM(q)=NALLOC(2); DN(q)=0;
! 360: PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;
! 361: } else if ( al ) {
! 362: UTOQ(al,q); SGN(q) = sgn;
! 363: } else
! 364: q = 0;
! 365: *rp = q;
1.1 noro 366: }
1.9 noro 367:
368: void Pdisnan(NODE arg,Q *rp)
369: {
1.13 ! noro 370: Real r;
! 371: double d;
1.12 fujimoto 372: #if defined(VISUAL) || defined(__MINGW32__)
1.10 ohara 373: int c;
374: #endif
1.9 noro 375:
1.13 ! noro 376: r = (Real)ARG0(arg);
! 377: if ( !r || !NUM(r) || !REAL(r) ) {
! 378: *rp = 0;
! 379: return;
! 380: }
! 381: d = ToReal(r);
1.12 fujimoto 382: #if defined(VISUAL) || defined(__MINGW32__)
1.13 ! noro 383: c = _fpclass(d);
! 384: if ( c == _FPCLASS_SNAN || c == _FPCLASS_QNAN ) *rp = ONE;
! 385: else if ( c == _FPCLASS_PINF || c == _FPCLASS_NINF ) STOQ(2,*rp);
1.10 ohara 386: #else
1.13 ! noro 387: if ( isnan(d) ) *rp = ONE;
! 388: else if ( isinf(d) ) STOQ(2,*rp);
1.10 ohara 389: #endif
1.13 ! noro 390: else *rp = 0;
1.9 noro 391: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>