![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to B1fun.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSD / contrib / B / src / bsmall|
Return to B1fun.c CVS log | Up to [CSRG BSD Unix] / 43BSD / contrib / B / src / bsmall |
1.1 root 1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1984. */
2: /* $Header: B1fun.c,v 1.1 84/06/28 00:48:54 timo Exp $ */
3:
4: /* B functions */
5: #include "b.h"
6: #include "b1obj.h"
7: #include "b2sem.h"
8: #include "B1num.h"
9:
10: #define Maxlen(x, y) (Length(x) > Length(y) ? Length(x) : Length(y))
11: #define Sumlen(x, y) (Length(x) + Length(y))
12:
13: value sum(x, y) register value x, y; {
14: value r, z;
15: Checknum(x); Checknum(y);
16: if (!Exact(x) || !Exact(y)) return mk_approx(Numval(x) + Numval(y));
17: z= mk_exact(Denominator(x), Denominator(y), 0);
18: r= mk_exact(Numerator(x)*Denominator(z)+Numerator(y)*Numerator(z),
19: Denominator(x)*Denominator(z), Maxlen(x, y));
20: release(z);
21: return r;
22: }
23:
24: value negated(x) register value x; {
25: Checknum(x);
26: if (!Exact(x)) return mk_approx(-Numval(x));
27: return mk_exact(-Numerator(x), Denominator(x), Length(x));
28: }
29:
30: value diff(x, y) register value x, y; {
31: value r, my;
32: r= sum(x, my= negated(y));
33: release(my);
34: return r;
35: }
36:
37: value inverse(x) value x;{
38: Checknum(x);
39: if (Numval(x) == 0) error("in x/y, y is zero");
40: if (!Exact(x)) return mk_approx(1.0/Numval(x));
41: return mk_exact(Denominator(x), Numerator(x), Length(x));
42: }
43:
44: value prod(x, y) register value x, y; {
45: value a, b, r;
46: Checknum(x); Checknum(y);
47: if (!Exact(x) || !Exact(y)) return mk_approx(Numval(x) * Numval(y));
48: a= mk_exact(Numerator(x), Denominator(y), 0);
49: b= mk_exact(Numerator(y), Denominator(x), 0);
50: r= mk_exact(Numerator(a)*Numerator(b), Denominator(a)*Denominator(b),
51: Sumlen(x, y));
52: release(a); release(b);
53: return r;
54: }
55:
56: value quot(x, y) register value x, y; {
57: value r, iy;
58: r= prod(x, iy= inverse(y));
59: release(iy);
60: return r;
61: }
62:
63: #define Even(x) ((x) == Two*floor((x)/Two))
64: value power(x, y) register value x, y; {
65: Checknum(x); Checknum(y);
66: if (Exact(y)) {
67: integer py= Numerator(y), qy= Denominator(y);
68: if (Integral(y) && Exact(x)) {
69: integer px, qx, ppx, pqx, Ppx, Pqx;
70: if (py == Zero) return mk_int(One);
71: if (py > Zero) {
72: px= Numerator(x);
73: qx= Denominator(x);
74: } else {
75: py= -py;
76: px= Denominator(x);
77: qx= Numerator(x);
78: }
79: ppx= pqx= One;
80: Ppx= px; Pqx= qx;
81: while (py >= Two) {
82: if (!Even(py)) {
83: ppx*= Ppx; pqx*= Pqx;
84: }
85: Ppx*= Ppx; Pqx*= Pqx;
86: py= floor(py/Two);
87: }
88: ppx*= Ppx; pqx*= Pqx;
89: return mk_exact(ppx, pqx, 0);
90: } /* END Integral(y) && Exact(x) */
91: else {
92: double vx= Numval(x);
93: short sx= vx < 0 ? -1 : vx == 0 ? 0 : 1;
94: if (sx < 0 && Even(qy))
95: error("in x**(p/q), x is negative and q is even");
96: if (sx == 0 && py < Zero)
97: error("0**y with negative y");
98: if (sx < 0 && Even(py)) sx= 1;
99: return mk_approx(sx * pow(fabs(vx), py/qy));
100: }
101: } /* END Exact(y) */
102: else {
103: double vx= Numval(x), vy= Approxval(y);
104: if (vy == 0) return mk_approx(1.0);
105: if (vx < 0)
106: error("in x**y, x is negative and y is not exact");
107: if (vx == 0 && vy < 0)
108: error("0E0**y with negative y");
109: return mk_approx(pow(vx, vy));
110: }
111: }
112:
113: value root2(n, x) register value n, x; {
114: value r, in;
115: Checknum(n);
116: if (Numval(n) == 0) error("in x root y, x is zero");
117: r= power(x, in= inverse(n));
118: release(in);
119: return r;
120: }
121:
122: value absval(x) register value x; {
123: Checknum(x);
124: if (!Exact(x)) return mk_approx(fabs(Numval(x)));
125: return mk_exact((integer) fabs((double) Numerator(x)), Denominator(x), Length(x));
126: }
127:
128: value signum(x) register value x; {
129: double v= numval(x);
130: return mk_int(v < 0 ? -One : v == 0 ? Zero : One);
131: }
132:
133: value floorf(x) register value x; {
134: return mk_int(floor(numval(x)));
135: }
136:
137: value ceilf(x) register value x; {
138: return mk_int(ceil(numval(x)));
139: }
140:
141: value round1(x) register value x; {
142: return mk_int(floor(numval(x) + .5));
143: }
144:
145: value round2(n, x) register value n, x; {
146: value ten, tenp, xtenp, r0, r;
147: Checknum(n);
148: if (!Integral(n)) error("in n round x, n is not an integer");
149: ten= mk_integer(10);
150: tenp= power(ten, n);
151: xtenp= prod(x, tenp);
152: r0= round1(xtenp);
153: r= mk_exact(Numerator(r0), Numerator(tenp), propintlet((int) Numerator(n)));
154: release(ten); release(tenp); release(xtenp); release(r0);
155: return r;
156: }
157:
158: value mod(a, n) register value a, n; {
159: value f, p, d;
160: Checknum(a); Checknum(n);
161: f= mk_int(floor(Numval(a) / Numval(n)));
162: p= prod(n, f);
163: d= diff(a, p);
164: release(f); release(p);
165: return d;
166: }
167:
168: double lastran;
169:
170: setran (seed) double seed;
171: {double x;
172: x= seed >= 0 ? seed : -seed;
173: while (x >= 1) x/= 10;
174: lastran= floor(67108864.0 * x);
175: }
176:
177: set_random(v) value v; {
178: setran((double) hash(v));
179: }
180:
181: value random() /* 0 <= r < 1 */
182: {double p;
183: p= 26353589.0 * lastran + 1;
184: lastran= p - 67108864.0 * floor (p / 67108864.0);
185: return mk_approx(lastran / 67108864.0);
186: }
187:
188: value root1(v) value v; {
189: value two= mk_integer(2);
190: v= root2(two, v);
191: release(two);
192: return(v);
193: }
194:
195: value pi() { return mk_approx(3.141592653589793238462); }
196: value e() { return mk_approx(exp(1.0)); }
197:
198: value sin1(v) value v; { return mk_approx(sin(numval(v))); }
199: value cos1(v) value v; { return mk_approx(cos(numval(v))); }
200: value tan1(v) value v; { return mk_approx(tan(numval(v))); }
201: value atn1(v) value v; { return mk_approx(atan(numval(v))); }
202: value exp1(v) value v; { return mk_approx(exp(numval(v))); }
203: value log1(v) value v; { return mk_approx(log(numval(v))); }
204:
205: value log2(u, v) value u, v;{
206: return mk_approx(log(numval(v)) / log(numval(u)));
207: }
208:
209: value atn2(u, v) value u, v; {
210: return mk_approx(atan2(numval(v), numval(u)));
211: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.