aboutsummaryrefslogtreecommitdiff
path: root/src/algebra/interval.spad.pamphlet
blob: 7fbd9fe8f321181c9a8e19d1ab5886782358eefb (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra interval.spad}
\author{Mike Dewar}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{category INTCAT IntervalCategory}
<<category INTCAT IntervalCategory>>=
)abbrev category INTCAT IntervalCategory
+++ Author: Mike Dewar
+++ Date Created: November 1996
+++ Date Last Updated:
+++ Basic Functions:
+++ Related Constructors: 
+++ Also See:
+++ AMS Classifications:
+++ Keywords:
+++ References:
+++ Description:
+++ This category implements of interval arithmetic and transcendental
+++ functions over intervals.
IntervalCategory(R: Join(FloatingPointSystem,TranscendentalFunctionCategory)):
 Category == Join(GcdDomain, OrderedSet, TranscendentalFunctionCategory, RadicalCategory, RetractableTo(Integer)) with
  approximate
  interval : (R,R) -> %
    ++ interval(inf,sup) creates a new interval, either \axiom{[inf,sup]} if
    ++ \axiom{inf <= sup} or \axiom{[sup,in]} otherwise.
  qinterval : (R,R) -> %
    ++ qinterval(inf,sup) creates a new interval \axiom{[inf,sup]}, without
    ++ checking the ordering on the elements.
  interval : R -> %
    ++ interval(f) creates a new interval around f.
  interval : Fraction Integer -> %
    ++ interval(f) creates a new interval around f.
  inf : % -> R
    ++ inf(u) returns the infinum of \axiom{u}.
  sup : % -> R
    ++ sup(u) returns the supremum of \axiom{u}.
  width : % -> R
    ++ width(u) returns \axiom{sup(u) - inf(u)}.
  positive? : % -> Boolean
    ++ positive?(u) returns \axiom{true} if every element of u is positive,
    ++ \axiom{false} otherwise.
  negative? : % -> Boolean
    ++ negative?(u) returns \axiom{true} if every element of u is negative,
    ++ \axiom{false} otherwise.
  contains? : (%,R) -> Boolean
    ++ contains?(i,f) returns true if \axiom{f} is contained within the interval
    ++ \axiom{i}, false otherwise.

@
\section{domain INTRVL Interval}
<<domain INTRVL Interval>>=
)abbrev domain INTRVL Interval
+++ Author: Mike Dewar
+++ Date Created: November 1996
+++ Date Last Updated:
+++ Basic Functions:
+++ Related Constructors: 
+++ Also See:
+++ AMS Classifications:
+++ Keywords:
+++ References:
+++ Description:
+++ This domain is an implementation of interval arithmetic and transcendental
+++ functions over intervals.
Interval(R:Join(FloatingPointSystem,TranscendentalFunctionCategory)): IntervalCategory(R) == add 

  import Integer
--  import from R

  Rep := Record(Inf:R, Sup:R)

  roundDown(u:R):R == 
    if zero?(u) then float(-1,-(bits() pretend Integer))
                else float(mantissa(u) - 1,exponent(u))

  roundUp(u:R):R   == 
    if zero?(u) then float(1, -(bits()) pretend Integer)
                else float(mantissa(u) + 1,exponent(u))

  -- Sometimes the float representation does not use all the bits (e.g. when
  -- representing an integer in software using arbitrary-length Integers as
  -- your mantissa it is convenient to keep them exact).  This function 
  -- normalises things so that rounding etc. works as expected.  It is only
  -- called when creating new intervals.
  normaliseFloat(u:R):R == 
    zero? u => u
    m : Integer := mantissa u
    b : Integer := bits() pretend Integer
    l : Integer := length(m)
    if l < b then 
      BASE : Integer := base()$R pretend Integer
      float(m*BASE**((b-l) pretend PositiveInteger),exponent(u)-b+l)
    else
      u

  interval(i:R,s:R):% == 
    i > s =>  [roundDown normaliseFloat s,roundUp normaliseFloat i]
    [roundDown normaliseFloat i,roundUp normaliseFloat s]

  interval(f:R):% ==  
    zero?(f) => 0
    one?(f)  => 1
    -- This next part is necessary to allow e.g. mapping between Expressions:
    -- AXIOM assumes that Integers stay as Integers!
--    import from Union(value1:Integer,failed:"failed")
    fnew : R := normaliseFloat f
    retractIfCan(f)@Union(Integer,"failed") case "failed" =>
      [roundDown fnew, roundUp fnew]
    [fnew,fnew]

  qinterval(i:R,s:R):% ==
    [roundDown normaliseFloat i,roundUp normaliseFloat s]

  exactInterval(i:R,s:R):% == [i,s]
  exactSupInterval(i:R,s:R):% == [roundDown i,s]
  exactInfInterval(i:R,s:R):% == [i,roundUp s]

  inf(u:%):R == u.Inf
  sup(u:%):R == u.Sup
  width(u:%):R == u.Sup - u.Inf

  contains?(u:%,f:R):Boolean == (f > inf(u)) and (f < sup(u))

  positive?(u:%):Boolean == inf(u) > 0
  negative?(u:%):Boolean == negative? sup(u)

  a:% < b:% ==
    if inf(a) < inf(b) then
      true
    else if inf(a) > inf(b) then
      false
    else
      sup(a) < sup(b)

  a:% + b:% == 
    -- A couple of blatent hacks to preserve the Ring Axioms!
    if zero?(a) then return(b) else if zero?(b) then return(a)
    if a = b then return qinterval(2*inf(a),2*sup(a))
    qinterval(inf(a) + inf(b), sup(a) + sup(b))


  a:% - b:% ==  
    if zero?(a) then return(-b) else if zero?(b) then return(a)
    if a = b then 0 else qinterval(inf(a) - sup(b), sup(a) - inf(b))


  a:% * b:% == 
    -- A couple of blatent hacks to preserve the Ring Axioms!
    if one?(a) then return(b) else if one?(b) then return(a)
    if zero?(a) then return(0) else if zero?(b) then return(0)
    prods : List R :=  sort [inf(a)*inf(b),sup(a)*sup(b),
                             inf(a)*sup(b),sup(a)*inf(b)]
    qinterval(first prods, last prods)


  a:Integer * b:% == 
    if (a > 0) then 
      qinterval(a*inf(b),a*sup(b))
    else if negative? a then
      qinterval(a*sup(b),a*inf(b))
    else
      0

  a:PositiveInteger * b:% == qinterval(a*inf(b),a*sup(b))

  a:% ** n:PositiveInteger == 
    contains?(a,0) and zero?((n pretend Integer) rem 2) =>
      interval(0,max(inf(a)**n,sup(a)**n)) 
    interval(inf(a)**n,sup(a)**n)


  -(a:%) == exactInterval(-sup(a),-inf(a))

  a:% = b:% == (inf(a)=inf(b)) and (sup(a)=sup(b))
  a:% ~= b:% == (inf(a)~=inf(b)) or (sup(a)~=sup(b))

  1 == 
	one : R := normaliseFloat 1
	[one,one]

  0 == [0,0]

  recip(u:%):Union(%,"failed") == 
   contains?(u,0) => "failed"
   vals:List R := sort [1/inf(u),1/sup(u)]$List(R)
   qinterval(first vals, last vals)


  unit?(u:%):Boolean == contains?(u,0)

  (u:% exquo v:%): Union(%,"failed") ==
   contains?(v,0) => "failed"
   one?(v) => u
   u=v => 1
   u=-v => -1
   vals:List R := sort [inf(u)/inf(v),inf(u)/sup(v),sup(u)/inf(v),sup(u)/sup(v)]$List(R)
   qinterval(first vals, last vals)


  gcd(u:%,v:%):% == 1

  coerce(u:Integer):% ==
    ur := normaliseFloat(u::R)
    exactInterval(ur,ur)


  interval(u:Fraction Integer):% == 
--    import   log2 : % -> %
--             coerce : Integer -> %
--             retractIfCan : % -> Union(value1:Integer,failed:"failed")
--    from Float
    flt := u::R

    -- Test if the representation in R is exact
    --den := denom(u)::Float
    bin : Union(Integer,"failed") := retractIfCan(log2(denom(u)::Float))
    bin case Integer and length(numer u)$Integer < (bits() pretend Integer) => 
      flt := normaliseFloat flt
      exactInterval(flt,flt)

    qinterval(flt,flt)


  retractIfCan(u:%):Union(Integer,"failed") == 
    not zero? width(u) => "failed"
    retractIfCan inf u
  

  retract(u:%):Integer == 
    not zero? width(u) =>
      error "attempt to retract a non-Integer interval to an Integer"
    retract inf u
  

  coerce(u:%):OutputForm ==
    bracket([coerce inf(u), coerce sup(u)]$List(OutputForm))

  characteristic:NonNegativeInteger == 0


  -- Explicit export from TranscendentalFunctionCategory
  pi():% == qinterval(pi(),pi())

  -- From ElementaryFunctionCategory
  log(u:%):% == 
    positive?(u) => qinterval(log inf u, log sup u)
    error "negative logs in interval"
  

  exp(u:%):% == qinterval(exp inf u, exp sup u)

  u:% ** v:% == 
    zero?(v) => if zero?(u) then error "0**0 is undefined" else 1
    one?(u)  => 1
    expts : List R :=  sort [inf(u)**inf(v),sup(u)**sup(v),
                             inf(u)**sup(v),sup(u)**inf(v)]
    qinterval(first expts, last expts)

  -- From TrigonometricFunctionCategory

  -- This function checks whether an interval contains a value of the form
  -- `offset + 2 n pi'.
  hasTwoPiMultiple(offset:R,ipi:R,i:%):Boolean == 
    next : Integer := retract ceiling( (inf(i) - offset)/(2*ipi) )
    contains?(i,offset+2*next*ipi)
  

  -- This function checks whether an interval contains a value of the form
  -- `offset + n pi'.
  hasPiMultiple(offset:R,ipi:R,i:%):Boolean == 
    next : Integer := retract ceiling( (inf(i) - offset)/ipi )
    contains?(i,offset+next*ipi)
  

  sin(u:%):% == 
    ipi : R := pi()$R
    hasOne? : Boolean := hasTwoPiMultiple(ipi/(2::R),ipi,u)
    hasMinusOne? : Boolean := hasTwoPiMultiple(3*ipi/(2::R),ipi,u)

    if hasOne? and hasMinusOne? then 
      exactInterval(-1,1)
    else 
      vals : List R := sort [sin inf u, sin sup u]
      if hasOne? then
        exactSupInterval(first vals, 1)
      else if hasMinusOne? then
        exactInfInterval(-1,last vals)
      else
        qinterval(first vals, last vals)
    
  

  cos(u:%):% == 
    ipi : R := pi()
    hasOne? : Boolean := hasTwoPiMultiple(0,ipi,u)
    hasMinusOne? : Boolean := hasTwoPiMultiple(ipi,ipi,u)

    if hasOne? and hasMinusOne? then 
      exactInterval(-1,1)
    else 
      vals : List R := sort [cos inf u, cos sup u]
      if hasOne? then
        exactSupInterval(first vals, 1)
      else if hasMinusOne? then
        exactInfInterval(-1,last vals)
      else
        qinterval(first vals, last vals)
    
  

  tan(u:%):% == 
    ipi : R := pi()
    if width(u) > ipi then
      error "Interval contains a singularity"
    else 
      -- Since we know the interval is less than pi wide, monotonicity implies
      -- that there is no singularity.  If there is a singularity on a endpoint
      -- of the interval the user will see the error generated by R.
      lo : R := tan inf u 
      hi : R := tan sup u

      lo > hi => error "Interval contains a singularity"
      qinterval(lo,hi)
    
  

  csc(u:%):% == 
    ipi : R := pi()
    if width(u) > ipi then
      error "Interval contains a singularity"
    else 
--      import from Integer
      -- singularities are at multiples of Pi
      if hasPiMultiple(0,ipi,u) then error "Interval contains a singularity"
      vals : List R := sort [csc inf u, csc sup u]
      if hasTwoPiMultiple(ipi/(2::R),ipi,u) then 
        exactInfInterval(1,last vals)
      else if hasTwoPiMultiple(3*ipi/(2::R),ipi,u) then
        exactSupInterval(first vals,-1)
      else
        qinterval(first vals, last vals)
    
  

  sec(u:%):% == 
    ipi : R := pi()
    if width(u) > ipi then
      error "Interval contains a singularity"
    else 
--      import from Integer
      -- singularities are at Pi/2 + n Pi
      if hasPiMultiple(ipi/(2::R),ipi,u) then
        error "Interval contains a singularity"
      vals : List R := sort [sec inf u, sec sup u]
      if hasTwoPiMultiple(0,ipi,u) then 
        exactInfInterval(1,last vals)
      else if hasTwoPiMultiple(ipi,ipi,u) then
        exactSupInterval(first vals,-1)
      else
        qinterval(first vals, last vals)
    
  


  cot(u:%):% == 
    ipi : R := pi()
    if width(u) > ipi then
      error "Interval contains a singularity"
    else 
      -- Since we know the interval is less than pi wide, monotonicity implies
      -- that there is no singularity.  If there is a singularity on a endpoint
      -- of the interval the user will see the error generated by R.
      hi : R := cot inf u 
      lo : R := cot sup u

      lo > hi => error "Interval contains a singularity"
      qinterval(lo,hi)
    
  

  -- From ArcTrigonometricFunctionCategory

  asin(u:%):% == 
    lo : R := inf(u)
    hi : R := sup(u)
    if (lo < -1) or (hi > 1) then error "asin only defined on the region -1..1"
    qinterval(asin lo,asin hi)
  

  acos(u:%):% == 
    lo : R := inf(u)
    hi : R := sup(u)
    if (lo < -1) or (hi > 1) then error "acos only defined on the region -1..1"
    qinterval(acos hi,acos lo)
  

  atan(u:%):% == qinterval(atan inf u, atan sup u)

  acot(u:%):% == qinterval(acot sup u, acot inf u)

  acsc(u:%):% == 
    lo : R := inf(u)
    hi : R := sup(u)
    if ((lo <= -1) and (hi >= -1)) or ((lo <= 1) and (hi >= 1)) then
      error "acsc not defined on the region -1..1"
    qinterval(acsc hi, acsc lo)
  

  asec(u:%):% == 
    lo : R := inf(u)
    hi : R := sup(u)
    if ((lo < -1) and (hi > -1)) or ((lo < 1) and (hi > 1)) then
      error "asec not defined on the region -1..1"
    qinterval(asec lo, asec hi)
  

  -- From HyperbolicFunctionCategory

  tanh(u:%):% == qinterval(tanh inf u, tanh sup u)

  sinh(u:%):% == qinterval(sinh inf u, sinh sup u)

  sech(u:%):% == 
    negative? u => qinterval(sech inf u, sech sup u)
    positive? u => qinterval(sech sup u, sech inf u)
    vals : List R := sort [sech inf u, sech sup u]
    exactSupInterval(first vals,1)
  

  cosh(u:%):% == 
    negative? u => qinterval(cosh sup u, cosh inf u)
    positive? u => qinterval(cosh inf u, cosh sup u)
    vals : List R := sort [cosh inf u, cosh sup u]
    exactInfInterval(1,last vals)
  

  csch(u:%):% == 
    contains?(u,0) => error "csch: singularity at zero"
    qinterval(csch sup u, csch inf u)
  

  coth(u:%):% == 
    contains?(u,0) => error "coth: singularity at zero"
    qinterval(coth sup u, coth inf u)
  

  -- From ArcHyperbolicFunctionCategory

  acosh(u:%):% == 
    inf(u)<1 => error "invalid argument: acosh only defined on the region 1.."
    qinterval(acosh inf u, acosh sup u)
  

  acoth(u:%):% == 
    lo : R := inf(u)
    hi : R := sup(u)
    if ((lo <= -1) and (hi >= -1)) or ((lo <= 1) and (hi >= 1)) then
      error "acoth not defined on the region -1..1"
    qinterval(acoth hi, acoth lo)
  

  acsch(u:%):% == 
    contains?(u,0) => error "acsch: singularity at zero"
    qinterval(acsch sup u, acsch inf u)
  

  asech(u:%):% == 
    lo : R := inf(u)
    hi : R := sup(u)
    if  (lo <= 0) or (hi > 1) then 
      error "asech only defined on the region 0 < x <= 1"
    qinterval(asech hi, asech lo)
  

  asinh(u:%):% == qinterval(asinh inf u, asinh sup u)

  atanh(u:%):% == 
    lo : R := inf(u)
    hi : R := sup(u)
    if  (lo <= -1) or (hi >= 1) then 
      error "atanh only defined on the region -1 < x < 1"
    qinterval(atanh lo, atanh hi)
  

  -- From RadicalCategory
  u:% ** n:Fraction Integer == interval(inf(u)**n,sup(u)**n)
  
@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
<<license>>

<<category INTCAT IntervalCategory>>
<<domain INTRVL Interval>>

@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}