From: ··········@YahooGroups.Com
Subject: I'm not the only one who values interval-arithmetic
Date: 
Message-ID: <REM-2003sep12-002@Yahoo.Com>
I've been advocating for interval arithmetic ever since the late 1960's
when I implemented a small bit of it on the IBM 1620. Recently I've
been advocating it for inclusion in Common Lisp. Yesterday while
catching up with back issues of SCIENCE, I found an article on the
computerized proof of Kepler's conjecture (now theorem) on optimal
close-packing of spheres, issue of 2003.Aug.29, page 1186. Here's the
key quote: "The numerical computation was all made mathematically
rigorous by using 'interval arithmetic' rather than the much faster,
but nonrigorous, floating-point arithmetic."

By the way, there's nothing wrong with floating-point representation of
the endpoints of an interval. So the quote is slightly misleading. What
they mean as fast but nonrigorous is the standard best-guess
(shot-in-dark) floating-point calculations which are built into most
modern CPUs. I'm curious what representation was used for the endpoints
of the intervals in the Kepler's-conjecture proof (fixed-point, exact
rationals, or in fact scientific notation i.e. floating-point after
all). I'm also curious what programming language they used, and whether
they used a standard interval-arithmetic package or one they
handcrafted specifically for the calculations for this one proof.

I continue to wish that CPU designers would implement high-speed
fixed-point and floating-point interval arithmetic, either directly or
via at least implementing high-speed floor/ceiling results which can
then be combined to yield the desired intervals. And regardless of any
hardware support for high speed, I wish programming languages such as
Lisp provided interval arithmetic as a data type and a package of
arithmetic and trigometric functions that used them. Has any Lisp
implementation expert created such a package, presumably more efficient
than anything I could hack together myself, and made it generally
available for use in CMUCL and other CLs?

From: Richard Fateman
Subject: intervals in Lisp, also extending rationals to include infinity
Date: 
Message-ID: <bjtfa5$1p52$1@agate.berkeley.edu>
Your belief that the current CPUs provide nonrigorous
floating point answers might have been correct in 1960, but
almost all CPUs sold today that do floating point, do it using
IEEE 754 binary floating point standard rules. These specify
all the bits in single and double formats, provide exactly the
rounding modes needed to do endpoint calculations, and also
provide representations for Not a Number and (signed) Infinities.
Also signed zeros, if you have an interest in them.

There are several ways of defining intervals including using
rational endpoints (probably a bad idea) and using two different
models of the real line, affine or projective.

Thus there is not just one prescription for what common lisp should do
to include intervals.  I am not familiar with what SBC / CMU
includes.  A selection of routines is in a directory
http://www.cs.berkeley.edu/~fateman/intervals
including some documentation, and a draft of a paper on the topic.

The work is somewhat sketchy, but includes, for example a sine routine.
This work was done circa 1991, and was intended for use as an lternative
arithmetic underneath a computer algebra system, "MockMMA" which has the 
same external syntax as Mathematica (as of version 2.0 or so).

There is also an extension to CL rationals which actually may make
implementations faster which allows you to divide by zero.
Most people would agree that  (+ 2  (/ 1  (/ 1 0)))   ie.  2 + 
1/infinity,  is 2.

So there is no reason to blow up at the (/ 1 0).



RJF




··········@YahooGroups.Com wrote:
> I've been advocating for interval arithmetic ever since the late 1960's
> when I implemented a small bit of it on the IBM 1620. Recently I've
> been advocating it for inclusion in Common Lisp. Yesterday while
> catching up with back issues of SCIENCE, I found an article on the
> computerized proof of Kepler's conjecture (now theorem) on optimal
> close-packing of spheres, issue of 2003.Aug.29, page 1186. Here's the
> key quote: "The numerical computation was all made mathematically
> rigorous by using 'interval arithmetic' rather than the much faster,
> but nonrigorous, floating-point arithmetic."
> 
> By the way, there's nothing wrong with floating-point representation of
> the endpoints of an interval. So the quote is slightly misleading. What
> they mean as fast but nonrigorous is the standard best-guess
> (shot-in-dark) floating-point calculations which are built into most
> modern CPUs. I'm curious what representation was used for the endpoints
> of the intervals in the Kepler's-conjecture proof (fixed-point, exact
> rationals, or in fact scientific notation i.e. floating-point after
> all). I'm also curious what programming language they used, and whether
> they used a standard interval-arithmetic package or one they
> handcrafted specifically for the calculations for this one proof.
> 
> I continue to wish that CPU designers would implement high-speed
> fixed-point and floating-point interval arithmetic, either directly or
> via at least implementing high-speed floor/ceiling results which can
> then be combined to yield the desired intervals. And regardless of any
> hardware support for high speed, I wish programming languages such as
> Lisp provided interval arithmetic as a data type and a package of
> arithmetic and trigometric functions that used them. Has any Lisp
> implementation expert created such a package, presumably more efficient
> than anything I could hack together myself, and made it generally
> available for use in CMUCL and other CLs?
From: Pascal Bourguignon
Subject: Re: intervals in Lisp, also extending rationals to include infinity
Date: 
Message-ID: <87llstmke8.fsf@thalassa.informatimago.com>
Richard Fateman <·······@cs.berkeley.edu> writes:
> There is also an extension to CL rationals which actually may make
> implementations faster which allows you to divide by zero.
> Most people would agree that  (+ 2  (/ 1  (/ 1 0)))   ie.  2 +
> 1/infinity,  is 2.
>
> So there is no reason to blow up at the (/ 1 0).

Perhaps not at (/ 1 0), but what about (/ 0 0)
or even (/ 0.000000000000000000000000000001 0)?

Do  you consider the  numbers in  computers are  fixed numbers,  or as
parts of functions?  

In math when we write 1/0=infinity, we are actually

writting:      lim  x/y = infinity
            x->1,y->0

or more precisely:  for all z in R, there is an epsilon in R, such as
                    for all x, 1<x<1+epsilon and 
                    for all y, 0<y<0+epsilon, x/y>z

So, when we write 1/0, we are actually denoting functions whose limits
are 1 and 0 (that is, functions with some specific properties).

The division operation is merely not defined for 0 as divisor.

In lisp, you can always rather write: 

    (lim (x 1) (y 0) (/ 2 (/ 1 (/ (lambda (x) x) (lambda (y) y)))))

and compute on that...

if you wanted to add an :infinity value to be able to compute (/ 1 0),
you would as well need an  infinity of such infinity values!  Think of
(+ 2 (/ x (/ 1 0))) with x -> infinity [or (+ 2 (/ :infinity (/ 1 0)))
vs. (+ 2 (/ (* 2 :infinity) (/ 1 0)))].


--
__Pascal_Bourguignon__
http://www.informatimago.com/
Do not adjust your mind, there is a fault in reality.
From: Wolfhard Buß
Subject: Re: intervals in Lisp, also extending rationals to include infinity
Date: 
Message-ID: <m3u17dhngu.fsf@buss-14250.user.cis.dfn.de>
Pascal Bourguignon:

>             lim  x/y = infinity
>             x->1,y->0
> 
> or more precisely:  for all z in R, there is an epsilon in R, such as
>                     for all x, 1<x<1+epsilon and 
>                     for all y, 0<y<0+epsilon, x/y>z

Modulo two typos.


-- 
"Hurry if you still want to see something. Everything is vanishing."
                                       --  Paul C�zanne (1839-1906)
From: Gorbag
Subject: Re: I'm not the only one who values interval-arithmetic
Date: 
Message-ID: <Aep8b.4061$c6.3386@bos-service2.ext.raytheon.com>
<··········@YahooGroups.Com> wrote in message
······················@Yahoo.Com...
> I've been advocating for interval arithmetic ever since the late 1960's
> when I implemented a small bit of it on the IBM 1620.
>...Has any Lisp
> implementation expert created such a package, presumably more efficient
> than anything I could hack together myself, and made it generally
> available for use in CMUCL and other CLs?

May I suggest you go ahead and hack something together yourself and see if
others start to use it? The popularity will drive more efficient
implementation.

Water's series package is a nice idea but it never caught on, thus it never
was implemented very efficiently.

PCL on the other hand DID catch on, and became CLOS.

The key is to implement and publish. If it's useful, the efficient
implementations will follow. Things that aren't conceptually useful won't
catch on to begin with (well, I suppose C++ might be a counterexample, but
here too, the inefficient form was published first and eventually more
efficient implementations were created).
From: Kent M Pitman
Subject: Re: I'm not the only one who values interval-arithmetic
Date: 
Message-ID: <sfwpti5zsve.fsf@shell01.TheWorld.com>
"Gorbag" <·············@nospam.mac.com> writes:

> <··········@YahooGroups.Com> wrote in message
> ······················@Yahoo.Com...
> > I've been advocating for interval arithmetic ever since the late 1960's
> > when I implemented a small bit of it on the IBM 1620.
> >...Has any Lisp
> > implementation expert created such a package, presumably more efficient
> > than anything I could hack together myself, and made it generally
> > available for use in CMUCL and other CLs?
> 
> May I suggest you go ahead and hack something together yourself and see if
> others start to use it? The popularity will drive more efficient
> implementation.
> 
> Water's series package is a nice idea but it never caught on, thus it never
> was implemented very efficiently.
> 
> PCL on the other hand DID catch on, and became CLOS.
> 
> The key is to implement and publish. If it's useful, the efficient
> implementations will follow. Things that aren't conceptually useful won't
> catch on to begin with (well, I suppose C++ might be a counterexample, but
> here too, the inefficient form was published first and eventually more
> efficient implementations were created).

[Mostly replying to Robert here, but also speaking somewhat generally.]

This is a good example of an issue illustrating the important shift that
has de facto occurred in the CL community over the past decade and that
has to be worked into our culture.

It used to be that people waited for a new rev of the standard to tell
them what tools are being used.  The "New Way" is for you to _first_ be 
successful with something and _then_ lobby for its acceptance by the 
community based on that success.

It was part of the early bootstrapping of languages that we just tried
to think up a set of base tools that had enough critical mass and
internal consistency that it would serve as a seed to get the
community started.  This was a gamble, but a necessary gamble.

But it is a practical fact that no central agency can be smart enough
to anticipate what the community will need on an ongoing basis, and to
pretend otherwise is to condemn the community to wait until a huge,
expensive, global synchronization occurs.  And it is no longer necessary
to make a gamble.  You have tools that are "good enough" that you can 
demonstrate need without just standing around saying "I need help".

Even just saying "I want interval arithmetic" is highly ill-specified.
It's probably the case that it really matters what the details of this
are, in terms of what the internal representation is, how printing works,
etc.  And that will probably require experimentation that is best, if
you take my point, not being done with "concrete", the only available
mechanism for standards.  That is, standards are things that are cast in
concrete for the long haul.  But experimentation such as is needed to 
get things right the first time should be done with lighter weight
materials.

What is required, IMO, is not a "lisp implementation expert" but rather
a "numerics expert".  

As an aside, it occurs to me that this would have been a good thing to
have in implementing BASIC.  I didn't realize until after my stint
programming in Basic that all apparent integers in BASIC are really
floating point, and that arithmetic to compute indices yields something
that has to be rounded to get an integer index (ick).  Most computations
don't widen the interval large enough for it to matter, but it might have
been useful to know if they did so that one didn't get the wrong index.
In THAT case, making it built into the language might have helped the
language, and also, there was no precise number format.  At least in
CL, you can use exact integers, conses, and maintain the decimal point
separately... at least until you figure out how to articulate why that's
not adequate, since it has the sound of being adequate to me.
I'm pretty sure this is what Macsyma did for bigfloats...
Maybe examining the way bigfloats worked in DOE Macsyma, which I 
understand to be available as open sourceware of some kind, would help
you get some ideas.
From: Rob Warnock
Subject: Re: I'm not the only one who values interval-arithmetic
Date: 
Message-ID: <l8ydnZsijJZzlP6iXTWc-w@speakeasy.net>
Kent M Pitman  <······@world.std.com> wrote:
+---------------
| It used to be that people waited for a new rev of the standard to tell
| them what tools are being used.  The "New Way" is for you to _first_ be 
| successful with something and _then_ lobby for its acceptance by the 
| community based on that success.
+---------------

Actually, that's the old, old way, too. It used to be the case that
standards were defined by codifying "best practice" (modulo some
compromise if there were multiple extant "best" practices), but then...

<FLAME>
Standards committees [not Lisp, particularly, but in general] started
going wild thinking that *they* were the architects of the future
(not the people who actually knew how stuff worked!) and we went into
a long period where little new happened until *after* some standards
committee had already frozen the new thing beyond the point of redeeming
any fatal flaws (e.g., FDDI, ATM, etc.).
</FLAME>

Fortunately, that seems to be changing somewhat, and we're fortunately
reverting (at least somewhat) to the old, old way of "demonstrate one
or more working implementations first" [which is the same as what Kent
calls the "New Way"].


-Rob

-----
Rob Warnock, PP-ASEL-IA		<····@rpw3.org>
627 26th Avenue			<URL:http://rpw3.org/>
San Mateo, CA 94403		(650)572-2607
From: ·············@comcast.net
Subject: Re: I'm not the only one who values interval-arithmetic
Date: 
Message-ID: <vfrwe35d.fsf@comcast.net>
"Gorbag" <·············@nospam.mac.com> writes:

> <··········@YahooGroups.Com> wrote in message
> ······················@Yahoo.Com...
>> I've been advocating for interval arithmetic ever since the late 1960's
>> when I implemented a small bit of it on the IBM 1620.
>>...Has any Lisp
>> implementation expert created such a package, presumably more efficient
>> than anything I could hack together myself, and made it generally
>> available for use in CMUCL and other CLs?
>
> May I suggest you go ahead and hack something together yourself and see if
> others start to use it? The popularity will drive more efficient
> implementation.
>
> Water's series package is a nice idea but it never caught on, thus it never
> was implemented very efficiently.

Huh?  It seems to be quite efficient.  I think the reason it never
caught on is because it has some non-intuitive quirks.  It takes a bit
of practice to use it.
From: David Golden
Subject: SERIES
Date: 
Message-ID: <235b265c.0309140234.4d68b9e3@posting.google.com>
·············@comcast.net wrote in message news:<············@comcast.net>...

> > Water's series package is a nice idea but it never caught on, thus it never
> > was implemented very efficiently.
> 
> Huh?  It seems to be quite efficient.  I think the reason it never
> caught on is because it has some non-intuitive quirks.  It takes a bit
> of practice to use it.

Perhaps one reason it hasn't caught on is because you can't currently
download tarballs of the thing from the obvious place to get them -
the http://series.sourceforge.net/ front page that most people link to has
a broken download link, you have to be sufficiently familiar
with sourceforge to know that you can go to 
http://sourceforge.net/projects/series/ to get to a project's files.

Don't underestimate how discouraging it is for the casually interested
when the big, obvious, download link doesn't work...

I've just emailed the series-bugs list pointing this out.
From: ··········@YahooGroups.Com
Subject: Re: I'm not the only one who values interval-arithmetic
Date: 
Message-ID: <REM-2004apr18-001@Yahoo.Com>
> Date: Fri, 12 Sep 2003 15:35:59 -0400
> From: "Gorbag" <·············@nospam.mac.com>
> May I suggest you go ahead and hack something together yourself and
> see if others start to use it? The popularity will drive more
> efficient implementation.

Unfortunately hardly anyone has been willing to even look at most of
the software I've written, and as I've gotten older I no longer have
the emotional energy to spend thousands of hours doing work that nobody
pays me for and nobody appreciates and nobody even looks at. Also I
couldn't decide how to juggle the various kinds of floating-point
arithmetic etc. to make a decent interval-arithmetic package ("package"
in the general software sense, not the CL package (TM) usage of the
word) or module ("module" in the general software sense, not in the
mathematics sense of a set with operations on it etc.), and I didn't
have the emotional energy to decide without any encouragement from
anone providing money or companionship. So I haven't worked on it
hardly at all since September.

But yesterday as I was lying in bed I got the urge to implement the
following function:
;Given two non-negative rational numbers, find the rational number
; between them (on the closed interval) which has the smallest
; numerator+denominator.
using this method:
;Generate continued-fractions of them in parallel so long as they
; are the same, then when the first difference is found, truncate
; or reduce a term to create a c.f. that is shortest while remaining
; within the interval.
or as the comment in my code says in different words:
;Given two ratios, with 0 < rlo <= rhi, compute continued fractions
; of each until there's a difference or until remainder is zero
; i.e. cf was exact integer, and return the simplest
; cf between them. cf = (n1 n2 ... ny nz) = nz + 1/(ny + 1/(...))
; which could be one or the other endpoint. if rlo=rhi this just
; computes the exact cf at twice the usual cost.
(defun ratios-to-cf-between (rlo rhi) ...)
then re-build a rational number from the resultant c.f.:
(defun ratio-interval-to-simplest-ratio-within (rats)
  (let* ((midcf (ratios-to-cf-between (car rats) (cadr rats)))
         (midrat (cf-to-ratio midcf)))
    midrat))

I imagine this is the same algorithm used inside the RATIONALIZE
function of Common LISP, where it passes the upper and lower bound of
the interval emulated by a single floating-point number. For example,
comparing RATIONALIZE with my new code:
* (setq x (sqrt 2))
1.4142135
* (setq xrat1 (rationalize x))
4756/3363
* (integer-decode-float x)
11863283
-23
1
* (setq lo (* (expt 2 -23) (- 11863283 1/2)))
23726565/16777216
* (setq hi (* (expt 2 -23) (+ 11863283 1/2)))
23726567/16777216
* (setq xrat2 (ratio-intervals-to-simplest-ratio-within (list lo hi)))
4756/3363
* (= xrat1 xrat2)
T
Yup, my code produces exactly the same result as RATIONALIZE when given
appropriate input.

So anyway, I decided not to use floating-point numbers at all, but
rather to use rational numbers consistently, so I can implement
interval-arithmetic once and for all to cover *all* levels of
precision, instead of having to switch from float to long float to
bigfloat and deal with IEEE rounding modes and whatever other mess.
Once I got started it was difficult to stop, so yesterday by bedtime,
and continuing today to noonish, I had written 26 new functions (and
now it's four hour later during which time I collected the demos
described below and composed this article to post), including the three
mentionned above and also these of interest:
(defun ratio-make-smaller+simpler (rat) ...)
(defun ratio-make-larger+simpler (rat) ...)
(defun ratio-to-simpler-interval (rat) ...)
(defun ratio-interval-make-wider+simpler (rats) ...)
(defun ratio-interval-simplify-expand-to-threshold (rats thresh) ...)
(defun intervals-divide (intnum intden) ...)
(defun ratio-interval-to-simplest-ratio-within (rats) ...)
(defun ratio-interval-format-f (rats) ...)
(defun show-interval-refinement (oldint newint) ...)
;Above of general use, below specific to Newton's method for square root:
(defun ratio-to-quick-sqrt-seed (rat) ...)
(defun interval-to-sqrt-firstguess+bounds (nums) ...)
(defun interval+guess-to-sqrt-bounds (nums guess) ...)
(defun interval-newton-sqrt (nums oldrats) ...)

For a text file with a demo of those general utilities, see:
http://www.rawbw.com/~rem/demo-intarith-util.txt
Ask me if you have any questions about what any of the functions do.

Now for the fun part, Newton's method for computing square root, but
using intervals of rationals for input (usually) and output (always)!

;Given a *nonnegative* rational number (ratio or integer),
; find what range the number is in by powers of two,
; and generate a very crude but quickly-determined approximate
; square root as a ratio suitable for starting point for Newton's method.
(defun ratio-to-quick-sqrt-seed (rat)
  (prog (numflr numlen halfnb oddbit)
    (case rat ((0 1) (return rat)))
    (setq numflr (floor rat))
    (setq numlen (integer-length numflr))
    (cond ((< 0 numlen)
           (multiple-value-setq (halfnb oddbit) (floor numlen 2))
;range flr len hnb odd
;[0,1)  0   0  (0 0) (don't use, return 0,1 exactly, or invert)
;[1,2)  1   1  (0 1) -- 1 to 1.414 -- use 1.2 = 6/5
;[2,4) 2..3 2  (1 0) -- 1.414 to 2 -- use 1.666... = 2**1 * 5/6
;[4,8) 4..7 3  (1 1) -- 2 to 2.828 -- use 2.4      = 2**1 * 6/5
           (return (* (expt 2 halfnb)
                      (ecase oddbit (0 5/6) (1 6/5))))))
    (return (/ 1 (ratio-to-quick-sqrt-seed (/ 1 rat))))
    ))

The idea is to jump immediately to a point moderately close to the
correct square root, within a range of appx. 1 : 1.4, i.e. with an
error of only about 25% worst-case. From that point, Newton's method
will start rapidly-converging with the very first actual iteration.
Note: This is used when we want to take sqrt of an exact rational .

;Given a single rational number or an interval of rational numbers,
; compute a very crude approximate square root as seed for Newton's
; method, which is first return value, and very crude absolute
; bounds on the correct answer, namely between 1 and the original number,
; second return value. Numbers must be strictly non-negative of course!!
(defun interval-to-sqrt-firstguess+bounds (nums)
  (prog (guesses guess oldrats)
    (or (consp nums) (setq nums (list nums nums)))
    (setq guesses (mapcar #'ratio-to-quick-sqrt-seed nums))
    (setq guess (ratio-interval-to-simplest-ratio-within guesses))
    (setq oldrats (sort (cons 1 (copy-list nums)) #'<))
    (setq oldrats (cons (first oldrats) (last oldrats)))
    (return (values guess oldrats))
    ))

For 'guess', the seed, the starting point for Newton's method, one end
of the consequent interval (the other end will be NUM/guess), the idea
is that, given an interval instead of just a single number, as input
for the number whose square-root is desired, compute the seed (from
previous function) at each endpoint, and if they are different then
pick the simplest rational number between them, else just use the seed
as-is.

For 'oldrats', a guaranteed interval within which the square root must
reside, well sqrt(n) is between 1 and n, so that's the initial interval
generated above.

(defun interval+guess-to-sqrt-bounds (nums guess)
  (prog (unguess rats allspan nbprec uspan spanratio nbextra oldcx thresh)
    (setq unguess (intervals-divide nums guess))
    (setq rats (intervals-convext-hull guess unguess))
    ...
    (return rats)))

Well, the part shown is the obvious part of course. guess and 1/guess
are bounds for the correct square root, except here if nums is an
interval then unguess = 1/guess will also be an interval.
intervals-convext-hull takes care of the task of generating the
smallest interval that includes both guess and 1/guess.

But if nums is **very** precise, a very very narrow interval, then
unguess copies that extreme and unnecessary precision at early stages
when Newton's method is nowhere near the root yet. The ... part omited
above compares how much significance already exists within the interval
from guess to 1/guess, against the significance within unguess itself,
knowing that Newton's method generally doubles the number of
significant digits with each iteration (squares the precision), so
there's no point in working with unguess much more precise than twice
the significance we already have, so unguess is widened until it's
reasonable, by calling ratio-interval-simplify-expand-to-threshold. I'm
not publishing that source here for these reasons:
-1- It's not fully tuned yet.
-2- I might want to apply for a patent on the algorithm.
-3- Lacking a patent, I want to keep this a trade secret, and maybe
sell it to a suitable employer/agent/whatever.

;Given the original number (as exact rational, or interval of rationals),
; and previously-known bounds on the square-root of that number,
; run one regular Newton's method interation followed by one
; simplify-the-representation iteration, show refinements of interval
; from original, and return the final interval (bounds on square-root).
(defun interval-newton-sqrt (nums oldrats)
  (prog (guess1 midrats guess2 newrats precrat)
    (show-interval oldrats)
    (setq guess1 (/ (+ (car oldrats) (cadr oldrats)) 2)) (float guess1)
    (setq midrats (interval+guess-to-sqrt-bounds nums guess1))
    (setq midrats (show-interval-refinement oldrats midrats))
    (show-interval midrats)
    (setq guess2 (ratio-intervals-to-simplest-ratio-within midrats)) (float guess2)
    (setq newrats (interval+guess-to-sqrt-bounds nums guess2))
    (show-interval-refinement oldrats newrats)
    (show-interval newrats)
    (float (setq precrat (oldint+newint-to-precision-ratio oldrats newrats)) 1d0)
    (or (< 2 precrat)
      (format t "Interval didn't shrink enough:~%Old: ~A~%New: ~A~%~
                 Refinement (i.e. shrinkage) factor = ~F~%** RETURN to go on:"
        (ratio-interval-format-f oldrats)
        (ratio-interval-format-f newrats) precrat) (read-line))
    (return newrats)))

That's the key idea I invented this weekend: Alternating one regular
Newton's method (taking arithmetic mean of endpoints, and inverse of
that mean, as bounds on geometric mean), and one simplification
iteration (to prevent numerators and denominators of endpoints from
growing unnecessarily beyond what's needed to represent a given amount
of precision. This is what makes rational-number interval-arithmetic
practical, as of this weekend when I invented this. Now that I've
disclosed that algorithm, I can't patent that particular part of the
overall algorithm internationally, but I could still patent it within
the next year within the USA and a few other countries, so does anybody
think it's worth patenting?

Anyway, text files with a demos of Newton's method for
rational-endpoint continued fractions, using my code described above,
are available here:

With exact 2 as input, compute bounds on sqrt(2):
http://www.rawbw.com/~rem/demo-intarith-sqrt-2.txt

With only bounds on 2 as input, compute bounds on sqrt(2), manually
narrowing the bounds on 2 whenever Newton's method has done as well as
possible with the earlier bounds:
http://www.rawbw.com/~rem/demo-intarith-sqrt-2b.txt

Combine the above two methods: With exact 2, compute bounds on sqrt(2),
but then use those bounds to compute sqrt(sqrt(2)), manually emulating
lazy-evaluation whereby sqrt(2) is made more precise *only* when needed
as more precise input to compute sqrt(sqrt(2)):
http://www.rawbw.com/~rem/demo-intarith-sqrt-sqrt-2.txt

Next I plan to automate the process of lazy evaluation, by having an
uninterned symbol for each number represented in the dataflow for
interval-arithmetic calculations. For example, to do that last
calculation I'd set up this dataflow:
2 --> two --[SQRT]-> sqrt2 --[SQRT]-> frrt2 -->[OutputDemand]
Then OutputDemand would request a particular precision, and that
request would work its way back through the other processes and then
intervals of the specified accuracy would work their way forward to
OutputDemand.

But the real test would be this dataflow:

2 -->two--[SQRT]->sqrt2
      |                 three<--3   sqrt2-->[NEGATE]-->nsqrt2
      |  sqrt2            |                              |
      V    V              V                              V
    [MULTIPLY]-->sqrt8->[ADD]->x1--[SQRT]->x2--------->[ADD]-->y

y-->[OutputDemand]

I know the exact value of y, do you?
It'll be fun watching interval arithmetic compute bounds on that
number, and see it guess the exact value correctly (thanks to
ratio-interval-to-simplest-ratio-within) after only a few iterations,
as a demo that all this code really works.

If enough people show an interest in this software, maybe I'll:
-1- Set up a Web-server-side application to demonstrate it live
-2- Finish the basic arithmetic functions for interval arithmetic
-3- Implement Taylor/MacLaurin series for transcendental functions
     (exp[e], log[e], sin, cos, etc.) using interval arithmetic
     on alternating series only so that the convex hull of any
     two adjacent partial sums strictly bound the correct answer
-4- Set up a HTTP-RPC-sexpr interface so that customers can use my
     software from within their own CL applications

I wish somebody would hire me to write software for them.

ObHackSighting: Giving GIFs reasonable names, so that people who whose
Web access is text only can guess what they're missing and thereby make
sense out of the document:
   Linkname: Taylor Series
        URL: http://www.efunda.com/math/taylor_series/taylor_series.cfm
   When this expansion converges over a certain range of [x.gif] , that
   is, [rngoesto0.gif] , then the expansion is called the Taylor Series
   of [fx.gif] expanded about [a.gif] .
   If [ais0.gif] the series is called the MacLaurin Series:
(Why was I looking at that WebPage at this monent? Because I needed to
know the correct spelling of 'MacLaurin' in my message above, so I did
a Google search for "Taylor series" and the first hit was this winner
that not only contained the spelling but had this GIF naming hack I
couldn't resist using for show-and-tell.)

Oh, almost forgot: The other new idea is using inscribed and circumscribed
2**n-sided polygons to bound area of unit circle which is PI. Computing
some amd c



xxx 
some amd
xxx
sine and cosine of half-angle involves square root, which I've now
implemented.
From: ··········@YahooGroups.Com
Subject: Re: I'm not the only one who values interval-arithmetic
Date: 
Message-ID: <REM-2004apr18-005@Yahoo.Com>
Update from what I posted this afternoon: I went ahead and manually set
up the dataflow for  sqrt(3 + 2*sqrt(2)) - sqrt(2)  using my new
rational-endpoints interval-arithmetic software, and collected a
transcript of the demo:
http://www.rawbw.com/~rem/demo-intarith-sqrtsurd.txt
where the correct answer is clearly bounded extremely closely to a
particular value which is the correct mathematical value. I also
include a demo of how my interval arithmetic is better than
double-float arithmetic when the correct value is very close to
something "obvious" but not actually equal to it. Double-float shows it
to possibly be exactly that wrong value, whereas my interval arithmetic
eventually clearly bounds the value away from the "obvious" guess.
From: Joe Marshall
Subject: Re: I'm not the only one who values interval-arithmetic
Date: 
Message-ID: <8ygqu5ww.fsf@ccs.neu.edu>
··········@YahooGroups.Com writes:

> Update from what I posted this afternoon: I went ahead and manually set
> up the dataflow for  sqrt(3 + 2*sqrt(2)) - sqrt(2)  using my new
> rational-endpoints interval-arithmetic software, and collected a
> transcript of the demo:
> http://www.rawbw.com/~rem/demo-intarith-sqrtsurd.txt
> where the correct answer is clearly bounded extremely closely to a
> particular value which is the correct mathematical value. 


For those who are interested, I have an implementation of exact real
arithmetic using continued fractions.

(cf:render (cf*cf:- (cf:sqrt (rat*cf:+ 3 (rat*cf:* 2 (cf:sqrt 2)))) 2))
1.000000000000000000000000000000000000000000000000000000000000000...

For the `Rump Polynomial' I get this:
(cf:render (f1 77617 33096))
-0.827396059946821368141165095479816291999033115784384819917814841...

More digits available upon request.
From: ··········@YahooGroups.Com
Subject: Re: I'm not the only one who values interval-arithmetic
Date: 
Message-ID: <REM-2004apr21-011@Yahoo.Com>
> Date: Tue, 20 Apr 2004 09:37:19 -0400
> From: Joe Marshall <···@ccs.neu.edu>
> For those who are interested, I have an implementation of exact real
> arithmetic using continued fractions.

I believe Bill Gosper did basically the same thing circa 1975 using
MacSyma and/or MacLISP. Did you basically re-invent the wheel, or get
the idea from him, or translate his existing code to CL, or ...?

(I probably also re-invented a wheel with my past four days of work,
but it was interesting to do.)

> (cf:render (cf*cf:- (cf:sqrt (rat*cf:+ 3 (rat*cf:* 2 (cf:sqrt 2)))) 2))
> 1.000000000000000000000000000000000000000000000000000000000000000...

I believe you have a typo there. That last term of 2 you're subtracting
should be (cf:sqrt 2), just like the deep-inner term just before it,
right?

Does your cf code have a way to automatically pipeline a long chain of
calculations such as a power series or nested square-roots as I'm doing
(see later below)?

In any case, do you have your code interface to CGI (or any other
WebServer application interface) as a demo for people to play with?
I'll probably put a demo of my stuff up on CGI sometime soon if a few
readers here show an interest in the form of desire to play with it.

By the way, tonight I finished implementing one more use of my Newton's
method interval-rationals SQRT algorithm: I programmed half-angle
formulas, for small angles these are numerically stable:
  cos(A/2) = sqrt((cos(A)+1)/2)
  sin(A/2) = sin(A)/(2*cos(A/2))
then by daisy-chaining these starting with cos(PI/2)=0 and sin(PI/2)=1
I compute sine and cosine of 2*pi/2**n, and from the last of those in
the chain I compute area of inscribed and circumscribed polygons which
are lower and upper bound on PI respectively. Allowing Newton's method
to run until each link in the chain has more than 50 digits of
accuracy, and running the links until the bounds-intervals overlap,
which happens at n=112 (I get exactly two bits more accuracy with each
incf n, except near the end of the chain where the inaccuracy in
endpoints slows down the convergence slightly), I get these results:
cosine:
0.99999999999999999999999999999999999999999999999999999999999999999926[2..B]
sine:
0.00000000000000000000000000000000121009747292310784037823488144581130995333
  364392349547007623490994[7..A]
PI:
3.14159265358979323846264338327950288419716939937510582097494459230[1..D]
I've only memorized the first fifteen digits after the decimal point,
but I can't believe my program would give those correctly and make a
mistake later, so I trust those numbers.

Next I guess I should program the general case of power series
restricted to alternating with strictly decreasing absolute value, then
apply that general mechanism to special cases of exp(x) and the
trigonometric functions. For example, exp of any negative argument is
alternating, and if the absolute value of the argument is less than 1
the terms are strictly decreasing, so if the argument as given is not
in that range (-1, 0) I do a transformation of the argument into that
range and a transformation of the result back, for example exp(1/2) =
1/exp(-1/2), and exp(-1.5) = exp(-0.75)**2, etc.
From: Joe Marshall
Subject: Re: I'm not the only one who values interval-arithmetic
Date: 
Message-ID: <4qrd8kky.fsf@ccs.neu.edu>
··········@YahooGroups.Com writes:

>> Date: Tue, 20 Apr 2004 09:37:19 -0400
>> From: Joe Marshall <···@ccs.neu.edu>
>> For those who are interested, I have an implementation of exact real
>> arithmetic using continued fractions.
>
> I believe Bill Gosper did basically the same thing circa 1975 using
> MacSyma and/or MacLISP. Did you basically re-invent the wheel, or get
> the idea from him, or translate his existing code to CL, or ...?

I was unable to find actual code, but I found a memo by Gosper at

   http://www.tweedledum.com/rwg/cfup.htm

that describes the technique, and Peter Potts work at

  http://www.purplefinder.com/~potts/exact.html

is also a good resource.  So I didn't invent anything new, I just
translated the algorithms into working code.

> (I probably also re-invented a wheel with my past four days of work,
> but it was interesting to do.)
>
>> (cf:render (cf*cf:- (cf:sqrt (rat*cf:+ 3 (rat*cf:* 2 (cf:sqrt 2)))) 2))
>> 1.000000000000000000000000000000000000000000000000000000000000000...
>
> I believe you have a typo there. That last term of 2 you're subtracting
> should be (cf:sqrt 2), just like the deep-inner term just before it,
> right?

Yes, that is a typo.

> Does your cf code have a way to automatically pipeline a long chain of
> calculations such as a power series or nested square-roots as I'm doing
> (see later below)?

Well, sort of.  The primitive functions pipeline the computation
(which is why you need cf:render to actually pump the pipeline to get
digits out).  If you have a recursive expansion for a power series or
square root where the recursive term is an argument to a primitive, it
should pipeline automatically.  Otherwise, you would have to be aware
of the pipeline.

> In any case, do you have your code interface to CGI (or any other
> WebServer application interface) as a demo for people to play with?

Nope, but you can just run the code.

> By the way, tonight I finished implementing one more use of my Newton's
> method interval-rationals SQRT algorithm: I programmed half-angle
> formulas, for small angles these are numerically stable:
>   cos(A/2) = sqrt((cos(A)+1)/2)
>   sin(A/2) = sin(A)/(2*cos(A/2))
> then by daisy-chaining these starting with cos(PI/2)=0 and sin(PI/2)=1
> I compute sine and cosine of 2*pi/2**n, and from the last of those in
> the chain I compute area of inscribed and circumscribed polygons which
> are lower and upper bound on PI respectively. 

I haven't added trig functions yet.  I was reaching the point of
diminishing return: the math gets really, really hard.  I kludged the
square-root by finding the fixed-point of (lambda (y) (/ x y)) where x
is the number to find the square root of.  I just terminated it after
20 or so iterations.

> PI:
> 3.14159265358979323846264338327950288419716939937510582097494459230[1..D]

Using this continued fraction:

       4              1
       --  =  1 + ---------
       pi                 4
                  3 + ---------
                              9
                      5 + ---------
                                  16
                          7 + ---------
                                        .
                                          .
                                            .

we get this code:

(define pi
  (let ()
    (define odds
      (let loop ((i 1))
        (cons-stream i (loop (+ i 2)))))

    (define squares
      (let loop ((i 1))
        (cons-stream (* i i) (loop (add1 i)))))

    (normalize-cf (cons-stream 1 squares) odds 0 4 1 0)))

The (cons-stream 1 squares) gives us our numerator stream and the 
0 4 1 0 are the coefficients a, b, c, and d of the formula  

   ax + b
  --------
   cx + d

where X is the continued fraction, so this undoes the 4/pi. 


=> (cf:render pi 10 66)
  3.14159265358979323846264338327950288419716939937510582097494459230...

> PI:
> 3.14159265358979323846264338327950288419716939937510582097494459230[1..D]


Try computing this function (from Muller):

a0 = 61/11

a1 = 11/2

                      3000
             1130 - --------
                      a
                        n-2
a  = 111 - ------------------
 n            a
               n-1

(imagine that the n-1 and n-2 are subscripts)

The answer converges to 6, but if you use floating point, it will
converge to 100.  (even if you use an enormous amount of precision)
From: Raymond Toy
Subject: Re: I'm not the only one who values interval-arithmetic
Date: 
Message-ID: <4nk78dn5pa.fsf@edgedsp4.rtp.ericsson.se>
>>>>> "RobertMaas" == RobertMaas  <··········@YahooGroups.Com> writes:

    RobertMaas> By the way, there's nothing wrong with floating-point representation of
    RobertMaas> the endpoints of an interval. So the quote is slightly misleading. What
    RobertMaas> they mean as fast but nonrigorous is the standard best-guess
    RobertMaas> (shot-in-dark) floating-point calculations which are built into most
    RobertMaas> modern CPUs. I'm curious what representation was used for the endpoints

I understand that the UltraSparc III CPU has support for interval
arithmetic.  But I don't really know what support means here.  Perhaps
faster or better directed rounding operations?

    RobertMaas> hardware support for high speed, I wish programming languages such as
    RobertMaas> Lisp provided interval arithmetic as a data type and a package of
    RobertMaas> arithmetic and trigometric functions that used them. Has any Lisp
    RobertMaas> implementation expert created such a package, presumably more efficient
    RobertMaas> than anything I could hack together myself, and made it generally
    RobertMaas> available for use in CMUCL and other CLs?

CMUCL (and SBCL) has interval arithmetic built-in and is used by the
compiler.  I used to have (some) of it as a separate package, but I
don't think I have that anymore.  Plus, it's not "real" interval
arithmetic in the sense that directed rounding is not used to compute
the bounds of an interval.  Doing that in the context of the compiler
is trickier since you can do mixed type operations so keeping track of
the epsilons is much harder.  I decided to give up on that issue.

Ray
From: ··········@YahooGroups.Com
Subject: Re: I'm not the only one who values interval-arithmetic
Date: 
Message-ID: <REM-2003sep14-001@Yahoo.Com>
{{Date: Fri, 12 Sep 2003 15:54:09 -0400
  From: Raymond Toy <···@rtp.ericsson.se>
  I understand that the UltraSparc III CPU has support for interval
  arithmetic.}}

Ah, thanks for the pointer. I looked on Google and found:
File that you are currently viewing
   Linkname: UltraSPARC Processors, Documentation
        URL: http://www.sun.com/processors/manuals/index.html
Link that you currently have selected
   Linkname: Interval Arithmetic in High performance Technical Computing
        URL: http://www.sun.com/processors/whitepapers/ia12_wp.pdf
Unfortunately I don't have any way to view PDF files over VT100 dialup
link into Unix shell, so maybe I'll have to take a look when I'm at the
public library next time.

{{CMUCL (and SBCL) has interval arithmetic built-in and is used by the
  compiler.}}

Hmm, when I was browsing the online HTML-format CMUCL manual
http://cvs2.cons.org/ftp-area/cmucl/doc/cmu-user/
I didn't happen to see anything about interval arithmetic. Looking
carefully just now, I see rounding mode discussed in regard to floating
point.
http://cvs2.cons.org/ftp-area/cmucl/doc/cmu-user/extensions.html#toc11
http://cvs2.cons.org/ftp-area/cmucl/doc/cmu-user/extensions.html#toc15
The concept index:
http://cvs2.cons.org/ftp-area/cmucl/doc/cmu-user/cmu-user.hcnd.html
doesn't mention the word 'interval'.

{{it's not "real" interval arithmetic in the sense that directed
  rounding is not used to compute the bounds of an interval.}}

You mean it doesn't compute the worst-case bounds correctly at all,
just uses floating-point shot-in-dark result instead of correct lower
or upper bound, so it's basically a sham?

Back to rounding mode (see URLs above): 
   IEEE floating point specifies four possible rounding modes: ...
   :positive-infinity ...
   :negative-infinity ...
those would be the two modes useful for interval arithmetic. But:
      Warning:
   Although the rounding mode can be changed with
   set-floating-point-modes, use of any value other than the default
   (:nearest) can cause unusual behavior, since it will affect rounding
   done by Common Lisp system code as well as rounding in user code.
So has anyone reading this thread actually tried using this feature of
CMUCL so they can advise me how to avoid trouble? For example, if I
effectively bind the rounding mode around each single arithmetic
operation, including no program logic or anything else within the
binding, would that avoid all such problems? Would it be too paranoid?
   extensions:set-floating-point-modes &key :traps :rounding-mode
   :fast-mode :accrued-exceptions
   :current-exceptions :precision-control
i.e. for my case:
(extensions:set-floating-point-modes :rounding-mode :negative-infinity)
(extensions:set-floating-point-modes :rounding-mode :positive-infinity)
So that would allow me to clobber the current mode and make it whatever
mode I wanted. But first saving the old mode seems more messy:
   extensions:get-floating-point-modes
(extensions:get-floating-point-modes) ==>
(:TRAPS (:OVERFLOW :INVALID :DIVIDE-BY-ZERO) :ROUNDING-MODE :NEAREST
 :CURRENT-EXCEPTIONS (:INEXACT) :ACCRUED-EXCEPTIONS (:INEXACT) :FAST-MODE NIL)
So does that cons up a new list each time it's called, which would make
it expensive relative to the single machine instruction implied by the
floating-point arithmetic operation? Or does it maintain this list
without re-building it, in which case it's probably dangerous for any
user program to modify it?
* (/ 1e0 3e0)
0.33333334
* (integer-decode-float *)
11184811
-25
1
* (extensions:set-floating-point-modes :rounding-mode :negative-infinity)
* (/ 1e0 3e0)
0.3333333
* (integer-decode-float *)
11184810
-25
1
Well, it does seem to basically work, so maybe I'll hack together a
proper (not sham) interval-arithmetic package for the elementary
operations (+ - * /) and SQRT too. For these low-level operations, I'll
have them take low and high inputs as separate arguments and return low
and high result as separate return values. These can then be
incorporated into single-handle representation however desired later.

I suppose, after initial experimentation, the first main thing to write
will be a macro WITH-FLOAT-ROUNDING-MODE analagous to the existing
   [Macro]
   ext:with-float-traps-masked traps &body body
From: Raymond Toy
Subject: Re: I'm not the only one who values interval-arithmetic
Date: 
Message-ID: <3F64D751.5040008@rtp.ericsson.se>
··········@YahooGroups.Com wrote:
> 
> {{CMUCL (and SBCL) has interval arithmetic built-in and is used by the
>   compiler.}}
> 
> Hmm, when I was browsing the online HTML-format CMUCL manual
> http://cvs2.cons.org/ftp-area/cmucl/doc/cmu-user/
> I didn't happen to see anything about interval arithmetic. Looking
> carefully just now, I see rounding mode discussed in regard to floating

Because it's not exposed to the user.  It's only used by the compiler 
for type derivation.

> {{it's not "real" interval arithmetic in the sense that directed
>   rounding is not used to compute the bounds of an interval.}}
> 
> You mean it doesn't compute the worst-case bounds correctly at all,
> just uses floating-point shot-in-dark result instead of correct lower
> or upper bound, so it's basically a sham?

You are correct---it does not compute worst-case bounds at all.  So far, 
it's been good enough for the compiler because people only roughly 
declare things.  For me, about I all ever declare is that a single or 
double float never contains a negative value or negative zero, so that 
sqrt and log and friends can be inlined when possible.

And it doesn't do worse-case bounds because I was basically lazy.  It's 
also complicated by the fact that I didn't want to handle the mixed 
float type issues.

> those would be the two modes useful for interval arithmetic. But:
>       Warning:
>    Although the rounding mode can be changed with
>    set-floating-point-modes, use of any value other than the default
>    (:nearest) can cause unusual behavior, since it will affect rounding
>    done by Common Lisp system code as well as rounding in user code.
> So has anyone reading this thread actually tried using this feature of
> CMUCL so they can advise me how to avoid trouble? For example, if I
> effectively bind the rounding mode around each single arithmetic
> operation, including no program logic or anything else within the
> binding, would that avoid all such problems? Would it be too paranoid?

I believe it's just a warning that you shouldn't set it to other values 
and expect everything else to work.  Things like special functions and 
such.  If you keep the rounding changes localized so only basic 
arithmetic functions are affected, you should be ok and get the desired 
results.

>    extensions:set-floating-point-modes &key :traps :rounding-mode
>    :fast-mode :accrued-exceptions
>    :current-exceptions :precision-control
> i.e. for my case:
> (extensions:set-floating-point-modes :rounding-mode :negative-infinity)
> (extensions:set-floating-point-modes :rounding-mode :positive-infinity)
> So that would allow me to clobber the current mode and make it whatever
> mode I wanted. But first saving the old mode seems more messy:
>    extensions:get-floating-point-modes
> (extensions:get-floating-point-modes) ==>
> (:TRAPS (:OVERFLOW :INVALID :DIVIDE-BY-ZERO) :ROUNDING-MODE :NEAREST
>  :CURRENT-EXCEPTIONS (:INEXACT) :ACCRUED-EXCEPTIONS (:INEXACT) :FAST-MODE NIL)
> So does that cons up a new list each time it's called, which would make

Yes, I think it does cons a new list everytime.  Obviously, there is a 
way to get at the machine bits that hold this info.  It's not exposed to 
the user, though.

> Well, it does seem to basically work, so maybe I'll hack together a
> proper (not sham) interval-arithmetic package for the elementary
> operations (+ - * /) and SQRT too. For these low-level operations, I'll
> have them take low and high inputs as separate arguments and return low
> and high result as separate return values. These can then be
> incorporated into single-handle representation however desired later.
> 
> I suppose, after initial experimentation, the first main thing to write
> will be a macro WITH-FLOAT-ROUNDING-MODE analagous to the existing
>    [Macro]
>    ext:with-float-traps-masked traps &body body
> 

If you can make the interval arithmetic package look somewhat like 
what's in CMUCL and SBCL, I'm sure they'd be grateful. :-)  These 
routines operate on an interval type containing the low and high bounds. 
  And it's complicated by supporting open and closed intervals too, and 
mixing operations of many different operand types.

Ray
From: ··········@YahooGroups.Com
Subject: Re: I'm not the only one who values interval-arithmetic
Date: 
Message-ID: <REM-2004apr19-002@Yahoo.Com>
> From: Raymond Toy <···@rtp.ericsson.se>
> > Hmm, when I was browsing the online HTML-format CMUCL manual
> > http://cvs2.cons.org/ftp-area/cmucl/doc/cmu-user/
> > I didn't happen to see anything about interval arithmetic. Looking
> > carefully just now, I see rounding mode discussed in regard to floating
> Because it's not exposed to the user.  It's only used by the compiler
> for type derivation.

But is it written in LISP, composed of regular functions in some
package? Or is it only assembly-language code not interfaced in any way
that LISP software can use it without special knowledge? If it's
written in LISP, is there some internal documentation that says the
name of the package and something about how it works?

> You are correct---it does not compute worst-case bounds at all.

Yuk. Does the compiler ever compile wrong code, or fail to or
spuriously issue a warning or error signal, because it slops the end of
an interval past a threshold for a decision?

> And it doesn't do worse-case bounds because I was basically lazy.

But there's a feeling of major dissatisfaction doing that, or at least
IMO there should be such a feeling. Like why bother to write the code
at all if the boundaries are wrong? Doing boundaries correctly
shouldn't be really difficult, should it? Using rational numbers, I
didn't have any problem doing boundaries correctly. (But see below
about open/closed intervals.)

> It's also complicated by the fact that I didn't want to handle the
> mixed float type issues.

Hmm, when I wrote an interval arithmetic package recently I avoid
floats totally, using rational numbers throughout, for basically the
same reason, not wanting to deal with different kinds of numbers as
special cases.

> If you can make the interval arithmetic package look somewhat like
> what's in CMUCL and SBCL, I'm sure they'd be grateful.

"they? ??

Would they be grateful enough to pay me real money for my work so I
won't become homeless in a few months when I max out my credit cards
paying the rent?

If what they did for compiler-use only isn't documented anywhere I can
read and isn't accessible from me for experimenting with it to get a
feel for how it really works, how could I know whether I'm doing
something that looks similar or not?

> These routines operate on an interval type containing the low and
> high bounds.

Yes, that's how I do it too. That makes a lot more sense (for internal
representation) than center and tolerance, because bounds of a result
can be computed directly from bounds of arguments to function.

> And it's complicated by supporting open and closed intervals too,

For my purpose, actual numeric results, I have no need for that. I use
closed intervals only. But I can see for type declarations, such as
saying a number is non-negative, you have to have an open bound at
infinity, so once you start down that path ...

If I get in the mood someday, I might generalize my software to include
both open and closed endpoints, just to see if it makes the software
seriously more complicated. (And of course if anybody ever agrees to
pay me for my time I can be motivated to do what my boss wants.)

> and mixing operations of many different operand types.

Yeah. I avoided that by using rational numbers exclusively.
The only type-special case in my code is for integers which must be
coerced to an interval by mapping NUM -> (NUM NUM) before some
operations that look at the individual bounds. If I have data in
floating point format, I know whether it's intended to be an exact
binary fraction or an interval between two binary fractions, and so I
manually coerce it to the appropriate single-number or range before
passing to my interval arithmetic code.

> Date: Sun, 14 Sep 2003 21:02:10 GMT
(I didn't see your article when it first appeared, because I didn't
have any efficient method for finding all followups (to stuff I had
posted) until just a few nights ago when I finally found your article
and put it in the queue to compose a followup myself. Sorry for very
belated response.)
From: Raymond Toy
Subject: Re: I'm not the only one who values interval-arithmetic
Date: 
Message-ID: <sxdad16262j.fsf@edgedsp4.rtp.ericsson.se>
>>>>> "RobertMaas" == RobertMaas  <··········@YahooGroups.Com> writes:

    >> From: Raymond Toy <···@rtp.ericsson.se>
    >> > Hmm, when I was browsing the online HTML-format CMUCL manual
    >> > http://cvs2.cons.org/ftp-area/cmucl/doc/cmu-user/
    >> > I didn't happen to see anything about interval arithmetic. Looking
    >> > carefully just now, I see rounding mode discussed in regard to floating
    >> Because it's not exposed to the user.  It's only used by the compiler
    >> for type derivation.

    RobertMaas> But is it written in LISP, composed of regular functions in some
    RobertMaas> package? Or is it only assembly-language code not interfaced in any way
    RobertMaas> that LISP software can use it without special knowledge? If it's
    RobertMaas> written in LISP, is there some internal documentation that says the
    RobertMaas> name of the package and something about how it works?

Yes, it's pure Lisp.  Some of the basic operations like interval
addition, subtraction, etc., are separate.  However, lots of it is
also tied to the compiler internals.

If you really want to know it's in src/code/srctran.lisp for the most
part.  Look for make-interval, interval-add, etc.  But you also need
to look at the deftransforms there.

Sorry, the only documentation is the comment in the code.  Also, note
that the intervals assume that -0.0 is different from 0.0, unlike most
of Lisp.

    >> You are correct---it does not compute worst-case bounds at all.

    RobertMaas> Yuk. Does the compiler ever compile wrong code, or fail to or
    RobertMaas> spuriously issue a warning or error signal, because it slops the end of
    RobertMaas> an interval past a threshold for a decision?

Not that I know of.  But most people don't write hairy declarations.
I think it's usually of the form (or (member 0d0) (double-float 0d0))
to tell the compiler the number cannot be negative, or (double-float
0d0 1d0) so things like (sqrt (- 1 x)) doesn't use generic sqrt.

    >> And it doesn't do worse-case bounds because I was basically lazy.

    RobertMaas> But there's a feeling of major dissatisfaction doing that, or at least
    RobertMaas> IMO there should be such a feeling. Like why bother to write the code
    RobertMaas> at all if the boundaries are wrong? Doing boundaries correctly
    RobertMaas> shouldn't be really difficult, should it? Using rational numbers, I
    RobertMaas> didn't have any problem doing boundaries correctly. (But see below
    RobertMaas> about open/closed intervals.)

I just didn't feel like it, it hasn't been a problem in practice, and
it's a lot of work for what would seem to be a small gain.  People
just don't write hairy declarations.

The problem is dealing with mixed types and keep track of it all.
Using rationals is one approach, but what do you do when you
(single-float 1.0 (2.0))?

    >> If you can make the interval arithmetic package look somewhat like
    >> what's in CMUCL and SBCL, I'm sure they'd be grateful.

    RobertMaas> "they? ??

    RobertMaas> Would they be grateful enough to pay me real money for my work so I
    RobertMaas> won't become homeless in a few months when I max out my credit cards
    RobertMaas> paying the rent?

If you use any free lisp compiler be sure to send them some money
too. 

    >> and mixing operations of many different operand types.

    RobertMaas> Yeah. I avoided that by using rational numbers exclusively.

Don't you have problems with rationals with numerators and
denominators growing very large very quickly?  Or do you simplify them
to something more manageable after each operation?  If so, doesn't
that get rather expensive?

Ray
From: John
Subject: Installing mod_lisp in win32
Date: 
Message-ID: <oprveog9crfhfgdd@news.chello.no>
I am having trouble installing mod_lisp on under Windows XP.
 -rwxrwx---+   1 Administ SYSTEM      32851 Jul  8 21:02 mod_isapi.so
-rwxrwx---    1 Administ SYSTEM       9216 Sep 12 21:14 mod_lisp.so
-rwxrwx---+   1 Administ SYSTEM      28757 Jul  8 21:02 mod_log_config.so
This is how the listing shows (Cygwin ls -l)
Note that the alternate access flag (+) s not set.
 When I try to load the module Apache reports:
 No Access to file. (in norwegian)
 I changed the group and user (chnmod, chgrp) but how do I set this 
alternate access?
(if this is inded the flaw)

 John

-- 
Using M2, Opera's revolutionary e-mail client: http://www.opera.com/m2/