1/* BEGIN LICENSE BLOCK 2 * Version: CMPL 1.1 3 * 4 * The contents of this file are subject to the Cisco-style Mozilla Public 5 * License Version 1.1 (the "License"); you may not use this file except 6 * in compliance with the License. You may obtain a copy of the License 7 * at www.eclipse-clp.org/license. 8 * 9 * Software distributed under the License is distributed on an "AS IS" 10 * basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See 11 * the License for the specific language governing rights and limitations 12 * under the License. 13 * 14 * The Original Code is The ECLiPSe Constraint Logic Programming System. 15 * The Initial Developer of the Original Code is Cisco Systems, Inc. 16 * Portions created by the Initial Developer are 17 * Copyright (C) 1996-2006 Cisco Systems, Inc. All Rights Reserved. 18 * 19 * Contributor(s): 20 * 21 * END LICENSE BLOCK */ 22 23/*---------------------------------------------------------------------- 24* System: ECLiPSe Constraint Logic Programming System 25* Version: $Id: intervals.c,v 1.13 2015/04/02 14:35:46 jschimpf Exp $ 26* 27 28Supported operations: 29 30+ E unary plus 31- E unary minus 32abs(E) absolute value 33floor(E) round down to integral value 34ceiling(E) round up to integral value 35truncate(E) round towards zero to integral value 36round(E) round to nearest integral value 37sgn(E) sign value (exception if zero-crossing) 38 39E1 + E2 addition 40E1 - E2 subtraction 41E1 * E2 multiplication 42E1 / E2 division 43E1 ^ E2 power operation 44min(E1,E2) minimum of 2 values 45max(E1,E2) maximum of 2 values 46 47sin(E) trigonometric function 48cos(E) trigonometric function 49tan(E) trigonometric function 50asin(E) trigonometric function 51acos(E) trigonometric function 52atan(E) trigonometric function 53atan(Y,X) trigonometric function 54exp(E) exponential function e^x 55ln(E) natural logarithm 56sqrt(E) square root 57 58Conversions: 59 60integer(E) convert to integer, value must be exact 61float(E) return "middle" value (logsplit) 62breal(E) convert to interval 63 rounds out floats 64breal_from_bounds(L, U) create from bounds 65breal_min(E) return lower bound of an interval 66breal_max(E) return upper bound of an interval 67 68breal_bounds(I, L, U) 69 returns upper and lower bounds of an interval I 70 71 Extra predicates for intervals: 72sqr(E) 73+-(E,Res) not a function, can be done with -/1 and union 74rsqr(E,Res) +- sqrt, not a function 75pow(E,Int,Res) 76pow(E,Int,Res) Z=pow(X,Y), reverse of X=pow(Z,N), N odd 77rpow_even(X,Y,Z,Res) Reverse of Z^N=X, N even. 78 Computes Res = intersect( X^Y, Z), Y=1/N 79 80accuracy(E,-Width:double) 81interval_union(X,Y,Z) 82linsplit(X,MinWidth,Ratio,-Split:double) 83logsplit(X,MinWidth,Ratio,-Split:double) 84mindeltas(+Xold, +Xnew, -LwdDelta:double, -UpbDelta:double) 85 86*----------------------------------------------------------------------*/ 87 88#include "config.h" 89#include "sepia.h" 90#include "types.h" 91#include "embed.h" 92#include "error.h" 93#include "mem.h" 94#include "dict.h" 95#include "emu_export.h" 96#include "opcode.h" 97#include "intervals.h" 98#include "rounding_control.h" 99 100/* 101 * Contents: Floating-Point Interval Constraint Solver 102 * (C part: primitives for interval arithmetic) 103 * 104 * Author: Joachim Schimpf, IC-Parc 105 * 106 * 107 * The following code relies heavily on IEEE-conformant floating-point 108 * arithmetic for rounding, infinity handling and undefined values (NaN). 109 * 110 * General hints about programming with IEEE-floats: 111 * Great care has to be taken when intermediate results are undefined, 112 * eg. after computations involving infinities (inf-inf etc). The result 113 * is then an IEEE-NaN. NaNs propagate through further computations 114 * and comparisons involving them do always fail. The latter means that 115 * (A>=B) and (!(A<B)) are not semantically equivalent!!! 116 * Note that we never return NaNs to the Prolog level (well, at least 117 * not intentionally...). 118 * Negative Zeros: A test like x < 0.0 does not succeed for -0.0! 119 * : The Test x == 0.0 will succeed for both 0.0 AND -0.0. 120 * : In order to assign a value of -0.0, one must use the 121 * expression -1*(0.0) as MSVC incorrectly handles -0.0 122 * in source code. 123 * 124 * About safe rounding: 125 * We use two approaches to ensure safe rounding. For basic operations 126 * (+, -, *, /) where we are sure it will work correctly, we use processor 127 * rounding modes to guarantee the safety of the results (see 128 * rounding_control.h). In other cases, where it's less clear how the 129 * rounding modes will affect the operations (sin, cos and others), we 130 * always round the results one step outwards. 131 */ 132 133 134/* Declarations of global variables used by the macros in rounding_control.h */ 135 136Declare_Rounding_Control_State 137; 138 139#define Return_Double(v, t, d) { \ 140 value dval; \ 141 Make_Double_Val(dval, d); \ 142 Bind_Var(v, t, dval.all, TDBL); \ 143} 144 145#ifdef IEEE_INEXACT 146int result_inexact; 147#endif 148 149#define min(x,y) ((x) < (y) ? (x) : (y)) 150#define max(x,y) ((x) > (y) ? (x) : (y)) 151 152#define NEGZERO (-1*(0.0)) 153 154#define IsExactIntIvl(pw,_i) \ 155 (IvlLwb(pw) == IvlUpb(pw) && (modf(IvlLwb(pw), _i) == 0.0)) 156 157 158/* ---------------------------------------------------------------------- 159 * Forward declarations 160 * ---------------------------------------------------------------------- */ 161static int 162linsplit(double minwidth, double ratio, double lwb, double upb, int upper, double *split); 163static int 164logsplit(double minwidth, double ratio, double orig_lwb, double orig_upb, int upper, double *split); 165 166static dident d_undecidable; 167 168/* ---------------------------------------------------------------------- 169 * Auxiliary functions 170 * ---------------------------------------------------------------------- */ 171 172static int 173samesign(double x, double y) 174{ 175 ieee_double dx, dy; 176 dx.as_dbl = x; 177 dy.as_dbl = y; 178#if (SIZEOF_WORD == 8) 179 return (dx.as_int & SIGN_BIT) == (dy.as_int & SIGN_BIT); 180#else /* assume SIZEOF_WORD == 4 */ 181 return (dx.as_struct.mant1 & SIGN_BIT) == (dy.as_struct.mant1 & SIGN_BIT); 182#endif 183} 184 185#if (SIZEOF_WORD == 8) 186 187double 188ec_ieee_up(double x) 189{ 190 ieee_double res; 191 res.as_dbl = x; 192 if ((res.as_int & ~SIGN_BIT) == 0) /* x == +-0 */ 193 { 194 return MINDOUBLE; 195 } 196 else if ((word) res.as_int >= 0) /* x > 0 */ 197 { 198 if (res.as_int < UWSUF(0x7ff0000000000000)) /* x < Inf */ 199 { 200 res.as_int += 1; 201 return res.as_dbl; 202 } 203 } 204 else if (res.as_int <= UWSUF(0xfff0000000000000)) /* -Inf <= x < 0 */ 205 { 206 res.as_int -= 1; 207 return res.as_dbl; 208 } 209 return x; /* NaN */ 210} 211 212 213double 214ec_ieee_down(double x) 215{ 216 ieee_double res; 217 res.as_dbl = x; 218 if ((res.as_int & ~SIGN_BIT) == 0) /* x == +-0 */ 219 { 220 return -MINDOUBLE; 221 } 222 else if ((word) res.as_int > 0) /* x > 0 */ 223 { 224 if (res.as_int <= UWSUF(0x7ff0000000000000)) /* x <= Inf */ 225 { 226 res.as_int -= 1; 227 return res.as_dbl; 228 } 229 } 230 else if (res.as_int < UWSUF(0xfff0000000000000)) /* -Inf < x < 0 */ 231 { 232 res.as_int += 1; 233 return res.as_dbl; 234 } 235 return x; /* NaN */ 236} 237 238#else /* assume SIZEOF_WORD == 4 */ 239 240#define Inc64(h,l) { \ 241 if ((l += 1) == 0) \ 242 h += 1; \ 243 } 244 245#define Dec64(h,l) { \ 246 if (l == 0) { \ 247 l -= 1; \ 248 h -= 1; \ 249 } else { \ 250 l -= 1; \ 251 } \ 252 } 253 254double 255ec_ieee_up(double x) 256{ 257 ieee_double res; 258 unsigned m1, m0; 259 res.as_dbl = x; 260 m1 = res.as_struct.mant1; 261 m0 = res.as_struct.mant0; 262 if ((m1 & ~SIGN_BIT | m0) == 0) /* x == +-0 */ 263 { 264 return MINDOUBLE; 265 } 266 else if (((word) m1) >= 0) /* x > 0 */ 267 { 268 if (m1 < 0x7ff00000) /* x < Inf */ 269 { 270 Inc64(m1,m0); 271 res.as_struct.mant1 = m1; 272 res.as_struct.mant0 = m0; 273 return res.as_dbl; 274 } 275 } 276 else if (m1 < 0xfff00000 277 || m1 == 0xfff00000 && m0 == 0) /* -Inf <= x < 0 */ 278 { 279 Dec64(m1,m0); 280 res.as_struct.mant1 = m1; 281 res.as_struct.mant0 = m0; 282 return res.as_dbl; 283 } 284 return x; /* NaN */ 285} 286 287 288double 289ec_ieee_down(double x) 290{ 291 ieee_double res; 292 unsigned m1, m0; 293 res.as_dbl = x; 294 m1 = res.as_struct.mant1; 295 m0 = res.as_struct.mant0; 296 if ((m1 & ~SIGN_BIT | m0) == 0) /* x == +-0 */ 297 { 298 return -MINDOUBLE; 299 } 300 else if ((word) m1 > 0 || m1 == 0 && m0 > 0) /* x > 0 */ 301 { 302 if (m1 < 0x7ff00000 303 || m1 == 0x7ff00000 && m0 == 0) /* x <= Inf */ 304 { 305 Dec64(m1,m0); 306 res.as_struct.mant1 = m1; 307 res.as_struct.mant0 = m0; 308 return res.as_dbl; 309 } 310 } 311 else if (m1 < 0xfff00000) /* -Inf < x < 0 */ 312 { 313 Inc64(m1,m0); 314 res.as_struct.mant1 = m1; 315 res.as_struct.mant0 = m0; 316 return res.as_dbl; 317 } 318 return x; /* NaN */ 319} 320 321#endif 322 323extern void 324ec_i_add( 325 volatile double xl, /* Make sure the additions don't get reordered. */ 326 volatile double xu, 327 double yl, double yu, 328 double *lwb, double *upb) 329{ 330 volatile double tmp_upb; 331 volatile double tmp_lwb; 332 333 set_round_up(); 334 tmp_upb = xu + yu; 335 set_round_down(); 336 tmp_lwb = xl + yl; 337 restore_round_mode(); 338 339 if (!GoodFloat(tmp_lwb) || !GoodFloat(tmp_upb)) { 340 /* One or other of the resulting bounds is a NaN, so the "true" */ 341 /* value could be anything */ 342 *lwb = -HUGE_VAL; 343 *upb = HUGE_VAL; 344 return; 345 } 346 347 /* For any breal with one bound set to zero, the sign of the zero 348 should be the same as the sign of the other bound*/ 349 if (tmp_lwb == 0.0 && tmp_upb > 0.0) { 350 /* ensure lwb is positive zero */ 351 *lwb = 0.0; 352 } else { 353 *lwb = tmp_lwb; 354 } 355 if (tmp_upb == 0.0 && tmp_lwb < 0.0) { 356 /* ensure upb is negative zero */ 357 *upb = NEGZERO; 358 } else { 359 *upb = tmp_upb; 360 } 361} 362 363extern void 364ec_i_sub(double xl, double xu, double yl, double yu, double *lwb, double *upb) 365{ 366 ec_i_add(xl, xu, -yu, -yl, lwb, upb); 367} 368 369/* Multiplies a bounded real by another bounded real. */ 370 371/* 372** The usual way to multiply two bounded numbers of unknown sign is to 373** compute the product of the 4 relevant bound pairs and select the smallest 374** and largest as the resulting bounds. If we're using processor rounding 375** modes we probably need to do the calculations twice: once rounding up and 376** once rounding down. 377** 378** Unfortunately, since zeros compare equal, it is hard to ensure that -0.0 379** is selected over 0.0 for a lower bound, and vice-versa for the upper 380** bound. As a result, we do something more complicated, based on the 381** following table. The table gives, for each possible combination of signs 382** of the input variables, the signs of each of the 4 products, and which 383** products will give the lower and upper bounds of the results: 384** 385** z = x * y 386** 387** sign of | sign of | min/max 388** inputs | products | products 389** ------------+-------------+------------- 390** xl xu yl yu | ll lu ul uu | zl zu 391** ------------+-------------+------------- 392** - - - - | + + + + | uu ll 393** - - - + | + - + - | lu ll 394** - - + + | - - - - | lu ul 395** - + - - | + + - - | ul ll 396** - + - + | + - - + | lu/ul ll/uu 397** - + + + | - - + + | lu uu 398** + + - - | - - - - | ul lu 399** + + - + | - + - + | ul uu 400** + + + + | + + + + | ll uu 401** 402** For computing the upper bound, we first determine whether xl and yl have 403** the same sign. Extracting the relevant subtables from the above we have: 404** 405** Case samesign(xl, yl): 406** 407** sign of | sign of | max 408** inputs | products | product 409** ------------+-------------+--------- 410** xl xu yl yu | ll lu ul uu | zu 411** ------------+-------------+--------- 412** - - - - | + + + + | ll 413** - - - + | + - + - | ll 414** - + - - | + + - - | ll 415** - + - + | + - - + | ll/uu 416** + + + + | + + + + | uu 417** 418** In this case the upper bound can only be ll or uu, and since ll 419** always has positive sign, it is safe wrt zeros to choose ll over uu 420** if they compare equal (if they're zero, it will yield the positive 421** one). 422** 423** Case !samesign(xl, yl): 424** 425** sign of | sign of | max 426** inputs | products | product 427** ------------+-------------+--------- 428** xl xu yl yu | ll lu ul uu | zu 429** ------------+-------------+--------- 430** - - + + | - - - - | ul 431** - + + + | - - + + | uu 432** + + - - | - - - - | lu 433** + + - + | - + - + | uu 434** 435** In this case we don't have to consider ll, and when products compare 436** equal it is safe to choose uu (if any of the other products have 437** positive sign, uu will too). 438** 439** For the lower bounds we can obtain similarly useful distinctions by 440** comparing the signs for xl and yu: 441** 442** Case samesign(xl, yu): 443** 444** sign of | sign of | min 445** inputs | products | product 446** ------------+-------------+--------- 447** xl xu yl yu | ll lu ul uu | zl 448** ------------+-------------+--------- 449** - - - - | + + + + | uu 450** - + - - | + + - - | ul 451** + + - + | - + - + | ul 452** + + + + | + + + + | ll 453** 454** In this case we don't have to consider lu, and when products compare 455** equal it is safe to choose ul (if any of the other products have 456** negative sign, ul will too). 457** 458** Case !samesign(xl, yu): 459** 460** sign of | sign of | min 461** inputs | products | product 462** ------------+-------------+--------- 463** xl xu yl yu | ll lu ul uu | zl 464** ------------+-------------+--------- 465** - - - + | + - + - | lu 466** - - + + | - - - - | lu 467** - + - + | + - - + | lu/ul 468** - + + + | - - + + | lu 469** + + - - | - - - - | ul 470** 471** In this case the lower bound can only be lu or ul, and since lu 472** always has negative sign, it is safe wrt zeros to choose lu over ul 473** if they compare equal (if they're zero, it will yield the negative 474** one). 475*/ 476 477extern void 478ec_i_mul( 479 volatile double xl, /* Make sure multiplications don't get cached. */ 480 volatile double xu, 481 double yl, double yu, 482 double *lwb, double *upb) 483{ 484 double ll; 485 double lu; 486 double ul; 487 double uu; 488 double ll_lu, ul_uu; 489 490 /* 491 ** Note that we catch NaNs in ll and uu when computing the upper bound 492 ** and NaNs in lu and ul when computing the lower bound --- we assume 493 ** the rounding mode does not affect a NaN result. 494 */ 495 496 /* Compute the upper bound. */ 497 498 set_round_up(); 499 500 ll = xl * yl; 501 if (!GoodFloat(ll)) goto I_MUL_BAD_FLOAT; 502 uu = xu * yu; 503 if (!GoodFloat(uu)) goto I_MUL_BAD_FLOAT; 504 505 if (samesign(xl, yl)) { 506 /* Upper bound is either ll or uu. */ 507 /* When equal, choosing ll over uu is safe wrt signed zeros. */ 508 *upb = ll >= uu ? ll : uu; 509 } else { 510 /* Upper bound can't be ll. */ 511 /* When equal, choosing uu is safe wrt signed zeros. */ 512 ul = xu * yl; 513 lu = xl * yu; 514 ul_uu = uu >= ul ? uu : ul; 515 *upb = ul_uu >= lu ? ul_uu : lu; 516 } 517 518 /* Compute the lower bound. */ 519 520 set_round_down(); 521 522 lu = xl * yu; 523 if (!GoodFloat(lu)) goto I_MUL_BAD_FLOAT; 524 ul = xu * yl; 525 if (!GoodFloat(ul)) goto I_MUL_BAD_FLOAT; 526 527 if (samesign(xl, yu)) { 528 /* Lower bound can't be lu. */ 529 /* When equal, choosing ul is safe wrt signed zeros. */ 530 ll = xl * yl; 531 uu = xu * yu; 532 ul_uu = ul <= uu ? ul : uu; 533 *lwb = ul_uu <= ll ? ul_uu : ll; 534 } else { 535 /* Lower bound is either lu or ul. */ 536 /* When equal, choosing lu over ul is safe wrt signed zeros. */ 537 *lwb = lu <= ul ? lu : ul; 538 } 539 540 restore_round_mode(); 541 return; 542 543I_MUL_BAD_FLOAT: 544 /* We got a NaN, so just return an infinite interval. */ 545 restore_round_mode(); 546 *lwb = -HUGE_VAL; 547 *upb = HUGE_VAL; 548 return; 549} 550 551/* Divides a bounded real by another bounded real. */ 552/* See the notes for multiplication above. */ 553/* Note that it might be possible to exploit the fact that yl and yu will */ 554/* have the same sign in the main body of the code. */ 555extern void 556ec_i_div( 557 volatile double xl, /* Make sure divisions don't get cached. */ 558 volatile double xu, 559 double yl, double yu, 560 double *lwb, double *upb) 561{ 562 double ll; 563 double lu; 564 double ul; 565 double uu; 566 double ll_lu, ul_uu; 567 568 /* Correctly handle division by negative zero */ 569 if (yl == 0.0) { 570 /* This will ensure that if yl == 0.0 or yl == -0.0 then 571 it will be set to 0.0 */ 572 yl = 0.0; 573 } 574 575 if (yu == 0.0) { 576 /* This will ensure that if yu == 0.0 or yu == -0.0 then 577 it will be set to -0.0 */ 578 yu = NEGZERO; 579 } 580 581 /* If divisor spans zero, nothing to be done */ 582 if (!samesign(yl, yu)) { 583 *lwb = -HUGE_VAL; 584 *upb = HUGE_VAL; 585 return; 586 } 587 588 /* 589 ** Note that we catch NaNs in ll and uu when computing the upper bound 590 ** and NaNs in lu and ul when computing the lower bound --- we assume 591 ** the rounding mode does not affect a NaN result. 592 */ 593 594 /* Compute the upper bound. */ 595 596 set_round_up(); 597 598 ll = xl / yu; 599 if (!GoodFloat(ll)) goto I_DIV_BAD_FLOAT; 600 uu = xu / yl; 601 if (!GoodFloat(uu)) goto I_DIV_BAD_FLOAT; 602 603 if (samesign(xl, yu)) { 604 /* Upper bound is either ll or uu. */ 605 /* When equal, choosing ll over uu is safe wrt signed zeros. */ 606 *upb = ll >= uu ? ll : uu; 607 } else { 608 /* Upper bound can't be ll. */ 609 /* When equal, choosing uu is safe wrt signed zeros. */ 610 ul = xu / yu; 611 lu = xl / yl; 612 ul_uu = uu >= ul ? uu : ul; 613 *upb = ul_uu >= lu ? ul_uu : lu; 614 } 615 616 /* Compute the lower bound. */ 617 618 set_round_down(); 619 620 lu = xl / yl; 621 if (!GoodFloat(lu)) goto I_DIV_BAD_FLOAT; 622 ul = xu / yu; 623 if (!GoodFloat(ul)) goto I_DIV_BAD_FLOAT; 624 625 if (samesign(xl, yl)) { 626 /* Lower bound can't be lu. */ 627 /* When equal, choosing ul is safe wrt signed zeros. */ 628 ll = xl / yu; 629 uu = xu / yl; 630 ul_uu = ul <= uu ? ul : uu; 631 *lwb = ul_uu <= ll ? ul_uu : ll; 632 } else { 633 /* Lower bound is either lu or ul. */ 634 /* When equal, choosing lu over ul is safe wrt signed zeros. */ 635 *lwb = lu <= ul ? lu : ul; 636 } 637 638 restore_round_mode(); 639 return; 640 641I_DIV_BAD_FLOAT: 642 /* We got a NaN, so just return an infinite interval. */ 643 restore_round_mode(); 644 *lwb = -HUGE_VAL; 645 *upb = HUGE_VAL; 646 return; 647} 648 649 650static void 651i_exp(double xl, double xu, double *lwb, double *upb) 652{ 653 /* Catch the special cases of raising 'e' to +Inf and -Inf as some 654 platforms give incorrect results w.r.t. IEEE754 specs */ 655 if ( xl == -HUGE_VAL ) { 656 *lwb = 0.0; 657 } else if ( xl == HUGE_VAL ) { 658 *lwb = HUGE_VAL; 659 } else { 660 *lwb = down(exp(xl)); 661 } 662 if ( xu == -HUGE_VAL ) { 663 *upb = 0.0; 664 } else if ( xu == HUGE_VAL ) { 665 *upb = HUGE_VAL; 666 } else { 667 *upb = up(exp(xu)); 668 } 669} 670 671static void 672i_sin(double xl, double xu, double *lwb, double *upb, int cosflag) 673{ 674 double width; 675 volatile double xl1 = xl; /* Hack to stop reordering */ 676 677 set_round_up(); 678 width = xu-xl1; 679 restore_round_mode(); 680 681 if (width >= 2*M_PI) 682 { 683 *lwb = -1.0; *upb = 1.0; 684 } 685 else 686 { 687 double ls, lc, us, uc; 688 689 if (cosflag) { 690 sincos(xl, &lc, &ls); 691 sincos(xu, &uc, &us); 692 lc = -lc; uc = -uc; 693 } else { 694 sincos(xl, &ls, &lc); 695 sincos(xu, &us, &uc); 696 } 697 698 if (lc >= 0) { 699 if (uc >= 0) { 700 if (width >= M_PI) { 701 *lwb = -1.0; *upb = 1.0; 702 } else { 703 *lwb = down(ls); *upb = up(us); 704 } 705 } else { 706 *lwb = down(min(ls, us)); *upb = 1.0; 707 } 708 } else { 709 if (uc >= 0) { 710 *lwb = -1.0; *upb = up(max(ls, us)); 711 } else { 712 if (width >= M_PI) { 713 *lwb = -1.0; *upb = 1.0; 714 } else { 715 *lwb = down(us); *upb = up(ls); 716 } 717 } 718 } 719 } 720} 721 722 723static double 724ipow(double x, int n, int roundup) /* currently n > 0 */ 725{ 726 double res; 727 int neg = 0; 728 if (x < 0) { 729 x = -x; 730 if (n&1) { 731 neg = !neg; 732 roundup = !roundup; 733 } 734 } 735 if (n < 5) { 736 /* faster, but less precise... */ 737 volatile double x1 = x; /* Hack to stop reordering */ 738 if (roundup) { 739 set_round_up(); 740 } else { 741 set_round_down(); 742 } 743 res = x; 744 while(--n > 0) res *= x1; 745 restore_round_mode(); 746 } else { 747 res = roundup ? up(Pow(x,(double)n)) : down(Pow(x,(double)n)); 748 } 749 return neg ? -res : res; 750} 751 752 753static int 754linsplit(double minwidth, double ratio, double lwb, double upb, int upper, double *split) 755{ 756 /* 757 * Splits an interval given by lwb - upb, at 1/Ratio if 758 * upper == FALSE or 1/1-Ratio if upper = TRUE giving a 759 * new split point *split. 760 * We fail if the input interval is already narrower than MinWidth. 761 * Or if the piece removed is smaller than minwidth and the new 762 * interval is greater than minwidth. 763 */ 764 765 double halfwidth,halfmin; /* never a NaN */ 766 double halfchunk; 767 double scale; 768 769 /* Quick check which ensures (among other things) that the lower */ 770 /* bound is not +inf and the upper bound is not -inf. */ 771 if (lwb == upb) { 772 if (minwidth == 0.0) { 773 *split = lwb; 774 return 0; 775 } else { 776 /* Input interval is too small. */ 777 return -1; 778 } 779 } 780 781 if (!finite(upb)) { 782 upb = MAXDOUBLE; 783 } 784 if (!finite(lwb)) { 785 lwb = -MAXDOUBLE; 786 } 787 788 scale = fabs(upper ? upb : lwb); 789 if (scale > 1.0) minwidth = minwidth * scale; 790 791 halfmin = minwidth/2; 792 793 halfwidth = upb/2.0 - lwb/2.0; /* avoiding inf */ 794 if (halfwidth < halfmin) { 795 return -1; 796 } 797 798 halfchunk = ratio * halfwidth; 799 if (halfchunk < halfmin && halfwidth-halfchunk > halfmin) { 800 return -1; 801 } 802 /* PROBLEM: we might cut off very small bits here... */ 803 if (upper) { 804 *split = upb - 2*halfchunk; 805 } else { 806 *split = lwb + 2*halfchunk; 807 } 808 809 /* 810 ** In extreme cases (such as the user specifying a precision too 811 ** small), it might be possible for the split to lie outside the 812 ** bounds; in such cases we round it inwards. 813 */ 814 if (*split > upb) { 815 *split = upb; 816 } else if (*split < lwb) { 817 *split = lwb; 818 } 819 820 return 0; 821} 822 823 824static int 825logsplit(double minwidth, double ratio, double orig_lwb, double orig_upb, int upper, double *split) 826{ 827 /* 828 ** Splits the interval in such a way that between -1 and 1 the split 829 ** is linear, and outside that it's logarithmic. The rationale 830 ** behind this choice of splitting is that between -1 and 1 the 831 ** absolute error is smaller than the relative error (and hence 832 ** should be used for deciding when to stop splitting), while 833 ** outside that region the reverse is true. 834 ** 835 ** The splitting is done by mapping values x > 1 to 1 + log(x) and 836 ** values x < -1 to -1 - log(-x), and then calling the linear 837 ** splitting routine. The result is then transformed back using the 838 ** reverse function: values y > 1 are mapped to exp(y - 1) and 839 ** values y < -1 to -exp(-y - 1). 840 */ 841 842#define convert_bound_log(x) \ 843 if (x > 1) { \ 844 x = 1 + log(x); \ 845 } else if (x < -1) { \ 846 x = -1 - log(-x); \ 847 } 848 849#define unconvert_bound_log(y) \ 850 if (y > 1) { \ 851 y = exp(y - 1); \ 852 } else if (y < -1) { \ 853 y = -exp(-y - 1); \ 854 } 855 856 int result; 857 double lwb = orig_lwb, upb = orig_upb; 858 859 /* Quick check which ensures (among other things) that the lower */ 860 /* bound is not +inf and the upper bound is not -inf. */ 861 if (lwb == upb) { 862 if (minwidth == 0.0) { 863 *split = lwb; 864 return 0; 865 } else { 866 /* Input interval is too small. */ 867 return -1; 868 } 869 } 870 871 if (!finite(upb)) { 872 upb = MAXDOUBLE; 873 } 874 if (!finite(lwb)) { 875 lwb = -MAXDOUBLE; 876 } 877 878 convert_bound_log(lwb); 879 convert_bound_log(upb); 880 881 result = linsplit(minwidth, ratio, lwb, upb, upper, split); 882 if (result != 0) { 883 return result; 884 } 885 886 unconvert_bound_log(*split); 887 888 /* 889 ** In extreme cases (such as the user specifying a precision too 890 ** small), the split may not lie between the bounds; in such cases 891 ** we round it inwards. 892 */ 893 if (*split > orig_upb) { 894 *split = orig_upb; 895 } else if (*split < orig_lwb) { 896 *split = orig_lwb; 897 } 898 899 return 0; 900} 901 902 903extern 904int 905ec_ria_unop(int op, double xl, double xu, double *zl_ptr, double *zu_ptr) 906{ 907 double lwb, upb; 908 909 switch (op) 910 { 911 case RIA_UN_SQR: /* sqr */ 912 { 913 volatile double xl1 = xl, xu1 = xu; /* hack to stop reordering */ 914 if (xl >= 0.0) { 915 set_round_down(); 916 lwb = xl1 * xl1; 917 set_round_up(); 918 upb = xu1 * xu1; 919 } else if (xu < 0.0) { 920 set_round_down(); 921 lwb = xu1 * xu1; 922 set_round_up(); 923 upb = xl1 * xl1; 924 } else { 925 set_round_up(); 926 lwb = 0.0; 927 upb = -xl > xu ? xl1 * xl1 : xu1 * xu1; 928 } 929 restore_round_mode(); 930 } 931 break; 932 case RIA_UN_SQRT: /* sqrt */ 933 if (xu >= 0.0) { 934 upb = up(sqrt(xu)); 935 if (xl < 0.0) 936 lwb = 0.0; 937 else { 938 lwb = sqrt(xl); 939 if (lwb > 0.0) lwb = down(lwb); 940 } 941 } else { 942 Fail_; 943 } 944 break; 945 case RIA_UN_SIN: /* sin */ 946 i_sin(xl, xu, &lwb, &upb, 0); 947 break; 948 case RIA_UN_COS: /* cos */ 949 i_sin(xl, xu, &lwb, &upb, 1); 950 break; 951 case RIA_UN_EXP: /* exp */ 952 i_exp(xl, xu, &lwb, &upb); 953 break; 954 case RIA_UN_LN: /* ln */ 955 if (xu >= 0.0) { 956 lwb = xl > 0.0 ? down(log(xl)) : -HUGE_VAL; 957 upb = up(log(xu)); 958 } else { 959 Fail_; 960 } 961 break; 962 case RIA_UN_ATAN: /* atan */ 963 lwb = down(atan(xl)); 964 upb = up(atan(xu)); 965 break; 966 case RIA_UN_PI: /* pi */ 967 /* argument is dummy */ 968 lwb = 3.141592653589793; 969 upb = 3.1415926535897935; 970 break; 971 case RIA_UN_ABS: /* abs */ 972 if (xl >= 0.0) { 973 lwb = xl; upb = xu; 974 } else if (xu < 0.0) { 975 lwb = -xu; upb = -xl; 976 } else { 977 lwb = 0.0; upb = max(xu, -xl); 978 } 979 break; 980 case RIA_UN_ROUNDOUT: /* roundout */ 981 lwb = down(xl); 982 upb = up(xu); 983 break; 984 case RIA_UN_NEG: /* - */ 985 lwb = -xu; 986 upb = -xl; 987 break; 988 case RIA_UN_WIDTH: /* width */ 989 lwb = upb = xu-xl; 990 break; 991 default: 992 Bip_Error(RANGE_ERROR); 993 } 994 995 *zl_ptr = lwb; 996 *zu_ptr = upb; 997 998 Succeed_; 999} 1000 1001 1002extern 1003int 1004ec_ria_binop(int op, double xl, double xu, double yl, double yu, double *zl_ptr, double *zu_ptr) 1005{ 1006 double lwb, upb; 1007 int result; 1008 1009 switch (op) 1010 { 1011 case RIA_BIN_ADD: /* + */ 1012 ec_i_add(xl, xu, yl, yu, &lwb, &upb); 1013 break; 1014 case RIA_BIN_SUB: /* - */ 1015 ec_i_sub(xl, xu, yl, yu, &lwb, &upb); 1016 break; 1017 case RIA_BIN_MULT: /* * */ 1018 ec_i_mul(xl, xu, yl, yu, &lwb, &upb); 1019 break; 1020 case RIA_BIN_DIV: /* / */ 1021 ec_i_div(xl, xu, yl, yu, &lwb, &upb); 1022 break; 1023 case RIA_BIN_RSQR: /* rsqr i.e. +-sqrt */ 1024 { 1025 /* second argument is used to check result range */ 1026 double rmin, rmax; 1027 if (xu >= 0.0) { 1028 rmax = up(sqrt(xu)); 1029 rmin = xl <= 0.0 ? 0.0 : down(sqrt(xl)); 1030 if (yl > -rmin) { 1031 upb = rmax; lwb = rmin; 1032 } else if (yu < rmin) { 1033 upb = -rmin; lwb = -rmax; 1034 } else { 1035 upb = rmax; lwb = -rmax; 1036 } 1037 } else { 1038 Fail_; 1039 } 1040 break; 1041 } 1042 case RIA_BIN_POW_INT: /* pow(X,int) */ 1043 { 1044 word n = (word) yl; 1045 if (n > 0) 1046 { 1047 if (n&1) { /* odd */ 1048 lwb = ipow(xl,n,DOWN); 1049 upb = ipow(xu,n,UP); 1050 } else { /* even */ 1051 if (xl >= 0.0) { 1052 lwb = ipow(xl,n,DOWN); 1053 upb = ipow(xu,n,UP); 1054 } else if (xu < 0.0) { 1055 lwb = ipow(xu,n,DOWN); 1056 upb = ipow(xl,n,UP); 1057 } else { /* 0-crossing */ 1058 lwb = 0.0; 1059 upb = xu > -xl ? ipow(xu,n,UP) 1060 : ipow(xl,n,UP); 1061 } 1062 } 1063 } 1064 else if (n < 0) 1065 { 1066 Bip_Error(RANGE_ERROR); 1067 } 1068 else 1069 { 1070 lwb = 1.0; upb = 1.0; 1071 } 1072 break; 1073 } 1074 case RIA_BIN_RPOW_ODD: /* Z=pow(X,Y), reverse of X=pow(Z,N), N odd */ 1075 lwb = xl < 0.0 ? 1076 -up(Pow(-xl, -xl>1? yu: yl)) : 1077 down(Pow(xl, xl>1? yl: yu)); 1078 upb = xu < 0.0 ? 1079 -down(Pow(-xu, -xu>1? yl: yu)) : 1080 up(Pow(xu, xu>1? yu: yl)); 1081 break; 1082#if 0 1083 case 7: /* pow(X,Y) */ 1084 if (xu < 0.0) { 1085 Fail_; 1086 } else { 1087 if (yl > 0.0) { 1088 lwb = xl <= 0.0 ? 0.0 : 1089 down(pow(xl, xl>1? yl: yu)); 1090 upb = up(pow(xu, xu>1? yu: yl)); 1091 } else if (yu < 0.0) { 1092 lwb = down(pow(xu, xu>1? yl: yu)); 1093 upb = xl <= 0.0 ? HUGE_VAL : 1094 up(pow(xl, xl>1? yu: yl)); 1095 } else if (xl <= 0.0) { /* X and Y-range include 0 */ 1096 lwb = 0.0; upb = HUGE_VAL; 1097 } else { /* Y-range includes 0 */ 1098 double ll = pow(xl,yl); 1099 double lu = pow(xl,yu); 1100 double ul = pow(xu,yl); 1101 double uu = pow(xu,yu); 1102 lwb = down(min(min(ll,lu),min(ul,uu))); 1103 upb = up(max(max(ll,lu),min(ul,uu))); 1104 } 1105 } 1106 break; 1107 case RIA_BIN_RELAX: /* relax(X,const) */ 1108 lwb = down(xl-yu); 1109 upb = up(xu+yu); 1110 break; 1111#endif 1112 case RIA_BIN_MIN: /* min(X,Y) */ 1113 lwb = min(xl,yl); 1114 upb = min(xu,yu); 1115 break; 1116 case RIA_BIN_MAX: /* max(X,Y) */ 1117 lwb = max(xl,yl); 1118 upb = max(xu,yu); 1119 break; 1120 case RIA_BIN_LOGSPLIT: /* logsplit(X,MinWidth,Ratio) */ 1121 if (0 != logsplit(yl, yu, xl, xu, 0, &lwb)) 1122 Fail_; 1123 break; 1124 case RIA_BIN_PLUSMINUS: /* +- */ 1125 { 1126 /* second argument is used to check result range */ 1127 double absmin, absmax; 1128 lwb = xl; 1129 upb = xu; 1130 if (lwb >= 0.0) { 1131 absmin = lwb; absmax = upb; 1132 } else if (upb <= 0.0) { 1133 absmin = -upb; absmax = -lwb; 1134 } else { 1135 absmin = 0.0; absmax = max(-lwb, upb); 1136 } 1137 if (yl > -absmin) { 1138 lwb = absmin; upb = absmax; 1139 } else if (yu < absmin) { 1140 lwb = -absmax; upb = -absmin; 1141 } else { 1142 lwb = -absmax; upb = absmax; 1143 } 1144 break; 1145 } 1146 case RIA_BIN_MIN_DELTA: /* compute min of relative and absolute delta */ 1147 { 1148 double lwb_delta = yl-xl; /* absolute delta */ 1149 double upb_delta = xu-yu; /* absolute delta */ 1150 1151#ifdef _WIN32 1152 /* Hack for Windows, where inf/inf is not a NAN, but a small float */ 1153 /* (meaning an infinite bound changes return a negligible delta). */ 1154 if (lwb_delta == HUGE_VAL) { 1155 lwb = lwb_delta; 1156 } else 1157#endif 1158 if (lwb_delta > 0.0) { /* lwb_delta may be NaN! */ 1159 lwb = lwb_delta/fabs(xl); /* relative delta */ 1160 if (!(lwb < lwb_delta)) /* lwb mab be NaN! */ 1161 lwb = lwb_delta; 1162 } else { 1163 lwb = 0.0; /* no change */ 1164 } 1165#ifdef _WIN32 1166 /* Hack for Windows, where inf/inf is not a NAN, but a small float */ 1167 /* (meaning an infinite bound changes return a negligible delta). */ 1168 if (upb_delta == HUGE_VAL) { 1169 upb = upb_delta; 1170 } else 1171#endif 1172 if (upb_delta > 0.0) { /* upb_delta may be NaN */ 1173 upb = upb_delta/fabs(xu); /* relative delta */ 1174 if (!(upb < upb_delta)) /* upb may be NaN */ 1175 upb = upb_delta; 1176 } else { 1177 upb = 0.0; /* no change */ 1178 } 1179 break; 1180 } 1181 case RIA_BIN_LINSPLIT: /* linsplit(X,MinWidth,Ratio) */ 1182 if (0 != linsplit(yl, yu, xl, xu, 0, &lwb)) 1183 Fail_; 1184 break; 1185 case RIA_BIN_LINSPLIT_UPPER: /* linsplit(X,MinWidth,Ratio) */ 1186 if (0 != linsplit(yl, yu, xl, xu, 1, &lwb)) 1187 Fail_; 1188 break; 1189 case RIA_BIN_LOGSPLIT_UPPER: /* logsplit(X,MinWidth,Ratio) */ 1190 if (0 != logsplit(yl, yu, xl, xu, 1, &lwb)) 1191 Fail_; 1192 break; 1193#if 0 1194 case 17: /* round(XL,XH,Precision) */ 1195 lwb = xl - remainder(xl,yl); 1196 upb = xu - remainder(xu,yu); 1197 break; 1198#endif 1199 default: 1200 Bip_Error(RANGE_ERROR); 1201 } 1202 1203 *zl_ptr = lwb; 1204 *zu_ptr = upb; 1205 1206 Succeed_; 1207} 1208 1209 1210extern 1211int 1212ec_ria_ternop(int op, double xl, double xu, double yl, double yu, double zl, double zu, double *rl_ptr, double *ru_ptr) 1213{ 1214 double lwb, upb; 1215 int result; 1216 switch (op) 1217 { 1218 case RIA_TERN_RPOW_EVEN: /* rpow_even */ 1219 /* Reverse of Z^N=X, N>=2 and even. Computes R = intersect( X^Y, Z), Y=1/N */ 1220 if (xu >= 0.0) { 1221 lwb = xl <= 0.0 ? 0.0 : down(Pow(xl, yl)); 1222 upb = up(SafePow(xu, yu)); 1223 if (zl > -lwb) { 1224 ; /* only positive solution */ 1225 } else if (zu < lwb) { 1226 double aux = upb; /* only negative solution */ 1227 upb = -lwb; 1228 lwb = -aux; 1229 } else { 1230 lwb = -upb; /* both solutions */ 1231 } 1232 } else { 1233 Fail_; 1234 } 1235 break; 1236 case RIA_TERN_UNION: /* union(X,Y,Z) */ 1237 if (xu < yl) { /* disjoint, X below Y */ 1238 lwb = xu < zl ? yl : xl; 1239 upb = zu < yl ? xu : yu; 1240 if (lwb > upb) { Fail_; } 1241 } else if (yu < xl) { /* disjoint, Y below X */ 1242 lwb = yu < zl ? xl : yl; 1243 upb = zu < xl ? yu : xu; 1244 if (lwb > upb) { Fail_; } 1245 } else { 1246 lwb = min(xl,yl); 1247 upb = max(xu,yu); 1248 } 1249 break; 1250 case RIA_TERN_DIV: /* / */ 1251 if (!samesign(yl, yu)) { 1252 volatile double xl1 = xl, xu1 = xu; /* hack to stop reordering */ 1253 /* 1254 ** Want to work out the union of the intervals: 1255 ** X / yu .. +inf 1256 ** -inf .. X / yl 1257 ** and match them against the interval for Z. 1258 */ 1259 set_round_down(); 1260 lwb = xl1 / yu; 1261 set_round_up(); 1262 upb = xu1 / yl; 1263 if (zl > upb) { 1264 /* Lower computed interval falls outside Z's range, */ 1265 /* so return upper computed interval. */ 1266 upb = HUGE_VAL; 1267 } else if (zu < lwb) { 1268 /* Upper computed interval falls outside Z's range, */ 1269 /* so return lower computed interval. */ 1270 lwb = -HUGE_VAL; 1271 } else { 1272 /* Z overlaps both intervals, so can't choose one yet. */ 1273 lwb = -HUGE_VAL; 1274 upb = HUGE_VAL; 1275 } 1276 restore_round_mode(); 1277 } else { 1278 ec_i_div(xl, xu, yl, yu, &lwb, &upb); 1279 } 1280 break; 1281 default: 1282 Bip_Error(RANGE_ERROR); 1283 } 1284 1285 *rl_ptr = lwb; 1286 *ru_ptr = upb; 1287 1288 Succeed_; 1289} 1290 1291 1292/* ---------------------------------------------------------------------- 1293 * Prolog Interface 1294 * ---------------------------------------------------------------------- */ 1295 1296int 1297p_breal_from_bounds(value vl, type tl, value vu, type tu, value vx, type tx) 1298{ 1299 double lo, hi; 1300 value ivl; 1301 int err; 1302 1303 Check_Number(tl); 1304 Check_Number(tu); 1305 Check_Output_Type(tx, TIVL); 1306 1307 /* Get the lower bound to use, coercing it (safely) if required. */ 1308 if (IsDouble(tl)) { 1309 lo = Dbl(vl); 1310 } else { 1311 if (IsInterval(tl)) { 1312 ivl = vl; 1313 } else { 1314 err = tag_desc[TagType(tl)].coerce_to[TIVL](vl, &ivl); 1315 if (err != PSUCCEED) return err; 1316 } 1317 lo = IvlLwb(ivl.ptr); 1318 } 1319 1320 /* Get the upper bound to use, coercing it (safely) if required. */ 1321 if (IsDouble(tu)) { 1322 hi = Dbl(vu); 1323 } else { 1324 value ivl; 1325 if (IsInterval(tu)) { 1326 ivl = vu; 1327 } else { 1328 err = tag_desc[TagType(tu)].coerce_to[TIVL](vu, &ivl); 1329 if (err != PSUCCEED) return err; 1330 } 1331 hi = IvlUpb(ivl.ptr); 1332 } 1333 1334 Return_Unify_Interval(vx, tx, lo, hi); 1335} 1336 1337 1338int 1339p_breal_min(value vx, type tx, value vmin, type tmin) 1340{ 1341 double min; 1342 1343 Check_Output_Float(tmin); 1344 1345 if (IsDouble(tx)) { 1346 Return_Unify_Pw(vmin, tmin, vx, tx); 1347 } else if (IsInterval(tx)) { 1348 Return_Unify_Double(vmin, tmin, IvlLwb(vx.ptr)); 1349 } else { 1350 value ivl; 1351 int result; 1352 1353 Check_Number(tx); 1354 result = tag_desc[TagType(tx)].coerce_to[TIVL](vx, &ivl); 1355 Return_If_Not_Success(result); 1356 1357 Return_Unify_Double(vmin, tmin, IvlLwb(ivl.ptr)); 1358 } 1359} 1360 1361 1362int 1363p_breal_max(value vx, type tx, value vmax, type tmax) 1364{ 1365 Check_Output_Float(tmax); 1366 1367 if (IsDouble(tx)) { 1368 Return_Unify_Pw(vmax, tmax, vx, tx); 1369 } else if (IsInterval(tx)) { 1370 Return_Unify_Double(vmax, tmax, IvlUpb(vx.ptr)); 1371 } else { 1372 value ivl; 1373 int result; 1374 1375 Check_Number(tx); 1376 result = tag_desc[TagType(tx)].coerce_to[TIVL](vx, &ivl); 1377 Return_If_Not_Success(result); 1378 1379 Return_Unify_Double(vmax, tmax, IvlUpb(ivl.ptr)); 1380 } 1381} 1382 1383 1384int 1385p_breal_bounds(value vx, type tx, value vmin, type tmin, value vmax, type tmax) 1386{ 1387 Prepare_Requests 1388 1389 Check_Output_Float(tmin); 1390 Check_Output_Float(tmax); 1391 1392 if (IsDouble(tx)) { 1393 Request_Unify_Pw(vmin, tmin, vx, tx); 1394 Request_Unify_Pw(vmax, tmax, vx, tx); 1395 } else if (IsInterval(tx)) { 1396 Request_Unify_Double(vmin, tmin, IvlLwb(vx.ptr)); 1397 Request_Unify_Double(vmax, tmax, IvlUpb(vx.ptr)); 1398 } else { 1399 value ivl; 1400 int result; 1401 1402 Check_Number(tx); 1403 result = tag_desc[TagType(tx)].coerce_to[TIVL](vx, &ivl); 1404 Return_If_Not_Success(result); 1405 1406 Request_Unify_Double(vmin, tmin, IvlLwb(ivl.ptr)); 1407 Request_Unify_Double(vmax, tmax, IvlUpb(ivl.ptr)); 1408 } 1409 1410 Return_Unify 1411} 1412 1413 1414/* 1415** NOTE: 1416** The following predicates (ec_ria_unop, ec_ria_binop, ec_ria_ternop) are 1417** intended to only be called with fresh new variables for the result 1418** bounds (zl/zu for unop/binop, rl/ru for ternop). Using old variables 1419** will likely do the wrong thing. 1420*/ 1421 1422int 1423p_ria_unop(value vop, type top, value v_xl, type t_xl, value v_xu, type t_xu, value v_zl, type t_zl, value v_zu, type t_zu) 1424{ 1425 double lwb, upb; 1426 int result; 1427 1428 Check_Double(t_xl); Check_Double(t_xu); 1429 Check_Ref(t_zl); Check_Ref(t_zu); 1430 1431 result = ec_ria_unop(vop.nint, Dbl(v_xl), Dbl(v_xu), &lwb, &upb); 1432 Return_If_Not_Success(result); 1433 1434 Return_Double(v_zl, t_zl, lwb); 1435 Return_Double(v_zu, t_zu, upb); 1436 Succeed_; 1437} 1438 1439 1440int 1441p_ria_binop(value vop, type top, value v_xl, type t_xl, value v_xu, type t_xu, value v_yl, type t_yl, value v_yu, type t_yu, value v_zl, type t_zl, value v_zu, type t_zu) 1442{ 1443 double lwb, upb; 1444 int result; 1445 1446 Check_Double(t_xl); Check_Double(t_xu); 1447 Check_Double(t_yl); Check_Double(t_yu); 1448 Check_Ref(t_zl); Check_Ref(t_zu); 1449 1450 result = ec_ria_binop(vop.nint, Dbl(v_xl), Dbl(v_xu), Dbl(v_yl), Dbl(v_yu), 1451 &lwb, &upb); 1452 Return_If_Not_Success(result); 1453 1454 Return_Double(v_zl, t_zl, lwb); 1455 Return_Double(v_zu, t_zu, upb); 1456 Succeed_; 1457} 1458 1459 1460int 1461p_ria_ternop(value vop, type top, value v_xl, type t_xl, value v_xu, type t_xu, value v_yl, type t_yl, value v_yu, type t_yu, value v_zl, type t_zl, value v_zu, type t_zu, value v_rl, type t_rl, value v_ru, type t_ru) 1462{ 1463 double lwb, upb; 1464 int result; 1465 1466 Check_Double(t_xl); Check_Double(t_xu); 1467 Check_Double(t_yl); Check_Double(t_yu); 1468 Check_Double(t_zl); Check_Double(t_zu); 1469 Check_Ref(t_rl); Check_Ref(t_ru); 1470 1471 result = ec_ria_ternop(vop.nint, Dbl(v_xl), Dbl(v_xu), Dbl(v_yl), Dbl(v_yu), 1472 Dbl(v_zl), Dbl(v_zu), &lwb, &upb); 1473 Return_If_Not_Success(result); 1474 1475 Return_Double(v_rl, t_rl, lwb); 1476 Return_Double(v_ru, t_ru, upb); 1477 Succeed_; 1478} 1479 1480 1481/*-------------------------------------------------------------------------- 1482 * Methods 1483 *--------------------------------------------------------------------------*/ 1484 1485#define IVL_STRING_SIZE 64 1486 1487 1488static int 1489_int_ivl(value in, value *out) /* CAUTION: we allow out == &in */ 1490{ 1491 pword *pw; 1492 1493#ifdef DOUBLE_INT_LIMIT 1494 /* We're on a machine where an integer may not be exactly representable 1495 * as a double. 1496 */ 1497 if (in.nint > DOUBLE_INT_LIMIT || in.nint < -DOUBLE_INT_LIMIT) { 1498 /* The integer is not exactly representable, so we need to widen. */ 1499 Push_Interval(pw, down((double) in.nint), up((double) in.nint)); 1500 } else 1501#endif 1502 Push_Interval(pw, (double) in.nint, (double) in.nint); 1503 1504 out->ptr = pw; 1505 Succeed_; 1506} 1507 1508 1509static int 1510_dbl_ivl(value in, value *out) /* CAUTION: we allow out == &in */ 1511{ 1512 pword *pw; 1513 double in_dbl = Dbl(in) ; 1514 Push_Interval(pw, in_dbl, in_dbl); 1515 out->ptr = pw; 1516 Succeed_; 1517} 1518 1519 1520static int 1521_big_ivl(value in, value *out) /* CAUTION: we allow out == &in */ 1522{ 1523 pword *pw; 1524 value dval; 1525 double d; 1526 int err; 1527 err = tag_desc[TBIG].coerce_to[TDBL](in, &dval); 1528 if (err != PSUCCEED) return(err); 1529 d = Dbl(dval); 1530 if (d >= DOUBLE_INT_LIMIT_AS_DOUBLE || d <= -DOUBLE_INT_LIMIT_AS_DOUBLE) { 1531 /* The integer is not exactly representable, so we need to widen. */ 1532 Push_Interval(pw, down(d), up(d)); 1533 } else { 1534 Push_Interval(pw, d, d); 1535 } 1536 out->ptr = pw; 1537 Succeed_; 1538} 1539 1540 1541/*ARGSUSED*/ 1542static int 1543_ivl_dbl(value in, value *out) /* CAUTION: we allow out == &in */ 1544{ 1545 double d; 1546 if (logsplit(0.0, 0.5, IvlLwb(in.ptr), IvlUpb(in.ptr), 0, &d) != 0) 1547 { 1548 Bip_Error(RANGE_ERROR); 1549 } 1550 Make_Double_Val(*out, d); 1551 Succeed_; 1552} 1553 1554 1555/*ARGSUSED*/ 1556static int 1557_ivl_string_size(value v, type t, int quoted) 1558{ 1559 return IVL_STRING_SIZE; 1560} 1561 1562 1563/*ARGSUSED*/ 1564static int 1565_ivl_to_string(value v, type t, char *buf, int quoted) 1566{ 1567 int len; 1568 type tdbl; 1569 value dval; 1570 tdbl.kernel = TDBL; 1571#ifdef UNBOXED_DOUBLES 1572 Make_Double_Val(dval, IvlLwb(v.ptr)); 1573 len = tag_desc[TDBL].to_string(dval, tdbl, buf, quoted); 1574 buf[len++] = '_'; 1575 buf[len++] = '_'; 1576 Make_Double_Val(dval, IvlUpb(v.ptr)); 1577 len += tag_desc[TDBL].to_string(dval, tdbl, buf+len, quoted); 1578#else 1579 dval.ptr = v.ptr; /* dirty */ 1580 len = tag_desc[TDBL].to_string(dval, tdbl, buf, quoted); 1581 buf[len++] = '_'; 1582 buf[len++] = '_'; 1583 dval.ptr = (pword*)((double*)(v.ptr)+1); /* dirty */ 1584 len += tag_desc[TDBL].to_string(dval, tdbl, buf+len, quoted); 1585#endif 1586 return len; 1587} 1588 1589 1590/*ARGSUSED*/ 1591static int 1592_write_ivl(int quoted, stream_id stream, value v, type t) 1593{ 1594 char buf[IVL_STRING_SIZE]; 1595 int len = _ivl_to_string(v, t, buf, quoted); 1596 return ec_outf(stream, buf, len); 1597} 1598 1599 1600static int 1601_ivl_from_string(char *s, pword *result, int base) 1602{ 1603 (void) string_to_number(s, result, (stream_id) 0, 0); 1604 if (IsTag(result->tag.kernel, TEND)) 1605 { Bip_Error(BAD_FORMAT_STRING) } 1606 Succeed_; 1607} 1608 1609 1610static int 1611_ivl_nop(value v1, pword *pres) 1612{ 1613 pres->tag.kernel = TIVL; 1614 pres->val.all = v1.all; 1615 Succeed_; 1616} 1617 1618 1619static int 1620_ivl_add(value v1, value v2, pword *pres) 1621{ 1622 double lwb, upb; 1623 double v1_lwb, v1_upb, v2_lwb, v2_upb ; 1624 v1_lwb = IvlLwb(v1.ptr); 1625 v1_upb = IvlUpb(v1.ptr); 1626 v2_lwb = IvlLwb(v2.ptr); 1627 v2_upb = IvlUpb(v2.ptr); 1628 ec_i_add(v1_lwb, v1_upb, v2_lwb, v2_upb, &lwb, &upb); 1629 Make_Interval(pres, lwb, upb); 1630 Succeed_; 1631} 1632 1633 1634static int 1635_ivl_sub(value v1, value v2, pword *pres) 1636{ 1637 double lwb, upb; 1638 double v1_lwb, v1_upb, v2_lwb, v2_upb ; 1639 v1_lwb = IvlLwb(v1.ptr); 1640 v1_upb = IvlUpb(v1.ptr); 1641 v2_lwb = IvlLwb(v2.ptr); 1642 v2_upb = IvlUpb(v2.ptr); 1643 ec_i_sub(v1_lwb, v1_upb, v2_lwb, v2_upb, &lwb, &upb ); 1644 Make_Interval(pres, lwb, upb); 1645 Succeed_; 1646} 1647 1648 1649static int 1650_ivl_mul(value v1, value v2, pword *pres) 1651{ 1652 double lwb, upb; 1653 double v1_lwb, v1_upb, v2_lwb, v2_upb ; 1654 v1_lwb = IvlLwb(v1.ptr); 1655 v1_upb = IvlUpb(v1.ptr); 1656 v2_lwb = IvlLwb(v2.ptr); 1657 v2_upb = IvlUpb(v2.ptr); 1658 ec_i_mul(v1_lwb, v1_upb, v2_lwb, v2_upb, &lwb, &upb); 1659 Make_Interval(pres, lwb, upb); 1660 Succeed_; 1661} 1662 1663 1664static int 1665_ivl_div(value v1, value v2, pword *pres) 1666{ 1667 double lwb, upb; 1668 double v1_lwb, v1_upb, v2_lwb, v2_upb ; 1669 v1_lwb = IvlLwb(v1.ptr); 1670 v1_upb = IvlUpb(v1.ptr); 1671 v2_lwb = IvlLwb(v2.ptr); 1672 v2_upb = IvlUpb(v2.ptr); 1673 ec_i_div(v1_lwb, v1_upb, v2_lwb, v2_upb, &lwb, &upb); 1674 Make_Interval(pres, lwb, upb); 1675 Succeed_; 1676} 1677 1678 1679static int 1680_ivl_min(value v1, value v2, pword *pres) 1681{ 1682 double t1, t2; 1683 double lwb, upb; 1684 t1 = IvlLwb(v1.ptr); 1685 t2 = IvlLwb(v2.ptr); 1686 lwb = min(t1,t2); 1687 t1 = IvlUpb(v1.ptr); 1688 t2 = IvlUpb(v2.ptr); 1689 upb = min(t1,t2); 1690 Make_Interval(pres, lwb, upb); 1691 Succeed_; 1692} 1693 1694 1695static int 1696_ivl_max(value v1, value v2, pword *pres) 1697{ 1698 double t1, t2; 1699 double lwb, upb; 1700 t1 = IvlLwb(v1.ptr); 1701 t2 = IvlLwb(v2.ptr); 1702 lwb = max(t1,t2); 1703 t1 = IvlUpb(v1.ptr); 1704 t2 = IvlUpb(v2.ptr); 1705 upb = max(t1,t2); 1706 Make_Interval(pres, lwb, upb); 1707 Succeed_; 1708} 1709 1710 1711static int 1712_ivl_neg(value v1, pword *pres) 1713{ 1714 Make_Interval(pres, -IvlUpb(v1.ptr), -IvlLwb(v1.ptr)); 1715 Succeed_; 1716} 1717 1718 1719/* 1720 * This function is only called from the parser and gets passed 1721 * intervals that have been constructed by the lexer. These might 1722 * be proper intervals, in which case they are normally negated. 1723 * 1724 * Otherwise they are "raw intervals". This happens because the lexer 1725 * returns "-1.1__-0.9" as two separate tokens, "-" and "1.1__-0.9", 1726 * The latter is then a "raw interval" which doesn't satisfy lwb=<upb. 1727 * The parser then calls this function to combine the sign and the raw 1728 * interval into a proper interval. This is done by just negating the 1729 * lower bound. The raw-flag is reset in the parser or read_token. 1730 */ 1731static int 1732_ivl_chgsign(value v1, pword *pres) 1733{ 1734 if (!RawInterval(v1.ptr)) 1735 return _ivl_neg(v1, pres); 1736 1737 Make_Interval(pres, -IvlLwb(v1.ptr), IvlUpb(v1.ptr)); 1738 Succeed_; 1739} 1740 1741 1742static int 1743_ivl_abs(value v1, pword *pres) 1744{ 1745 double lwb, upb; 1746 if (IvlLwb(v1.ptr) >= 0.0) { 1747 lwb = IvlLwb(v1.ptr); upb = IvlUpb(v1.ptr); 1748 } else if (IvlUpb(v1.ptr) < 0.0) { 1749 lwb = -IvlUpb(v1.ptr); upb = -IvlLwb(v1.ptr); 1750 } else { 1751 double t1 = IvlUpb(v1.ptr); 1752 double t2 = -IvlLwb(v1.ptr); 1753 lwb = 0.0; upb = max(t1, t2); 1754 } 1755 Make_Interval(pres, lwb, upb); 1756 Succeed_; 1757} 1758 1759 1760static int 1761_ivl_sqrt(value v1, pword *pres) 1762{ 1763 double lwb, upb; 1764 if (IvlLwb(v1.ptr) < 0.0) { 1765 /* the result could be imaginary and is thus not representable */ 1766 Bip_Error(ARITH_EXCEPTION); 1767 } 1768 upb = up(sqrt(IvlUpb(v1.ptr))); 1769 lwb = sqrt(IvlLwb(v1.ptr)); 1770 /* don't go below zero, but preserve negative zeros from sqrt() */ 1771 if (lwb > 0.0) lwb = down(lwb); 1772 Make_Interval(pres, lwb, upb); 1773 Succeed_; 1774} 1775 1776 1777static int 1778_ivl_sin(value v1, pword *pres) 1779{ 1780 double lwb, upb; 1781 i_sin(IvlLwb(v1.ptr), IvlUpb(v1.ptr), &lwb, &upb, 0); 1782 Make_Interval(pres, lwb, upb); 1783 Succeed_; 1784} 1785 1786 1787static int 1788_ivl_cos(value v1, pword *pres) 1789{ 1790 double lwb, upb; 1791 i_sin(IvlLwb(v1.ptr), IvlUpb(v1.ptr), &lwb, &upb, 1); 1792 Make_Interval(pres, lwb, upb); 1793 Succeed_; 1794} 1795 1796static int 1797_ivl_exp(value v1, pword *pres) 1798{ 1799 double lwb, upb; 1800 i_exp(IvlLwb(v1.ptr), IvlUpb(v1.ptr), &lwb, &upb); 1801 Make_Interval(pres, lwb, upb); 1802 Succeed_; 1803} 1804 1805 1806static int 1807_ivl_ln(value v1, pword *pres) 1808{ 1809 double lwb, upb; 1810 lwb = IvlLwb(v1.ptr); 1811 if (lwb < 0.0) { 1812 /* result could be undefined, thus not representable as an interval */ 1813 Bip_Error(ARITH_EXCEPTION); 1814 } 1815 lwb = log(lwb); 1816 if (lwb > 0.0) lwb = down(lwb); /* don't go below zero */ 1817 upb = up(log(IvlUpb(v1.ptr))); 1818 Make_Interval(pres, lwb, upb); 1819 Succeed_; 1820} 1821 1822 1823static int 1824_ivl_atan(value v1, pword *pres) 1825{ 1826 double lwb, upb; 1827 lwb = down(atan(IvlLwb(v1.ptr))); 1828 upb = up(atan(IvlUpb(v1.ptr))); 1829 Make_Interval(pres, lwb, upb); 1830 Succeed_; 1831} 1832 1833static int 1834_ivl_atan2(value vy, value vx, pword *pres) 1835{ 1836 double lwb, upb; 1837 double xl, xu, yl, yu; 1838 yl = IvlLwb(vy.ptr); 1839 yu = IvlUpb(vy.ptr); 1840 xl = IvlLwb(vx.ptr); 1841 xu = IvlUpb(vx.ptr); 1842 if (samesign(xl,-1.0) && !samesign(yl, yu)) { 1843 lwb = -3.1415926535897935; 1844 upb = 3.1415926535897935; 1845 } else { 1846 double t1, t2; 1847 double ll = Atan2(yl,xl); 1848 double lu = Atan2(yl,xu); 1849 double ul = Atan2(yu,xl); 1850 double uu = Atan2(yu,xu); 1851 t1 = min(ll, lu); 1852 t2 = min(ul, uu); 1853 lwb = down(min(t1, t2)); 1854 t1 = max(ll, lu); 1855 t2 = max(ul, uu); 1856 upb = up(max(t1, t2)); 1857 } 1858 Make_Interval(pres, lwb, upb); 1859 Succeed_; 1860} 1861 1862 1863/* TODO: rewrite this with set_round_xxx */ 1864 1865static int 1866_ivl_pow(value v1, value v2, pword *pres) 1867{ 1868 double xl, xu, yl, yu; 1869 double lwb, upb, yi; 1870 1871 xl = IvlLwb(v1.ptr); 1872 xu = IvlUpb(v1.ptr); 1873 yl = IvlLwb(v2.ptr); 1874 yu = IvlUpb(v2.ptr); 1875 1876 if (xl >= 0.0) { 1877 /* base is nonnegative, result as well */ 1878 if (yl > 0.0) { 1879 /* increasing */ 1880 lwb = Pow(xl, xl>1? yl: yu); 1881 upb = Pow(xu, xu>1? yu: yl); 1882 } else if (yu < 0.0) { 1883 /* decreasing */ 1884 lwb = SafePow(xu, xu>1? yl: yu); 1885 upb = SafePow(xl, xl>1? yu: yl); 1886 } else { /* Y-range includes 0 */ 1887 double t1, t2; 1888 double ll = SafePow(xl,yl); 1889 double lu = Pow(xl,yu); 1890 double ul = SafePow(xu,yl); 1891 double uu = Pow(xu,yu); 1892 t1 = min(ll, lu); 1893 t2 = min(ul, uu); 1894 lwb = min(t1, t2); 1895 t1 = max(ll, lu); 1896 t2 = max(ul, uu); 1897 upb = max(t1, t2); 1898 } 1899 if (lwb > 0.0) lwb = down(lwb); /* don't go below 0.0 */ 1900 upb = up(upb); 1901 1902 } else if (IsExactIntIvl(v2.ptr, &yi)) { /* xl < 0.0 */ 1903 /* base may be negative, exponent must be integral */ 1904 double l = Pow(xl,yi); 1905 double u = SafePow(xu,yi); 1906 lwb = down(min(l, u)); 1907 upb = up(max(l, u)); 1908 1909 if (xu >= 0.0) { /* X-range includes 0 */ 1910 l = SafePow(NEGZERO,yi); 1911 u = SafePow(0.0,yi); 1912 /* l,u are +-0, 1, or +-inf, don't round */ 1913 lwb = min(lwb,l); 1914 upb = max(upb,l); 1915 lwb = min(lwb,u); 1916 upb = max(upb,u); 1917 } 1918 1919 } else { 1920 /* If the base lower bound is negative the result could be complex*/ 1921 Bip_Error(ARITH_EXCEPTION); 1922 } 1923 Make_Interval(pres, lwb, upb); 1924 Succeed_; 1925} 1926 1927 1928static int 1929_ivl_floor(value v1, pword *pres) 1930{ 1931 double lwb, upb; 1932 lwb = floor(IvlLwb(v1.ptr)); /* integers, no inaccuracy */ 1933 upb = floor(IvlUpb(v1.ptr)); 1934 Make_Interval(pres, lwb, upb); 1935 Succeed_; 1936} 1937 1938 1939static int 1940_ivl_ceil(value v1, pword *pres) 1941{ 1942 double lwb, upb; 1943 lwb = Ceil(IvlLwb(v1.ptr)); /* integers, no inaccuracy */ 1944 upb = Ceil(IvlUpb(v1.ptr)); 1945 Make_Interval(pres, lwb, upb); 1946 Succeed_; 1947} 1948 1949 1950static int 1951_ivl_truncate(value v1, pword *pres) 1952{ 1953 double lwb, upb; 1954#ifdef HAVE_TRUNC 1955 lwb = trunc(IvlLwb(v1.ptr)); /* integers, no inaccuracy */ 1956 upb = trunc(IvlUpb(v1.ptr)); 1957#else 1958 lwb = IvlLwb(v1.ptr); 1959 lwb = lwb < 0 ? Ceil(lwb) : floor(lwb); 1960 upb = IvlUpb(v1.ptr); 1961 upb = upb < 0 ? Ceil(upb) : floor(upb); 1962#endif 1963 Make_Interval(pres, lwb, upb); 1964 Succeed_; 1965} 1966 1967 1968static int 1969_ivl_sgn(value v1, pword *pres) 1970{ 1971 int res; 1972 if (IvlLwb(v1.ptr) > 0.0) 1973 res = 1; 1974 else if (IvlUpb(v1.ptr) < 0.0) 1975 res = -1; 1976 else if (IvlLwb(v1.ptr) == 0.0 && IvlUpb(v1.ptr) == 0.0) 1977 res = 0; 1978 else { Bip_Error(ARITH_EXCEPTION); } 1979 Make_Integer(pres, res); 1980 Succeed_; 1981} 1982 1983 1984static int 1985_ivl_round(value v1, pword *pres) 1986{ 1987 double lwb, upb; 1988#if defined(HAVE_RINT) && !defined(HP_RINT) && !defined(HAVE_RINT_BUG) 1989 lwb = rint(IvlLwb(v1.ptr)); 1990 upb = rint(IvlUpb(v1.ptr)); 1991#else 1992 /* 1993 * Round to even number if we are exactly in the middle. 1994 * Make sure we round to -0.0 if between -0.5 and -0.0 1995 */ 1996 lwb = Ceil(IvlLwb(v1.ptr)); 1997 if (lwb - IvlLwb(v1.ptr) > 0.5 || (lwb - IvlLwb(v1.ptr) == 0.5 && ((word)lwb & 1))) 1998 lwb -= 1.0; 1999 upb = Ceil(IvlUpb(v1.ptr)); 2000 if (upb - IvlUpb(v1.ptr) > 0.5 || (upb - IvlUpb(v1.ptr) == 0.5 && ((word)upb & 1))) 2001 upb -= 1.0; 2002#endif /* rint */ 2003 Make_Interval(pres, lwb, upb); 2004 Succeed_; 2005} 2006 2007 2008static int 2009_ivl_int2(value v1, pword *pres) 2010{ 2011 double ipart; 2012 if (IsExactIntIvl(v1.ptr, &ipart)) 2013 { 2014 value dval; 2015 Make_Double_Val(dval, ipart); 2016 return tag_desc[TDBL].arith_op[ARITH_FIX](dval, pres); 2017 } 2018 else 2019 { 2020 Bip_Error(ARITH_EXCEPTION); 2021 } 2022} 2023 2024 2025/* 2026 * This is term comparison, which must be totally defined. 2027 * Much of this ordering is arbitrary, but the following invariants hold: 2028 * 2029 * Where ==/2 succeeds 2030 * term compare must give = (compatible with _equal_ivl()) 2031 * Where arith compare yields < or > 2032 * term compare should give the same. 2033 * Where arith compare yields = 2034 * term compare will give = 2035 * except for different zeros, where it can give <, > or = 2036 * Where arith compare delays (or ==/2 is undecidable) 2037 * term compare can give >, < or = 2038 */ 2039 2040static int 2041_compare_ivl(value v1, value v2) 2042{ 2043 double lwb1, upb1; 2044 double lwb2, upb2; 2045 2046 if (v1.ptr == v2.ptr) 2047 return 0; 2048 2049 if ((lwb1 = IvlLwb(v1.ptr)) > (upb2 = IvlUpb(v2.ptr))) 2050 return 1; /* unambiguous */ 2051 2052 if ((upb1 = IvlUpb(v1.ptr)) < (lwb2 = IvlLwb(v2.ptr))) 2053 return -1; /* unambiguous */ 2054 2055 /* overlap, but possibly just at -0.0,0.0 */ 2056 return PedanticGreater(lwb1,lwb2) ? 1 /* arbitrary (overlap) */ 2057 : PedanticLess(lwb1,lwb2) ? -1 /* arbitrary (overlap) */ 2058 : PedanticGreater(upb1,upb2) ? 1 /* arbitrary (overlap) */ 2059 : PedanticLess(upb1,upb2) ? -1 /* arbitrary (overlap) */ 2060 /* both lower bounds and both upper bounds are exactly equal */ 2061 : !(GlobalFlags & BREAL_EXCEPTIONS) ? 0 2062 : lwb1 == upb2 ? 0 /* zero width */ 2063 : v1.ptr > v2.ptr ? 1 : -1; /* arbitrary (same bounds) */ 2064} 2065 2066 2067/* 2068 * This is arithmetic comparison, which may delay if it is not decidable. 2069 * Needs the kind of comparison in *relation (BILe, etc) as input! 2070 * Returns PDELAY if undecidable, otherwise PSUCCEED and -1,0,1 in *relation. 2071 * The exact meaning of the result depends on the kind of comparison: 2072 * 2073 * Lt,Ge Le,Gt Eq,Ne LeGe 2074 * -1 < =< < =< 2075 * 0 = = = = 2076 * 1 >= > > >= 2077 */ 2078 2079static int 2080_arith_compare_ivl(value v1, value v2, int *relation) 2081{ 2082 /* We want identical terms to compare equal regardless of anything else. */ 2083 if (v1.ptr == v2.ptr) { 2084 *relation = 0; 2085 return PSUCCEED; 2086 } 2087 2088 /* 2089 * Note that in some of the cases below, we might return 1 or -1 when 2090 * the values could actually be equal. This works out OK though because 2091 * the test in the calling code will give the right answer. 2092 */ 2093 switch(*relation) { 2094 case BILe: /* -1 (=<), 1 (>) or undecidable */ 2095 case BIGt: 2096 if (IvlUpb(v1.ptr) <= IvlLwb(v2.ptr)) { 2097 *relation = -1; /* Means -1 or 0 */ 2098 return PSUCCEED; 2099 } else if (IvlLwb(v1.ptr) > IvlUpb(v2.ptr)) { 2100 *relation = 1; 2101 return PSUCCEED; 2102 } else { 2103 *relation = 0; 2104 return PDELAY; 2105 } 2106 break; 2107 2108 case BIGe: /* -1 (<), 1 (>=) or undecidable */ 2109 case BILt: 2110 if (IvlLwb(v1.ptr) >= IvlUpb(v2.ptr)) { 2111 *relation = 1; /* Means 1 or 0 */ 2112 return PSUCCEED; 2113 } else if (IvlUpb(v1.ptr) < IvlLwb(v2.ptr)) { 2114 *relation = -1; 2115 return PSUCCEED; 2116 } else { 2117 *relation = 0; 2118 return PDELAY; 2119 } 2120 break; 2121 2122 case BIEq: /* -1 (<), 0 (=), 1 (>) or undecidable */ 2123 case BINe: 2124 if (IvlUpb(v1.ptr) < IvlLwb(v2.ptr)) { 2125 *relation = -1; 2126 return PSUCCEED; 2127 } else if (IvlLwb(v1.ptr) > IvlUpb(v2.ptr)) { 2128 *relation = 1; 2129 return PSUCCEED; 2130 } else if (IvlLwb(v1.ptr) == IvlUpb(v2.ptr) && 2131 IvlUpb(v1.ptr) == IvlLwb(v2.ptr)) { 2132 *relation = 0; 2133 return PSUCCEED; 2134 } else { 2135 *relation = 0; 2136 return PDELAY; 2137 } 2138 break; 2139 2140 case BILeGe: /* -1 (=<), 0 (=), 1 (>=) or undecidable */ 2141 if (IvlUpb(v1.ptr) <= IvlLwb(v2.ptr)) { 2142 if (IvlLwb(v1.ptr) >= IvlUpb(v2.ptr)) 2143 *relation = 0; 2144 else 2145 *relation = -1; /* Means -1 or 0 */ 2146 return PSUCCEED; 2147 } else if (IvlLwb(v1.ptr) >= IvlUpb(v2.ptr)) { 2148 if (IvlUpb(v1.ptr) <= IvlLwb(v2.ptr)) 2149 *relation = 0; 2150 else 2151 *relation = 1; /* Means 1 or 0 */ 2152 return PSUCCEED; 2153 } else { 2154 *relation = 0; 2155 return PDELAY; 2156 } 2157 break; 2158 } 2159 assert(0); 2160} 2161 2162 2163static int 2164_equal_ivl(pword *pw1, pword *pw2) 2165{ 2166 if (pw1 == pw2) 2167 return 1; /* pointers identical */ 2168 2169 if (GlobalFlags & BREAL_EXCEPTIONS) 2170 { 2171 if (IvlUpb(pw1) < IvlLwb(pw2) || IvlUpb(pw2) < IvlLwb(pw1)) 2172 return 0; /* no overlap */ 2173 2174 /* overlap: check for zero width */ 2175 if (IvlLwb(pw1) == IvlUpb(pw2) && IvlUpb(pw1) == IvlLwb(pw2)) 2176 { 2177 if (IvlLwb(pw1) != 0) 2178 return 1; /* nonzero and zero width: identical */ 2179 2180 if (PedanticZeroEq(IvlLwb(pw1), IvlLwb(pw2)) 2181 && PedanticZeroEq(IvlUpb(pw1), IvlUpb(pw2))) 2182 return 1; /* identical zeros */ 2183 } 2184 2185 /* nonzero width and overlap: undecidable */ 2186 { 2187 value v; 2188 v.did = d_undecidable; 2189 Exit_Block(v, tdict); 2190 } 2191 } 2192 else 2193 { 2194#if SIZEOF_DOUBLE < SIZEOF_WORD || SIZEOF_DOUBLE > 2*SIZEOF_WORD 2195 PROBLEM: Cannot deal with word size SIZEOF_WORD. 2196#else 2197 return ((uword*)BufferStart(pw1))[0] == ((uword*)BufferStart(pw2))[0] 2198 && ((uword*)BufferStart(pw1))[1] == ((uword*)BufferStart(pw2))[1] 2199#if SIZEOF_DOUBLE > SIZEOF_WORD 2200 && ((uword*)BufferStart(pw1))[2] == ((uword*)BufferStart(pw2))[2] 2201 && ((uword*)BufferStart(pw1))[3] == ((uword*)BufferStart(pw2))[3] 2202#endif 2203#endif 2204 ; 2205 } 2206} 2207 2208 2209/*ARGSUSED*/ 2210static int 2211_unimp_err(void) 2212{ 2213 return UNIMPLEMENTED; 2214} 2215 2216 2217/*-------------------------------------------------------------------------- 2218 * Initialize intervals 2219 *--------------------------------------------------------------------------*/ 2220 2221void 2222ec_intervals_init(void) 2223{ 2224 2225#ifdef SAFE_ROUNDING 2226#ifdef IEEE_ROUND_DOWN 2227 char *res; 2228 ieee_flags("set", "direction", "negative", &res); 2229#endif 2230#ifdef IEEE_INEXACT 2231 ieee_handler("set", "inexact", inexact_handler); 2232#endif 2233#endif 2234 2235 tag_desc[TIVL].tag_name = in_dict("breal", 0); 2236 tag_desc[TIVL].type_name = in_dict("breal", 0); 2237 2238 tag_desc[TINT].coerce_to[TIVL] = _int_ivl; /* coerce to interval */ 2239 tag_desc[TDBL].coerce_to[TIVL] = _dbl_ivl; 2240 tag_desc[TBIG].coerce_to[TIVL] = _big_ivl; 2241 2242 tag_desc[TIVL].coerce_to[TDBL] = _ivl_dbl; /* coerce from interval */ 2243 2244 tag_desc[TIVL].string_size = _ivl_string_size; 2245 tag_desc[TIVL].to_string = _ivl_to_string; 2246 tag_desc[TIVL].write = _write_ivl; 2247 tag_desc[TIVL].compare = _compare_ivl; 2248 tag_desc[TIVL].arith_compare = _arith_compare_ivl; 2249 tag_desc[TIVL].equal = _equal_ivl; 2250 2251 tag_desc[TIVL].arith_op[ARITH_PLUS] = _ivl_nop; 2252 tag_desc[TIVL].arith_op[ARITH_NEG] = _ivl_neg; 2253 tag_desc[TIVL].arith_op[ARITH_ABS] = _ivl_abs; 2254 tag_desc[TIVL].arith_op[ARITH_ADD] = _ivl_add; 2255 tag_desc[TIVL].arith_op[ARITH_SUB] = _ivl_sub; 2256 tag_desc[TIVL].arith_op[ARITH_MUL] = _ivl_mul; 2257 tag_desc[TIVL].arith_op[ARITH_DIV] = _ivl_div; 2258 tag_desc[TIVL].arith_op[ARITH_MIN] = _ivl_min; 2259 tag_desc[TIVL].arith_op[ARITH_MAX] = _ivl_max; 2260 tag_desc[TIVL].arith_op[ARITH_SIN] = _ivl_sin; 2261 tag_desc[TIVL].arith_op[ARITH_COS] = _ivl_cos; 2262 tag_desc[TIVL].arith_op[ARITH_TAN] = _unimp_err; 2263 tag_desc[TIVL].arith_op[ARITH_EXP] = _ivl_exp; 2264 tag_desc[TIVL].arith_op[ARITH_LN] = _ivl_ln; 2265 tag_desc[TIVL].arith_op[ARITH_ASIN] = _unimp_err; 2266 tag_desc[TIVL].arith_op[ARITH_ACOS] = _unimp_err; 2267 tag_desc[TIVL].arith_op[ARITH_ATAN] = _ivl_atan; 2268 tag_desc[TIVL].arith_op[ARITH_ATAN2] = _ivl_atan2; 2269 tag_desc[TIVL].arith_op[ARITH_SQRT] = _ivl_sqrt; 2270 tag_desc[TIVL].arith_op[ARITH_POW] = _ivl_pow; 2271 tag_desc[TIVL].arith_op[ARITH_FLOOR] = _ivl_floor; 2272 tag_desc[TIVL].arith_op[ARITH_CEIL] = _ivl_ceil; 2273 tag_desc[TIVL].arith_op[ARITH_TRUNCATE] = _ivl_truncate; 2274 tag_desc[TIVL].arith_op[ARITH_ROUND] = _ivl_round; 2275 tag_desc[TIVL].arith_op[ARITH_SGN] = _ivl_sgn; 2276 tag_desc[TIVL].arith_op[ARITH_NICERAT] = _unimp_err; 2277 tag_desc[TIVL].arith_op[ARITH_INT] = _ivl_int2; 2278 2279 tag_desc[TIVL].arith_op[ARITH_CHGSIGN] = _ivl_chgsign; 2280 2281 tag_desc[TIVL].from_string = _ivl_from_string; 2282 2283 built_in(in_dict("ria_unop", 5), p_ria_unop, B_UNSAFE|U_GROUND) 2284 -> mode = BoundArg(4, CONSTANT) | BoundArg(5, CONSTANT); 2285 built_in(in_dict("ria_binop", 7), p_ria_binop, B_UNSAFE|U_GROUND) 2286 -> mode = BoundArg(6, CONSTANT) | BoundArg(7, CONSTANT); 2287 (void) built_in(in_dict("ria_ternop", 9), p_ria_ternop, B_UNSAFE|U_GROUND); 2288 /* no space in mode mask to store this information (arity too high): */ 2289 /* -> mode = BoundArg(8, CONSTANT) | BoundArg(9, CONSTANT); */ 2290 (void) exported_built_in(in_dict("breal_from_bounds", 3), p_breal_from_bounds, B_UNSAFE|U_SIMPLE); 2291 (void) exported_built_in(in_dict("breal_min", 2), p_breal_min, B_UNSAFE|U_SIMPLE); 2292 (void) exported_built_in(in_dict("breal_max", 2), p_breal_max, B_UNSAFE|U_SIMPLE); 2293 (void) exported_built_in(in_dict("breal_bounds", 3), p_breal_bounds, B_UNSAFE|U_SIMPLE); 2294 2295 d_undecidable = in_dict("undecidable comparison of bounded reals", 0); 2296} 2297