Changeset 39960 in project for wiki/eggref/5/math


Ignore:
Timestamp:
04/11/21 01:04:05 (2 months ago)
Author:
Diego
Message:

Add (math flonum) documentation, re-order sections

File:
1 edited

Legend:

Unmodified
Added
Removed
  • wiki/eggref/5/math

    r39757 r39960  
    55=== Introduction
    66
    7 math is a chicken port of racket's math library.
    8 
    9 The following documentation is largely taken directly from the racket
    10 documentation, but tweaked for the CHICKEN implementation.
     7math is a CHICKEN port of racket's [[https://docs.racket-lang.org/math/|math]]
     8library.
     9
     10The following documentation is largely a direct copy of the racket
     11documentation, re-formatted for svnwiki and tweaked for the CHICKEN
     12implementation.
     13
     14=== Development status
     15
     16This egg is still largely under active development. There may be missing
     17modules and features.
     18
     19=== Implementation caveats
     20
     21* It's possible some undefined behavior may occur with arguments of the wrong
     22  type, since a good amount of the original functions were originally defined
     23  in typed racket, which AFAIK would catch those and throw an exception.
     24
     25* In some places the original implementation references {{unsafe-}} {{fx}} and
     26  {{fl}} operators, but these are actually just aliased to safe operators. This
     27  implementation just uses CHICKEN's {{chicken.fixnum}} module, which is
     28  unsafe. This may also lead to undefined behavior in some cases.
    1129
    1230=== Modules
     
    1533
    1634Constants and elementary functions
     35
     36'''Note:''' This is currently missing procedures from {{racket/math}} that the
     37original racket library re-exports.
    1738
    1839===== Constants
     
    5475
    5576Returns {{#t}} when {{v}} is an inexact complex number. Analogous to
    56 [[https://wiki.call-cc.org/man/5/Module%20(chicken%20base)#flonum|flonum?]]
     77{{flonum?}}.
    5778
    5879<procedure>(number->float-complex x) -> cplxnum</procedure>
     
    7596> (power-of-two? 1/2)
    7697#t
    77 > (power-of-two? (flnext 2.0))
     98> (power-of-two? (fpnext 2.0))
    7899#f
    79100</enscript>
     
    119140sampling.
    120141
    121 ===== Measuring error
     142===== Measuring Error
    122143
    123144<procedure>(absolute-error x r) -> number</procedure>
     
    141162> (absolute-error 1e-20 0.0)
    1421631e-20
    143 > (absolute-error (- 1.0 (fl 4999999/5000000)) 1/5000000)
     164> (absolute-error (- 1.0 (fp 4999999/5000000)) 1/5000000)
    1441655.751132903242251e-18
    145166</enscript>
     
    167188> (relative-error 1e-20 0.0)
    168189+inf.0
    169 > (relative-error (- 1.0 (fl 4999999/5000000)) 1/5000000)
     190> (relative-error (- 1.0 (fp 4999999/5000000)) 1/5000000)
    1701912.8755664516211255e-11
    171192</enscript>
     
    176197error for measuring the error in a flonum approximation. An even better one is
    177198error in ulps; see {{fpulp}} and {{fpulp-error}}.
     199
     200==== (math flonum)
     201
     202Flonum functions, including high-accuracy support
     203
     204'''Note:''' This library uses the {{fp}} prefix rather than the original
     205library's {{fl}} prefix for consistency with {{(chicken flonum)}}.
     206
     207'''Note:''' This is currently missing procedures from {{racket/flonum}} that
     208the original racket library re-exports.
     209
     210'''Note:''' Currently, this module is not very well tested.
     211
     212===== Additional Flonum Fnuctions
     213
     214<procedure>(fp x) -> flonum</procedure>
     215
     216; x : number
     217
     218Equivalent to {{exact->inexact}}, but much easier to read and write.
     219
     220Examples:
     221<enscript highlight="scheme">
     222> (fp 1/2)
     2230.5
     224> (fp 0.5)
     2250.5
     226> (fp #i0.5)
     2270.5
     228</enscript>
     229
     230<procedure>(fpsgn x) -> flonum</procedure>
     231<procedure>(fpeven? x) -> flonum</procedure>
     232<procedure>(fpodd? x) -> flonum</procedure>
     233
     234; x : flonum
     235
     236Like {{sgn}}, {{even?}}, and {{odd?}}, but restricted to flonum input.
     237
     238Examples:
     239
     240<enscript highlight="scheme">
     241> (map fpsgn '(-2.0 -0.0 0.0 2.0))
     242'(-1.0 0.0 0.0 1.0)
     243> (map fpeven? '(2.0 1.0 0.5))
     244'(#t #f #f)
     245> (map fpodd? '(2.0 1.0 0.5))
     246'(#f #t #f)
     247</enscript>
     248
     249<procedure>(fprational? x) -> boolean</procedure>
     250<procedure>(fpinfinite? x) -> boolean</procedure>
     251<procedure>(fpnan? x) -> boolean</procedure>
     252<procedure>(fpinteger? x) -> boolean</procedure>
     253
     254; x : flonum
     255
     256Like {{rational?}}, {{infinite?}}, {{nan?}}, and {{integer?}}, but restricted
     257to flonum input.
     258
     259<procedure>(fphypot x y) -> flonum</procedure>
     260
     261; x : flonum
     262; y : flonum
     263
     264Computes {{(fpsqrt (+ (* x x) (* y y)))}} in way that overflows only when the
     265answer is too large.
     266
     267Examples:
     268<enscript highlight="scheme">
     269> (fpsqrt (+ (* 1e+200 1e+200) (* 1e+199 1e+199)))
     270+inf.0
     271> (fphypot 1e+200 1e+199)
     2721.0049875621120889e+200
     273</enscript>
     274
     275<procedure>(fpsum xs) -> flonum</procedure>
     276
     277; xs : (list-of flonum)
     278
     279Like {{(apply + xs)}}, but incurs rounding error only once.
     280
     281Examples:
     282<enscript highlight="scheme">
     283> (+ 1.0 1e-16)
     2841.0
     285> (+ (+ 1.0 1e-16) 1e-16)
     2861.0
     287> (fpsum '(1.0 1e-16 1e-16))
     2881.0000000000000002
     289</enscript>
     290
     291The sum function does the same for heterogenous lists of reals.
     292
     293Worst-case time complexity is O(''n''^2), though the pathological inputs needed
     294to observe quadratic time are exponentially improbable and are hard to generate
     295purposely. Expected time complexity is O(''n'' log(''n'')).
     296
     297See {{fpvector-sums}} for a variant that computes all the partial sums in
     298{{xs}}.
     299
     300<procedure>(fpsinh x) -> flonum</procedure>
     301<procedure>(fpcosh x) -> flonum</procedure>
     302<procedure>(fptanh x) -> flonum</procedure>
     303
     304; x : flonum
     305
     306Return the [[https://en.wikipedia.org/wiki/Hyperbolic_function|hyperbolic sine,
     307cosine, and tangent]] of {{x}}, respectively.
     308
     309Maximum observed error is 2 ulps (see {{fpulp}}), making these functions
     310(currently) much more accurate than their {{(math base)}} counterparts
     311('''Note:''' currently missing). They also return sensible values on the
     312largest possible domain.
     313
     314<procedure>(fpsinh x) -> flonum</procedure>
     315<procedure>(fpcosh x) -> flonum</procedure>
     316<procedure>(fptanh x) -> flonum</procedure>
     317
     318; x : flonum
     319
     320Return the [[https://en.wikipedia.org/wiki/Inverse_hyperbolic_function|inverse
     321hyperbolic sine, cosine, and tangent]] of {{y}}, respectively.
     322
     323These functions are as robust and accurate as their corresponding inverses.
     324
     325<procedure>(fpfactorial n) -> flonum</procedure>
     326<procedure>(fpbinomial n k) -> flonum</procedure>
     327<procedure>(fppermutations n k) -> flonum</procedure>
     328<procedure>(fpmultinomial n ks) -> flonum</procedure>
     329
     330; n : flonum
     331; k : flonum
     332; ks : (list-of flonum)
     333
     334Like {{(fp (factorial (fp->exact-integer n)))}} and so on, but computed in
     335constant time. Also, these return {{+nan.0}} instead of raising exceptions.
     336
     337For factorial-like functions that return sensible values for non-integers, see
     338{{log-gamma}} and {{log-beta}} ('''Note:''' currently missing).
     339
     340<procedure>(fplog1p x) -> flonum</procedure>
     341<procedure>(fpexpm1 x) -> flonum</procedure>
     342
     343; x : flonum
     344
     345Like {{(fplog (+ 1.0 x))}} and {{(- (fpexp x) 1.0)}}, but accurate when {{x}}
     346is small (within 1 ulp - see {{fpulp}}).
     347
     348For example, one difficult input for {{(fplog (+ 1.0 x))}} and {{(- (fpexp x) 1.0)}} is
     349{{x = 1e-14}}, which {{fplog1p}} and {{fpexpm1}} compute correctly:
     350
     351<enscript highlight="scheme">
     352> (fplog (+ 1.0 1e-14))
     3539.992007221626358e-15
     354> (fplog1p 1e-14)
     3559.99999999999995e-15
     356> (- (fpexp 1e-14) 1.0)
     3579.992007221626409e-15
     358> (fpexpm1 1e-14)
     3591.0000000000000049e-14
     360</enscript>
     361
     362These functions are mutual inverses:
     363
     364[[image:https://docs.racket-lang.org/math/pict_2.png]]
     365
     366Notice that both graphs pass through the origin. Thus, inputs close to 0.0,
     367around which flonums are particularly dense, result in outputs that are also
     368close to 0.0. Further, both functions are approximately the identity function
     369near 0.0, so the output density is approximately the same.
     370
     371Many flonum functions defined in terms of {{fplog}} and {{fpexp}} become much
     372more accurate when their defining expressions are put in terms of {{fplog1p}}
     373and {{fpexpm1}}. The functions exported by this module and by
     374math/special-functions use them extensively.
     375
     376One notorious culprit is {{(fpexpt (- 1.0 x) y)}}, when {{x}} is near 0.0.
     377Computing it directly too often results in the wrong answer:
     378
     379<enscript highlight="scheme">
     380> (fpexpt (- 1.0 1e-20) 1e+20)
     3811.0
     382</enscript>
     383
     384We should expect that multiplying a number just less than 1.0 by itself that
     385many times would result in something less than 1.0. The problem comes from
     386subtracting such a small number from 1.0 in the first place:
     387
     388<enscript highlight="scheme">
     389> (- 1.0 1e-20)
     3901.0
     391</enscript>
     392
     393Fortunately, we can compute this correctly by putting the expression in terms
     394of {{fplog1p}}, which avoids the error-prone subtraction:
     395
     396<enscript highlight="scheme">
     397> (fpexp (* 1e+20 (fplog1p (- 1e-20))))
     3980.36787944117144233
     399</enscript>
     400
     401See {{fpexpt1p}}, which is more accurate still.
     402
     403<procedure>(fpexpt1p x y) -> flonum</procedure>
     404
     405; x : flonum
     406; y : flonum
     407
     408Like {{(fpexpt (+ 1.0 x) y)}}, but accurate for any {{x}} and {{y}}.
     409
     410<procedure>(fpexpt+ x1 x2 y) -> flonum</procedure>
     411
     412; x1 : flonum
     413; x2 : flonum
     414; y : flonum
     415
     416Like {{(fpexpt (+ x1 x2) y)}}, but more accurate.
     417
     418<procedure>(fpexp2 x) -> flonum</procedure>
     419
     420; x : flonum
     421
     422Equivalent to {{fpexpt 2.0 x}}, but faster when {{x}} is an integer.
     423
     424<procedure>(fplog2 x) -> flonum</procedure>
     425
     426; x : flonum
     427
     428Computes the base-2 log of {{x}} more accurately than {{(/ (fplog x) (fplog
     4292.0))}}. In particular, {{(fplog2 x)}} is correct for any power of two {{x}}.
     430
     431Examples:
     432<enscript highlight="scheme">
     433> (fplog2 4.5)
     4342.169925001442312
     435> (/ (fplog (fpexp2 -1066.0)) (fplog 2.0))
     436-1066.0000000000002
     437> (fplog2 (fpexp2 -1066.0))
     438-1066.0
     439</enscript>
     440
     441Maximum observed error is 0.5006 ulps (see {{fpulp}}), but is almost always no
     442more than 0.5 (i.e. it is almost always correct).
     443
     444<procedure>(fplogb b x) -> flonum</procedure>
     445
     446; b : flonum
     447; x : flonum
     448
     449Computes the base-b log of {{x}} more accurately than {{(/ (fplog x) (fplog
     450b))}}, and handles limit values correctly.
     451
     452Maximum observed error is 2.1 ulps (see {{fpulp}}), but is usually less than
     4530.7 (i.e. near rounding error).
     454
     455Except possibly at limit values (such as 0.0 and {{+inf.0}}, and {{b = 1.0}})
     456and except when the inner expression underflows or overflows, {{fplogb}}
     457approximately meets these identities for {{b > 0.0}}:
     458
     459* Left inverse: {{(fplogb b (fpexpt b y)) = y}}
     460* Right inverse: {{(fpexpt b (fplogb b x)) = x}} when {{x > 0.0}}
     461
     462Unlike with {{fpexpt}}, there is no standard for {{fplogb}}’s behavior at limit
     463values. Fortunately, deriving the following rules (applied in order) is not
     464prohibitively difficult.
     465
     466<table>
     467<tr> <th>Case                            </th> <th style="text-align:center">Condition      </th> <th style="text-align:right">Value     </th> </tr>
     468<tr> <td>{{(fplogb b 1.0)}}              </td> <td align="center">                          </td> <td align="right">{{0.0}}   </td> </tr>
     469<tr> <td>{{(fplogb 1.0 x)}}              </td> <td align="center">                          </td> <td align="right">{{+nan.0}}</td> </tr>
     470<tr> <td>{{(fplogb b x)}}                </td> <td align="center">{{b < 0.0}} or {{x < 0.0}}</td> <td align="right">{{+nan.0}}</td> </tr>
     471<tr> <td>''Double limits''               </td> <td align="center">                          </td> <td align="right">          </td> </tr>
     472<tr> <td>{{(fplogb 0.0 0.0)}}            </td> <td align="center">                          </td> <td align="right">{{+inf.0}}</td> </tr>
     473<tr> <td>{{(fplogb 0.0 +inf.0)}}         </td> <td align="center">                          </td> <td align="right">{{-inf.0}}</td> </tr>
     474<tr> <td>{{(fplogb +inf 0.0)}}           </td> <td align="center">                          </td> <td align="right">{{-inf.0}}</td> </tr>
     475<tr> <td>{{(fplogb +inf +inf.0)}}        </td> <td align="center">                          </td> <td align="right">{{+inf.0}}</td> </tr>
     476<tr> <td>''Limits with respect to {{b}}''</td> <td align="center">                          </td> <td align="right">          </td> </tr>
     477<tr> <td>{{(fplogb 0.0 x)}}              </td> <td align="center">{{x < 1.0}}               </td> <td align="right">{{0.0}}   </td> </tr>
     478<tr> <td>{{(fplogb 0.0 x)}}              </td> <td align="center">{{x > 1.0}}               </td> <td align="right">{{-0.0}}  </td> </tr>
     479<tr> <td>{{(fplogb +inf.0 x)}}           </td> <td align="center">{{x > 1.0}}               </td> <td align="right">{{0.0}}   </td> </tr>
     480<tr> <td>{{(fplogb +inf.0 x)}}           </td> <td align="center">{{x < 1.0}}               </td> <td align="right">{{-0.0}}  </td> </tr>
     481<tr> <td>''Limits with respect to {{x}}''</td> <td align="center">                          </td> <td align="right">          </td> </tr>
     482<tr> <td>{{(fplogb b 0.0)}}              </td> <td align="center">{{x < 1.0}}               </td> <td align="right">{{+inf.0}}</td> </tr>
     483<tr> <td>{{(fplogb b 0.0)}}              </td> <td align="center">{{x > 1.0}}               </td> <td align="right">{{-inf.0}}</td> </tr>
     484<tr> <td>{{(fplogb b +inf.0)}}           </td> <td align="center">{{x > 1.0}}               </td> <td align="right">{{+inf.0}}</td> </tr>
     485<tr> <td>{{(fplogb b +inf.0)}}           </td> <td align="center">{{x < 1.0}}               </td> <td align="right">{{-inf.0}}</td> </tr>
     486</table>
     487
     488Most of these rules are derived by taking limits of the mathematical base-b log
     489function. Except for {{(fplogb 1.0 x)}}, when doing so gives rise to
     490ambiguities, they are resolved using {{fpexpt}}’s behavior, which follows the
     491IEEE 754 and C99 standards for {{pow}.
     492
     493For example, consider {{(fplogb 0.0 0.0)}}. Taking an interated limit, we get ∞
     494if the outer limit is with respect to {{x}}, or 0 if the outer limit is with
     495respect to {{b}}. This would normally mean {{(fplogb 0.0 0.0) = +nan.0}}.
     496
     497However, choosing {{+inf.0}} ensures that these additional left-inverse and
     498right-inverse identities hold:
     499
     500<enscript highlight="scheme">
     501(fplogb 0.0 (fpexpt 0.0 +inf.0)) = +inf.0
     502(fpexpt 0.0 (fplogb 0.0 0.0)) = 0.0
     503</enscript>
     504
     505Further, choosing 0.0 does not ensure that any additional identities hold.
     506
     507<procedure>(fpbracketed-root f a b) -> flonum</procedure>
     508
     509; f : (flonum -> flonum)
     510; a : flonum
     511; b : flonum
     512
     513Uses the [[https://en.wikipedia.org/wiki/Brent%27s_method|Brent-Dekker method]]
     514to find a floating-point root of {{f}} (a flonum {{x}} for which {{(f x)}} is very
     515near a zero crossing) between {{a}} and {{b}}. The values {{(f a)}} and {{(f
     516b)}} must have opposite signs, but a and b may be in any order.
     517
     518Examples:
     519<enscript highlight="scheme">
     520> (define (f x) (+ 1.0 (* (+ x 3.0) (sqr (- x 1.0)))))
     521> (define x0 (fpbracketed-root f -4.0 2.0))
     522> (f (fpprev x0))
     523-7.105427357601002e-15
     524> (f x0)
     5256.661338147750939e-16
     526> (fpbracketed-root f -1.0 2.0)
     527+nan.0
     528</enscript>
     529
     530Caveats:
     531
     532* There is no guarantee that {{fpbracketed-root}} will find any particular
     533  root. Moreover, future updates to its implementation could make it find
     534  different ones.
     535
     536* There is currently no guarantee that it will find the closest {{x}} to an
     537  exact root.
     538
     539* It currently runs for at most 5000 iterations.
     540
     541It usually requires far fewer iterations, especially if the initial bounds
     542{{a}} and {{b}} are tight.
     543
     544<procedure>(make-fpexpt x) -> (flonum -> flonum)</procedure>
     545
     546; x : number
     547
     548Equivalent to {{(λ (y) (fpexpt x y))}} when {{x}} is a flonum, but much more accurate
     549for large {{y}} when {{x}} cannot be exactly represented by a flonum.
     550
     551Suppose we want to compute {{π^y}}, where {{y}} is a flonum. If we use flexpt
     552with an approximation of the irrational base {{π}}, the error is low near zero,
     553but grows with distance from the origin. Using {{make-fpexpt}}, the error is near
     554rounding error everywhere.
     555
     556<procedure>(fpsqrt1pm1 x) -> flonum</procedure>
     557
     558; x : flonum
     559
     560Like {{(- (fpsqrt (+ 1.0 x)) 1.0)}}, but accurate when {{x}} is small.
     561
     562<procedure>(fplog1pmx x) -> flonum</procedure>
     563
     564; x : flonum
     565
     566Like {{(- (fplog1p x) x)}}, but accurate when {{x}} is small.
     567
     568<procedure>(fpexpsqr x) -> flonum</procedure>
     569
     570; x : flonum
     571
     572Like {{(fpexp (* x x))}}, but accurate when {{x}} is large.
     573
     574<procedure>(fpgauss x) -> flonum</procedure>
     575
     576; x : flonum
     577
     578Like {{(fpexp (- (* x x)))}}, but accurate when {{x}} is large.
     579
     580<procedure>(fpexp1p x) -> flonum</procedure>
     581
     582; x : flonum
     583
     584Like {{(fpexp (+ 1.0 x))}}, but accurate when {{x}} is near a power of 2.
     585
     586<procedure>(fpsinpix x) -> flonum</procedure>
     587<procedure>(fpcospix x) -> flonum</procedure>
     588<procedure>(fptanpix x) -> flonum</procedure>
     589
     590; x : flonum
     591
     592Like {{(fpsin (* pi x))}}, {{(fpcos (* pi x))}} and {{(fptan (* pi x))}}, respectively, but
     593accurate near roots and singularities. When {{x = (+ n 0.5)}} for some integer {{n}},
     594{{(fptanpix x) = +nan.0}}.
     595
     596<procedure>(fpcscpix x) -> flonum</procedure>
     597<procedure>(fpsecpix x) -> flonum</procedure>
     598<procedure>(fpcotpix x) -> flonum</procedure>
     599
     600; x : flonum
     601
     602Like {{(/ 1.0 (fpsinpix x))}}, {{(/ 1.0 (fpcospix x))}} and {{(/ 1.0 (fptanpix
     603x))}}, respectively, but the first two return {{+nan.0}} at singularities and
     604{{fpcotpix}} avoids a double reciprocal.
     605
     606
     607===== Log-Space Arithmetic
     608
     609It is often useful, especially when working with probabilities and probability
     610densities, to represent nonnegative numbers in ''log space'', or as the natural
     611logs of their true values. Generally, the reason is that the ''smallest''
     612positive flonum is ''too'' large.
     613
     614For example, say we want the probability density of the standard normal
     615distribution (the bell curve) at 50 standard deviations from zero:
     616
     617('''Note''': {{(math distributions)}} is still un-implemented in CHICKEN, but
     618should arrive in a later release)
     619
     620<enscript highlight="scheme">
     621> (import (math distributions))
     622> (pdf (normal-dist) 50.0)
     6230.0
     624</enscript>
     625
     626Mathematically, the density is nonzero everywhere, but the density at 50 is
     627less than {{+min.0}}. However, its density in log space, or its log-density, is
     628representable:
     629
     630<enscript highlight="scheme">
     631> (pdf (normal-dist) 50.0 #t)
     632-1250.9189385332047
     633</enscript>
     634
     635While this example may seem contrived, it is very common, when computing the
     636density of a ''vector'' of data, for the product of the densities to be too
     637small to represent directly.
     638
     639In log space, exponentiation becomes multiplication, multiplication becomes
     640addition, and addition becomes tricky. See {{lg+}} and {{lgsum}} for solutions.
     641
     642<procedure>(lg* logx logy) -> flonum</procedure>
     643<procedure>(lg/ logx logy) -> flonum</procedure>
     644<procedure>(lgprod logxs) -> flonum</procedure>
     645
     646; logx : flonum
     647; logy : flonum
     648; logxs : (list-of flonum)
     649
     650Equivalent to {{(fp+ logx logy)}}, {{(fp- logx logy)}} and {{(fpsum logxs)}},
     651respectively.
     652
     653<procedure>(lg+ logx logy) -> flonum</procedure>
     654<procedure>(lg- logx logy) -> flonum</procedure>
     655
     656; logx : flonum
     657; logy : flonum
     658
     659Like {{(fplog (+ (fpexp logx) (fpexp logy)))}} and {{(fplog (- (fpexp logx)
     660(fpexp logy)))}}, respectively, but more accurate and less prone to overflow
     661and underflow.
     662
     663When {{logy > logx}}, {{lg-}} returns {{+nan.0}}. Both functions correctly
     664treat {{-inf.0}} as log-space 0.0.
     665
     666To add more than two log-space numbers with the same guarantees, use {{lgsum}}.
     667
     668Examples:
     669<enscript highlight="scheme">
     670> (lg+ (fplog 0.5) (fplog 0.2))
     671-0.35667494393873234
     672> (fpexp (lg+ (fplog 0.5) (fplog 0.2)))
     6730.7000000000000001
     674> (lg- (fplog 0.5) (fplog 0.2))
     675-1.203972804325936
     676> (fpexp (lg- (fplog 0.5) (fplog 0.2)))
     6770.30000000000000004
     678> (lg- (fplog 0.2) (fplog 0.5))
     679+nan.0
     680</enscript>
     681
     682Though more accurate than a naive implementation, both functions are prone to
     683catastrophic cancellation (see {{fpulp-error}}) in regions where they output a
     684value close to 0.0 (or log-space 1.0). While these outputs have high relative
     685error, their absolute error is very low, and when exponentiated, nearly have
     686just rounding error. Further, catastrophic cancellation is unavoidable when
     687{{logx}} and {{logy}} themselves have error, which is by far the common case.
     688
     689These are, of course, excuses—but for floating-point research generally. There
     690are currently no reasonably fast algorithms for computing {{lg+}} and {{lg-}}
     691with low relative error. For now, if you need that kind of accuracy, use
     692{{(math bigfloat)}} ('''Note''': still unimplemented in CHICKEN).
     693
     694<procedure>(lgsum logxs) -> flonum</procedure>
     695
     696; logxs : flonum
     697
     698Like folding {{lg+}} over {{logxs}}, but more accurate. Analogous to {{fpsum}}.
     699
     700<procedure>(lg1+ logx) -> flonum</procedure>
     701<procedure>(lg1- logx) -> flonum</procedure>
     702
     703; logx : flonum
     704
     705Equivalent to {{(lg+ (fplog 1.0) logx)}} and {{(lg- (fplog 1.0) logx)}},
     706respectively, but faster.
     707
     708<procedure>(fpprobability? x [log?])</procedure>
     709
     710; x : flonum
     711; log? : boolean
     712
     713When {{log?}} is {{#f}}, returns {{#t}} when {{(<= 0.0 x 1.0)}}. When {{log?}}
     714is {{#t}}, returns {{#t}} when {{(<= -inf.0 x 0.0)}}
     715
     716Examples:
     717<enscript highlight="scheme">
     718> (fpprobability? -0.1)
     719#f
     720> (fpprobability? 0.5)
     721#t
     722> (fpprobability? +nan.0 #t)
     723#f
     724</enscript>
     725
     726===== Debugging Flonum Functions
     727
     728The following functions and constants are useful in authoring and debugging
     729flonum functions that must be accurate on the largest possible domain.
     730
     731====== Measuring Floating-Point Error
     732
     733<procedure>(fpulp x) -> flonum</procedure>
     734
     735; x : flonum
     736
     737Returns {{x}}'s ''ulp'' or '''u'''nit in '''l'''ast '''p'''lace: the magnitude
     738of the least significant bit in {{x}}.
     739
     740Examples:
     741<enscript highlight="scheme">
     742> (fpulp 1.0)
     7432.220446049250313e-16
     744> (fpulp 1e-100)
     7451.2689709186578246e-116
     746> (fpulp 1e+200)
     7471.6996415770136547e+184
     748</enscript>
     749
     750<procedure>(fpulp-error x r) -> flonum</procedure>
     751
     752; x : flonum
     753; r : number
     754
     755Returns the absolute number of ulps difference between {{x}} and {{r}}.
     756
     757For non-rational arguments such as {{+nan.0}}, {{fpulp-error}} returns 0.0 if
     758{{(eqv? x r)}}; otherwise it returns {{+inf.0}}.
     759
     760A flonum function with maximum error 0.5 ulps exhibits only rounding error; it
     761is ''correct''. A flonum function with maximum error no greater than a few ulps
     762is ''accurate''. Most moderately complicated flonum functions, when implemented
     763directly, seem to have over a hundred thousand ulps maximum error.
     764
     765Examples:
     766<enscript highlight="scheme">
     767> (fpulp-error 0.5 1/2)
     7680.0
     769> (fpulp-error 0.14285714285714285 1/7)
     7700.2857142857142857
     771> (fpulp-error +inf.0 +inf.0)
     7720.0
     773> (fpulp-error +inf.0 +nan.0)
     774+inf.0
     775> (fpulp-error 1e-20 0.0)
     776+inf.0
     777> (fpulp-error (- 1.0 (fl 4999999/5000000)) 1/5000000)
     778217271.6580864
     779</enscript>
     780
     781he last example subtracts two nearby flonums, the second of which had already
     782been rounded, resulting in horrendous error. This is an example of ''catastrophic
     783cancellation'. Avoid subtracting nearby flonums whenever possible. [1]
     784
     785See {{relative-error}} for a similar way to measure approximation error when
     786the approximation is not necessarily represented by a flonum.
     787
     788[1] You can make an exception if the result is to be exponentiated. If {{x}}
     789has small {{absolute-error}}, then {{(exp x)}} has small {{relative-error}} and
     790small {{fpulp-error}}
     791
     792====== Flonum Constants
     793
     794<constant>-max.0</constant>
     795<constant>-min.0</constant>
     796<constant>+min.0</constant>
     797<constant>+max.0</constant>
     798
     799The nonzero, rational flonums with maximum and minimum magnitude.
     800
     801Example:
     802<enscript highlight="scheme">
     803> (list -max.0 -min.0 +min.0 +max.0)
     804(-1.7976931348623157e+308 -5e-324 5e-324 1.7976931348623157e+308)
     805</enscript>
     806
     807<constant>epsilon.0</constant>
     808
     809The smallest flonum that can be added to 1.0 to yield a larger number, or the
     810magnitude of the least significant bit in 1.0.
     811
     812Examples:
     813<enscript highlight="scheme">
     814> epsilon.0
     8152.220446049250313e-16
     816> (fpulp 1.0)
     8172.220446049250313e-16
     818</enscript>
     819
     820Epsilon is often used in stopping conditions for iterative or additive
     821approximation methods. For example, the following function uses it to stop
     822Newton’s method to compute square roots. (Please do not assume this example is
     823robust.)
     824
     825<enscript highlight="scheme">
     826(define (newton-sqrt x)
     827  (let loop ([y  (* 0.5 x)])
     828    (define dy (/ (- x (sqr y)) (* 2.0 y)))
     829    (if (<= (abs dy) (abs (* 0.5 epsilon.0 y)))
     830        (+ y dy)
     831        (loop (+ y dy)))))
     832</enscript>
     833
     834When {{(<= (abs dy) (abs (* 0.5 epsilon.0 y)))}}, adding dy to y rarely results
     835in a different flonum. The value 0.5 can be changed to allow looser
     836approximations. This is a good idea when the approximation does not have to be
     837as close as possible (e.g. it is only a starting point for another
     838approximation method), or when the computation of dy is known to be inaccurate.
     839
     840Approximation error is often understood in terms of relative error in epsilons.
     841Number of epsilons relative error roughly corresponds with error in ulps,
     842except when the approximation is subnormal.
     843
     844====== Low-Level Flonum Operations
     845
     846<procedure>(flonum->bit-field x) -> integer</procedure>
     847
     848; x : flonum
     849
     850Returns the bits comprising {{x}} as an integer. A convenient shortcut for
     851composing {{integer-bytes->integer}} with {{real->floating-point-bytes}}
     852('''Note''': neither of these is in CHICKEN).
     853
     854Examples:
     855<enscript highlight="scheme">
     856> (number->string (flonum->bit-field -inf.0) 16)
     857"fff0000000000000"
     858> (number->string (flonum->bit-field +inf.0) 16)
     859"7ff0000000000000"
     860> (number->string (flonum->bit-field -0.0) 16)
     861"8000000000000000"
     862> (number->string (flonum->bit-field 0.0) 16)
     863"0"
     864> (number->string (flonum->bit-field -1.0) 16)
     865"bff0000000000000"
     866> (number->string (flonum->bit-field 1.0) 16)
     867"3ff0000000000000"
     868> (number->string (flonum->bit-field +nan.0) 16)
     869"7ff8000000000000"
     870</enscript>
     871
     872<procedure>(bit-field->flonum i) -> flonum</procedure>
     873
     874; i : integer
     875
     876The inverse of {{flonum->bit-field}}.
     877
     878<procedure>(flonum->ordinal x) -> integer</procedure>
     879
     880; x : flonum
     881
     882Returns the signed ordinal index of {{x}} in a total order over flonums.
     883
     884When inputs are not {{+nan.0}}, this function is monotone and symmetric; i.e.
     885if {{(fp<= x y)}} then {{(<= (flonum->ordinal x) (flonum->ordinal y))}}, and
     886{{(= (flonum->ordinal (- x)) (- (flonum->ordinal x)))}}.
     887
     888Examples:
     889<enscript highlight="scheme">
     890> (flonum->ordinal -inf.0)
     891-9218868437227405312
     892> (flonum->ordinal +inf.0)
     8939218868437227405312
     894> (flonum->ordinal -0.0)
     8950
     896> (flonum->ordinal 0.0)
     8970
     898> (flonum->ordinal -1.0)
     899-4607182418800017408
     900> (flonum->ordinal 1.0)
     9014607182418800017408
     902> (flonum->ordinal +nan.0)
     9039221120237041090560
     904</enscript>
     905
     906These properties mean that {{flonum->ordinal}} does not distinguish -0.0 and 0.0.
     907
     908<procedure>(ordinal->flonum i) -> flonum</procedure>
     909
     910- i : integer
     911
     912The inverse of {{flonum->ordinal}}.
     913
     914<procedure>(flonums-between x y) -> integer</procedure>
     915
     916; x : flonum
     917; y : flonum
     918
     919Returns the number of flonums between x and y, excluding one endpoint.
     920Equivalent to {{(- (flonum->ordinal y) (flonum->ordinal x))}}.
     921
     922Examples:
     923<enscript highlight="scheme">
     924> (flonums-between 0.0 1.0)
     9254607182418800017408
     926> (flonums-between 1.0 2.0)
     9274503599627370496
     928> (flonums-between 2.0 3.0)
     9292251799813685248
     930> (flonums-between 1.0 +inf.0)
     9314611686018427387904
     932</enscript>
     933
     934<procedure>(fpstep x n) -> flonum</procedure>
     935
     936; x : flonum
     937; n : integer
     938
     939Returns the flonum {{n}} flonums away from {{x}}, according to
     940{{flonum->ordinal}}. If {{x}} is {{+nan.0}}, returns {{+nan.0}}.
     941
     942Examples:
     943<enscript highlight="scheme">
     944> (fpstep 0.0 1)
     9455e-324
     946> (fpstep (fpstep 0.0 1) -1)
     9470.0
     948> (fpstep 0.0 -1)
     949-5e-324
     950> (fpstep +inf.0 1)
     951+inf.0
     952> (fpstep +inf.0 -1)
     9531.7976931348623157e+308
     954> (fpstep -inf.0 -1)
     955-inf.0
     956> (fpstep -inf.0 1)
     957-1.7976931348623157e+308
     958> (fpstep +nan.0 1000)
     959+nan.0
     960</enscript>
     961
     962<procedure>(fpnext x) -> flonum</procedure>
     963<procedure>(fpprev x) -> flonum</procedure>
     964
     965; x : flonum
     966
     967Equivalent to {{(flstep x 1)}} and {{(flstep x -1)}}, respectively.
     968
     969<procedure>(fpsubnormal? x) -> boolean</procedure>
     970
     971; x : flonum
     972
     973Returns {{#t}} when {{x}} is a
     974[[https://en.wikipedia.org/wiki/Denormal_number|subnormal number]].
     975
     976Though flonum operations on subnormal numbers are still often implemented by
     977software exception handling, the situation is improving. Robust flonum
     978functions should handle subnormal inputs correctly, and reduce error in outputs
     979as close to zero ulps as possible (see {{fpulp}}).
     980
     981<constant>-max-subnormal.0</constant>
     982<constant>+max-subnormal.0</constant>
     983
     984The maximum positive and negative subnormal flonums. A flonum {{x}} is
     985subnormal when it is not zero and {{(<= (abs x) +max-subnormal.0)}}.
     986
     987Example:
     988<enscript highlight="scheme">
     989> +max-subnormal.0
     9902.225073858507201e-308
     991</enscript>
     992
     993===== Double-Double Operations
     994
     995For extra precision, floating-point computations may use two nonoverlapping
     996flonums to represent a single number. Such pairs are often called
     997''double-double'' numbers. The exact sum of the pair is the number it
     998represents. (Because they are nonoverlapping, the floating-point sum is equal
     999to the largest.)
     1000
     1001For speed, especially with arithmetic operations, there is no data type for
     1002double-double numbers. They are always unboxed: given as two arguments, and
     1003received as two values. In both cases, the number with higher magnitude is
     1004first.
     1005
     1006Inputs are never checked to ensure they are sorted and nonoverlapping, but
     1007outputs are guaranteed to be sorted and nonoverlapping if inputs are.
     1008
     1009<procedure>(fp2 x [y]) -> flonum flonum</procedure>
     1010
     1011; x : number (or flonum if y passed)
     1012; y : flonum
     1013
     1014Converts a real number or the sum of two flonums into a double-double.
     1015
     1016<enscript highlight="scheme">
     1017> (fp 1/7)
     10180.14285714285714285
     1019> (relative-error (fp 1/7) 1/7)
     10205.551115123125783e-17
     1021> (define-values (x2 x1) (fp2 1/7))
     1022> (list x2 x1)
     1023(0.14285714285714285 7.93016446160826e-18)
     1024> (fp (relative-error (+ (inexact->exact x2)
     1025                         (inexact->exact x1))
     1026                      1/7))
     10273.0814879110195774e-33
     1028</enscript>
     1029
     1030Notice that the exact sum of {{x2}} and {{x1}} in the preceeding example has very low
     1031relative error.
     1032
     1033If {{x}} is not rational, {{fp2}} returns (values x 0.0).
     1034
     1035<procedure>(fp2->real x2 x1) -> number</procedure>
     1036
     1037; x2 : flonum
     1038; x1 : flonum
     1039
     1040Returns the exact sum of {{x2}} and {{x1}} if {{x2}} is rational, {{x2}}
     1041otherwise.
     1042
     1043<enscript highlight="scheme">
     1044> (define-values (x2 x1) (fp2 1/7))
     1045> (fp2->real x2 x1)
     104646359793379775246683308002939465/324518553658426726783156020576256
     1047</enscript>
     1048
     1049<procedure>(fp2? x2 x1) -> boolean</procedure>
     1050
     1051When {{x2}} is rational, returns {{#t}} when {{(fpabs x2) > (fpabs x1)}} and
     1052{{x2}} and {{x1}} are nonoverlapping. When {{x2}} is not rational, returns
     1053{{(fp= x1 0.0)}}.
     1054
     1055Examples:
     1056<enscript highlight="scheme">
     1057> (define-values (x2 x1) (fl2 1/7))
     1058> (fl2? x2 x1)
     1059#t
     1060> (fl2? 0.14285714285714285 0.07692307692307693)
     1061#f
     1062> (fl2? +inf.0 0.0001)
     1063#f
     1064</enscript>
     1065
     1066This function is quite slow, so it is used only for testing.
     1067
     1068<procedure>(fp+/error x y) -> flonum flonum</procedure>
     1069<procedure>(fp-/error x y) -> flonum flonum</procedure>
     1070<procedure>(fp*/error x y) -> flonum flonum</procedure>
     1071<procedure>(fp//error x y) -> flonum flonum</procedure>
     1072<procedure>(fpsqr/error x) -> flonum flonum</procedure>
     1073<procedure>(fpsqrt/error x) -> flonum flonum</procedure>
     1074<procedure>(fpexp/error x) -> flonum flonum</procedure>
     1075<procedure>(fpexpm1/error x) -> flonum flonum</procedure>
     1076
     1077Compute the same values as {{(fp+ x y)}}, {{(fp- x y)}}, {{(fp* x y)}}, {{(fp/
     1078x y)}}, {{(fp* x x)}}, {{(fpsqrt x)}}, {{(fpexp x)}} and {{(fpexpm1 x)}}, but
     1079return the normally rounded-off low-order bits as the second value. The result
     1080is an unboxed double-double.
     1081
     1082Use these functions to generate double-double numbers directly from the results
     1083of floating-point operations.
     1084
     1085For {{fpexp/error}} and {{fpexpm1/error}}, the largest observed error is 3
     1086ulps. (See {{fp2ulp}}.) For the rest, the largest observed error is 0.5 ulps.
     1087
     1088
     1089<procedure>(fp2zero? x2 x1) -> boolean</procedure>
     1090<procedure>(fp2rational? x2 x1) -> boolean</procedure>
     1091<procedure>(fp2positive? x2 x1) -> boolean</procedure>
     1092<procedure>(fp2negative? x2 x1) -> boolean</procedure>
     1093<procedure>(fp2infinite? x2 x1) -> boolean</procedure>
     1094<procedure>(fp2nan? x2 x1) -> boolean</procedure>
     1095
     1096; x2 : flonum
     1097; x1 : flonum
     1098
     1099Like {{zero?}}, {{rational?}}, {{positive?}}, {{negative?}}, {{infinite?}} and
     1100{{nan?}}, but for double-double flonums.
     1101
     1102
     1103<procedure>(fp2+ x2 x1 y2 [y1]) -> flonum flonum</procedure>
     1104<procedure>(fp2- x2 x1 y2 [y1]) -> flonum flonum</procedure>
     1105<procedure>(fp2* x2 x1 y2 [y1]) -> flonum flonum</procedure>
     1106<procedure>(fp2/ x2 x1 y2 [y1]) -> flonum flonum</procedure>
     1107<procedure>(fp2abs x2 [x1]) -> flonum flonum</procedure>
     1108<procedure>(fp2sqr x2 [x1]) -> flonum flonum</procedure>
     1109<procedure>(fp2sqrt x2 [x1]) -> flonum flonum</procedure>
     1110
     1111; x2 : flonum
     1112; x1 : flonum
     1113; y2 : flonum
     1114; y1 : flonum
     1115
     1116Arithmetic and square root for double-double flonums. For arithmetic, error is
     1117less than 8 ulps. (See {{fp2ulp}}.) For {{fp2sqr}} and {{fp2sqrt}}, error is
     1118less than 1 ulp, and {{fp2abs}} is exact.
     1119
     1120<procedure>(fp2= x2 x1 y2 y1)</procedure>
     1121<procedure>(fp2> x2 x1 y2 y1)</procedure>
     1122<procedure>(fp2< x2 x1 y2 y1)</procedure>
     1123<procedure>(fp2>= x2 x1 y2 y1)</procedure>
     1124<procedure>(fp2<= x2 x1 y2 y1)</procedure>
     1125
     1126; x2 : flonum
     1127; x1 : flonum
     1128; y2 : flonum
     1129; y1 : flonum
     1130
     1131Comparison functions for double-double flonums.
     1132
     1133
     1134<procedure>(fp2exp x2 x1)</procedure>
     1135<procedure>(fp2log x2 x1)</procedure>
     1136<procedure>(fp2expm1 x2 x1)</procedure>
     1137<procedure>(fp2log1p x2 x1)</procedure>
     1138
     1139Like {{fpexp}}, {{fplog}}, {{fpexpm1}} and {{fplog1p}}, but for double-double
     1140flonums.
     1141
     1142For {{fp2exp}} and {{fp2expm1}}, error is less than 3 ulps. (See {{flpulp}}.)
     1143For {{fp2log}} and {{}fp2log1p}}, error is less than 2 ulps.
     1144
     1145====== Debugging Double-Double Functions
     1146
     1147<procedure>(fp2ulp x2 x1) -> flonum</procedure>
     1148<procedure>(fp2ulp-error x2 x1 r) -> flonum</procedure>
     1149
     1150; x2 : flonum
     1151; x1 : flonum
     1152; r : flonum
     1153
     1154Like {{fpulp}} and {{fpulp-error}}, but for double-double flonums.
     1155
     1156The unit in last place of a double-double is that of the higher-order of the
     1157pair, shifted 52 bits right.
     1158
     1159Examples:
     1160
     1161<enscript highlight="scheme">
     1162> (fp2ulp 1.0 0.0)
     11634.930380657631324e-32
     1164> (let-values ([(x2 x1)  (fl2 1/7)])
     1165    (fp2ulp-error x2 x1 1/7))
     11660.07142857142857142
     1167</enscript>
     1168
     1169<constant>+max.hi</constant>
     1170<constant>+max.lo</constant>
     1171<constant>-max.hi</constant>
     1172<constant>-max.lo</constant>
     1173
     1174The maximum-magnitude, unboxed double-double flonums.
     1175
     1176<constant>+max-subnormal.hi</constant>
     1177<constant>-max-subnormal.hi</constant>
     1178
     1179The high-order flonum of the maximum-magnitude, subnormal double-double
     1180flonums.
     1181
     1182<enscript highlight="scheme">
     1183> +max-subnormal.0
     11842.225073858507201e-308
     1185> +max-subnormal.hi
     11861.0020841800044864e-292
     1187</enscript>
     1188
     1189Try to avoid computing with double-doubles in the subnormal range in
     1190intermediate computations.
     1191
     1192====== Low-level Double-Double Operations
     1193
     1194The following syntactic forms are fast versions of functions like
     1195{{fp+/error}}. They are fast because they make assumptions about the magnitudes
     1196of and relationships between their arguments, and do not handle non-rational
     1197double-double flonums properly.
     1198
     1199<syntax>(fast-mono-fp+/error x y)</syntax>
     1200<syntax>(fast-mono-fp-/error x y)</syntax>
     1201
     1202Return two values: {{(fp+ x y)}} or {{(fp- x y)}}, and its rounding error. Both
     1203assume {{(fpabs x) > (fpabs y)}}. The values are unspecified when {{x}} or
     1204{{y}} is not rational.
     1205
     1206
     1207<syntax>(fast-fp+/error x y)</syntax>
     1208<syntax>(fast-fp-/error x y)</syntax>
     1209
     1210Like {{fast-mono-fp+/error}} and {{fast-mono-fp-/error}}, but do not assume
     1211{{(fpabs x) > (fpabs y)}}.
     1212
     1213<syntax>(fast-fl*/error x y)</syntax>
     1214<syntax>(fast-fl//error x y)</syntax>
     1215<syntax>(fast-flsqr/error x)</syntax>
     1216
     1217Like {{fp*/error}}, {{fp//error}} and {{fpsqr/error}}, but faster, and may
     1218return garbage when an argument is subnormal or nearly infinite.
     1219
     1220<syntax>(fpsplit x)</syntax>
     1221
     1222Returns nonoverlapping {{(values y2 y1)}}, each with 26 bits precision, with
     1223{{(fpabs y2) > (fpabs y1)}}, such that {{(fp+ y2 y1) = x}}. For {{(flabs x) >
     12241.3393857490036326e+300}}, returns {{(values +nan.0 +nan.0)}}.
     1225
     1226Used to implement double-double multiplication.
     1227
     1228===== Additional Flonum Vector Functions
     1229
     1230<procedure>(build-fpvector n proc) -> f64vector</procedure>
     1231
     1232; n : integer
     1233; proc : (fixnum -> flonum)
     1234
     1235Creates a length-n flonum vector by applying {{proc}} to the indexes from 0 to
     1236{{(- n 1)}}. Analogous to build-vector.
     1237
     1238Example:
     1239<enscript highlight="scheme">
     1240> (build-fpvector 10 fl)
     1241#f64(0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0)
     1242</enscript>
     1243
     1244<syntax>(inline-build-fpvector n proc)</syntax>
     1245
     1246; n : integer
     1247; proc : (fixnum -> flonum)
     1248
     1249Like {{build-flvector}}, but always inlined. This increases speed at the expense of
     1250code size.
     1251
     1252<syntax>(fpvector-map proc xs xss ...) -> f64vector</syntax>
     1253
     1254; proc : (flonum flonum ... -> flonum)
     1255; xs : f64vector
     1256; xss : f64vector
     1257
     1258Applies {{proc}} to the corresponding elements of {{xs}} and {{xss}}. Analogous
     1259to {{vector-map}}.
     1260
     1261The {{proc}} is meant to accept the same number of arguments as the number of its
     1262following flonum vector arguments.
     1263
     1264<syntax>(inline-fpvector-map proc xs xss ...)</syntax>
     1265
     1266; proc : (flonum flonum ... -> flonum)
     1267; xs : f64vector
     1268; xss : f64vector
     1269
     1270Like {{flvector-map}}, but always inlined.
     1271
     1272<procedure>(fpvector-copy! dest dest-start src [src-start src-end]) -> void</procedure>
     1273
     1274; dest : f64vector
     1275; dest-start : integer
     1276; src : f64vector
     1277; src-start : integer
     1278; src-end : integer
     1279
     1280Like {{vector-copy!}} but for flonum vectors.
     1281
     1282<procedure>(list->fpvector vs) -> f64vector</procedure>
     1283<procedure>(fpvector->list xs) -> (list-of flonum)</procedure>
     1284<procedure>(vector->fpvector vs) -> f64vector</procedure>
     1285<procedure>(fpvector->vector xs) -> (vector-of flonum)</procedure>
     1286
     1287; vs : (list-of number)
     1288; xs : f64vector
     1289
     1290Convert between lists and flonum vectors, and between vectors and flonum
     1291vectors.
     1292
     1293<procedure>(fpvector+ xs ys) -> f64vector</procedure>
     1294<procedure>(fpvector* xs ys) -> f64vector</procedure>
     1295<procedure>(fpvector- xs [ys]) -> f64vector</procedure>
     1296<procedure>(fpvector/ xs [ys]) -> f64vector</procedure>
     1297<procedure>(fpvector-scale xs y) -> f64vector</procedure>
     1298<procedure>(fpvector-abs xs) -> f64vector</procedure>
     1299<procedure>(fpvector-sqr xs) -> f64vector</procedure>
     1300<procedure>(fpvector-sqrt xs) -> f64vector</procedure>
     1301<procedure>(fpvector-min xs) -> f64vector</procedure>
     1302<procedure>(fpvector-max xs) -> f64vector</procedure>
     1303
     1304; xs : f64vector
     1305; ys : f64vector
     1306; y : flonum
     1307
     1308Arithmetic lifted to operate on flonum vectors.
     1309
     1310<procedure>(fpvector-sum xs) -> flonum</procedure>
     1311
     1312; xs : f64vector
     1313
     1314Like {{fpsum}}, but operates on flonum vectors. In fact, {{fpsum}} is defined in terms
     1315of {{fpvector-sum}}.
     1316
     1317<procedure>(fpvector-sums xs) -> f64vector</procedure>
     1318
     1319; xs : f64vector
     1320
     1321Computes the partial sums of the elements in {{xs}} in a way that incurs rounding
     1322error only once for each partial sum.
     1323
     1324Example:
     1325<enscript highlight="scheme">
     1326> (fpvector-sums
     1327   #f64(1.0 1e-16 1e-16 1e-16 1e-16 1e+100 -1e+100))
     1328#f64(1.0 1.0 1.0 1.0 1.0 1e+100 1.0)
     1329</enscript>
     1330
     1331Compare the same example computed by direct summation:
     1332
     1333<enscript highlight="scheme">
     1334> (import srfi-1)
     1335> (cdr
     1336   (reverse
     1337    (fold (lambda (x xs) (cons (+ x (first xs)) xs))
     1338          (list 0.0)
     1339          '(1.0 1e-16 1e-16 1e-16 1e-16 1e+100 -1e+100))))
     1340'(1.0 1.0 1.0 1.0 1.0 1e+100 0.0)
     1341</enscript>
    1781342
    1791343==== (math number-theory)
     
    4281592only be set using with-modulus. (The basic modular operators cache parameter
    4291593reads, and this restriction guarantees that the cached values are current.
    430 '''NOTE:''' I'm not entirely sure this is true for the chicken port, as a
     1594'''NOTE:''' I'm not entirely sure this is true for the CHICKEN port, as a
    4311595slightly more complicated racket syntax-case has been turned into a simple
    4321596syntax-rule for {{(parameterize ...)}})
     
    14702634[[https://docs.racket-lang.org/math/]]
    14712635
    1472 === Development status
    1473 
    1474 This egg is still largely under active development, with the exception of the
    1475 {{math.number-theory}} module, which is finished and ready for use.
    1476 
    1477 
    1478 === Implementation caveats
    1479 
    1480 * It's possible some undefined behavior may occur with arguments of the wrong
    1481   type, since a good amount of the original functions were originally defined
    1482   in typed racket, which AFAIK would catch those and throw an exception.
    1483 
    1484 * In some places the original implementation references {{unsafe-}} {{fx}} and
    1485   {{fl}} operators, but these are actually just aliased to safe operators. This
    1486   implementation just uses chicken's {{chicken.fixnum}} module, which is
    1487   unsafe. This may also lead to undefined behavior in some cases.
    1488 
    14892636=== Author
    14902637
     
    15052652=== Version History
    15062653
     2654; 0.3.1 : Ensure to export all api procedures, correct module dependencies for build
    15072655; 0.3.0 : Finish (math base) and (math flonum), bug and typo fixes, credit original authors
    15082656; 0.2.3 : Fix broken .egg file
Note: See TracChangeset for help on using the changeset viewer.