![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to util.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / cmd / lfactor|
Return to util.c CVS log | Up to [Research Unix] / researchv10no / cmd / lfactor |
1.1 root 1: #define testing 0
2: /*
3: * compos performs the composition of two reduced_______ forms.
4: * Thus we are entitled to assume that
5: * |b| <= a <= c
6: * and
7: * a,b < 1e16
8: */
9: #include "mp.h"
10: #include "form.h"
11: form
12: compos(f1,f2)
13: form f1, f2;
14: {
15: form f3;
16: mint m, r, y;
17: mint yq, yr;
18: mint v1, v2;
19: mint d, d1, s, n;
20: mint x1, x2, y1, y2;
21: mint temp, temp1, temp2;
22: mint zero;
23: double dinv();
24: double modf();
25:
26: zero = itom(0);
27: n.high = s.high = 0;
28: s.low = (f1.b.low + f2.b.low)*0.5;
29: n.low = s.low - f2.b.low;
30: #if testing
31: temp1 = msub(mult(f1.c,f1.a),mult(f2.c,f2.a));
32: if(mcmp(temp1,mult(n,s)) != 0){
33: printf("compos: unequal determinants!\n");
34: tobort(5);
35: }
36: #endif
37: d.high = 0;
38: x1.high = y1.high = 0;
39: d.low = dinv(f1.a.low,f2.a.low,&x1.low,&y1.low);
40: d1.high = 0;
41: x2.high = 0;
42: y2.high = 0;
43: v1.high = v2.high = 0;
44: if(d.low == 1){
45: d1.low = 1;
46: v1.low = f1.a.low;
47: v2.low = f2.a.low;
48: m = mult(y1,n);
49: }else{
50: d1.low = dinv(s.low,d.low,&x2.low,&y2.low);
51: v1.low = f1.a.low/d1.low;
52: v2.low = f2.a.low/d1.low;
53: mdiv(mult(y1,y2), v1, &temp);
54: temp1 = mult(temp, n);
55: temp2 = mult(x2, f2.c);
56: m = msub(temp1, temp2);
57: }
58: mdiv(m, v1, &r);
59: if((r.low + r.low) > v1.low)
60: r.low -= v1.low;
61: else if((r.low+r.low) < -v1.low)
62: r.low += v1.low;
63: yq = mdiv(mult(r,r), v1, &yr);
64: temp = mult(f2.c, d1);
65: temp = madd(temp, mult(r,f2.b));
66: y = madd(temp,mult(v2,yr));
67: temp = mdiv(y, v1, &temp2);
68: if(mcmp(temp2,zero) != 0) tobort(10);
69: f3.a = mult(v1, v2);
70: f3.b = madd(f2.b, smult(mult(v2,r),2));
71: f3.c = madd(temp, mult(yq,v2));
72: reduce(&f3);
73: return(f3);;
74: }
75: void
76: reduce(f)
77: form *f;
78: {
79: mint babs;
80: mint zero, one;
81: mint mtemp, k;
82: mint ka;
83: double dt1, dt2, dtemp;
84:
85: zero = itom(0);
86: one = itom(1);
87: back:
88: if(mcmp(f->a,f->c) > 0){
89: mtemp = f->a;
90: f->a = f->c;
91: f->c = mtemp;
92: f->b = mchs(f->b);
93: }
94: babs = f->b;
95: if(mcmp(f->b,zero) < 0)
96: babs = mchs(babs);
97: if(mcmp(f->a,babs) >= 0){
98: if((mcmp(f->a,f->c)==0) || (mcmp(f->a,babs)==0)){
99: f->b = babs;
100: }
101: return;
102: }
103: dt1 = f->b.high*1e16 + f->b.low;
104: dt2 = f->a.high*1e16 + f->a.low;
105: modf(dt1/dt2, &dtemp);
106: if(dtemp>0) dtemp += 1;
107: else if(dtemp<0) dtemp -= 1;
108: modf(dtemp/2, &dtemp);
109: modf(dtemp/1e16, &k.high);
110: k.low = dtemp - 1e16*k.high;
111: ka = mult(k,f->a);
112: f->b = msub(f->b,madd(ka,ka));
113: f->c = msub(f->c,mult(k,madd(ka,f->b)));
114: goto back;
115: }
116: tobort(n)
117: int n;
118: {
119: printf("Error %d\n", n);
120: abort();
121: }
122:
123: form
124: formpow(f,h)
125: form f;
126: double h;
127: {
128: form ftemp;
129: int i,n;
130: int bit[120];
131: double junk;
132: double temph;
133: double modf();
134:
135: if(h == 0.){
136: ftemp = f;
137: ftemp.b = mchs(ftemp.b);
138: ftemp = compos(f, ftemp);
139: return(ftemp);
140: }
141:
142: i = 0;
143: temph = h;
144: while(temph != 0){
145: modf(temph/2, &junk);
146: if(2.*junk != temph)
147: bit[i++] = 1;
148: else
149: bit[i++] = 0;
150: temph = junk;
151: }
152: n = i-1;
153:
154: ftemp = f;
155: for(i=n; i>0; i--){
156: ftemp = compos(ftemp, ftemp);
157: if(bit[i-1] != 0){
158: ftemp = compos(f, ftemp);
159: }
160: }
161: return(ftemp);
162: }
163: void
164: formout(f)
165: form f;
166: {
167: printf("(");
168: mout1(f.a);
169: printf(",");
170: mout1(f.b);
171: printf(",");
172: mout1(f.c);
173: printf(")");
174: return;
175: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.