![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to d-itv.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / engine|
Return to d-itv.c CVS log | Up to [local] / OpenXM_contrib2 / asir2018 / engine |
1.1 noro 1: /*
1.5 ! kondoh 2: * $OpenXM: OpenXM_contrib2/asir2018/engine/d-itv.c,v 1.4 2019/10/17 03:03:12 kondoh Exp $
1.1 noro 3: */
4: #if defined(INTERVAL)
5: #include <float.h>
6: #include "ca.h"
7: #include "base.h"
1.3 kondoh 8: #if 0
1.1 noro 9: #if defined(PARI)
10: #include "genpari.h"
11: #endif
1.3 kondoh 12: #endif
1.1 noro 13:
1.4 kondoh 14: double addulpd(double);
15: double subulpd(double);
16: extern int mpfr_roundmode;
17: Num tobf(Num,int);
18:
19:
1.1 noro 20: #if defined(ITVDEBUG)
21: void printbinint(int d)
22: {
23: int i, j, mask;
24: union {
25: int x;
26: char c[sizeof(int)];
27: } a;
28:
29: a.x = d;
30: for(i=sizeof(int)-1;i>=0;i--) {
31: mask = 0x80;
32: for(j=0;j<8;j++) {
33: if (a.c[i] & mask) fprintf(stderr,"1");
34: else fprintf(stderr,"0");
35: mask >>= 1;
36: }
37: }
38: fprintf(stderr,"\n");
39: }
40: #endif
41:
1.4 kondoh 42: #if 0
1.1 noro 43: double NatToRealUp(N a, int *expo)
44: {
45: int *p;
46: int l,all,i,j,s,tb,top,tail;
47: unsigned int t,m[2];
48: N b,c;
49: int kk, carry, rem;
50:
51: p = BD(a); l = PL(a) - 1;
52: for ( top = 0, t = p[l]; t; t >>= 1, top++ );
53: all = top + BSH*l; tail = (53-top)%BSH; j = l-(53-top)/BSH-1;
54:
55:
56: if ( j >= 0 ) {
57:
58: #if defined(ITVDEBUG)
59: fprintf(stderr," %d : tail = %d\n", j, tail);
60: printbinint(p[j]);
61: #endif
62: kk = (1<< (BSH - tail)) - 1;
63: rem = p[j] & kk;
64: if ( !rem )
65: for (i=0;i<j;i++)
66: if ( p[j] ) {
67: carry = 1;
68: break;
69: }
70: else carry = 1;
71: if ( carry ) {
72: b = NALLOC(j+1+1);
73: PL(b) = j+1;
74: p = BD(b);
75: for(i=0;i<j;i++) p[i] = 0;
76: p[j]=(1<< (BSH - tail));
77:
78: addn(a,b,&c);
79:
80: p = BD(c); l = PL(c) - 1;
81: for ( top = 0, t = p[l]; t; t >>= 1, top++ );
82: all = top + BSH*l; tail = (53-top)%BSH; i = l-(53-top)/BSH-1;
83: } else i = j;
84: } else i = j;
85:
86:
87: m[1] = i < 0 ? 0 : p[i]>>(BSH-tail);
88: for ( j = 1, i++, tb = tail; i <= l; i++ ) {
89: s = 32-tb; t = i < 0 ? 0 : p[i];
90: if ( BSH > s ) {
91: m[j] |= ((t&((1<<s)-1))<<tb);
92: if ( !j )
93: break;
94: else {
95: j--; m[j] = t>>s; tb = BSH-s;
96: }
97: } else {
98: m[j] |= (t<<tb); tb += BSH;
99: }
100: }
101: s = (all-1)+1023;
102: m[0] = (m[0]&((1<<20)-1))|(MIN(2046,s)<<20); *expo = MAX(s-2046,0);
103: #ifdef vax
104: t = m[0]; m[0] = m[1]; m[1] = t; itod(m);
105: #endif
106: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) || defined(ANDROID)
107: t = m[0]; m[0] = m[1]; m[1] = t;
108: #endif
109: return *((double *)m);
110: }
111:
112: static double Q2doubleDown(Q a)
113: {
114: double nm,dn,t,snm,rd;
115: int enm,edn,e;
116: unsigned int *p,s;
117:
118: nm = (SGN(a)>0)? NatToReal(NM(a),&enm): NatToRealUp(NM(a),&enm);
119:
120: if ( INT(a) )
121: if ( enm >= 1 )
122: error("Q2doubleDown : Overflow");
123: else
124: return SGN(a)>0 ? nm : -nm;
125: else {
126: dn = (SGN(a)>0)? NatToReal(DN(a),&edn): NatToRealUp(DN(a),&edn);
127:
128: if ( ((e = enm - edn) >= 1024) || (e <= -1023) )
129: error("Q2doubleDown : Overflow"); /* XXX */
130: else {
131: t = 0.0; p = (unsigned int *)&t; /* Machine */
132: *p = ((e+1023)<<20);
133: #ifdef vax
134: s = p[0]; p[0] = p[1]; p[1] = s; itod(p);
135: #endif
136: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) | defined(ANDROID)
137: s = p[0]; p[0] = p[1]; p[1] = s;
138: #endif
139: FPMINUSINF
140: snm = (SGN(a)>0) ? nm : -nm;
141: rd = snm / dn * t;
142: FPNEAREST
143: return rd;
144: }
145: }
146: }
147:
148:
149: static double Q2doubleUp(Q a)
150: {
151: double nm,dn,t,snm,rd;
152: int enm,edn,e;
153: unsigned int *p,s;
154:
155: nm = (SGN(a)>0)? NatToRealUp(NM(a),&enm): NatToReal(NM(a),&enm);
156:
157: if ( INT(a) )
158: if ( enm >= 1 )
159: error("Q2doubleUp : Overflow");
160: else
161: return SGN(a)>0 ? nm : -nm;
162: else {
163: dn = (SGN(a)>0)? NatToRealUp(DN(a),&edn): NatToReal(DN(a),&edn);
164:
165: if ( ((e = enm - edn) >= 1024) || (e <= -1023) )
166: error("Q2doubleUp : Overflow"); /* XXX */
167: else {
168: t = 0.0; p = (unsigned int *)&t; /* Machine */
169: *p = ((e+1023)<<20);
170: #ifdef vax
171: s = p[0]; p[0] = p[1]; p[1] = s; itod(p);
172: #endif
173: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) || defined(ANDROID)
174: s = p[0]; p[0] = p[1]; p[1] = s;
175: #endif
176: #if 0
177: FPPLUSINF
178: snm = (SGN(a)>0) ? nm : -nm;
179: rd = snm / dn * t;
180: #endif
181:
182: FPMINUSINF
183: snm = (SGN(a)>0) ? -nm : nm;
184: rd = snm / dn * t;
185: FPNEAREST
186: return (-rd);
187: }
188: }
189: }
190:
191: static double PARI2doubleDown(BF a)
192: {
1.3 kondoh 193: //GEN p;
194: Num p;
1.1 noro 195: double d;
196:
197: ritopa(a, &p);
198: d = rtodbl(p);
199: cgiv(p);
200: return d;
201: }
202:
203: static double PARI2doubleUp(BF a)
204: {
205: return PARI2doubleDown(a);
206: }
1.3 kondoh 207: #endif
1.1 noro 208:
1.4 kondoh 209: static double Q2doubleDown(Q a)
210: {
211: Num b;
212: double rd;
213: int current_roundmode=0;
214:
215: current_roundmode = mpfr_roundmode;
216: mpfr_roundmode = MPFR_RNDD;
217: b = tobf((Num)a,60);
218: mpfr_roundmode = current_roundmode;
219: rd = mpfr_get_d(BDY((BF)b),MPFR_RNDD);
220: return rd;
221: }
222:
223: static double Q2doubleUp(Q a)
224: {
225: Num b;
226: double rd;
227: int current_roundmode=0;
228:
229: current_roundmode = mpfr_roundmode;
230: mpfr_roundmode = MPFR_RNDU;
231: b = tobf((Num)a,60);
232: mpfr_roundmode = current_roundmode;
233: rd = mpfr_get_d(BDY((BF)b),MPFR_RNDU);
234: return rd;
235: }
236:
1.1 noro 237: double subulpd(double d)
238: {
239: double inf;
240:
241: FPMINUSINF
242: inf = d - DBL_MIN;
243: FPNEAREST
244: return inf;
245: }
246:
247: double addulpd(double d)
248: {
249: double md, sup;
250:
251: md = -d;
252: FPMINUSINF
253: sup = md - DBL_MIN;
254: FPNEAREST
255: return (-sup);
256: }
257:
1.5 ! kondoh 258: double toRealDown(Num a)
1.1 noro 259: {
260: double inf;
261:
262: if ( ! a ) return 0.0;
263: switch ( NID(a) ) {
264: case N_Q:
265: inf = Q2doubleDown((Q)a); break;
266: case N_R:
267: inf = subulpd(BDY((Real)a)); break;
268: case N_B:
1.4 kondoh 269: inf = mpfr_get_d(BDY((BF)a),MPFR_RNDD);
270: break;
1.1 noro 271: case N_IP:
1.5 ! kondoh 272: inf = toRealDown(INF((Itv)a));
1.1 noro 273: break;
274: case N_IntervalDouble:
275: inf = INF((IntervalDouble)a); break;
1.5 ! kondoh 276: case N_IntervalBigFloat:
! 277: inf = mpfr_get_d(BDY((BF)INF((IntervalBigFloat)a)),MPFR_RNDD);
! 278: break;
1.1 noro 279: case N_A:
280: default:
281: inf = 0.0;
1.5 ! kondoh 282: error("toRealDown: not supported operands.");
1.1 noro 283: break;
284: }
285: return inf;
286: }
287:
1.5 ! kondoh 288: double toRealUp(Num a)
1.1 noro 289: {
290: double sup;
291:
292: if ( ! a ) return 0.0;
293: switch ( NID(a) ) {
294: case N_Q:
295: sup = Q2doubleUp((Q)a); break;
296: case N_R:
297: sup = addulpd(BDY((Real)a)); break;
298: case N_B:
1.4 kondoh 299: sup = mpfr_get_d(BDY((BF)a),MPFR_RNDU);
300: break;
1.1 noro 301: case N_IP:
1.5 ! kondoh 302: sup = toRealUp(SUP((Itv)a)); break;
1.1 noro 303: case N_IntervalDouble:
304: sup = SUP((IntervalDouble)a); break;
1.5 ! kondoh 305: case N_IntervalBigFloat:
! 306: sup = mpfr_get_d(BDY((BF)INF((IntervalBigFloat)a)),MPFR_RNDU);
! 307: break;
1.1 noro 308: case N_A:
309: default:
310: sup = 0.0;
1.5 ! kondoh 311: error("toRealUp: not supported operands.");
1.1 noro 312: break;
313: }
314: return sup;
315: }
316:
317:
318: void Num2double(Num a, double *inf, double *sup)
319: {
1.5 ! kondoh 320: *inf = 0.0;
! 321: *sup = 0.0;
! 322: if (a && NUM(a) )
1.1 noro 323: switch ( NID(a) ) {
324: case N_Q:
325: *inf = Q2doubleDown((Q)a);
326: *sup = Q2doubleUp((Q)a);
327: break;
328: case N_R:
329: *inf = BDY((Real)a);
330: *sup = BDY((Real)a);
331: break;
332: case N_B:
1.3 kondoh 333: *inf = mpfr_get_d(BDY((BF)a), MPFR_RNDD);
334: *sup = mpfr_get_d(BDY((BF)a), MPFR_RNDU);
1.1 noro 335: break;
336: case N_IP:
1.5 ! kondoh 337: *inf = toRealDown(INF((Itv)a));
! 338: *sup = toRealUp(SUP((Itv)a));
1.1 noro 339: break;
340: case N_IntervalDouble:
341: *inf = INF((IntervalDouble)a);
342: *sup = SUP((IntervalDouble)a);
343: break;
344: case N_A:
345: default:
346: error("Num2double: not supported operands.");
347: break;
348: }
349: }
350:
351:
352: void additvd(Num a, Num b, IntervalDouble *c)
353: {
354: double ai,as,mas, bi,bs;
355: double inf,sup;
356:
357: if ( !a ) {
358: *c = (IntervalDouble)b;
359: } else if ( !b ) {
360: *c = (IntervalDouble)a;
361: #if 0
362: } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
363: addnum(0,a,b,c);
364: #endif
365: } else {
366:
367: Num2double((Num)a,&ai,&as);
368: Num2double((Num)b,&bi,&bs);
369: mas = -as;
370: FPMINUSINF
371: inf = ai + bi;
372: sup = mas - bs; /* as + bs = ( - ( (-as) - bs ) ) */
373: FPNEAREST
374: MKIntervalDouble(inf,(-sup),*c);
375: }
376: }
377:
378: void subitvd(Num a, Num b, IntervalDouble *c)
379: {
380: double ai,as,mas, bi,bs;
381: double inf,sup;
382:
383: if ( !a ) {
384: *c = (IntervalDouble)b;
385: } else if ( !b ) {
386: *c = (IntervalDouble)a;
387: #if 0
388: } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
389: addnum(0,a,b,c);
390: #endif
391: } else {
392: Num2double(a,&ai,&as);
393: Num2double(b,&bi,&bs);
394: FPMINUSINF
395: inf = ai - bs;
396: sup = ( bi - as ); /* as - bi = ( - ( bi - as ) ) */
397: FPNEAREST
398: MKIntervalDouble(inf,(-sup),*c);
399: }
400: }
401:
402: void mulitvd(Num a, Num b, IntervalDouble *c)
403: {
404: double ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p;
405: double inf, sup;
406: int ba,bb;
407:
408: if ( !a || !b ) {
409: *c = 0;
410: #if 0
411: } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
412: mulnum(0,a,b,c);
413: #endif
414: } else {
415: Num2double(a,&ai,&as);
416: Num2double(b,&bi,&bs);
417:
418: FPMINUSINF
419:
420: a2 = -as;
421: b2 = -bs;
422:
423: if ( ba = ( a2 < 0.0 ) ) {
424: a1 = ai;
425: } else {
426: a1 = a2;
427: a2 = ai;
428: }
429: if ( bb = ( b2 < 0.0 ) ) {
430: b1 = bi;
431: } else {
432: b1 = b2;
433: b2 = bi;
434: }
435:
436: c2 = - a2 * b2;
437: if ( b1 < 0.0 ) {
438: c1 = - a2 * b1;
439: if ( a1 < 0.0 ) {
440: p = - a1 * b2;
441: if ( p < c1 ) c1 = p;
442: p = - a1 * b1;
443: if ( p < c2 ) c2 = p;
444: }
445: } else {
446: c1 = (a1<0.0)?( - a1 * b2 ):( a1 * b1 );
447: }
448: if ( ba == bb ) {
449: inf = c1;
450: sup = - c2;
451: } else {
452: inf = c2;
453: sup = - c1;
454: }
455: FPNEAREST
456: MKIntervalDouble(inf,sup,*c);
457: }
458: }
459:
460: int initvd(Num a, IntervalDouble b)
461: {
462: Real inf, sup;
463:
464: if ( NID(b) == N_IntervalDouble ) {
465: MKReal(INF(b), inf);
466: MKReal(SUP(b), sup);
467: } else return 0;
468: if ( compnum(0,(Num)inf,a) > 0 ) return 0;
469: else if ( compnum(0,a,(Num)sup) > 0 ) return 0;
470: else return 1;
471: }
472:
473: void divitvd(Num a, Num b, IntervalDouble *c)
474: {
475: double ai,as,bi,bs,a1,a2,b1,b2,c1,c2;
476: double inf, sup;
477: int ba,bb;
478:
479: if ( !b ) {
480: *c = 0;
481: error("divitvd : division by 0");
482: } else if ( !a ) {
483: *c = 0;
484: #if 0
485: } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
486: divnum(0,a,b,c);
487: #endif
488: } else {
489: Num2double(a,&ai,&as);
490: Num2double(b,&bi,&bs);
491:
492: FPMINUSINF
493:
494: a2 = -as;
495: b2 = -bs;
496:
497: if ( ba = ( a2 < 0.0 ) ) {
498: a1 = ai;
499: } else {
500: a1 = a2;
501: a2 = ai;
502: }
503: if ( bb = ( b2 < 0.0 ) ) {
504: b1 = bi;
505: } else {
506: b1 = b2;
507: b2 = bi;
508: }
509:
510: if ( b1 <= 0.0 ) {
511: *c = 0;
512: error("divitvd : division by 0 interval");
513: } else {
514: c2 = a2 / b1;
515: c1 = (a1<0.0)?( a1 / b1 ):( - a1 / b2 );
516: }
517: if ( ba == bb ) {
518: inf = c1;
519: sup = - c2;
520: } else {
521: inf = c2;
522: sup = - c1;
523: }
524: FPNEAREST
525: MKIntervalDouble(inf,sup,*c);
526: }
527: }
528:
529: void pwritvd(Num a, Num e, IntervalDouble *c)
530: {
1.4 kondoh 531: long ei;
1.1 noro 532: IntervalDouble t;
533:
534: if ( !e ) {
535: *c = (IntervalDouble)ONE;
536: } else if ( !a ) {
537: *c = 0;
538: #if 0
539: } else if ( NID(a) <= N_IP ) {
540: pwrnum(0,a,e,c);
541: #endif
542: } else if ( !INT(e) ) {
543: #if defined(PARI) && 0
544: gpui_ri((Obj)a,(Obj)c,(Obj *)c);
545: #else
546: error("pwritvd : can't calculate a fractional power");
547: #endif
548: } else {
1.4 kondoh 549: //ei = QTOS((Q)e);
550: ei = mpz_get_si(BDY((Z)e));
1.3 kondoh 551: if (ei<0) ei = -ei;
1.1 noro 552: pwritv0d((IntervalDouble)a,ei,&t);
1.4 kondoh 553: // if ( SGN((Q)e) < 0 )
554: if ( sgnq((Q)e) < 0 )
1.1 noro 555: divnum(0,(Num)ONE,(Num)t,(Num *)c);
556: else
557: *c = t;
558: }
559: }
560:
1.4 kondoh 561: void pwritv0d(IntervalDouble a, long e, IntervalDouble *c)
1.1 noro 562: {
563: double inf, sup;
564: double t, Xmin, Xmax;
565: int i;
566:
567: if ( e == 1 )
568: *c = a;
569: else {
570: if ( !(e%2) ) {
571: if ( initvd(0,a) ) {
572: Xmin = 0.0;
573: t = - INF(a);
574: Xmax = SUP(a);
575: if ( t > Xmax ) Xmax = t;
576: } else {
577: if ( INF(a) > 0.0 ) {
578: Xmin = INF(a);
579: Xmax = SUP(a);
580: } else {
581: Xmin = - SUP(a);
582: Xmax = - INF(a);
583: }
584: }
585: } else {
586: Xmin = INF(a);
587: Xmax = SUP(a);
588: }
589: FPPLUSINF
590: sup = (Xmax!=0.0)?pwrreal0(Xmax,e):0.0;
591: FPMINUSINF
592: inf = (Xmin!=0.0)?pwrreal0(Xmin,e):0.0;
593: FPNEAREST
594: MKIntervalDouble(inf,sup,*c);
595: }
596: }
597:
598: void chsgnitvd(IntervalDouble a, IntervalDouble *c)
599: {
600: double inf,sup;
601:
602: if ( !a )
603: *c = 0;
604: #if 0
605: else if ( NID(a) <= N_IP )
606: chsgnnum(a,c);
607: #endif
608: else {
609: inf = - SUP(a);
610: sup = - INF(a);
611: MKIntervalDouble(inf,sup,*c);
612: }
613: }
614:
615: int cmpitvd(IntervalDouble a, IntervalDouble b)
616: /* a > b: 1 */
617: /* a = b: 0 */
618: /* a < b: -1 */
619: {
620: double ai,as,bi,bs;
621:
622: if ( !a ) {
623: if ( !b || (NID(b)<=N_IP) ) {
624: return compnum(0,(Num)a,(Num)b);
625: } else {
626: bi = INF(b);
627: bs = SUP(b);
628: if ( bi > 0.0 ) return -1;
629: else if ( bs < 0.0 ) return 1;
630: else return 0;
631: }
632: } else if ( !b ) {
633: if ( !a || (NID(a)<=N_IP) ) {
634: return compnum(0,(Num)a,(Num)b);
635: } else {
636: ai = INF(a);
637: as = SUP(a);
638: if ( ai > 0.0 ) return 1;
639: else if ( as < 0.0 ) return -1;
640: else return 0;
641: }
642: } else {
643: bi = INF(b);
644: bs = SUP(b);
645: ai = INF(a);
646: as = SUP(a);
647: if ( ai > bs ) return 1;
648: else if ( bi > as ) return -1;
649: else return 0;
650: }
651: }
652:
653: void miditvd(IntervalDouble a, Num *b)
654: {
655: double t;
656: Real rp;
657:
658: if ( ! a ) *b = 0;
659: else if ( (NID(a) <= N_IP) )
660: *b = (Num)a;
661: else {
662: t = (INF(a) + SUP(a))/2.0;
663: MKReal(t,rp);
664: *b = (Num)rp;
665: }
666: }
667:
668: void cupitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c)
669: {
670: double ai,as,bi,bs;
671: double inf,sup;
672:
673: ai = INF(a);
674: as = SUP(a);
675: bi = INF(b);
676: bs = SUP(b);
677: inf = MIN(ai,bi);
678: sup = MAX(as,bs);
679: MKIntervalDouble(inf,sup,*c);
680: }
681:
682: void capitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c)
683: {
684: double ai,as,bi,bs;
685: double inf,sup;
686:
687: ai = INF(a);
688: as = SUP(a);
689: bi = INF(b);
690: bs = SUP(b);
691: inf = MAX(ai,bi);
692: sup = MIN(as,bs);
693: if ( inf > sup ) *c = 0;
694: else MKIntervalDouble(inf,sup,*c);
695: }
696:
697:
698: void widthitvd(IntervalDouble a, Num *b)
699: {
700: double t;
701: Real rp;
702:
703: if ( ! a ) *b = 0;
704: else {
705: t = SUP(a) - INF(a);
706: MKReal(t,rp);
707: *b = (Num)rp;
708: }
709: }
710:
711: void absitvd(IntervalDouble a, Num *b)
712: {
713: double ai,as,bi,bs;
714: double t;
715: Real rp;
716:
717: if ( ! a ) *b = 0;
718: else {
719: ai = INF(a);
720: as = SUP(a);
721: if (ai < 0) bi = -ai; else bi = ai;
722: if (as < 0) bs = -as; else bs = as;
723: if ( bi > bs ) t = bi; else t = bs;
724: MKReal(t,rp);
725: *b = (Num)rp;
726: }
727: }
1.5 ! kondoh 728: void absintvald(IntervalDouble a, IntervalDouble *b)
! 729: {
! 730: double ai,as,inf,sup;
! 731:
! 732: ai = INF(a);
! 733: as = SUP(a);
! 734: if ( as < 0 ) {
! 735: sup = -ai;
! 736: inf = -as;
! 737: } else if (ai < 0) {
! 738: inf = 0.0;
! 739: sup = MAX(as, -ai);
! 740: } else {
! 741: inf = ai;
! 742: sup = as;
! 743: }
! 744: MKIntervalDouble(inf,sup,*b);
! 745: }
1.1 noro 746:
747: void distanceitvd(IntervalDouble a, IntervalDouble b, Num *c)
748: {
749: double ai,as,bi,bs;
750: double di,ds;
751: double t;
752: Real rp;
753:
754: ai = INF(a);
755: as = SUP(a);
756: bi = INF(b);
757: bs = SUP(b);
758: di = bi - ai;
759: ds = bs - as;
760:
761: if (di < 0) di = -di;
762: if (ds < 0) ds = -ds;
763: if (di > ds) t = di; else t = ds;
764: MKReal(t,rp);
765: *c = (Num)rp;
766: }
767:
768: #endif