Changeset 39960 in project
 Timestamp:
 04/11/21 01:04:05 (5 weeks ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

wiki/eggref/5/math
r39757 r39960 5 5 === Introduction 6 6 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. 7 math is a CHICKEN port of racket's [[https://docs.racketlang.org/math/math]] 8 library. 9 10 The following documentation is largely a direct copy of the racket 11 documentation, reformatted for svnwiki and tweaked for the CHICKEN 12 implementation. 13 14 === Development status 15 16 This egg is still largely under active development. There may be missing 17 modules 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. 11 29 12 30 === Modules … … 15 33 16 34 Constants and elementary functions 35 36 '''Note:''' This is currently missing procedures from {{racket/math}} that the 37 original racket library reexports. 17 38 18 39 ===== Constants … … 54 75 55 76 Returns {{#t}} when {{v}} is an inexact complex number. Analogous to 56 [[https://wiki.callcc.org/man/5/Module%20(chicken%20base)#flonumflonum?]] 77 {{flonum?}}. 57 78 58 79 <procedure>(number>floatcomplex x) > cplxnum</procedure> … … 75 96 > (poweroftwo? 1/2) 76 97 #t 77 > (poweroftwo? (f lnext 2.0))98 > (poweroftwo? (fpnext 2.0)) 78 99 #f 79 100 </enscript> … … 119 140 sampling. 120 141 121 ===== Measuring error142 ===== Measuring Error 122 143 123 144 <procedure>(absoluteerror x r) > number</procedure> … … 141 162 > (absoluteerror 1e20 0.0) 142 163 1e20 143 > (absoluteerror ( 1.0 (f l4999999/5000000)) 1/5000000)164 > (absoluteerror ( 1.0 (fp 4999999/5000000)) 1/5000000) 144 165 5.751132903242251e18 145 166 </enscript> … … 167 188 > (relativeerror 1e20 0.0) 168 189 +inf.0 169 > (relativeerror ( 1.0 (f l4999999/5000000)) 1/5000000)190 > (relativeerror ( 1.0 (fp 4999999/5000000)) 1/5000000) 170 191 2.8755664516211255e11 171 192 </enscript> … … 176 197 error for measuring the error in a flonum approximation. An even better one is 177 198 error in ulps; see {{fpulp}} and {{fpulperror}}. 199 200 ==== (math flonum) 201 202 Flonum functions, including highaccuracy support 203 204 '''Note:''' This library uses the {{fp}} prefix rather than the original 205 library's {{fl}} prefix for consistency with {{(chicken flonum)}}. 206 207 '''Note:''' This is currently missing procedures from {{racket/flonum}} that 208 the original racket library reexports. 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 218 Equivalent to {{exact>inexact}}, but much easier to read and write. 219 220 Examples: 221 <enscript highlight="scheme"> 222 > (fp 1/2) 223 0.5 224 > (fp 0.5) 225 0.5 226 > (fp #i0.5) 227 0.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 236 Like {{sgn}}, {{even?}}, and {{odd?}}, but restricted to flonum input. 237 238 Examples: 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 256 Like {{rational?}}, {{infinite?}}, {{nan?}}, and {{integer?}}, but restricted 257 to flonum input. 258 259 <procedure>(fphypot x y) > flonum</procedure> 260 261 ; x : flonum 262 ; y : flonum 263 264 Computes {{(fpsqrt (+ (* x x) (* y y)))}} in way that overflows only when the 265 answer is too large. 266 267 Examples: 268 <enscript highlight="scheme"> 269 > (fpsqrt (+ (* 1e+200 1e+200) (* 1e+199 1e+199))) 270 +inf.0 271 > (fphypot 1e+200 1e+199) 272 1.0049875621120889e+200 273 </enscript> 274 275 <procedure>(fpsum xs) > flonum</procedure> 276 277 ; xs : (listof flonum) 278 279 Like {{(apply + xs)}}, but incurs rounding error only once. 280 281 Examples: 282 <enscript highlight="scheme"> 283 > (+ 1.0 1e16) 284 1.0 285 > (+ (+ 1.0 1e16) 1e16) 286 1.0 287 > (fpsum '(1.0 1e16 1e16)) 288 1.0000000000000002 289 </enscript> 290 291 The sum function does the same for heterogenous lists of reals. 292 293 Worstcase time complexity is O(''n''^2), though the pathological inputs needed 294 to observe quadratic time are exponentially improbable and are hard to generate 295 purposely. Expected time complexity is O(''n'' log(''n'')). 296 297 See {{fpvectorsums}} 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 306 Return the [[https://en.wikipedia.org/wiki/Hyperbolic_functionhyperbolic sine, 307 cosine, and tangent]] of {{x}}, respectively. 308 309 Maximum 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 312 largest 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 320 Return the [[https://en.wikipedia.org/wiki/Inverse_hyperbolic_functioninverse 321 hyperbolic sine, cosine, and tangent]] of {{y}}, respectively. 322 323 These 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 : (listof flonum) 333 334 Like {{(fp (factorial (fp>exactinteger n)))}} and so on, but computed in 335 constant time. Also, these return {{+nan.0}} instead of raising exceptions. 336 337 For factoriallike functions that return sensible values for nonintegers, see 338 {{loggamma}} and {{logbeta}} ('''Note:''' currently missing). 339 340 <procedure>(fplog1p x) > flonum</procedure> 341 <procedure>(fpexpm1 x) > flonum</procedure> 342 343 ; x : flonum 344 345 Like {{(fplog (+ 1.0 x))}} and {{( (fpexp x) 1.0)}}, but accurate when {{x}} 346 is small (within 1 ulp  see {{fpulp}}). 347 348 For example, one difficult input for {{(fplog (+ 1.0 x))}} and {{( (fpexp x) 1.0)}} is 349 {{x = 1e14}}, which {{fplog1p}} and {{fpexpm1}} compute correctly: 350 351 <enscript highlight="scheme"> 352 > (fplog (+ 1.0 1e14)) 353 9.992007221626358e15 354 > (fplog1p 1e14) 355 9.99999999999995e15 356 > ( (fpexp 1e14) 1.0) 357 9.992007221626409e15 358 > (fpexpm1 1e14) 359 1.0000000000000049e14 360 </enscript> 361 362 These functions are mutual inverses: 363 364 [[image:https://docs.racketlang.org/math/pict_2.png]] 365 366 Notice that both graphs pass through the origin. Thus, inputs close to 0.0, 367 around which flonums are particularly dense, result in outputs that are also 368 close to 0.0. Further, both functions are approximately the identity function 369 near 0.0, so the output density is approximately the same. 370 371 Many flonum functions defined in terms of {{fplog}} and {{fpexp}} become much 372 more accurate when their defining expressions are put in terms of {{fplog1p}} 373 and {{fpexpm1}}. The functions exported by this module and by 374 math/specialfunctions use them extensively. 375 376 One notorious culprit is {{(fpexpt ( 1.0 x) y)}}, when {{x}} is near 0.0. 377 Computing it directly too often results in the wrong answer: 378 379 <enscript highlight="scheme"> 380 > (fpexpt ( 1.0 1e20) 1e+20) 381 1.0 382 </enscript> 383 384 We should expect that multiplying a number just less than 1.0 by itself that 385 many times would result in something less than 1.0. The problem comes from 386 subtracting such a small number from 1.0 in the first place: 387 388 <enscript highlight="scheme"> 389 > ( 1.0 1e20) 390 1.0 391 </enscript> 392 393 Fortunately, we can compute this correctly by putting the expression in terms 394 of {{fplog1p}}, which avoids the errorprone subtraction: 395 396 <enscript highlight="scheme"> 397 > (fpexp (* 1e+20 (fplog1p ( 1e20)))) 398 0.36787944117144233 399 </enscript> 400 401 See {{fpexpt1p}}, which is more accurate still. 402 403 <procedure>(fpexpt1p x y) > flonum</procedure> 404 405 ; x : flonum 406 ; y : flonum 407 408 Like {{(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 416 Like {{(fpexpt (+ x1 x2) y)}}, but more accurate. 417 418 <procedure>(fpexp2 x) > flonum</procedure> 419 420 ; x : flonum 421 422 Equivalent to {{fpexpt 2.0 x}}, but faster when {{x}} is an integer. 423 424 <procedure>(fplog2 x) > flonum</procedure> 425 426 ; x : flonum 427 428 Computes the base2 log of {{x}} more accurately than {{(/ (fplog x) (fplog 429 2.0))}}. In particular, {{(fplog2 x)}} is correct for any power of two {{x}}. 430 431 Examples: 432 <enscript highlight="scheme"> 433 > (fplog2 4.5) 434 2.169925001442312 435 > (/ (fplog (fpexp2 1066.0)) (fplog 2.0)) 436 1066.0000000000002 437 > (fplog2 (fpexp2 1066.0)) 438 1066.0 439 </enscript> 440 441 Maximum observed error is 0.5006 ulps (see {{fpulp}}), but is almost always no 442 more 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 449 Computes the baseb log of {{x}} more accurately than {{(/ (fplog x) (fplog 450 b))}}, and handles limit values correctly. 451 452 Maximum observed error is 2.1 ulps (see {{fpulp}}), but is usually less than 453 0.7 (i.e. near rounding error). 454 455 Except possibly at limit values (such as 0.0 and {{+inf.0}}, and {{b = 1.0}}) 456 and except when the inner expression underflows or overflows, {{fplogb}} 457 approximately 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 462 Unlike with {{fpexpt}}, there is no standard for {{fplogb}}âs behavior at limit 463 values. Fortunately, deriving the following rules (applied in order) is not 464 prohibitively difficult. 465 466 <table> 467 <tr> <th>Case </th> <th style="textalign:center">Condition </th> <th style="textalign: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 488 Most of these rules are derived by taking limits of the mathematical baseb log 489 function. Except for {{(fplogb 1.0 x)}}, when doing so gives rise to 490 ambiguities, they are resolved using {{fpexpt}}âs behavior, which follows the 491 IEEE 754 and C99 standards for {{pow}. 492 493 For example, consider {{(fplogb 0.0 0.0)}}. Taking an interated limit, we get â 494 if the outer limit is with respect to {{x}}, or 0 if the outer limit is with 495 respect to {{b}}. This would normally mean {{(fplogb 0.0 0.0) = +nan.0}}. 496 497 However, choosing {{+inf.0}} ensures that these additional leftinverse and 498 rightinverse 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 505 Further, choosing 0.0 does not ensure that any additional identities hold. 506 507 <procedure>(fpbracketedroot f a b) > flonum</procedure> 508 509 ; f : (flonum > flonum) 510 ; a : flonum 511 ; b : flonum 512 513 Uses the [[https://en.wikipedia.org/wiki/Brent%27s_methodBrentDekker method]] 514 to find a floatingpoint root of {{f}} (a flonum {{x}} for which {{(f x)}} is very 515 near a zero crossing) between {{a}} and {{b}}. The values {{(f a)}} and {{(f 516 b)}} must have opposite signs, but a and b may be in any order. 517 518 Examples: 519 <enscript highlight="scheme"> 520 > (define (f x) (+ 1.0 (* (+ x 3.0) (sqr ( x 1.0))))) 521 > (define x0 (fpbracketedroot f 4.0 2.0)) 522 > (f (fpprev x0)) 523 7.105427357601002e15 524 > (f x0) 525 6.661338147750939e16 526 > (fpbracketedroot f 1.0 2.0) 527 +nan.0 528 </enscript> 529 530 Caveats: 531 532 * There is no guarantee that {{fpbracketedroot}} 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 541 It usually requires far fewer iterations, especially if the initial bounds 542 {{a}} and {{b}} are tight. 543 544 <procedure>(makefpexpt x) > (flonum > flonum)</procedure> 545 546 ; x : number 547 548 Equivalent to {{(Î» (y) (fpexpt x y))}} when {{x}} is a flonum, but much more accurate 549 for large {{y}} when {{x}} cannot be exactly represented by a flonum. 550 551 Suppose we want to compute {{Ï^y}}, where {{y}} is a flonum. If we use flexpt 552 with an approximation of the irrational base {{Ï}}, the error is low near zero, 553 but grows with distance from the origin. Using {{makefpexpt}}, the error is near 554 rounding error everywhere. 555 556 <procedure>(fpsqrt1pm1 x) > flonum</procedure> 557 558 ; x : flonum 559 560 Like {{( (fpsqrt (+ 1.0 x)) 1.0)}}, but accurate when {{x}} is small. 561 562 <procedure>(fplog1pmx x) > flonum</procedure> 563 564 ; x : flonum 565 566 Like {{( (fplog1p x) x)}}, but accurate when {{x}} is small. 567 568 <procedure>(fpexpsqr x) > flonum</procedure> 569 570 ; x : flonum 571 572 Like {{(fpexp (* x x))}}, but accurate when {{x}} is large. 573 574 <procedure>(fpgauss x) > flonum</procedure> 575 576 ; x : flonum 577 578 Like {{(fpexp ( (* x x)))}}, but accurate when {{x}} is large. 579 580 <procedure>(fpexp1p x) > flonum</procedure> 581 582 ; x : flonum 583 584 Like {{(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 592 Like {{(fpsin (* pi x))}}, {{(fpcos (* pi x))}} and {{(fptan (* pi x))}}, respectively, but 593 accurate 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 602 Like {{(/ 1.0 (fpsinpix x))}}, {{(/ 1.0 (fpcospix x))}} and {{(/ 1.0 (fptanpix 603 x))}}, respectively, but the first two return {{+nan.0}} at singularities and 604 {{fpcotpix}} avoids a double reciprocal. 605 606 607 ===== LogSpace Arithmetic 608 609 It is often useful, especially when working with probabilities and probability 610 densities, to represent nonnegative numbers in ''log space'', or as the natural 611 logs of their true values. Generally, the reason is that the ''smallest'' 612 positive flonum is ''too'' large. 613 614 For example, say we want the probability density of the standard normal 615 distribution (the bell curve) at 50 standard deviations from zero: 616 617 ('''Note''': {{(math distributions)}} is still unimplemented in CHICKEN, but 618 should arrive in a later release) 619 620 <enscript highlight="scheme"> 621 > (import (math distributions)) 622 > (pdf (normaldist) 50.0) 623 0.0 624 </enscript> 625 626 Mathematically, the density is nonzero everywhere, but the density at 50 is 627 less than {{+min.0}}. However, its density in log space, or its logdensity, is 628 representable: 629 630 <enscript highlight="scheme"> 631 > (pdf (normaldist) 50.0 #t) 632 1250.9189385332047 633 </enscript> 634 635 While this example may seem contrived, it is very common, when computing the 636 density of a ''vector'' of data, for the product of the densities to be too 637 small to represent directly. 638 639 In log space, exponentiation becomes multiplication, multiplication becomes 640 addition, 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 : (listof flonum) 649 650 Equivalent to {{(fp+ logx logy)}}, {{(fp logx logy)}} and {{(fpsum logxs)}}, 651 respectively. 652 653 <procedure>(lg+ logx logy) > flonum</procedure> 654 <procedure>(lg logx logy) > flonum</procedure> 655 656 ; logx : flonum 657 ; logy : flonum 658 659 Like {{(fplog (+ (fpexp logx) (fpexp logy)))}} and {{(fplog ( (fpexp logx) 660 (fpexp logy)))}}, respectively, but more accurate and less prone to overflow 661 and underflow. 662 663 When {{logy > logx}}, {{lg}} returns {{+nan.0}}. Both functions correctly 664 treat {{inf.0}} as logspace 0.0. 665 666 To add more than two logspace numbers with the same guarantees, use {{lgsum}}. 667 668 Examples: 669 <enscript highlight="scheme"> 670 > (lg+ (fplog 0.5) (fplog 0.2)) 671 0.35667494393873234 672 > (fpexp (lg+ (fplog 0.5) (fplog 0.2))) 673 0.7000000000000001 674 > (lg (fplog 0.5) (fplog 0.2)) 675 1.203972804325936 676 > (fpexp (lg (fplog 0.5) (fplog 0.2))) 677 0.30000000000000004 678 > (lg (fplog 0.2) (fplog 0.5)) 679 +nan.0 680 </enscript> 681 682 Though more accurate than a naive implementation, both functions are prone to 683 catastrophic cancellation (see {{fpulperror}}) in regions where they output a 684 value close to 0.0 (or logspace 1.0). While these outputs have high relative 685 error, their absolute error is very low, and when exponentiated, nearly have 686 just rounding error. Further, catastrophic cancellation is unavoidable when 687 {{logx}} and {{logy}} themselves have error, which is by far the common case. 688 689 These are, of course, excusesâbut for floatingpoint research generally. There 690 are currently no reasonably fast algorithms for computing {{lg+}} and {{lg}} 691 with 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 698 Like 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 705 Equivalent to {{(lg+ (fplog 1.0) logx)}} and {{(lg (fplog 1.0) logx)}}, 706 respectively, but faster. 707 708 <procedure>(fpprobability? x [log?])</procedure> 709 710 ; x : flonum 711 ; log? : boolean 712 713 When {{log?}} is {{#f}}, returns {{#t}} when {{(<= 0.0 x 1.0)}}. When {{log?}} 714 is {{#t}}, returns {{#t}} when {{(<= inf.0 x 0.0)}} 715 716 Examples: 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 728 The following functions and constants are useful in authoring and debugging 729 flonum functions that must be accurate on the largest possible domain. 730 731 ====== Measuring FloatingPoint Error 732 733 <procedure>(fpulp x) > flonum</procedure> 734 735 ; x : flonum 736 737 Returns {{x}}'s ''ulp'' or '''u'''nit in '''l'''ast '''p'''lace: the magnitude 738 of the least significant bit in {{x}}. 739 740 Examples: 741 <enscript highlight="scheme"> 742 > (fpulp 1.0) 743 2.220446049250313e16 744 > (fpulp 1e100) 745 1.2689709186578246e116 746 > (fpulp 1e+200) 747 1.6996415770136547e+184 748 </enscript> 749 750 <procedure>(fpulperror x r) > flonum</procedure> 751 752 ; x : flonum 753 ; r : number 754 755 Returns the absolute number of ulps difference between {{x}} and {{r}}. 756 757 For nonrational arguments such as {{+nan.0}}, {{fpulperror}} returns 0.0 if 758 {{(eqv? x r)}}; otherwise it returns {{+inf.0}}. 759 760 A flonum function with maximum error 0.5 ulps exhibits only rounding error; it 761 is ''correct''. A flonum function with maximum error no greater than a few ulps 762 is ''accurate''. Most moderately complicated flonum functions, when implemented 763 directly, seem to have over a hundred thousand ulps maximum error. 764 765 Examples: 766 <enscript highlight="scheme"> 767 > (fpulperror 0.5 1/2) 768 0.0 769 > (fpulperror 0.14285714285714285 1/7) 770 0.2857142857142857 771 > (fpulperror +inf.0 +inf.0) 772 0.0 773 > (fpulperror +inf.0 +nan.0) 774 +inf.0 775 > (fpulperror 1e20 0.0) 776 +inf.0 777 > (fpulperror ( 1.0 (fl 4999999/5000000)) 1/5000000) 778 217271.6580864 779 </enscript> 780 781 he last example subtracts two nearby flonums, the second of which had already 782 been rounded, resulting in horrendous error. This is an example of ''catastrophic 783 cancellation'. Avoid subtracting nearby flonums whenever possible. [1] 784 785 See {{relativeerror}} for a similar way to measure approximation error when 786 the 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}} 789 has small {{absoluteerror}}, then {{(exp x)}} has small {{relativeerror}} and 790 small {{fpulperror}} 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 799 The nonzero, rational flonums with maximum and minimum magnitude. 800 801 Example: 802 <enscript highlight="scheme"> 803 > (list max.0 min.0 +min.0 +max.0) 804 (1.7976931348623157e+308 5e324 5e324 1.7976931348623157e+308) 805 </enscript> 806 807 <constant>epsilon.0</constant> 808 809 The smallest flonum that can be added to 1.0 to yield a larger number, or the 810 magnitude of the least significant bit in 1.0. 811 812 Examples: 813 <enscript highlight="scheme"> 814 > epsilon.0 815 2.220446049250313e16 816 > (fpulp 1.0) 817 2.220446049250313e16 818 </enscript> 819 820 Epsilon is often used in stopping conditions for iterative or additive 821 approximation methods. For example, the following function uses it to stop 822 Newtonâs method to compute square roots. (Please do not assume this example is 823 robust.) 824 825 <enscript highlight="scheme"> 826 (define (newtonsqrt 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 834 When {{(<= (abs dy) (abs (* 0.5 epsilon.0 y)))}}, adding dy to y rarely results 835 in a different flonum. The value 0.5 can be changed to allow looser 836 approximations. This is a good idea when the approximation does not have to be 837 as close as possible (e.g. it is only a starting point for another 838 approximation method), or when the computation of dy is known to be inaccurate. 839 840 Approximation error is often understood in terms of relative error in epsilons. 841 Number of epsilons relative error roughly corresponds with error in ulps, 842 except when the approximation is subnormal. 843 844 ====== LowLevel Flonum Operations 845 846 <procedure>(flonum>bitfield x) > integer</procedure> 847 848 ; x : flonum 849 850 Returns the bits comprising {{x}} as an integer. A convenient shortcut for 851 composing {{integerbytes>integer}} with {{real>floatingpointbytes}} 852 ('''Note''': neither of these is in CHICKEN). 853 854 Examples: 855 <enscript highlight="scheme"> 856 > (number>string (flonum>bitfield inf.0) 16) 857 "fff0000000000000" 858 > (number>string (flonum>bitfield +inf.0) 16) 859 "7ff0000000000000" 860 > (number>string (flonum>bitfield 0.0) 16) 861 "8000000000000000" 862 > (number>string (flonum>bitfield 0.0) 16) 863 "0" 864 > (number>string (flonum>bitfield 1.0) 16) 865 "bff0000000000000" 866 > (number>string (flonum>bitfield 1.0) 16) 867 "3ff0000000000000" 868 > (number>string (flonum>bitfield +nan.0) 16) 869 "7ff8000000000000" 870 </enscript> 871 872 <procedure>(bitfield>flonum i) > flonum</procedure> 873 874 ; i : integer 875 876 The inverse of {{flonum>bitfield}}. 877 878 <procedure>(flonum>ordinal x) > integer</procedure> 879 880 ; x : flonum 881 882 Returns the signed ordinal index of {{x}} in a total order over flonums. 883 884 When inputs are not {{+nan.0}}, this function is monotone and symmetric; i.e. 885 if {{(fp<= x y)}} then {{(<= (flonum>ordinal x) (flonum>ordinal y))}}, and 886 {{(= (flonum>ordinal ( x)) ( (flonum>ordinal x)))}}. 887 888 Examples: 889 <enscript highlight="scheme"> 890 > (flonum>ordinal inf.0) 891 9218868437227405312 892 > (flonum>ordinal +inf.0) 893 9218868437227405312 894 > (flonum>ordinal 0.0) 895 0 896 > (flonum>ordinal 0.0) 897 0 898 > (flonum>ordinal 1.0) 899 4607182418800017408 900 > (flonum>ordinal 1.0) 901 4607182418800017408 902 > (flonum>ordinal +nan.0) 903 9221120237041090560 904 </enscript> 905 906 These 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 912 The inverse of {{flonum>ordinal}}. 913 914 <procedure>(flonumsbetween x y) > integer</procedure> 915 916 ; x : flonum 917 ; y : flonum 918 919 Returns the number of flonums between x and y, excluding one endpoint. 920 Equivalent to {{( (flonum>ordinal y) (flonum>ordinal x))}}. 921 922 Examples: 923 <enscript highlight="scheme"> 924 > (flonumsbetween 0.0 1.0) 925 4607182418800017408 926 > (flonumsbetween 1.0 2.0) 927 4503599627370496 928 > (flonumsbetween 2.0 3.0) 929 2251799813685248 930 > (flonumsbetween 1.0 +inf.0) 931 4611686018427387904 932 </enscript> 933 934 <procedure>(fpstep x n) > flonum</procedure> 935 936 ; x : flonum 937 ; n : integer 938 939 Returns the flonum {{n}} flonums away from {{x}}, according to 940 {{flonum>ordinal}}. If {{x}} is {{+nan.0}}, returns {{+nan.0}}. 941 942 Examples: 943 <enscript highlight="scheme"> 944 > (fpstep 0.0 1) 945 5e324 946 > (fpstep (fpstep 0.0 1) 1) 947 0.0 948 > (fpstep 0.0 1) 949 5e324 950 > (fpstep +inf.0 1) 951 +inf.0 952 > (fpstep +inf.0 1) 953 1.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 967 Equivalent to {{(flstep x 1)}} and {{(flstep x 1)}}, respectively. 968 969 <procedure>(fpsubnormal? x) > boolean</procedure> 970 971 ; x : flonum 972 973 Returns {{#t}} when {{x}} is a 974 [[https://en.wikipedia.org/wiki/Denormal_numbersubnormal number]]. 975 976 Though flonum operations on subnormal numbers are still often implemented by 977 software exception handling, the situation is improving. Robust flonum 978 functions should handle subnormal inputs correctly, and reduce error in outputs 979 as close to zero ulps as possible (see {{fpulp}}). 980 981 <constant>maxsubnormal.0</constant> 982 <constant>+maxsubnormal.0</constant> 983 984 The maximum positive and negative subnormal flonums. A flonum {{x}} is 985 subnormal when it is not zero and {{(<= (abs x) +maxsubnormal.0)}}. 986 987 Example: 988 <enscript highlight="scheme"> 989 > +maxsubnormal.0 990 2.225073858507201e308 991 </enscript> 992 993 ===== DoubleDouble Operations 994 995 For extra precision, floatingpoint computations may use two nonoverlapping 996 flonums to represent a single number. Such pairs are often called 997 ''doubledouble'' numbers. The exact sum of the pair is the number it 998 represents. (Because they are nonoverlapping, the floatingpoint sum is equal 999 to the largest.) 1000 1001 For speed, especially with arithmetic operations, there is no data type for 1002 doubledouble numbers. They are always unboxed: given as two arguments, and 1003 received as two values. In both cases, the number with higher magnitude is 1004 first. 1005 1006 Inputs are never checked to ensure they are sorted and nonoverlapping, but 1007 outputs 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 1014 Converts a real number or the sum of two flonums into a doubledouble. 1015 1016 <enscript highlight="scheme"> 1017 > (fp 1/7) 1018 0.14285714285714285 1019 > (relativeerror (fp 1/7) 1/7) 1020 5.551115123125783e17 1021 > (definevalues (x2 x1) (fp2 1/7)) 1022 > (list x2 x1) 1023 (0.14285714285714285 7.93016446160826e18) 1024 > (fp (relativeerror (+ (inexact>exact x2) 1025 (inexact>exact x1)) 1026 1/7)) 1027 3.0814879110195774e33 1028 </enscript> 1029 1030 Notice that the exact sum of {{x2}} and {{x1}} in the preceeding example has very low 1031 relative error. 1032 1033 If {{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 1040 Returns the exact sum of {{x2}} and {{x1}} if {{x2}} is rational, {{x2}} 1041 otherwise. 1042 1043 <enscript highlight="scheme"> 1044 > (definevalues (x2 x1) (fp2 1/7)) 1045 > (fp2>real x2 x1) 1046 46359793379775246683308002939465/324518553658426726783156020576256 1047 </enscript> 1048 1049 <procedure>(fp2? x2 x1) > boolean</procedure> 1050 1051 When {{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 1055 Examples: 1056 <enscript highlight="scheme"> 1057 > (definevalues (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 1066 This 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 1077 Compute the same values as {{(fp+ x y)}}, {{(fp x y)}}, {{(fp* x y)}}, {{(fp/ 1078 x y)}}, {{(fp* x x)}}, {{(fpsqrt x)}}, {{(fpexp x)}} and {{(fpexpm1 x)}}, but 1079 return the normally roundedoff loworder bits as the second value. The result 1080 is an unboxed doubledouble. 1081 1082 Use these functions to generate doubledouble numbers directly from the results 1083 of floatingpoint operations. 1084 1085 For {{fpexp/error}} and {{fpexpm1/error}}, the largest observed error is 3 1086 ulps. (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 1099 Like {{zero?}}, {{rational?}}, {{positive?}}, {{negative?}}, {{infinite?}} and 1100 {{nan?}}, but for doubledouble 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 1116 Arithmetic and square root for doubledouble flonums. For arithmetic, error is 1117 less than 8 ulps. (See {{fp2ulp}}.) For {{fp2sqr}} and {{fp2sqrt}}, error is 1118 less 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 1131 Comparison functions for doubledouble 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 1139 Like {{fpexp}}, {{fplog}}, {{fpexpm1}} and {{fplog1p}}, but for doubledouble 1140 flonums. 1141 1142 For {{fp2exp}} and {{fp2expm1}}, error is less than 3 ulps. (See {{flpulp}}.) 1143 For {{fp2log}} and {{}fp2log1p}}, error is less than 2 ulps. 1144 1145 ====== Debugging DoubleDouble Functions 1146 1147 <procedure>(fp2ulp x2 x1) > flonum</procedure> 1148 <procedure>(fp2ulperror x2 x1 r) > flonum</procedure> 1149 1150 ; x2 : flonum 1151 ; x1 : flonum 1152 ; r : flonum 1153 1154 Like {{fpulp}} and {{fpulperror}}, but for doubledouble flonums. 1155 1156 The unit in last place of a doubledouble is that of the higherorder of the 1157 pair, shifted 52 bits right. 1158 1159 Examples: 1160 1161 <enscript highlight="scheme"> 1162 > (fp2ulp 1.0 0.0) 1163 4.930380657631324e32 1164 > (letvalues ([(x2 x1) (fl2 1/7)]) 1165 (fp2ulperror x2 x1 1/7)) 1166 0.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 1174 The maximummagnitude, unboxed doubledouble flonums. 1175 1176 <constant>+maxsubnormal.hi</constant> 1177 <constant>maxsubnormal.hi</constant> 1178 1179 The highorder flonum of the maximummagnitude, subnormal doubledouble 1180 flonums. 1181 1182 <enscript highlight="scheme"> 1183 > +maxsubnormal.0 1184 2.225073858507201e308 1185 > +maxsubnormal.hi 1186 1.0020841800044864e292 1187 </enscript> 1188 1189 Try to avoid computing with doubledoubles in the subnormal range in 1190 intermediate computations. 1191 1192 ====== Lowlevel DoubleDouble Operations 1193 1194 The following syntactic forms are fast versions of functions like 1195 {{fp+/error}}. They are fast because they make assumptions about the magnitudes 1196 of and relationships between their arguments, and do not handle nonrational 1197 doubledouble flonums properly. 1198 1199 <syntax>(fastmonofp+/error x y)</syntax> 1200 <syntax>(fastmonofp/error x y)</syntax> 1201 1202 Return two values: {{(fp+ x y)}} or {{(fp x y)}}, and its rounding error. Both 1203 assume {{(fpabs x) > (fpabs y)}}. The values are unspecified when {{x}} or 1204 {{y}} is not rational. 1205 1206 1207 <syntax>(fastfp+/error x y)</syntax> 1208 <syntax>(fastfp/error x y)</syntax> 1209 1210 Like {{fastmonofp+/error}} and {{fastmonofp/error}}, but do not assume 1211 {{(fpabs x) > (fpabs y)}}. 1212 1213 <syntax>(fastfl*/error x y)</syntax> 1214 <syntax>(fastfl//error x y)</syntax> 1215 <syntax>(fastflsqr/error x)</syntax> 1216 1217 Like {{fp*/error}}, {{fp//error}} and {{fpsqr/error}}, but faster, and may 1218 return garbage when an argument is subnormal or nearly infinite. 1219 1220 <syntax>(fpsplit x)</syntax> 1221 1222 Returns nonoverlapping {{(values y2 y1)}}, each with 26 bits precision, with 1223 {{(fpabs y2) > (fpabs y1)}}, such that {{(fp+ y2 y1) = x}}. For {{(flabs x) > 1224 1.3393857490036326e+300}}, returns {{(values +nan.0 +nan.0)}}. 1225 1226 Used to implement doubledouble multiplication. 1227 1228 ===== Additional Flonum Vector Functions 1229 1230 <procedure>(buildfpvector n proc) > f64vector</procedure> 1231 1232 ; n : integer 1233 ; proc : (fixnum > flonum) 1234 1235 Creates a lengthn flonum vector by applying {{proc}} to the indexes from 0 to 1236 {{( n 1)}}. Analogous to buildvector. 1237 1238 Example: 1239 <enscript highlight="scheme"> 1240 > (buildfpvector 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>(inlinebuildfpvector n proc)</syntax> 1245 1246 ; n : integer 1247 ; proc : (fixnum > flonum) 1248 1249 Like {{buildflvector}}, but always inlined. This increases speed at the expense of 1250 code size. 1251 1252 <syntax>(fpvectormap proc xs xss ...) > f64vector</syntax> 1253 1254 ; proc : (flonum flonum ... > flonum) 1255 ; xs : f64vector 1256 ; xss : f64vector 1257 1258 Applies {{proc}} to the corresponding elements of {{xs}} and {{xss}}. Analogous 1259 to {{vectormap}}. 1260 1261 The {{proc}} is meant to accept the same number of arguments as the number of its 1262 following flonum vector arguments. 1263 1264 <syntax>(inlinefpvectormap proc xs xss ...)</syntax> 1265 1266 ; proc : (flonum flonum ... > flonum) 1267 ; xs : f64vector 1268 ; xss : f64vector 1269 1270 Like {{flvectormap}}, but always inlined. 1271 1272 <procedure>(fpvectorcopy! dest deststart src [srcstart srcend]) > void</procedure> 1273 1274 ; dest : f64vector 1275 ; deststart : integer 1276 ; src : f64vector 1277 ; srcstart : integer 1278 ; srcend : integer 1279 1280 Like {{vectorcopy!}} but for flonum vectors. 1281 1282 <procedure>(list>fpvector vs) > f64vector</procedure> 1283 <procedure>(fpvector>list xs) > (listof flonum)</procedure> 1284 <procedure>(vector>fpvector vs) > f64vector</procedure> 1285 <procedure>(fpvector>vector xs) > (vectorof flonum)</procedure> 1286 1287 ; vs : (listof number) 1288 ; xs : f64vector 1289 1290 Convert between lists and flonum vectors, and between vectors and flonum 1291 vectors. 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>(fpvectorscale xs y) > f64vector</procedure> 1298 <procedure>(fpvectorabs xs) > f64vector</procedure> 1299 <procedure>(fpvectorsqr xs) > f64vector</procedure> 1300 <procedure>(fpvectorsqrt xs) > f64vector</procedure> 1301 <procedure>(fpvectormin xs) > f64vector</procedure> 1302 <procedure>(fpvectormax xs) > f64vector</procedure> 1303 1304 ; xs : f64vector 1305 ; ys : f64vector 1306 ; y : flonum 1307 1308 Arithmetic lifted to operate on flonum vectors. 1309 1310 <procedure>(fpvectorsum xs) > flonum</procedure> 1311 1312 ; xs : f64vector 1313 1314 Like {{fpsum}}, but operates on flonum vectors. In fact, {{fpsum}} is defined in terms 1315 of {{fpvectorsum}}. 1316 1317 <procedure>(fpvectorsums xs) > f64vector</procedure> 1318 1319 ; xs : f64vector 1320 1321 Computes the partial sums of the elements in {{xs}} in a way that incurs rounding 1322 error only once for each partial sum. 1323 1324 Example: 1325 <enscript highlight="scheme"> 1326 > (fpvectorsums 1327 #f64(1.0 1e16 1e16 1e16 1e16 1e+100 1e+100)) 1328 #f64(1.0 1.0 1.0 1.0 1.0 1e+100 1.0) 1329 </enscript> 1330 1331 Compare the same example computed by direct summation: 1332 1333 <enscript highlight="scheme"> 1334 > (import srfi1) 1335 > (cdr 1336 (reverse 1337 (fold (lambda (x xs) (cons (+ x (first xs)) xs)) 1338 (list 0.0) 1339 '(1.0 1e16 1e16 1e16 1e16 1e+100 1e+100)))) 1340 '(1.0 1.0 1.0 1.0 1.0 1e+100 0.0) 1341 </enscript> 178 1342 179 1343 ==== (math numbertheory) … … 428 1592 only be set using withmodulus. (The basic modular operators cache parameter 429 1593 reads, and this restriction guarantees that the cached values are current. 430 '''NOTE:''' I'm not entirely sure this is true for the chickenport, as a1594 '''NOTE:''' I'm not entirely sure this is true for the CHICKEN port, as a 431 1595 slightly more complicated racket syntaxcase has been turned into a simple 432 1596 syntaxrule for {{(parameterize ...)}}) … … 1470 2634 [[https://docs.racketlang.org/math/]] 1471 2635 1472 === Development status1473 1474 This egg is still largely under active development, with the exception of the1475 {{math.numbertheory}} module, which is finished and ready for use.1476 1477 1478 === Implementation caveats1479 1480 * It's possible some undefined behavior may occur with arguments of the wrong1481 type, since a good amount of the original functions were originally defined1482 in typed racket, which AFAIK would catch those and throw an exception.1483 1484 * In some places the original implementation references {{unsafe}} {{fx}} and1485 {{fl}} operators, but these are actually just aliased to safe operators. This1486 implementation just uses chicken's {{chicken.fixnum}} module, which is1487 unsafe. This may also lead to undefined behavior in some cases.1488 1489 2636 === Author 1490 2637 … … 1505 2652 === Version History 1506 2653 2654 ; 0.3.1 : Ensure to export all api procedures, correct module dependencies for build 1507 2655 ; 0.3.0 : Finish (math base) and (math flonum), bug and typo fixes, credit original authors 1508 2656 ; 0.2.3 : Fix broken .egg file
Note: See TracChangeset
for help on using the changeset viewer.