Magazyn Amiga 13/MA_Cover_13.bin/source/c/apl/ao.c

C/C++ Source or Header | 1999-11-23 | 5.0 KB | 443 lines

#include "apl.h"

data
ex_add(d1, d2)
data d1, d2;
{
    if(fuzz(d1, -d2) == 0) return(zero);
    return(d1 + d2);
}

data
ex_sub(d1, d2)
data d1, d2;
{
    if(fuzz(d1, d2) == 0) return(zero);
    return(d1 - d2);
}

data
ex_mul(d1, d2)
data d1, d2;
{
    return(d1 * d2);
}

data
ex_div(d1, d2)
data d1, d2;
{
    /* 0 div 0 is 1 */
    if(d2 == zero) {
        if (d1 == zero) return(one);
        else error("div D");
    }
    return(d1 / d2);
}

data
ex_mod(d1, d2)
data d1, d2;
{
    double f;

    /* x mod 0 is defined to be x */
    if (d1 == zero) return(d2);
    if(d1 < zero) d1 = -d1;
    f = d2;
    d2 = d2 - d1 * floor(f/d1);
    return(d2);
}

data
ex_min(d1, d2)
data d1, d2;
{
    if(d1 < d2) return(d1);
    return(d2);
}

data
ex_max(d1, d2)
data d1, d2;
{
    if(d1 > d2) return(d1);
    return(d2);
}

data
ex_pwr(d1, d2)
data d1, d2;
{
    int s;
    double f1, f2;

    s = 0;
    f1 = d1;
    if(f1 > 0.) {
        f1 = d2 * log(f1);
        goto chk;
    }
    if(f1 == 0.) return(d2 == zero ? (data)1.0 : zero);

    /* check for integer exponent */
    f2 = floor(d2);
    if(fabs(d2 - f2) < thread.fuzz){
        s = (int)f2%2;
        f1 = d2*log(fabs(f1));
        goto chk;
    }
    /* should check rational d2 here */
    goto bad;

chk:
    if(f1 < maxexp) {
        d1 = exp(f1);
        if(s) d1 = -d1;
        return(d1);
    }
bad:
    error("pwr D");
}

data
ex_log(d1, d2)
data d1, d2;
{
    double f1, f2;

    f1 = d1;
    f2 = d2;
    if(f1 <= 0. || f2 <= 0.) error("log D");
    d1 = log(f2)/log(f1);
    return(d1);
}

data
ex_cir(d1, d2)
data d1, d2;
{
    double f, f1;

    switch(fix(d1)) {

    default:
        error("crc D");

    case -7:
        /* arc tanh */
        f = d2;
        if(f <= -1. || f >= 1.) error("tanh D");
        f = log((1.+f)/(1.-f))/2.;
        goto ret;

    case -6:
        /* arc cosh */
        f = d2;
        if(f < 1.) error("cosh D");
        f = log(f+sqrt(f*f-1.));
        goto ret;

    case -5:
        /* arc sinh */
        f = d2;
        f = log(f+sqrt(f*f+1.));
        goto ret;

    case -4:
        /* sqrt(-1 + x*x) */
        f = -one + d2*d2;
        goto sq;

    case -3:
        /* arc tan */
        f = d2;
        f = atan(f);
        goto ret;

    case -2:
        /* arc cos */
        f = d2;
        if(f < -1. || f > 1.) error("arc cos D");
        f = atan2(sqrt(1.-f*f), f);
        goto ret;

    case -1:
        /* arc sin */
        f = d2;
        if(f < -1. || f > 1.) error("arc sin D");
        f = atan2(f, sqrt(1.-f*f));
        goto ret;

    case 0:
        /* sqrt(1 - x*x) */
        f = one - d2*d2;
        goto sq;

    case 1:
        /* sin */
        f = d2;
        f = sin(f);
        goto ret;

    case 2:
        /* cos */
        f = d2;
        f = cos(f);
        goto ret;

    case 3:
        /* tan */
        f = d2;
        f1 = cos(f);
        if(f1 == 0.) f = exp(maxexp);
        else f = sin(f)/f1;
        goto ret;

    case 4:
        /* sqrt(1 + x*x) */
        f = one + d2*d2;
        goto sq;

    case 5:
        /* sinh */
        f = d2;
        if(f < -maxexp || f > maxexp) error("sinh D");
        f = exp(f);
        f = (f-1./f)/2.;
        goto ret;

    case 6:
        /* cosh */
        f = d2;
        if(f < -maxexp || f > maxexp) error("cosh D");
        f = exp(f);
        f = (f+1./f)/2.;
        goto ret;

    case 7:
        /* tanh */
        f = d2;
        if(f > maxexp) {
            f = 1.;
            goto ret;
        }
        if(f < -maxexp) {
            f = -1.;
            goto ret;
        }
        f = exp(f);
        f = (f-1./f)/(f+1./f);
        goto ret;
    }

sq:
    if(f < 0.) error("sqrt D");
    f = sqrt(f);

ret:
    d1 = f;
    return(d1);
}

data
ex_comb(d1, d2)
data d1, d2;
{
    double f;

    if(d1 > d2) return(zero);
    f = gamma(d2+1.) - gamma(d1+1.) - gamma(d2-d1+1.);
    if(f > maxexp) error("comb D");
    d1 = exp(f);
    return(d1);
}

data
ex_and(d1, d2)
data d1, d2;
{
    if(d1!=zero && d2!=zero) return(one);
    return(zero);
}

data
ex_or(d1, d2)
data d1, d2;
{
    if(d1!=zero || d2!=zero) return(one);
    return(zero);
}

data
ex_nand(d1, d2)
data d1, d2;
{
    if(d1!=zero && d2!=zero) return(zero);
    return(one);
}

data
ex_nor(d1, d2)
data d1, d2;
{
    if(d1!=zero || d2!=zero) return(zero);
    return(one);
}

data
ex_lt(d1, d2)
data d1, d2;
{
    if(fuzz(d1, d2) < 0) return(one);
    return(zero);
}

data
ex_le(d1, d2)
data d1, d2;
{
    if(fuzz(d1, d2) <= 0) return(one);
    return(zero);
}

data
ex_eq(d1, d2)
data d1, d2;
{
    if(fuzz(d1, d2) == 0) return(one);
    return(zero);
}

data
ex_ge(d1, d2)
data d1, d2;
{
    if(fuzz(d1, d2) >= 0) return(one);
    return(zero);
}

data
ex_gt(d1, d2)
data d1, d2;
{
    if(fuzz(d1, d2) > 0) return(one);
    return(zero);
}

data
ex_ne(d1, d2)
data d1, d2;
{
    if(fuzz(d1, d2) != 0) return(one);
    return(zero);
}

data
ex_plus(d)
data d;
{
    return(d);
}

data
ex_minus(d)
data d;
{
    return(-d);
}

data
ex_sgn(d)
data d;
{
    if(d == zero) return(zero);
    if(d < zero) return(-one);
    return(one);
}

data
ex_recip(d)
data d;
{
    if(d == zero) error("recip D");
    return(one/d);
}

data
ex_abs(d)
data d;
{
    if(d < zero) return(-d);
    return(d);
}

data
ex_floor(d)
data d;
{
    d = floor(d + thread.fuzz);
    return(d);
}

data
ex_ceil(d)
data d;
{
    d = ceil(d - thread.fuzz);
    return(d);
}

data
ex_exp(d)
data d;
{
    double f;

    f = d;
    if(f > maxexp) error ("exp D");
    d = exp(f);
    return(d);
}

data
ex_loge(d)
data d;
{
    double f;

    f = d;
    if(f <= 0.) error("log D");
    d = log(f);
    return(d);
}

data
ex_pi(d)
data d;
{
    d = pi * d;
    return(d);
}

data
ex_fac(d)
data d;
{
    double f;

    f = gamma(d+1.);
    if(f > maxexp) error("fac D");
    d = exp(f);
    if(signgam < 0) d = -d;            /* if (signgam) in version 6 */
    return(d);
}

data
ex_not(d)
data d;
{
    if(d == zero) return(one);
    return(zero);
}