4.3BSD - /usr/src/new/B/src/bint/b1nuR.c

   1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
   2: 
   3: /*
   4:   $Header: b1nuR.c,v 1.4 85/08/22 16:51:49 timo Exp $
   5: */
   6: 
   7: /* Rational arithmetic */
   8: 
   9: #include "b.h"
  10: #include "b0con.h"
  11: #include "b1obj.h"
  12: #include "b1num.h"
  13: #include "b3err.h"
  14: 
  15: /* Length calculations used for fraction sizes: */
  16: #define Maxlen(u, v) \
  17:     (Roundsize(u) > Roundsize(v) ? Roundsize(u) : Roundsize(v))
  18: #define Sumlen(u, v) (Roundsize(u)+Roundsize(v))
  19: #define Difflen(u, v) (Roundsize(u)-Roundsize(v))
  20: 
  21: /* To shut off lint and other warnings: */
  22: #undef Copy
  23: #define Copy(x) ((integer)copy((value)(x)))
  24: 
  25: /* Globally used constants */
  26: 
  27: rational rat_zero;
  28: rational rat_half;
  29: 
  30: /* Make a normalized rational from two integers */
  31: 
  32: Visible rational mk_rat(x, y, len) integer x, y; int len; {
  33:     rational a;
  34:     integer u,v;
  35: 
  36:     if (y == int_0) {
  37:         if (interrupted)
  38:             return rat_zero;
  39:         syserr(MESS(1200, "mk_rat(x, y) with y=0"));
  40:     }
  41: 
  42:     if (x == int_0 && len <= 0) return (rational) Copy(rat_zero);
  43: 
  44:     if (Msd(y) < 0) {   /* interchange signs */
  45:         u = int_neg(x);
  46:         v = int_neg(y);
  47:     } else {
  48:         u = Copy(x);
  49:         v = Copy(y);
  50:     }
  51: 
  52:     a = (rational) grab_rat();
  53:     if (len > 0 && len+2 <= Maxintlet) Length(a) = -2 - len;
  54: 
  55:     if (u == int_0 || v == int_1) {
  56:         /* No simplification possible */
  57:         Numerator(a) = Copy(u);
  58:         Denominator(a) = int_1;
  59:     } else {
  60:         integer g, abs_u;
  61: 
  62:         if (Msd(u) < 0) abs_u = int_neg(u);
  63:         else abs_u = Copy(u);
  64:         g = int_gcd(abs_u, v);
  65:         release((value) abs_u);
  66: 
  67:         if (g != int_1) {
  68:             Numerator(a) = int_quot(u, g);
  69:             Denominator(a) = int_quot(v, g);
  70:         } else {
  71:             Numerator(a) = Copy(u);
  72:             Denominator(a) = Copy(v);
  73:         }
  74:         release((value) g);
  75:     }
  76: 
  77:     release((value) u); release((value) v);
  78: 
  79:     return a;
  80: }
  81: 
  82: 
  83: /* Arithmetic on rational numbers */
  84: 
  85: /* Shorthands: */
  86: #define N(u) Numerator(u)
  87: #define D(u) Denominator(u)
  88: 
  89: Visible rational rat_sum(u, v) register rational u, v; {
  90:     integer t1, t2, t3, t4;
  91:     rational a;
  92: 
  93:     t2= int_prod(N(u), D(v));
  94:     t3= int_prod(N(v), D(u));
  95:     t1= int_sum(t2, t3);
  96:     t4= int_prod(D(u), D(v));
  97:     a= mk_rat(t1, t4, Maxlen(u, v));
  98:     release((value) t1); release((value) t2);
  99:     release((value) t3); release((value) t4);
 100: 
 101:     return a;
 102: }
 103: 
 104: 
 105: Visible rational rat_diff(u, v) register rational u, v; {
 106:     integer t1, t2, t3, t4;
 107:     rational a;
 108: 
 109:     t2= int_prod(N(u), D(v));
 110:     t3= int_prod(N(v), D(u));
 111:     t1= int_diff(t2, t3);
 112:     t4= int_prod(D(u), D(v));
 113:     a= mk_rat(t1, t4, Maxlen(u, v));
 114:     release((value) t1); release((value) t2);
 115:     release((value) t3); release((value) t4);
 116: 
 117:     return a;
 118: }
 119: 
 120: 
 121: Visible rational rat_prod(u, v) register rational u, v; {
 122:     integer t1, t2;
 123:     rational a;
 124: 
 125:     t1= int_prod(N(u), N(v));
 126:     t2= int_prod(D(u), D(v));
 127:     a= mk_rat(t1, t2, Sumlen(u, v));
 128:     release((value) t1); release((value) t2);
 129: 
 130:     return a;
 131: }
 132: 
 133: 
 134: Visible rational rat_quot(u, v) register rational u, v; {
 135:     integer t1, t2;
 136:     rational a;
 137: 
 138:     if (Numerator(v) == int_0) {
 139:         error(MESS(1201, "in u/v, v is zero"));
 140:         return (rational) Copy(rat_zero);
 141:     }
 142: 
 143:     t1= int_prod(N(u), D(v));
 144:     t2= int_prod(D(u), N(v));
 145:     a= mk_rat(t1, t2, Difflen(u, v));
 146:     release((value) t1); release((value) t2);
 147: 
 148:     return a;
 149: }
 150: 
 151: 
 152: Visible rational rat_neg(u) register rational u; {
 153:     register rational a;
 154: 
 155:     /* Avoid a real subtraction from zero */
 156: 
 157:     if (Numerator(u) == int_0) return (rational) Copy(u);
 158: 
 159:     a = (rational) grab_rat();
 160:     N(a) = int_neg(N(u));
 161:     D(a) = Copy(D(u));
 162:     Length(a) = Length(u);
 163: 
 164:     return a;
 165: }
 166: 
 167: 
 168: /* Rational number to the integral power */
 169: 
 170: Visible rational rat_power(a, n) rational a; integer n; {
 171:     integer u, v, tu, tv, temp;
 172: 
 173:     if (n == int_0) return mk_rat(int_1, int_1, 0);
 174: 
 175:     if (Msd(n) < 0) {
 176:         if (Numerator(a) == int_0) {
 177:             error(MESS(1202, "in 0**n, n is negative"));
 178:             return (rational) Copy(a);
 179:         }
 180:         if (Msd(Numerator(a)) < 0) {
 181:             u= int_neg(Denominator(a));
 182:             v = int_neg(Numerator(a));
 183:         }
 184:         else {
 185:             u = Copy(Denominator(a));
 186:             v = Copy(Numerator(a));
 187:         }
 188:         n = int_neg(n);
 189:     } else {
 190:         if (Numerator(a) == int_0) return (rational) Copy(a);
 191:             /* To avoid necessary simplification later on */
 192:         u = Copy(Numerator(a));
 193:         v = Copy(Denominator(a));
 194:         n = Copy(n);
 195:     }
 196: 
 197:     tu = int_1;
 198:     tv = int_1;
 199: 
 200:     while (n != int_0 && !interrupted) {
 201:         if (Odd(Lsd(n))) {
 202:             if (u != int_1) {
 203:                 temp = tu;
 204:                 tu = int_prod(u, tu);
 205:                 release((value) temp);
 206:             }
 207:             if (v != int_1) {
 208:                 temp = tv;
 209:                 tv = int_prod(v, tv);
 210:                 release((value) temp);
 211:             }
 212:             if (n == int_1)
 213:                 break; /* Avoid useless last squaring */
 214:         }
 215: 
 216:         /* Square u, v */
 217: 
 218:         if (u != int_1) {
 219:             temp = u;
 220:             u = int_prod(u, u);
 221:             release((value) temp);
 222:         }
 223:         if (v != int_1) {
 224:             temp = v;
 225:             v = int_prod(v, v);
 226:             release((value) temp);
 227:         }
 228: 
 229:         n = int_half(n);
 230:     } /* while (n!=0) */
 231: 
 232:     release((value) n);
 233:     release((value) u);
 234:     release((value) v);
 235:     a = (rational) grab_rat();
 236:     Numerator(a) = tu;
 237:     Denominator(a) = tv;
 238: 
 239:     return a;
 240: }
 241: 
 242: 
 243: /* Compare two rational numbers */
 244: 
 245: Visible relation rat_comp(u, v) register rational u, v; {
 246:     int sd, su, sv;
 247:     integer nu, nv;
 248: 
 249:     /* 1. Compare pointers */
 250:     if (u == v || N(u) == N(v) && D(u) == D(v)) return 0;
 251: 
 252:     /* 2. Either zero? */
 253:     if (N(u) == int_0) return int_comp(int_0, N(v));
 254:     if (N(v) == int_0) return int_comp(N(u), int_0);
 255: 
 256:     /* 3. Compare signs */
 257:     su = Msd(N(u));
 258:     sv = Msd(N(v));
 259:     su = (su>0) - (su<0);
 260:     sv = (sv>0) - (sv<0);
 261:     if (su != sv) return su > sv ? 1 : -1;
 262: 
 263:     /* 4. Compute numerator of difference and return sign */
 264:     nu= int_prod(N(u), D(v));
 265:     nv= int_prod(N(v), D(u));
 266:     sd= int_comp(nu, nv);
 267:     release((value) nu); release((value) nv);
 268:     return sd;
 269: }
 270: 
 271: Visible Procedure rat_init() {
 272:     rat_zero = (rational) grab_rat();
 273:     Numerator(rat_zero) = int_0;
 274:     Denominator(rat_zero) = int_1;
 275: 
 276:     rat_half = (rational) grab_rat();
 277:     Numerator(rat_half) = int_1;
 278:     Denominator(rat_half) = int_2;
 279: }
 280: 
 281: Visible Procedure endrat() {
 282:     release((value) rat_zero);
 283:     release((value) rat_half);
 284: }

Defined functions

Defined macros