4.3BSD - /usr/src/new/apl/src/ao.c

   1: static char Sccsid[] = "ao.c @(#)ao.c	1.1	10/1/82 Berkeley ";
   2: #include "apl.h"
   3: 
   4: data
   5: ex_add(d1, d2)
   6: data d1, d2;
   7: {
   8: 
   9:     if(fuzz(d1, -d2) == 0)
  10:         return(zero);
  11:     return(d1 + d2);
  12: }
  13: 
  14: data
  15: ex_sub(d1, d2)
  16: data d1, d2;
  17: {
  18: 
  19:     if(fuzz(d1, d2) == 0)
  20:         return(zero);
  21:     return(d1 - d2);
  22: }
  23: 
  24: data
  25: ex_mul(d1, d2)
  26: data d1, d2;
  27: {
  28: 
  29:     return(d1 * d2);
  30: }
  31: 
  32: data
  33: ex_div(d1, d2)
  34: data d1, d2;
  35: {
  36: 
  37:     /* 0 div 0 is 1 */
  38:     if(d2 == zero)
  39:         if (d1 == zero) return(one);
  40:             else error("div D");
  41:     return(d1 / d2);
  42: }
  43: 
  44: data
  45: ex_mod(d1, d2)
  46: data d1, d2;
  47: {
  48:     double f;
  49: 
  50:     /* x mod 0 is defined to be x */
  51:     if (d1 == zero)
  52:         return(d2);
  53:     if(d1 < zero)
  54:         d1 = -d1;
  55:     f = d2;
  56:     d2 = d2 - d1 * floor(f/d1);
  57:     return(d2);
  58: }
  59: 
  60: data
  61: ex_min(d1, d2)
  62: data d1, d2;
  63: {
  64: 
  65:     if(d1 < d2)
  66:         return(d1);
  67:     return(d2);
  68: }
  69: 
  70: data
  71: ex_max(d1, d2)
  72: data d1, d2;
  73: {
  74: 
  75:     if(d1 > d2)
  76:         return(d1);
  77:     return(d2);
  78: }
  79: 
  80: data
  81: ex_pwr(d1, d2)
  82: data d1, d2;
  83: {
  84:     register s;
  85:     double f1, f2;
  86: 
  87:     s = 0;
  88:     f1 = d1;
  89:     if(f1 > 0.) {
  90:         f1 = d2 * log(f1);
  91:         goto chk;
  92:     }
  93:     if(f1 == 0.) {
  94:         return(d2 == zero ? (data)1.0 : zero);
  95:     }
  96:     /* check for integer exponent */
  97:     f2 = floor(d2);
  98:     if(fabs(d2 - f2) < thread.fuzz){
  99:         s = (int)f2%2;
 100:         f1 = d2*log(fabs(f1));
 101:         goto chk;
 102:     }
 103:     /* should check rational d2 here */
 104:     goto bad;
 105: 
 106: chk:
 107:     if(f1 < maxexp) {
 108:         d1 = exp(f1);
 109:         if(s)
 110:             d1 = -d1;
 111:         return(d1);
 112:     }
 113: bad:
 114:     error("pwr D");
 115: }
 116: 
 117: data
 118: ex_log(d1, d2)
 119: data d1, d2;
 120: {
 121:     double f1, f2;
 122: 
 123:     f1 = d1;
 124:     f2 = d2;
 125:     if(f1 <= 0. || f2 <= 0.)
 126:         error("log D");
 127:     d1 = log(f2)/log(f1);
 128:     return(d1);
 129: }
 130: 
 131: data
 132: ex_cir(d1, d2)
 133: data d1, d2;
 134: {
 135:     double f, f1;
 136: 
 137:     switch(fix(d1)) {
 138: 
 139:     default:
 140:         error("crc D");
 141: 
 142:     case -7:
 143:         /* arc tanh */
 144:         f = d2;
 145:         if(f <= -1. || f >= 1.)
 146:             error("tanh D");
 147:         f = log((1.+f)/(1.-f))/2.;
 148:         goto ret;
 149: 
 150:     case -6:
 151:         /* arc cosh */
 152:         f = d2;
 153:         if(f < 1.)
 154:             error("cosh D");
 155:         f = log(f+sqrt(f*f-1.));
 156:         goto ret;
 157: 
 158:     case -5:
 159:         /* arc sinh */
 160:         f = d2;
 161:         f = log(f+sqrt(f*f+1.));
 162:         goto ret;
 163: 
 164:     case -4:
 165:         /* sqrt(-1 + x*x) */
 166:         f = -one + d2*d2;
 167:         goto sq;
 168: 
 169:     case -3:
 170:         /* arc tan */
 171:         f = d2;
 172:         f = atan(f);
 173:         goto ret;
 174: 
 175:     case -2:
 176:         /* arc cos */
 177:         f = d2;
 178:         if(f < -1. || f > 1.)
 179:             error("arc cos D");
 180:         f = atan2(sqrt(1.-f*f), f);
 181:         goto ret;
 182: 
 183:     case -1:
 184:         /* arc sin */
 185:         f = d2;
 186:         if(f < -1. || f > 1.)
 187:             error("arc sin D");
 188:         f = atan2(f, sqrt(1.-f*f));
 189:         goto ret;
 190: 
 191:     case 0:
 192:         /* sqrt(1 - x*x) */
 193:         f = one - d2*d2;
 194:         goto sq;
 195: 
 196:     case 1:
 197:         /* sin */
 198:         f = d2;
 199:         f = sin(f);
 200:         goto ret;
 201: 
 202:     case 2:
 203:         /* cos */
 204:         f = d2;
 205:         f = cos(f);
 206:         goto ret;
 207: 
 208:     case 3:
 209:         /* tan */
 210:         f = d2;
 211:         f1 = cos(f);
 212:         if(f1 == 0.)
 213:             f = exp(maxexp); else
 214:             f = sin(f)/f1;
 215:         goto ret;
 216: 
 217:     case 4:
 218:         /* sqrt(1 + x*x) */
 219:         f = one + d2*d2;
 220:         goto sq;
 221: 
 222:     case 5:
 223:         /* sinh */
 224:         f = d2;
 225:         if(f < -maxexp || f > maxexp)
 226:             error("sinh D");
 227:         f = exp(f);
 228:         f = (f-1./f)/2.;
 229:         goto ret;
 230: 
 231:     case 6:
 232:         /* cosh */
 233:         f = d2;
 234:         if(f < -maxexp || f > maxexp)
 235:             error("cosh D");
 236:         f = exp(f);
 237:         f = (f+1./f)/2.;
 238:         goto ret;
 239: 
 240:     case 7:
 241:         /* tanh */
 242:         f = d2;
 243:         if(f > maxexp) {
 244:             f = 1.;
 245:             goto ret;
 246:         }
 247:         if(f < -maxexp) {
 248:             f = -1.;
 249:             goto ret;
 250:         }
 251:         f = exp(f);
 252:         f = (f-1./f)/(f+1./f);
 253:         goto ret;
 254:     }
 255: 
 256: sq:
 257:     if(f < 0.)
 258:         error("sqrt D");
 259:     f = sqrt(f);
 260: 
 261: ret:
 262:     d1 = f;
 263:     return(d1);
 264: }
 265: 
 266: data
 267: ex_comb(d1, d2)
 268: data d1, d2;
 269: {
 270:     double f;
 271: 
 272:     if(d1 > d2)
 273:         return(zero);
 274:     f = gamma(d2+1.) - gamma(d1+1.) - gamma(d2-d1+1.);
 275:     if(f > maxexp)
 276:         error("comb D");
 277:     d1 = exp(f);
 278:     return(d1);
 279: }
 280: 
 281: data
 282: ex_and(d1, d2)
 283: data d1, d2;
 284: {
 285: 
 286:     if(d1!=zero && d2!=zero)
 287:         return(one);
 288:     return(zero);
 289: }
 290: 
 291: data
 292: ex_or(d1, d2)
 293: data d1, d2;
 294: {
 295: 
 296:     if(d1!=zero || d2!=zero)
 297:         return(one);
 298:     return(zero);
 299: }
 300: 
 301: data
 302: ex_nand(d1, d2)
 303: data d1, d2;
 304: {
 305: 
 306:     if(d1!=zero && d2!=zero)
 307:         return(zero);
 308:     return(one);
 309: }
 310: 
 311: data
 312: ex_nor(d1, d2)
 313: data d1, d2;
 314: {
 315: 
 316:     if(d1!=zero || d2!=zero)
 317:         return(zero);
 318:     return(one);
 319: }
 320: 
 321: data
 322: ex_lt(d1, d2)
 323: data d1, d2;
 324: {
 325: 
 326:     if(fuzz(d1, d2) < 0)
 327:         return(one);
 328:     return(zero);
 329: }
 330: 
 331: data
 332: ex_le(d1, d2)
 333: data d1, d2;
 334: {
 335: 
 336:     if(fuzz(d1, d2) <= 0)
 337:         return(one);
 338:     return(zero);
 339: }
 340: 
 341: data
 342: ex_eq(d1, d2)
 343: data d1, d2;
 344: {
 345: 
 346:     if(fuzz(d1, d2) == 0)
 347:         return(one);
 348:     return(zero);
 349: }
 350: 
 351: data
 352: ex_ge(d1, d2)
 353: data d1, d2;
 354: {
 355: 
 356:     if(fuzz(d1, d2) >= 0)
 357:         return(one);
 358:     return(zero);
 359: }
 360: 
 361: data
 362: ex_gt(d1, d2)
 363: data d1, d2;
 364: {
 365: 
 366:     if(fuzz(d1, d2) > 0)
 367:         return(one);
 368:     return(zero);
 369: }
 370: 
 371: data
 372: ex_ne(d1, d2)
 373: data d1, d2;
 374: {
 375: 
 376:     if(fuzz(d1, d2) != 0)
 377:         return(one);
 378:     return(zero);
 379: }
 380: 
 381: data
 382: ex_plus(d)
 383: data d;
 384: {
 385: 
 386:     return(d);
 387: }
 388: 
 389: data
 390: ex_minus(d)
 391: data d;
 392: {
 393: 
 394:     return(-d);
 395: }
 396: 
 397: data
 398: ex_sgn(d)
 399: data d;
 400: {
 401: 
 402:     if(d == zero)
 403:         return(zero);
 404:     if(d < zero)
 405:         return(-one);
 406:     return(one);
 407: }
 408: 
 409: data
 410: ex_recip(d)
 411: data d;
 412: {
 413: 
 414:     if(d == zero)
 415:         error("recip D");
 416:     return(one/d);
 417: }
 418: 
 419: data
 420: ex_abs(d)
 421: data d;
 422: {
 423: 
 424:     if(d < zero)
 425:         return(-d);
 426:     return(d);
 427: }
 428: 
 429: data
 430: ex_floor(d)
 431: data d;
 432: {
 433: 
 434:     d = floor(d + thread.fuzz);
 435:     return(d);
 436: }
 437: 
 438: data
 439: ex_ceil(d)
 440: data d;
 441: {
 442: 
 443:     d = ceil(d - thread.fuzz);
 444:     return(d);
 445: }
 446: 
 447: data
 448: ex_exp(d)
 449: data d;
 450: {
 451:     double f;
 452: 
 453:     f = d;
 454:     if(f > maxexp)
 455:         error ("exp D");
 456:     d = exp(f);
 457:     return(d);
 458: }
 459: 
 460: data
 461: ex_loge(d)
 462: data d;
 463: {
 464:     double f;
 465: 
 466:     f = d;
 467:     if(f <= 0.)
 468:         error("log D");
 469:     d = log(f);
 470:     return(d);
 471: }
 472: 
 473: data
 474: ex_pi(d)
 475: data d;
 476: {
 477: 
 478:     d = pi * d;
 479:     return(d);
 480: }
 481: 
 482: /* ex_rand has been moved to af.c (with ex_deal) */
 483: 
 484: data
 485: ex_fac(d)
 486: data d;
 487: {
 488:     double f;
 489: 
 490:     f = gamma(d+1.);
 491:     if(f > maxexp)
 492:         error("fac D");
 493:     d = exp(f);
 494:     if(signgam < 0)         /* if (signgam) in version 6 */
 495:         d = -d;
 496:     return(d);
 497: }
 498: 
 499: data
 500: ex_not(d)
 501: data d;
 502: {
 503: 
 504:     if(d == zero)
 505:         return(one);
 506:     return(zero);
 507: }

Defined functions

Defined variables