source: project/wiki/eggref/5/math @ 39960

Last change on this file since 39960 was 39960, checked in by Diego, 2 months ago

Add (math flonum) documentation, re-order sections

File size: 71.6 KB
Line 
1[[toc:]]
2
3== math
4
5=== Introduction
6
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.
29
30=== Modules
31
32==== (math base)
33
34Constants and elementary functions
35
36'''Note:''' This is currently missing procedures from {{racket/math}} that the
37original racket library re-exports.
38
39===== Constants
40
41<constant>phi.0</constant>
42
43An approximation of φ, the [[https://en.wikipedia.org/wiki/Golden_ratio|golden ratio]].
44
45<enscript highlight="scheme">
46> phi.0
471.618033988749895
48</enscript>
49
50<constant>euler.0</constant>
51
52An approximation of ''e'', or [[https://en.wikipedia.org/wiki/E_(mathematical_constant)|Euler's number]].
53
54<enscript highlight="scheme">
55> euler.0
562.718281828459045
57> (exp 1)
582.718281828459045
59</enscript>
60
61<constant>catalan.0</constant>
62
63An approximation of ''G'', or [[https://en.wikipedia.org/wiki/Catalan's_constant|Catalan's constant]].
64
65<enscript highlight="scheme">
66> catalan.0
670.915965594177219
68</enscript>
69
70===== Functions
71
72<procedure>(float-complex? v) -> boolean</procedure>
73
74; v : any
75
76Returns {{#t}} when {{v}} is an inexact complex number. Analogous to
77{{flonum?}}.
78
79<procedure>(number->float-complex x) -> cplxnum</procedure>
80
81; x : number
82
83Returns a new complex number with a flonum real part and a flonum imaginary
84part. Analogous to {{exact->inexact}}.
85
86<procedure>(power-of-two? x) -> boolean</procedure>
87
88; x : number
89
90Returns {{#t}} when {{x}} is an integer power of 2.
91
92Examples:
93<enscript highlight="scheme">
94> (power-of-two? 1.0)
95#t
96> (power-of-two? 1/2)
97#t
98> (power-of-two? (fpnext 2.0))
99#f
100</enscript>
101
102<procedure>(asinh z) -> number</procedure>
103<procedure>(acosh z) -> number</procedure>
104<procedure>(atanh z) -> number</procedure>
105
106; z : number
107
108The inverses of {{sinh}}, {{cosh}}, and {{tanh}}.
109
110<procedure>(sum xs) -> number</procedure>
111
112; xs : (list-of number)
113
114Like {{(apply + xs)}}, but incurs rounding error only once when adding inexact
115numbers. (In fact, the inexact numbers in {{xs}} are summed separately using
116{{fpsum}}).
117
118===== Random Number Generation
119
120<procedure>(random-natural k) -> integer</procedure>
121
122; k : integer
123
124Returns a random natural number less than {{k}}, which must be positive.
125
126<procedure>(random-integer a b) -> integer</procedure>
127
128; a : integer
129; b : integer
130
131Returns a random integer n such that {{(<= a n)}} and {{(< n b)}}.
132
133<procedure>(random-bits num)</procedure>
134
135; num : integer
136
137Returns a random natural smaller than {{(expt 2 num)}}; {{num}} must be
138positive. For powers of two, this is faster than using {{random-natural}},
139which is implemented in terms of {{random-bits}}, using biased rejection
140sampling.
141
142===== Measuring Error
143
144<procedure>(absolute-error x r) -> number</procedure>
145
146; x : number
147; r : number
148
149Usually computes {{(abs (- x r))}} using exact rationals, but handles
150non-rational reals such as {{+inf.0}} specially.
151
152Examples:
153<enscript highlight="scheme">
154> (absolute-error 1/2 1/2)
1550
156> (absolute-error 0.14285714285714285 1/7)
1577.93016446160826e-18
158> (absolute-error +inf.0 +inf.0)
1590.0
160> (absolute-error +inf.0 +nan.0)
161+inf.0
162> (absolute-error 1e-20 0.0)
1631e-20
164> (absolute-error (- 1.0 (fp 4999999/5000000)) 1/5000000)
1655.751132903242251e-18
166</enscript>
167
168<procedure>(relative-error x r) -> number</procedure>
169
170; x : number
171; r : number
172
173Measures how close an approximation {{x}} is to the correct value {{r}},
174relative to the magnitude of {{r}}.
175
176This function usually computes {{(abs (/ (- x r) r))}} using exact rationals, but
177handles non-rational reals such as {{+inf.0}} specially, as well as {{r = 0}}.
178
179<enscript highlight="scheme">
180> (relative-error 1/2 1/2)
1810
182> (relative-error 0.14285714285714285 1/7)
1835.551115123125783e-17
184> (relative-error +inf.0 +inf.0)
1850.0
186> (relative-error +inf.0 +nan.0)
187+inf.0
188> (relative-error 1e-20 0.0)
189+inf.0
190> (relative-error (- 1.0 (fp 4999999/5000000)) 1/5000000)
1912.8755664516211255e-11
192</enscript>
193
194In the last two examples, relative error is high because the result is near
195zero. (Compare the same examples with {{absolute-error}}.) Because flonums are
196particularly dense near zero, this makes relative error better than absolute
197error for measuring the error in a flonum approximation. An even better one is
198error 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>
1342
1343==== (math number-theory)
1344
1345Number-theoretic functions
1346
1347===== Congruences and modular arithmetic
1348<procedure>(divides? m n) -> boolean</procedure>
1349
1350; m : integer
1351; n : integer
1352
1353Returns {{#t}} if {{m}} divides {{n}}, {{#f}} otherwise.
1354
1355Examples:
1356<enscript highlight="scheme">
1357> (divides? 2 9)
1358#f
1359> (divides? 2 8)
1360#t
1361</enscript>
1362
1363Note that 0 cannot divide anything:
1364
1365<enscript highlight="scheme">
1366> (divides? 0 5)
1367#f
1368> (divides? 0 0)
1369#f
1370</enscript>
1371
1372Practically, if {{(divides? m n)}} is {{#t}}, then {{(/ n m)}} will return an integer.
1373
1374Wikipedia: [[https://en.wikipedia.org/wiki/Divisor|Divisor]]
1375
1376<procedure>(bezout a b c ...) -> (list-of integer)</procedure>
1377
1378; a : integer
1379; b : integer
1380; c : integer
1381
1382Given integers {{a b c ...}} returns a list of integers {{(list u v w ...)}}
1383such that {{(gcd a b c ...)}} = {{(+ (* a u) (* b v) (* c w) ...)}}.
1384
1385Examples:
1386<enscript highlight="scheme">
1387> (bezout 6 15)
1388(-2 1)
1389> (+ (* -2 6) (* 1 15))
13903
1391> (gcd 6 15)
13923
1393</enscript>
1394
1395Wikipedia: [[https://en.wikipedia.org/wiki/B%C3%A9zout's_identity|Bézout's Identity]]
1396
1397<procedure>(coprime? a b ...) -> boolean</procedure>
1398
1399; a : integer
1400; b : integer
1401
1402Returns {{#t}} if the integers {{a b ...}} are coprime. Formally, a set of
1403integers is considered coprime (also called relatively prime) if their greatest
1404common divisor is 1.
1405
1406Example:
1407<enscript highlight="scheme">
1408> (coprime? 2 6 15)
1409#t
1410</enscript>
1411
1412Wikipedia: [[https://en.wikipedia.org/wiki/Coprime|Coprime]]
1413
1414<procedure>(pairwise-coprime? a b ...) -> boolean</procedure>
1415
1416; a : integer
1417; b : integer
1418
1419Returns {{#t}} if the integers {{a b ...}} are ''pairwise'' coprime, meaning that each
1420pair of integers is coprime.
1421
1422The numbers 2, 6 and 15 are coprime, but not ''pairwise'' coprime, because 6 and 15
1423share the factor 3:
1424<enscript highlight="scheme">
1425> (pairwise-coprime? 2 6 15)
1426#f
1427</enscript>
1428
1429Wikipedia:[[https://en.wikipedia.org/wiki/Pairwise_coprime|Pairwise Coprime]]
1430
1431<procedure>(solve-chinese as ns) -> integer</procedure>
1432
1433; as : (list-of integer)
1434; ns : (list-of integer)
1435
1436Given a length-k list of integers as and a length-k list of coprime moduli ns,
1437(solve-chinese as ns) returns the least natural number x that is a solution to
1438the equations
1439
1440 x = a₁ (mod n₁)
1441  ...
1442 x = aₖ (mod nₖ)
1443
1444The solution {{x}} is less than {{(* n1 ... nk)}}.
1445
1446The moduli {{ns}} must all be positive.
1447
1448What is the least number {{x}} that when divided by 3 leaves a remainder of 2, when
1449divided by 5 leaves a remainder of 3, and when divided by 7 leaves a remainder
1450of 2?
1451<enscript highlight="scheme">
1452> (solve-chinese '(2 3 2) '(3 5 7))
145323
1454</enscript>
1455
1456Wikipedia: [[https://en.wikipedia.org/wiki/Chinese_remainder_theorem|Chinese Remainder Theorem]]
1457
1458<procedure>(quadratic-residue? a n) -> boolean</procedure>
1459
1460; a : integer
1461; n : integer
1462
1463Returns {{#t}} if {{a}} is a quadratic residue modulo {{n}}, otherwise {{#f}}. The modulus
1464{{n}} must be positive, and a must be nonnegative.
1465
1466Formally, {{a}} is a quadratic residue modulo {{n}} if there exists a number
1467{{x}} such that {{(* x x)}} = {{a}} (mod {{n}}). In other words,
1468{{(quadratic-residue? a n)}} is {{#t}} when {{a}} is a perfect square modulo
1469{{n}}.
1470
1471Examples:
1472<enscript highlight="scheme">
1473> (quadratic-residue? 0 4)
1474#f
1475> (quadratic-residue? 1 4)
1476#t
1477> (quadratic-residue? 2 4)
1478#f
1479> (quadratic-residue? 3 4)
1480#f
1481</enscript>
1482
1483Wikipedia: [[https://en.wikipedia.org/wiki/Quadratic_residue|Quadratic Residue]]
1484
1485<procedure>(quadratic-character a p) -> integer</procedure>
1486
1487; a : integer
1488; b : integer
1489
1490Returns the value of the quadratic character modulo the prime {{p}}. That is,
1491for a non-zero {{a}} the number 1 is returned when {{a}} is a quadratic
1492residue, and -1 is returned when {{a}} is a non-residue. If {{a}} is zero, then
14930 is returned.
1494
1495If {{a}} is negative or {{p}} is not positive, quadratic-character raises an
1496error. If {{p}} is not prime, (quadratic-character a p) is indeterminate.
1497
1498This function is also known as the ''Legendre symbol''.
1499
1500<enscript highlight="scheme">
1501> (quadratic-character 0 5)
15020
1503> (quadratic-character 1 5)
15041
1505> (quadratic-character 2 5)
1506-1
1507> (quadratic-character 3 5)
1508-1
1509</enscript>
1510
1511Wikipedia: [[https://en.wikipedia.org/wiki/Legendre_symbol|Legendre Symbol]]
1512
1513<procedure>(jacobi-symbol a n) -> integer</procedure>
1514
1515; a : integer
1516; n : integer
1517
1518Computes the Jacobi symbol for any nonnegative integer {{a}} and any positive
1519odd integer {{n}}.
1520
1521If {{n}} is not an odd positive integer, {{(jacobi-symbol a n)}} throws an exception.
1522
1523<enscript highlight="scheme">
1524> (jacobi-symbol 1 1)
15251
1526> (jacobi-symbol 8 11)
1527-1
1528> (jacobi-symbol 39 27)
15290
1530> (jacobi-symbol 22 59)
15311
1532> (jacobi-symbol 32 8)
1533Error: (jacobi-symbol) bad argument type - not an odd integer: 8
1534</enscript>
1535
1536Wikipedia: [[https://en.wikipedia.org/wiki/Jacobi_symbol|Jacobi Symbol]]
1537
1538<procedure>(modular-inverse a n) -> integer</procedure>
1539
1540; a : integer
1541; b : integer
1542
1543Returns the inverse of a modulo {{n}} if {{a}} and {{n}} are coprime, otherwise raises an
1544error. The modulus {{n}} must be positive, and {{a}} must be nonzero.
1545
1546Formally, if {{a}} and {{n}} are coprime, {{b}} = {{(modular-inverse a n)}} is the unique
1547natural number less than {{n}} such that {{(* a b)}} = {{1}} (mod {{n}}).
1548
1549<enscript highlight="scheme">
1550> (modular-inverse 2 5)
15513
1552> (modulo (* 2 3) 5)
15531
1554</enscript>
1555
1556Wikipedia: [[https://en.wikipedia.org/wiki/Modular_multiplicative_inverse|Multiplicative Inverse]]
1557
1558<procedure>(modular-expt a b n) -> integer</procedure>
1559
1560; a : integer
1561; b : integer
1562; n : integer
1563
1564Computes {{(modulo (expt a b) n)}}, but much more efficiently. The modulus {{n}} must
1565be positive, and the exponent {{b}} must be nonnegative.
1566
1567<enscript highlight="scheme">
1568> (modulo (expt -6 523) 19)
156913
1570> (modular-expt -6 523 19)
157113
1572> (modular-expt 9 158235208 19)
15734
1574> ; don't try this at home!
1575  (modulo (expt 9 158235208) 19)
15764
1577</enscript>
1578
1579===== Parameterized Modular Arithmetic
1580
1581The {{math.number-theory}} library supports modular arithmetic parameterized on
1582a current modulus. For example, the code
1583
1584<enscript highlight="scheme">
1585(with-modulus n
1586  (mod= (modexpt a b) c))
1587</enscript>
1588
1589corresponds with the mathematical statement ''a^b'' = ''c'' (mod ''n'').
1590
1591The current modulus is stored in a parameter that, for performance reasons, can
1592only be set using with-modulus. (The basic modular operators cache parameter
1593reads, and this restriction guarantees that the cached values are current.
1594'''NOTE:''' I'm not entirely sure this is true for the CHICKEN port, as a
1595slightly more complicated racket syntax-case has been turned into a simple
1596syntax-rule for {{(parameterize ...)}})
1597
1598Wikipedia: [[https://en.wikipedia.org/wiki/Modular_arithmetic|Modular Arithmetic]]
1599
1600<syntax>(with-modulus n body ...)</syntax>
1601
1602; n : integer
1603
1604Alters the current modulus within the dynamic extent of {{body}}. The
1605expression {{n}} must evaluate to a positive integer.
1606
1607By default, the current modulus is 1, meaning that every modular arithmetic
1608expression that does not raise an error returns 0.
1609
1610<procedure>(current-modulus) -> integer</procedure>
1611
1612Returns the current modulus.
1613
1614Examples:
1615
1616<enscript highlight="scheme">
1617> (current-modulus)
16181
1619> (with-modulus 5 (current-modulus))
16205
1621</enscript>
1622
1623<procedure>(mod x) -> integer</procedure>
1624
1625; x : exact rational
1626
1627Converts a rational number {{x}} to a natural number less than the current
1628modulus.
1629
1630If {{x}} is an integer, this is equivalent to {{(modulo x n)}}. If {{x}} is a
1631fraction, an integer input is generated by multiplying its numerator by its
1632denominator’s modular inverse.
1633
1634Examples:
1635
1636<enscript highlight="scheme">
1637> (with-modulus 7 (mod (* 218 7)))
16380
1639> (with-modulus 7 (mod 3/2))
16405
1641> (with-modulus 7 (mod/ 3 2))
16425
1643> (with-modulus 7 (mod 3/7))
1644Error: (modular-inverse) bad argument type - not coprime to modulus 7: 7
1645</enscript>
1646
1647<procedure>(mod+ a ...) -> integer</procedure>
1648<procedure>(mod* a ...) -> integer</procedure>
1649
1650; a : integer
1651
1652Equivalent to {{(modulo (+ a ...) (current-modulus))}} and {{(modulo (* a ...) (current-modulus))}},
1653respectively, but generate smaller intermediate values.
1654
1655<procedure>(modsqr a) -> integer</procedure>
1656<procedure>(modexpt a b) -> integer</procedure>
1657
1658; a : integer
1659; b : integer
1660
1661Equivalent to {{(mod* a a)}} and {{(modular-expt a b (current-modulus))}}, respectively.
1662
1663<procedure>(mod- a b ...) -> integer</procedure>
1664
1665; a : integer
1666; b : integer
1667
1668Equivalent to {{(modulo (- a b ...) (current-modulus))}}, but generates smaller
1669intermediate values. Note that {{(mod- a)}} = {{(mod (- a))}}.
1670
1671<procedure>(mod/ a b ...) -> integer</procedure>
1672
1673; a : integer
1674; b : integer
1675
1676Divides a by {{(* b ...)}}, by multiplying {{a}} by the multiplicative inverse of
1677{{(* b ...)}}. The one-argument variant returns the modular inverse of {{a}}.
1678
1679Note that {{(mod/ a b ...)}} is '''not''' equivalent to
1680{{(modulo (/ a b ...) (current-modulus))}}; see {{mod=}} for a demonstration.
1681
1682<procedure>(mod= a b ...) -> boolean</procedure>
1683<procedure>(mod< a b ...) -> boolean</procedure>
1684<procedure>(mod<= a b ...) -> boolean</procedure>
1685<procedure>(mod> a b ...) -> boolean</procedure>
1686<procedure>(mod>= a b ...) -> boolean</procedure>
1687
1688; a : integer
1689; b : integer
1690
1691Each of these is equivalent to {{(op (mod a) (mod b) ...)}}, where op is the
1692corresponding numeric comparison function. Additionally, when given one
1693argument, the inequality tests always return {{#t}}.
1694
1695Suppose we wanted to know why 17/4 = 8 (mod 15), but 51/12 (mod 15) is
1696undefined, even though normally 51/12 = 17/4. In code,
1697
1698<enscript highlight="scheme">
1699> (with-modulus 15 (mod/ 17 4))
17008
1701> (/ 51 12)
170217/4
1703> (with-modulus 15 (mod/ 51 12))
1704Error: (modular-inverse) bad argument type - not coprime to modulus 15: 12
1705</enscript>
1706
1707We could try to divide by brute force: find, modulo 15, all the numbers {{a}} for
1708which {{(mod* a 4)}} is 17, then find all the numbers {{b}} for which {{(mod* a 12)}} is
170951.
1710<enscript highlight="scheme">
1711(import srfi-42)
1712> (with-modulus 15
1713    (list-ec (:range a 15)
1714             (if (mod= (mod* a 4) 17))
1715      a))
1716(8)
1717> (with-modulus 15
1718    (list-ec (:range a 15)
1719             (if (mod= (mod* a 12) 51))
1720      a))
1721(3 8 13)
1722</enscript>
1723
1724So the problem isn't that {{b}} doesn't exist, it's that {{b}} isn't
1725''unique''.
1726
1727===== Primes
1728<procedure>(prime? z) -> boolean</procedure>
1729
1730; z : integer
1731
1732Returns {{#t}} if {{z}} is a prime, {{#f}} otherwise.
1733
1734Formally, an integer {{z}} is prime when the only positive divisors of {{z}}
1735are 1 and {{(abs z)}}.
1736
1737The positive primes below 20 are:
1738<enscript highlight="scheme">
1739> (filter prime? (iota 20 1))
1740(2 3 5 7 11 13 17 19)
1741</enscript>
1742
1743The corresponding negative primes are:
1744<enscript highlight="scheme">
1745> (filter prime? (iota 20 0 -1))
1746(-2 -3 -5 -7 -11 -13 -17 -19)
1747</enscript>
1748
1749Wikipedia: [[https://en.wikipedia.org/wiki/Prime_number|Prime Number]]
1750
1751<procedure>(odd-prime? z) -> boolean</procedure>
1752
1753; z : integer
1754
1755Returns {{#t}} if {{z}} is a odd prime, {{#f}} otherwise.
1756
1757<enscript highlight="scheme">
1758> (odd-prime? 2)
1759#f
1760> (odd-prime? 3)
1761#t
1762</enscript>
1763
1764<procedure>(nth-prime n) -> integer</procedure>
1765
1766; n : integer
1767
1768Returns the {{n}}th positive prime; {{n}} must be nonnegative.
1769
1770<enscript highlight="scheme">
1771> (nth-prime 0)
17722
1773> (nth-prime 1)
17743
1775> (nth-prime 2)
17765
1777</enscript>
1778
1779<procedure>(random-prime n) -> integer</procedure>
1780
1781; n : integer
1782
1783Returns a random prime smaller than {{n}}, which must be greater than 2.
1784
1785The function {{random-prime}} picks random numbers below {{n}} until a prime is
1786found.
1787
1788<enscript highlight="scheme">
1789> (random-prime 10)
17903
1791> (random-prime 10)
17922
1793> (random-prime 10)
17945
1795</enscript>
1796
1797<procedure>(next-prime z) -> integer</procedure>
1798
1799; z : integer
1800
1801Returns the first prime larger than {{z}}.
1802
1803<enscript highlight="scheme">
1804> (next-prime 4)
18055
1806> (next-prime 5)
18077
1808</enscript>
1809
1810<procedure>(prev-prime z) -> integer</procedure>
1811
1812Returns the first prime smaller than {{z}}.
1813
1814<enscript highlight="scheme">
1815> (prev-prime 4)
18163
1817> (prev-prime 5)
18183
1819</enscript>
1820
1821<procedure>(next-primes z n) -> (list-of integer)</procedure>
1822
1823; z : integer
1824; n : integer
1825
1826Returns list of the next {{n}} primes larger than {{z}}; {{n}} must be
1827nonnegative.
1828
1829<enscript highlight="scheme">
1830> (next-primes 2 4)
1831(3 5 7 11)
1832</enscript>
1833
1834<procedure>(prev-primes z n) -> (list-of integer)</procedure>
1835
1836; z : integer
1837; n : integer
1838
1839Returns list of the next {{n}} primes smaller than {{z}}; {{n}} must be
1840nonnegative.
1841
1842<enscript highlight="scheme">
1843> (prev-primes 13 4)
1844(11 7 5 3)
1845</enscript>
1846
1847<procedure>(factorize n) -> (list-of (list integer integer))</procedure>
1848
1849; n : integer
1850
1851Returns the factorization of a natural number {{n}}. The factorization consists of
1852a list of corresponding primes and exponents. The primes will be in ascending
1853order.
1854
1855The prime factorization of 600 = 2^3 * 3^1 * 5^2:
1856
1857<enscript highlight="scheme">
1858> (factorize 600)
1859((2 3) (3 1) (5 2))
1860</enscript>
1861
1862<procedure>(defactorize f) -> integer</procedure>
1863
1864; f : (list-of (list integer integer))
1865
1866Returns the natural number, whose factorization is given by {{f}}. The
1867factorization {{f}} is represented as described in {{factorize}}.
1868
1869<enscript highlight="scheme">
1870> (defactorize '((2 3) (3 1) (5 2)))
1871600
1872</enscript>
1873
1874Wikipedia: [[https://en.wikipedia.org/wiki/Integer_factorization|Integer Factorization]]
1875
1876<procedure>(divisors z) -> (list-of integer)</procedure>
1877
1878; z : integer
1879
1880Returns a list of all positive divisors of the integer {{z}}. The divisors appear
1881in ascending order.
1882
1883<enscript highlight="scheme">
1884> (divisors 120)
1885(1 2 3 4 5 6 8 10 12 15 20 24 30 40 60 120)
1886> (divisors -120)
1887(1 2 3 4 5 6 8 10 12 15 20 24 30 40 60 120)
1888</enscript>
1889
1890<procedure>(prime-divisors z) -> (list-of integer)</procedure>
1891
1892; z : integer
1893
1894Returns a list of all positive prime divisors of the integer {{z}}. The divisors
1895appear in ascending order.
1896
1897<enscript highlight="scheme">
1898> (prime-divisors 120)
1899'(2 3 5)
1900</enscript>
1901
1902<procedure>(prime-exponents z) -> (list-of integer)</procedure>
1903
1904; z : integer
1905
1906Returns a list of the exponents of in a factorization of the integer {{z}}.
1907
1908<enscript highlight="scheme">
1909> (define z (* 2 2 2 3 5 5))
1910> (prime-divisors z)
1911(2 3 5)
1912> (prime-exponents z)
1913(3 1 2)
1914</enscript>
1915
1916===== Roots
1917
1918<procedure>(integer-root n m) -> integer</procedure>
1919
1920; n : integer
1921; m : integer
1922
1923Returns the mth integer root of {{n}}. This is the largest integer {{r}} such that
1924{{(expt r m)}} <= {{n}}.
1925
1926<enscript highlight="scheme">
1927> (integer-root (expt 3 4) 4)
19283
1929> (integer-root (+ (expt 3 4) 1) 4)
19303
1931</enscript>
1932
1933<procedure>(integer-root/remainder n m) -> integer integer</procedure>
1934
1935; n : integer
1936; m : integer
1937
1938Returns two values. The first, {{r}}, is the {{m}}th integer root of {{n}}. The
1939second is {{n}}-{{r}}^{{m}}.
1940
1941<enscript highlight="scheme">
1942> (integer-root/remainder (expt 3 4) 4)
19433
19440
1945> (integer-root/remainder (+ (expt 3 4) 1) 4)
19463
19471
1948</enscript>
1949
1950===== Powers
1951
1952<procedure>(max-dividing-power a b) -> integer</procedure>
1953
1954; a : integer
1955; b : integre
1956
1957Returns the largest exponent, {{n}}, of a power with base {{a}} that divides
1958{{b}}.
1959
1960That is, {{(expt a n)}} divides {{b}} but {{(expt a (+ n 1))}} does not divide {{b}}.
1961
1962<enscript highlight="scheme">
1963> (max-dividing-power 3 (expt 3 4))
19644
1965> (max-dividing-power 3 5)
19660
1967</enscript>
1968
1969<procedure>(perfect-power m) -> (or (list integer integer) #f)</procedure>
1970
1971; m : integer
1972
1973If {{m}} is a perfect power, a list with two elements {{b}} and {{n}} such that
1974{{(expt b n)}} = {{m}} is returned, otherwise {{#f}} is returned.
1975
1976<enscript highlight="scheme">
1977> (perfect-power (expt 3 4))
1978(3 4)
1979> (perfect-power (+ (expt 3 4) 1))
1980#f
1981</enscript>
1982
1983<procedure>(perfect-power? m) -> boolean</procedure>
1984
1985; m : integer
1986
1987Returns {{#t}} if {{m}} is a perfect power, otherwise {{#f}}.
1988
1989<enscript highlight="scheme">
1990> (perfect-power? (expt 3 4))
1991#t
1992> (perfect-power? (+ (expt 3 4) 1))
1993#f
1994</enscript>
1995
1996Wikipedia: [[https://en.wikipedia.org/wiki/Perfect_power|Perfect Power]]
1997
1998<procedure>(prime-power m) -> (or (list integer integer) #f)</procedure>
1999
2000; m : integer
2001
2002If {{M}} is a power of the form {{(expt p n)}} where p is prime, then a list with the
2003prime and the exponent is returned, otherwise {{#f}} is returned.
2004
2005<enscript highlight="scheme">
2006> (prime-power (expt 3 4))
2007(3 4)
2008> (prime-power (expt 6 4))
2009#f
2010</enscript>
2011
2012<procedure>(prime-power? m) -> boolean</procedure>
2013
2014; m : integer
2015
2016Returns {{#t}} if {{m}} is a prime power, otherwise {{#f}}.
2017
2018<enscript highlight="scheme">
2019> (prime-power? (expt 3 4))
2020#t
2021> (prime-power? (expt 6 4))
2022#f
2023> (prime-power? 1)
2024#f
2025> (prime-power? 0)
2026#f
2027</enscript>
2028
2029<procedure>(odd-prime-power? m) -> boolean</procedure>
2030
2031; m : integer
2032
2033Returns {{#t}} if {{m}} is a power of an odd prime, otherwise {{#f}}.
2034
2035<enscript highlight="scheme">
2036> (odd-prime-power? (expt 2 4))
2037#f
2038> (odd-prime-power? (expt 3 4))
2039#t
2040> (odd-prime-power? (expt 15 4))
2041#f
2042</enscript>
2043
2044<procedure>(as-power m) -> integer integer</procedure>
2045
2046; m : integer
2047
2048Returns two values {{b}} and {{n}} such that {{m}} = {{(expt b n)}} and {{n}}
2049is maximal.
2050
2051<enscript highlight="scheme">
2052> (as-power (* (expt 2 4) (expt 3 4)))
20536
20544
2055> (expt 6 4)
20561296
2057> (* (expt 2 4) (expt 3 4))
20581296
2059> (as-power (* (expt 2 4) (expt 3 5)))
20603888
20611
2062</enscript>
2063
2064<procedure>(prefect-square m) -> (or integer #f)</procedure>
2065
2066Returns {{(sqrt m)}} if {{m}} is perfect square, otherwise {{#f}}.
2067
2068<enscript highlight="scheme">
2069> (perfect-square 9)
20703
2071> (perfect-square 10)
2072#f
2073</enscript>
2074
2075===== Multiplicative and Arithmetic Functions
2076
2077The functions in this section are ''multiplicative'' (with exception of the Von
2078Mangoldt function). In number theory, a multiplicative function is a function
2079{{f}} such that {{(f (* a b))}} = {{(* (f a) (f b))}} for all coprime natural
2080numbers {{a}} and {{b}}.
2081
2082<procedure>(totient n) -> integer</procedure>
2083
2084; n : integer
2085
2086Returns the number of integers from 1 to {{n}} that are coprime with {{n}}.
2087
2088This function is known as Euler's totient or phi function.
2089
2090<enscript highlight="scheme">
2091> (totient 9)
20926
2093> (length (filter (curry coprime? 9) (range 10)))
20946
2095</enscript>
2096
2097Wikipedia: [[https://en.wikipedia.org/wiki/Euler%27s_totient_function|Euler's Totient]]
2098
2099<procedure>(moebius-mu n) -> integer</procedure>
2100
2101; n : integer
2102
2103Returns:
2104
2105* 1 if {{n}} is a square-free product of an even number of primes
2106* -1 if {{n}} is a square-free product of an odd number of primes
2107* 0 if {{n}} has a multiple prime factor
2108
2109<enscript highlight="scheme">
2110> (moebius-mu (* 2 3 5))
2111-1
2112> (moebius-mu (* 2 3 5 7))
21131
2114> (moebius-mu (* 2 2 3 5 7))
21150
2116</enscript>
2117
2118Wikipedia: [[https://en.wikipedia.org/wiki/M%C3%B6bius_function|Moebius Function]]
2119
2120<procedure>(divisor-sum n k) -> integer</procedure>
2121
2122; n : integer
2123; k : integer
2124
2125Returns sum of the {{k}}th powers of all divisors of {{n}}.
2126
2127<enscript highlight="scheme">
2128> (divisor-sum 12 2)
2129210
2130> (apply + (map sqr (divisors 12)))
2131210
2132</enscript>
2133
2134Wikipedia: [[https://en.wikipedia.org/wiki/Divisor_function|Divisor Function]]
2135
2136<procedure>(prime-omega n) -> integer</procedure>
2137
2138; n : integer
2139
2140Counting multiplicities the number of prime factors of {{n}} is returned.
2141
2142<enscript highlight="scheme">
2143> (prime-omega (* 2 2 2 3 3 5))
21446
2145</enscript>
2146
2147OEIS: [[https://oeis.org/A001222|Big Omega]], [[https://oeis.org/wiki/Omega(n),_number_of_prime_factors_of_n_(with_multiplicity)|Big Omega]]
2148
2149<procedure>(mangoldt-lambda n) -> number</procedure>
2150
2151; n : integer
2152
2153The von Mangoldt function. If n=p^k for a prime {{p}} and an integer {{k>=1}}
2154then {{(log n)}} is returned. Otherwise {{0}} is returned.
2155
2156Note: The Von Mangoldt function is not multiplicative.
2157
2158<enscript highlight="scheme">
2159> (mangoldt-lambda (* 3 3))
21601.0986122886681098
2161> (log 3)
21621.0986122886681098
2163</enscript>
2164
2165Wikipedia: [[https://en.wikipedia.org/wiki/Von_Mangoldt_function|Von Mangoldt Function]]
2166
2167===== Number Sequences
2168
2169<procedure>(bernoulli-number n) -> ratnum</procedure>
2170
2171; n : integer
2172
2173Returns the {{n}}th Bernoulli number; {{n}} must be nonnegative.
2174
2175<enscript highlight="scheme">
2176> (map bernoulli-number (iota 9))
2177(1 -1/2 1/6 0 -1/30 0 1/42 0 -1/30)
2178</enscript>
2179
2180Note that these are the ''first'' Bernoulli numbers, since
2181{{(bernoulli-number 1)}} = {{-1/2}}.
2182
2183Wikipedia: [[https://en.wikipedia.org/wiki/Bernoulli_number|Bernoulli Number]]
2184
2185<procedure>(eulerian-number n k) -> integer</procedure>
2186
2187; n : integer
2188; k : integer
2189
2190Returns the Eulerian number {{<n,k>}}; both arguments must be nonnegative.
2191
2192<enscript highlight="scheme">
2193> (eulerian-number 5 2)
219466
2195</enscript>
2196
2197Wikipedia: [[http://mathworld.wolfram.com/EulerianNumber.html|Eulerian Number]]
2198
2199<procedure>(fibonacci n) -> integer</procedure>
2200
2201; n : integer
2202
2203Returns the {{n}}th Fibonacci number; {{n}} must be nonnegative.
2204
2205The ten first Fibonacci numbers:
2206<enscript highlight="scheme">
2207> (map fibonacci (iota 10))
2208'(0 1 1 2 3 5 8 13 21 34)
2209</enscript>
2210
2211Wikipedia: [[https://en.wikipedia.org/wiki/Fibonacci_number|Fibonacci Number]]
2212
2213<procedure>(make-fibonaci a b) -> (integer -> integer)</procedure>
2214
2215; a : integer
2216; b : integer
2217
2218Returns a function representing a Fibonacci sequence with the first two numbers
2219{{a}} and {{b}}. The {{fibonacci}} function is defined as {{(make-fibonacci 0 1)}}.
2220
2221The Lucas numbers are defined as a Fibonacci sequence starting with 2 and 1:
2222
2223<enscript highlight="scheme">
2224> (map (make-fibonacci 2 1) (iota 10))
2225(2 1 3 4 7 11 18 29 47 76)
2226</enscript>
2227
2228Wikipedia: [[https://wikipedia.org/wiki/Lucas_number|Lucas Number]]
2229
2230<procedure>(modular-fibonacci n m) -> integer</procedure>
2231
2232; n : integer
2233; m : integer
2234
2235Returns the {{n}}th Fibonacci number modulo {{m}}; {{n}} must be nonnegative
2236and {{m}} must be positive.
2237
2238The ten first Fibonacci numbers modulo 5:
2239<enscript highlight="scheme">
2240> (map (lambda (n) (modular-fibonacci n 5)) (range 10))
2241(0 1 1 2 3 0 3 3 1 4)
2242</enscript>
2243
2244<procedure>(make-modular-fibonacci a b) -> (integer integer -> integer)</procedure>
2245
2246; a : integer
2247; b : integer
2248
2249Like {{make-fibonacci}}, but makes a modular fibonacci sequence.
2250
2251<procedure>(farey-sequence n) -> (list-of ratnum)</procedure>
2252
2253; n : integer
2254
2255Returns a list of the numbers in the {{n}}th Farey sequence; {{n}} must be positive.
2256
2257The {{n}}th Farey sequence is the sequence of all completely reduced rational
2258numbers from 0 to 1 which denominators are less than or equal to {{n}}.
2259
2260<enscript highlight="scheme">
2261> (farey-sequence 1)
2262(0 1)
2263> (farey-sequence 2)
2264(0 1/2 1)
2265> (farey-sequence 3)
2266(0 1/3 1/2 2/3 1)
2267</enscript>
2268
2269Wikipedia: [[https://en.wikipedia.org/wiki/Farey_sequence|Farey Sequence]]
2270
2271<procedure>(tangent-number n) -> integer</procedure>
2272
2273; n : integer
2274
2275Returns the {{n}}th tangent number; {{n}} must be nonnegative.
2276
2277<enscript highlight="scheme">
2278> (tangent-number 1)
22791
2280> (tangent-number 2)
22810
2282> (tangent-number 3)
22832
2284</enscript>
2285
2286MathWorld: [[http://mathworld.wolfram.com/TangentNumber.html|Tangent Number]]
2287
2288===== Combinatorics
2289
2290<procedure>(factorial n) -> integer</procedure>
2291
2292; n : integer
2293
2294Returns the factorial of {{n}}, which must be nonnegative. The factorial of
2295{{n}} is the number {{(* n (- n 1) (- n 2) ... 1)}}.
2296
2297<enscript highlight="scheme">
2298> (factorial 3)
22996
2300> (factorial 0)
23011
2302</enscript>
2303
2304Wikipedia: [[https://en.wikipedia.org/wiki/Factorial|Factorial]]
2305
2306<procedure>(binomial n k) -> integer</procedure>
2307
2308; n : integer
2309; k : integer
2310
2311Returns the number of ways to choose a set of k items from a set of n items;
2312i.e. the order of the k items is not significant. Both arguments must be
2313nonnegative.
2314
2315When {{k > n}}, {{(binomial n k) = 0}}. Otherwise, {{(binomial n k)}} is
2316equivalent to {{(/ (factorial n) (factorial k) (factorial (- n k)))}},
2317but computed more quickly.
2318
2319<enscript highlight="scheme">
2320> (binomial 5 3)
232110
2322</enscript>
2323
2324Wikipedia: [[https://en.wikipedia.org/wiki/Binomial_coefficient|Binomial Coefficient]]
2325
2326<procedure>(permutations n k) -> integer</procedure>
2327
2328; n : integer
2329; k : integer
2330
2331Returns the number of ways to choose a sequence of {{k}} items from a set of n
2332items; i.e. the order of the {{k}} items is significant. Both arguments must be
2333nonnegative.
2334
2335When {{k > n}}, {{(permutations n k) = 0}}. Otherwise, {{(permutations n k)}}
2336is equivalent to {{(/ (factorial n) (factorial (- n k)))}}.
2337
2338<enscript highlight="scheme">
2339> (permutations 5 3)
234060
2341</enscript>
2342
2343Wikipedia: [[https://en.wikipedia.org/wiki/Permutation#Permutations_in_combinatorics|Permutations]]
2344
2345<procedure>(multinomial n ks) -> integer</procedure>
2346
2347; n : integer
2348; ks : (list-of integer)
2349
2350A generalization of binomial to multiple sets of choices; e.g.
2351{{(multinomial n (list k0 k1 k2))}} is the number of ways to choose a set of
2352{{k0}} items, a set of {{k1}} items, and a set of {{k2}} items from a set of {{n}}
2353items. All arguments must be nonnegative.
2354
2355When {{(apply + ks) = n}}, this is equivalent to
2356{{(apply / (factorial n) (map factorial ks))}}. Otherwise, multinomial returns 0.
2357
2358<enscript highlight="scheme">
2359> (multinomial 5 '(3 2))
236010
2361> (= (multinomial 8 '(5 3))
2362     (binomial 8 5)
2363     (binomial 8 3))
2364#t
2365> (multinomial 10 '(5 3 2))
23662520
2367> (multinomial 0 '())
23681
2369> (multinomial 4 '(1 1))
23700
2371</enscript>
2372
2373Wikipedia: [[https://en.wikipedia.org/wiki/Multinomial_theorem#Multinomial_coefficients|Multinomial Coefficient]]
2374
2375<procedure>(partitions n) -> integer</procedure>
2376
2377; n : integer
2378
2379Returns the number of partitions of {{n}}, which must be nonnegative. A partition
2380of a positive integer {{n}} is a way of writing {{n}} as a sum of positive integers.
2381The number 3 has the partitions {{(+ 1 1 1)}}, {{(+ 1 2)}} and {{(+ 3)}}.
2382
2383<enscript highlight="scheme">
2384> (partitions 3)
23853
2386> (partitions 4)
23875
2388</enscript>
2389
2390Wikipedia: [[https://en.wikipedia.org/wiki/Partition_(number_theory)|Partition]]
2391
2392===== Polygonal numbers
2393
2394<procedure>(triangle-number? n) -> boolean</procedure>
2395<procedure>(square-number? n) -> boolean</procedure>
2396<procedure>(pentagonal-number? n) -> boolean</procedure>
2397<procedure>(hexagonal-number? n) -> boolean</procedure>
2398<procedure>(heptagonal-number? n) -> boolean</procedure>
2399<procedure>(octagonal-number? n) -> boolean</procedure>
2400
2401; n : integer
2402
2403These functions check whether the input is a polygonal number of the types
2404triangle, square, pentagonal, hexagonal, heptagonal and octogonal respectively.
2405
2406Wikipedia: [[https://en.wikipedia.org/wiki/Polygonal_number|Polygonal Number]]
2407
2408
2409<procedure>(triangle-number n) -> boolean</procedure>
2410<procedure>(sqr n) -> boolean</procedure>
2411<procedure>(pentagonal-number n) -> boolean</procedure>
2412<procedure>(hexagonal-number n) -> boolean</procedure>
2413<procedure>(heptagonal-number n) -> boolean</procedure>
2414<procedure>(octagonal-number n) -> boolean</procedure>
2415
2416; n : integer
2417
2418These functions return the {{n}}th polygonal number of the corresponding type
2419of polygonal number.
2420
2421===== Fractions
2422
2423<procedure>(mediant x y) -> ratnum</procedure>
2424
2425; x : ratnum
2426; y : ratnum
2427
2428Computes the mediant of the numbers {{x}} and {{y}}. The mediant of two
2429fractions {{p/q}} and {{r/s}} in their lowest term is the number
2430{{(p+r)/(q+s)}}.
2431
2432<enscript highlight="scheme">
2433> (mediant 1/2 5/6)
24343/4
2435</enscript>
2436
2437Wikipedia: [[https://en.wikipedia.org/wiki/Mediant_(mathematics)|Mediant]]
2438
2439===== The Quadratic Equation
2440
2441<procedure>(quadratic-solutions a b c) -> (list-of number)</procedure>
2442  a : number
2443  b : number
2444  c : number
2445
2446Returns a list of all real solutions to the equation a {{x^2 + bx + c = 0}}
2447with real coefficients.
2448
2449<enscript highlight="scheme">
2450> (quadratic-solutions 1 0 -1)
2451'(-1 1)
2452> (quadratic-solutions 1 2 1)
2453'(-1)
2454> (quadratic-solutions 1 0 1)
2455'()
2456</enscript>
2457
2458<procedure>(quadratic-integer-solutions a b c) -> (list-of integer)</procedure>
2459
2460; a : number
2461; b : number
2462; c : number
2463
2464Returns a list of all integer solutions to the equation a {{x^2 + bx + c = 0}}
2465with real coefficients.
2466
2467<enscript highlight="scheme">
2468> (quadratic-integer-solutions 1 0 -1)
2469'(-1 1)
2470> (quadratic-integer-solutions 1 0 -2)
2471'()
2472</enscript>
2473
2474<procedure>(quadratic-natural-solutions a b c) -> (list-of integer)</procedure>
2475
2476; a : number
2477; b : number
2478; c : number
2479
2480Returns a list of all natural solutions to the equation a {{x^2 + bx + c = 0}}
2481with real coefficients.
2482
2483<enscript highlight="scheme">
2484> (quadratic-natural-solutions 1 0 -1)
2485'(1)
2486> (quadratic-natural-solutions 1 0 -2)
2487'()
2488</enscript>
2489
2490<procedure>(complex-quadratic-solutions a b c) -> (list-of number)</procedure>
2491
2492; a : number
2493; b : number
2494; c : number
2495
2496Returns a list of all complex solutions to the equation a {{x^2 + bx + c = 0}}.
2497This function allows complex coeffecients.
2498
2499<enscript highlight="scheme">
2500> (complex-quadratic-solutions 1 0 1)
2501(0-1i 0+1i)
2502> (complex-quadratic-solutions 1 0 (sqrt -1))
2503(-0.7071067811865476+0.7071067811865475i 0.7071067811865476-0.7071067811865475i)
2504> (complex-quadratic-solutions 1 0 1)
2505(0-1i 0+1i)
2506</enscript>
2507
2508===== The group Zn and Primitive Roots
2509
2510The numbers 0, 1, ..., n-1 with addition and multiplication modulo {{n}} is a ring
2511called {{Zn}}.
2512
2513The group of units in {{Zn}} with respect to multiplication modulo {{n}} is
2514called {{Un}}.
2515
2516The order of an element {{x}} in {{Un}} is the least {{k>0}} such that {{x^k=1 mod n}}.
2517
2518A generator the group {{Un}} is called a primitive root modulo {{n}}. Note that {{g}} is a
2519primitive root if and only if {{order(g)=totient(n)}}. A group with a generator is
2520called cyclic.
2521
2522Wikipedia: [[https://en.wikipedia.org/wiki/Multiplicative_group_of_integers_modulo_n|The Group Zn]]
2523
2524<procedure>(unit-group n) -> (list-of integer)</procedure>
2525
2526; n : integer
2527
2528Returns a list of all elements of {{Un}}, the unit group modulo {{n}}. The
2529modulus {{n}} must be positive.
2530
2531<enscript highlight="scheme">
2532> (unit-group 5)
2533(1 2 3 4)
2534> (unit-group 6)
2535(1 5)
2536</enscript>
2537
2538<procedure>(unit-group-order x n) -> integer</procedure>
2539
2540; x : integer
2541; n : integer
2542
2543Returns the order of {{x}} in the group {{Un}}; both arguments must be
2544positive. If {{x}} and n are not coprime, {{(unit-group-order x n)}} raises an
2545error.
2546
2547<enscript highlight="scheme">
2548> (unit-group-order 2 5)
25494
2550> (unit-group-order 2 6)
2551Error: (unit-group-order) expected coprime arguments; given 2 and 6
2552</enscript>
2553
2554<procedure>(unit-group-orders n) -> (list-of positive-integer)</procedure>
2555
2556; n : integer
2557
2558Returns a list {{(list (unit-group-order x0 n) (unit-group-order x1 n) ...)}}
2559where {{x0}}, {{x1}}, {{...}} are the elements of {{Un}}. The modulus {{n}}
2560must be positive.
2561
2562<enscript highlight="scheme">
2563> (unit-group-orders 5)
2564(1 4 4 2)
2565> (map (cut unit-group-order <> 5) (unit-group 5))
2566(1 4 4 2)
2567</enscript>
2568
2569<procedure>(primitive-root? x n) -> boolean</procedure>
2570
2571; x : integer
2572; n : integer
2573
2574Returns {{#t}} if the element {{x}} in {{Un}} is a primitive root modulo {{n}},
2575otherwise {{#f}} is returned. An error is signaled if {{x}} is not a member of
2576{{Un}}. Both arguments must be positive.
2577
2578<enscript highlight="scheme">
2579> (primitive-root? 1 5)
2580#f
2581> (primitive-root? 2 5)
2582#t
2583> (primitive-root? 5 5)
2584Error: (primitive-root?) expected coprime arguments; given 5 and 5
2585</enscript>
2586
2587<procedure>(exists-primitive-root? n) -> boolean</procedure>
2588
2589; n : integer
2590
2591Returns {{#t}} if the group {{Un}} has a primitive root (i.e. it is cyclic),
2592otherwise {{#f}} is returned. In other words, {{#t}} is returned if {{n}} is
2593one of {{1}}, {{2}}, {{4}}, {{p^e}}, {{2*p^e}} where {{p}} is an odd prime, and
2594{{#f}} otherwise. The modulus {{n}} must be positive.
2595
2596<enscript highlight="scheme">
2597> (exists-primitive-root? 5)
2598#t
2599> (exists-primitive-root? 6)
2600#t
2601> (exists-primitive-root? 12)
2602#f
2603</enscript>
2604
2605<procedure>(primitive-root n) -> (or integer false)</procedure>
2606
2607Returns a primitive root of {{Un}} if one exists, otherwise {{#f}} is returned.
2608The modulus {{n}} must be positive.
2609
2610<enscript highlight="scheme">
2611> (primitive-root 5)
26122
2613> (primitive-root 6)
26145
2615</enscript>
2616
2617<procedure>(primitive-roots n) -> (list-of integer)</procedure>
2618
2619; n : integer
2620
2621Returns a list of all primitive roots of {Un}. The modulus {{n}} must be positive.
2622
2623<enscript highlight="scheme">
2624> (primitive-roots 3)
2625(2)
2626> (primitive-roots 5)
2627(2 3)
2628> (primitive-roots 6)
2629(5)
2630</enscript>
2631
2632=== Original documentation
2633
2634[[https://docs.racket-lang.org/math/]]
2635
2636=== Author
2637
2638Neil Toronto and Jens Axel SÞgaard for Racket
2639
2640=== Maintainer
2641
2642[[/users/diego-mundo|Diego A. Mundo]]
2643
2644=== Repository
2645
2646[[https://github.com/dieggsy/chicken-math]]
2647
2648=== License
2649
2650GPL-3.0
2651
2652=== Version History
2653
2654; 0.3.1 : Ensure to export all api procedures, correct module dependencies for build
2655; 0.3.0 : Finish (math base) and (math flonum), bug and typo fixes, credit original authors
2656; 0.2.3 : Fix broken .egg file
2657; 0.2.2 : Re-organize internals, add (math base constants)
2658; 0.2.1 : Minor bug fixes
2659; 0.2.0 : Update (math number-theory quadratic) to reflect upstream
2660; 0.1.0 : Initial release
Note: See TracBrowser for help on using the repository browser.