suzanne.soy/racket
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Renamed functions

Browse Source 
* bernoulli -> bernoulli-number
 * farey -> farey-sequence
 * fibonacci/mod -> modular-fibonacci
 * order -> unit-group-order
 * orders -> unit-group-orders

Documented `make-fibonacci' and `make-modular-fibonacci'

Reworked text about loading external libraries in docs for `math/bigfloat'

Removed type aliases like Z, Q, Prime (I like them, but TR was printing them
in unexpected places like array return types)
This commit is contained in:
Neil Toronto 2012-11-27 22:23:05 -07:00
parent a5961dcf36
commit 96d1400654
12 changed files with 193 additions and 203 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
collects/math/private/bigfloat/bigfloat-hurwitz-zeta.rkt
View File

@ -41,7 +41,7 @@
[numer : Bigfloat s]
[denom : Bigfloat (bf/ fn n+q)]
[k : Nonnegative-Fixnum 0])
(define ck (bf (/ (bernoulli (* 2 (fx+ k 1))) (factorial (* 2 (fx+ k 1))))))
(define ck (bf (/ (bernoulli-number (* 2 (fx+ k 1))) (factorial (* 2 (fx+ k 1))))))
(define dy (bf* (bf* numer denom) ck))
(define new-y (bf+ y dy))
(cond [((bfabs dy) . bf<= . (bf* eps (bfabs new-y)))

4
collects/math/private/functions/zeta.rkt
View File

@ -184,7 +184,7 @@ P. Borwein. An Efficient Algorithm for the Riemann Zeta Function.
[(integer? s)
(cond [(zero? s) 1/2]
[(negative? s) (define k (- 1 s))
(* (/ (bernoulli k) k) (- (expt 2 k) 1))]
(* (/ (bernoulli-number k) k) (- (expt 2 k) 1))]
[else (fleta (fl s))])]
[else
(fleta (fl s))]))
@ -198,7 +198,7 @@ P. Borwein. An Efficient Algorithm for the Riemann Zeta Function.
[(integer? s)
(cond [(zero? s) -1/2]
[(negative? s) (define k (- 1 s))
(- (/ (bernoulli k) k))]
(- (/ (bernoulli-number k) k))]
[(eqv? s 1) (raise-argument-error 'zeta "Real, not One" s)]
[else (flzeta (fl s))])]
[else

8
collects/math/private/number-theory/bernoulli.rkt
View File

@ -5,7 +5,7 @@
"factorial.rkt"
"binomial.rkt")
(provide bernoulli)
(provide bernoulli-number)
;; Number of globally memoized Bernoulli numbers
(define num-global-bs 200)
@ -76,7 +76,7 @@
[else (error 'unreachable-code)])]))))
(bern n))
(: bernoulli (Integer -> Exact-Rational))
(define (bernoulli n)
(cond [(n . < . 0) (raise-argument-error 'bernoulli "Natural" n)]
(: bernoulli-number (Integer -> Exact-Rational))
(define (bernoulli-number n)
(cond [(n . < . 0) (raise-argument-error 'bernoulli-number "Natural" n)]
[else (bernoulli* n)]))

8
collects/math/private/number-theory/farey.rkt
View File

@ -2,16 +2,16 @@
(require "types.rkt")
(provide farey mediant)
(provide farey-sequence mediant)
(: mediant : Exact-Rational Exact-Rational -> Exact-Rational)
(define (mediant x y)
(/ (+ (numerator x) (numerator y))
(+ (denominator x) (denominator y))))
(: farey : Integer -> (Listof Exact-Rational))
(define (farey n)
(cond [(n . <= . 0) (raise-argument-error 'farey "Positive-Integer" n)]
(: farey-sequence : Integer -> (Listof Exact-Rational))
(define (farey-sequence n)
(cond [(n . <= . 0) (raise-argument-error 'farey-sequence "Positive-Integer" n)]
[else
(let loop ([a 1] [b 1] [c (sub1 n)] [d n] [#{fs : (Listof Exact-Rational)} '()])
(let ([fs (cons (/ a b) fs)])

44
collects/math/private/number-theory/fibonacci.rkt
View File

@ -2,8 +2,8 @@
(provide make-fibonacci
fibonacci
make-fibonacci/mod
fibonacci/mod)
make-modular-fibonacci
modular-fibonacci)
(: generator : Integer Integer Natural Natural Natural -> Integer)
(define (generator a b p q count)
@ -29,28 +29,30 @@
(define fibonacci (make-fibonacci 0 1))
(: generator/mod : Integer Integer Natural Natural Natural Positive-Integer -> Integer)
(define (generator/mod a b p q count mod)
(: modular-generator : Integer Integer Natural Natural Natural Positive-Integer -> Integer)
(define (modular-generator a b p q count mod)
(cond
[(zero? count) (modulo b mod)]
[(even? count)
(generator/mod a b
(modulo (+ (* p p) (* q q)) mod)
(modulo (+ (* 2 p q) (* q q)) mod)
(quotient count 2)
mod)]
(modular-generator
a b
(modulo (+ (* p p) (* q q)) mod)
(modulo (+ (* 2 p q) (* q q)) mod)
(quotient count 2)
mod)]
[else
(generator/mod (modulo (+ (* b q) (* a q) (* a p)) mod)
(modulo (+ (* b p) (* a q)) mod)
p
q
(- count 1)
mod)]))
(modular-generator
(modulo (+ (* b q) (* a q) (* a p)) mod)
(modulo (+ (* b p) (* a q)) mod)
p
q
(- count 1)
mod)]))
(: make-fibonacci/mod (Integer Integer -> (Integer Integer -> Integer)))
(define ((make-fibonacci/mod a b) n mod)
(cond [(n . < . 0) (raise-argument-error 'fibonacci "Natural" 0 n mod)]
[(mod . <= . 0) (raise-argument-error 'fibonacci "Positive-Integer" 1 n mod)]
[else (generator/mod b a 0 1 n mod)]))
(: make-modular-fibonacci (Integer Integer -> (Integer Integer -> Integer)))
(define ((make-modular-fibonacci a b) n mod)
(cond [(n . < . 0) (raise-argument-error 'modular-fibonacci "Natural" 0 n mod)]
[(mod . <= . 0) (raise-argument-error 'modular-fibonacci "Positive-Integer" 1 n mod)]
[else (modular-generator b a 0 1 n mod)]))
(define fibonacci/mod (make-fibonacci/mod 0 1))
(define modular-fibonacci (make-modular-fibonacci 0 1))

170
collects/math/private/number-theory/number-theory.rkt
View File

@ -66,11 +66,11 @@
;;; Powers
;;;
(: max-dividing-power : Z Z -> N)
(: max-dividing-power : Integer Integer -> Natural)
; (max-dividing-power p n) = m <=> p^m | n and p^(m+1) doesn't divide n
; In Mathematica this one is called IntegerExponent
(define (max-dividing-power p n)
(: find-start : Z Z -> Z)
(: find-start : Integer Integer -> Integer)
(define (find-start p-to-e e)
;(display (list 'fs 'p-to-e p-to-e 'e e)) (newline)
; p-to-e divides n and p-to-e = p^e
@ -81,7 +81,7 @@
(find-start p-to-e2 (* 2 e))
(* 2 e))]
[else (find-power p-to-e e)])))
(: find-power : Z Z -> Z)
(: find-power : Integer Integer -> Integer)
(define (find-power p-to-e e)
;(display (list 'fp 'p-to-e p-to-e 'e e)) (newline)
; p-to-e <= n < (square p-to-e)
@ -90,10 +90,10 @@
[(not (divides? p n)) 0]
[else (assert (find-start p 1) natural?)]))
(: max-dividing-power-naive : Z Z -> N)
(: max-dividing-power-naive : Integer Integer -> Natural)
(define (max-dividing-power-naive p n)
; sames as max-dividing-power but using naive algorithm
(: loop : Z Z -> Z)
(: loop : Integer Integer -> Integer)
(define (loop p-to-e e)
(if (divides? p-to-e n)
(loop (* p p-to-e) (+ e 1))
@ -110,15 +110,15 @@
; Example : (solve-chinese '(2 3 2) '(3 5 7)) = 23
(: solve-chinese : Zs (Listof Z) -> N)
(: solve-chinese : (Listof Integer) (Listof Integer) -> Natural)
(define (solve-chinese as ns)
(unless (andmap positive? ns)
(raise-argument-error 'solve-chinese "(Listof Positive-Integer)" 1 as ns))
; the ns should be coprime
(let* ([n (apply * ns)]
[cs (map (λ: ([ni : Z]) (quotient n ni)) ns)]
[cs (map (λ: ([ni : Integer]) (quotient n ni)) ns)]
[ds (map modular-inverse cs ns)]
[es (cast ds integers?)])
[es (cast ds (make-predicate (Listof Integer)))])
(cast (modulo (apply + (map * as cs es)) n) natural?)))
;;;
@ -137,7 +137,7 @@
; 'composite (with at least probability 1/2) if n is a composite non-Carmichael number
; a proper divisor of n (with at least probability 1/2) if n is a Carmichael number
; [MCA, p.509 - Algorithm 18.5]
(: prime-strong-pseudo-single? : Integer -> (U 'probably-prime 'composite N))
(: prime-strong-pseudo-single? : Integer -> (U 'probably-prime 'composite Natural))
(define (prime-strong-pseudo-single? n)
(cond
[(n . <= . 0) (raise-argument-error 'prime-strong-pseudo-single? "Positive-Integer" n)]
@ -170,19 +170,19 @@
[(= n 1) 'composite]
[else 'probably-prime]))
(define-type Strong-Test-Result (U 'very-probably-prime 'composite N))
(define-type Strong-Test-Result (U 'very-probably-prime 'composite Natural))
(: prime-strong-pseudo/explanation : N -> Strong-Test-Result)
(: prime-strong-pseudo/explanation : Natural -> Strong-Test-Result)
(define (prime-strong-pseudo/explanation n)
; run the strong test several times to improve probability
(: loop : Z (U Strong-Test-Result 'probably-prime) -> Strong-Test-Result)
(: loop : Integer (U Strong-Test-Result 'probably-prime) -> Strong-Test-Result)
(define (loop trials result)
(cond [(= trials 0) 'very-probably-prime]
[(eq? result 'probably-prime) (loop (sub1 trials) (prime-strong-pseudo-single? n))]
[else result]))
(loop prime-strong-pseudo-trials (prime-strong-pseudo-single? n)))
(: prime-strong-pseudo? : N -> Boolean)
(: prime-strong-pseudo? : Natural -> Boolean)
(define (prime-strong-pseudo? n)
(let ([explanation (prime-strong-pseudo/explanation n)])
(or (eq? explanation 'very-probably-prime)
@ -209,7 +209,8 @@
(prime-strong-pseudo? n))))))
(: next-prime : (case-> (N -> N) (Z -> Z)) )
(: next-prime : (case-> (Natural -> Natural)
(Integer -> Integer)))
(define (next-prime n)
(cond
[(negative? n) (- (prev-prime (abs n)))]
@ -225,16 +226,16 @@
n+2
(next-prime n+2)))]))
(: untyped-next-prime : Z -> Z)
(: untyped-next-prime : Integer -> Integer)
(define (untyped-next-prime z)
(next-prime z))
(: untyped-prev-prime : Z -> Z)
(: untyped-prev-prime : Integer -> Integer)
(define (untyped-prev-prime z)
(prev-prime z))
(: prev-prime : Z -> Z)
(: prev-prime : Integer -> Integer)
(define (prev-prime n)
(cond
[(negative? n) (- (next-prime (abs n)))]
@ -250,12 +251,12 @@
(prev-prime n-2)))]))
(: next-primes : Z Z -> Zs)
(: next-primes : Integer Integer -> (Listof Integer))
(define (next-primes m primes-wanted)
(cond
[(primes-wanted . < . 0) (raise-argument-error 'next-primes "Natural" 1 m primes-wanted)]
[else
(: loop : Z Z -> Zs)
(: loop : Integer Integer -> (Listof Integer))
(define (loop n primes-wanted)
(if (= primes-wanted 0)
'()
@ -265,12 +266,12 @@
'()))))
(loop m primes-wanted)]))
(: prev-primes : Z Z -> Zs)
(: prev-primes : Integer Integer -> (Listof Integer))
(define (prev-primes m primes-wanted)
(cond
[(primes-wanted . < . 0) (raise-argument-error 'prev-primes "Natural" 1 m primes-wanted)]
[else
(: loop : Z Z -> Zs)
(: loop : Integer Integer -> (Listof Integer))
(define (loop n primes-wanted)
(if (= primes-wanted 0)
'()
@ -281,14 +282,14 @@
(loop m primes-wanted)]))
(: nth-prime : Z -> Prime)
(: nth-prime : Integer -> Natural)
(define (nth-prime n)
(cond [(n . < . 0) (raise-argument-error 'nth-prime "Natural" n)]
[else
(for/fold: ([p : Prime 2]) ([m (in-range n)])
(for/fold: ([p : Natural 2]) ([m (in-range n)])
(next-prime p))]))
(: random-prime : Z -> Prime)
(: random-prime : Integer -> Natural)
(define (random-prime n)
(when (<= n 2)
(raise-argument-error 'random-prime "Natural > 2" n))
@ -301,25 +302,25 @@
;;; FACTORIZATION
;;;
(: factorize : N -> (Listof (List N N)))
(: factorize : Natural -> (Listof (List Natural Natural)))
(define (factorize n)
(if (< n *SMALL-PRIME-LIMIT*) ; NOTE: Do measurement of best cut
(factorize-small n)
(factorize-large n)))
(: defactorize : (Listof (List N N)) -> N)
(: defactorize : (Listof (List Natural Natural)) -> Natural)
(define (defactorize bes)
(cond [(empty? bes) 1]
[else (define be (first bes))
(* (expt (first be) (second be))
(defactorize (rest bes)))]))
(: factorize-small : N -> (Listof (List N N)))
(: factorize-small : Natural -> (Listof (List Natural Natural)))
(define (factorize-small n)
; fast for small n, but works correctly for large n too
(small-prime-factors-over n 2))
(: small-prime-factors-over : N Prime -> (Listof (List N N)))
(: small-prime-factors-over : Natural Natural -> (Listof (List Natural Natural)))
; Factor a number n without prime factors below the prime p.
(define (small-prime-factors-over n p) ; p prime
(cond
@ -338,7 +339,7 @@
;;; ALGORITHM 19.8 Pollard's rho method
; INPUT n>=3 neither a prime nor a perfect power
; OUTPUT Either a proper divisor of n or #f
(: pollard : N -> (U N False))
(: pollard : Natural -> (U Natural False))
(define (pollard n)
(let ([x0 (random-natural n)])
(do ([xi x0 (remainder (+ (* xi xi) 1) n)]
@ -360,7 +361,7 @@
[(even? n) `((2 1) ,@(pollard-factorize (quotient n 2)))]
[(divides? 3 n) `((3 1) ,@(pollard-factorize (quotient n 3)))]
[(simple-perfect-power n)
=> (λ: ([base-and-exp : (List N N)])
=> (λ: ([base-and-exp : (List Natural Natural)])
(cond
[(prime? (car base-and-exp)) (list base-and-exp)]
[else (map (λ: ([b-and-e : (List Natural Natural)])
@ -374,19 +375,20 @@
(pollard-factorize (quotient n divisor)))
(loop (pollard n))))])))
(: factorize-large : N -> (Listof (List N N)))
(: factorize-large : Natural -> (Listof (List Natural Natural)))
(define (factorize-large n)
(combine-same-base
(sort (pollard-factorize n) base-and-exponent<?)))
(: base-and-exponent<? : (U N (List N N)) (U N (List N N)) -> Boolean)
(: base-and-exponent<? ((U Natural (List Natural Natural)) (U Natural (List Natural Natural))
-> Boolean))
(define (base-and-exponent<? x y)
(let ([id-or-first
(λ: ([x : (U Integer (List Integer Integer))])
(if (number? x) x (first x)))])
(<= (id-or-first x) (id-or-first y))))
(: combine-same-base : (Listof (List N N)) -> (Listof (List N N)))
(: combine-same-base : (Listof (List Natural Natural)) -> (Listof (List Natural Natural)))
(define (combine-same-base list-of-base-and-exponents)
; list-of-base-and-exponents must be sorted
(let ([l list-of-base-and-exponents])
@ -408,7 +410,7 @@
; find-tail pred clist -> pair or false
; Return the first pair of clist whose car satisfies pred. If no pair does, return false.
(: find-tail : (Z -> Boolean) Zs -> (U False Zs))
(: find-tail : (Integer -> Boolean) (Listof Integer) -> (U False (Listof Integer)))
(define (find-tail pred xs)
(cond [(empty? xs) #f]
[(pred (car xs)) xs]
@ -419,14 +421,14 @@
;;; Powers
;;;
(: as-power : N+ -> (Values N N))
(: as-power : Exact-Positive-Integer -> (Values Natural Natural))
; Write a>0 as b^r with r maximal. Return b and r.
(define (as-power a)
(let ([r (apply gcd ((inst map N (List N N)) second (factorize a)))])
(let ([r (apply gcd ((inst map Natural (List Natural Natural)) second (factorize a)))])
(values (integer-root a r) r)))
(: prime-power : N -> (U (List Prime N) False))
(: prime-power : Natural -> (U (List Natural Natural) False))
; if n is a prime power, return list of prime and exponent in question,
; otherwise return #f
(define (prime-power n)
@ -435,24 +437,24 @@
(first (prime-divisors/exponents n))
#f)))
(: prime-power? : N -> Boolean)
(: prime-power? : Natural -> Boolean)
; Is n of the form p^m, with p is prime?
(define (prime-power? n)
(and (prime-power n) #t))
(: odd-prime-power? : N -> Boolean)
(: odd-prime-power? : Natural -> Boolean)
(define (odd-prime-power? n)
(let ([p/e (prime-power n)])
(and p/e
(odd? (first p/e)))))
(: perfect-power? : N -> Boolean)
(: perfect-power? : Natural -> Boolean)
(define (perfect-power? a)
(and (not (zero? a))
(let-values ([(base n) (as-power a)])
(and (> n 1) (> a 1)))))
(: simple-perfect-power : N -> (U (List N N) False))
(: simple-perfect-power : Natural -> (U (List Natural Natural) False))
(define (simple-perfect-power a)
; simple-perfect-power is used by pollard-fatorize
(and (not (zero? a))
@ -461,7 +463,7 @@
(list base n)
#f))))
(: perfect-power : N -> (U (List N N) False))
(: perfect-power : Natural -> (U (List Natural Natural) False))
; if a = b^n with b>1 and n>1
(define (perfect-power a)
(and (not (zero? a))
@ -470,50 +472,50 @@
(list base n)
#f))))
(: perfect-square : N -> (U N False))
(: perfect-square : Natural -> (U Natural False))
(define (perfect-square n)
(let ([sqrt-n (integer-sqrt n)])
(if (= (* sqrt-n sqrt-n) n)
sqrt-n
#f)))
(: powers-of : N N -> (Listof N))
(: powers-of : Natural Natural -> (Listof Natural))
; returns a list of numbers: a^0, ..., a^n
(define (powers-of a n)
(let: loop : (Listof N)
([i : N 0]
[a^i : N 1])
(let: loop : (Listof Natural)
([i : Natural 0]
[a^i : Natural 1])
(if (<= i n)
(cons a^i (loop (+ i 1) (* a^i a)))
'())))
(define prime-divisors/exponents factorize)
(: prime-divisors : N -> (Listof Prime))
(: prime-divisors : Natural -> (Listof Natural))
; return list of primes in a factorization of n
(define (prime-divisors n)
(map (inst car N (Listof N))
(map (inst car Natural (Listof Natural))
(prime-divisors/exponents n)))
(: prime-exponents : N -> (Listof N))
(: prime-exponents : Natural -> (Listof Natural))
; return list of exponents in a factorization of n
(define (prime-exponents n)
(map (inst cadr N N (Listof N))
(map (inst cadr Natural Natural (Listof Natural))
(prime-divisors/exponents n)))
(: prime-omega : N -> N)
(: prime-omega : Natural -> Natural)
; http://reference.wolfram.com/mathematica/ref/PrimeOmega.html
(define (prime-omega n)
(for/fold: ([sum : Natural 0]) ([e (in-list (prime-exponents n))])
(+ sum e)))
(: integer-root/remainder : N N -> (Values N N))
(: integer-root/remainder : Natural Natural -> (Values Natural Natural))
(define (integer-root/remainder a n)
(let ([i (integer-root a n)])
(values i (assert (- a (expt i n)) natural?))))
(: integer-root : N N -> N)
(: integer-root : Natural Natural -> Natural)
(define (integer-root x y)
; y'th root of x
(cond
@ -540,7 +542,7 @@
(let* ([top-bits (arithmetic-shift x (- (* length/y/2 y)))]
[nth-root-top-bits (integer-root top-bits y)])
(arithmetic-shift (+ nth-root-top-bits 1) length/y/2))])
(let: loop : Z ([g : Z init-g])
(let: loop : Integer ([g : Integer init-g])
(let* ([a (expt g (assert (- y 1) natural?))]
[b (* a y)]
[c (* a (- y 1))]
@ -560,14 +562,14 @@
natural?)]))
(: simple-as-power : N+ -> (Values N N))
(: simple-as-power : Exact-Positive-Integer -> (Values Natural Natural))
; For a>0 write it as a = b^r where r maximal
; return (values b r)
(define (simple-as-power a)
; (displayln (list 'simple-as-power a))
; Note: The simple version is used by pollard-factorize
(let: loop : (Values N N)
([n : N (integer-length a)])
(let: loop : (Values Natural Natural)
([n : Natural (integer-length a)])
(let-values ([(root rem) (integer-root/remainder a (add1 n))])
(if (zero? rem)
(values root (assert (add1 n) natural?))
@ -575,20 +577,20 @@
(loop (sub1 n))
(error 'simple-as-power "internal error"))))))
(: prime-power? : N -> Boolean)
(: prime-power? : Natural -> Boolean)
;;;
;;; DIVISORS
;;;
(: divisors : Z -> (Listof N))
(: divisors : Integer -> (Listof Natural))
; return the positive divisorts of n
(define (divisors n)
(cond [(zero? n) '()]
[else (define n+ (if (positive? n) n (- n)))
(sort (factorization->divisors (factorize n+)) <)]))
(: factorization->divisors : (Listof (List N N)) -> (Listof N))
(: factorization->divisors : (Listof (List Natural Natural)) -> (Listof Natural))
(define (factorization->divisors f)
(cond
[(null? f) '(1)]
@ -598,8 +600,8 @@
; f = p^n * g
(let ([divisors-of-g (factorization->divisors g)])
(apply append
((inst map (Listof N) N)
(λ: ([p^i : N]) (map (λ: ([d : N]) (* p^i d)) divisors-of-g))
((inst map (Listof Natural) Natural)
(λ: ([p^i : Natural]) (map (λ: ([d : Natural]) (* p^i d)) divisors-of-g))
(powers-of p n)))))]))
;;;
@ -622,11 +624,11 @@
; phi(n) = n * product (1 - ---- )
; i=1 pi
(: totient : N -> N)
(: totient : Natural -> Natural)
(define (totient n)
(let ((ps (prime-divisors n)))
(assert (* (quotient n (apply * ps))
(apply * (map (λ: ([p : N]) (sub1 p)) ps)))
(apply * (map (λ: ([p : Natural]) (sub1 p)) ps)))
natural?)))
(: every : (All (A) (A -> Boolean) (Listof A) -> Boolean))
@ -640,47 +642,47 @@
; mu(n) = 1 if n is a product of an even number of primes
; = -1 if n is a product of an odd number of primes
; = 0 if n has a multiple prime factor
(: moebius-mu : N -> (U -1 0 1))
(: moebius-mu : Natural -> (U -1 0 1))
(define (moebius-mu n)
(: one? : Z -> Boolean)
(: one? : Integer -> Boolean)
(define (one? x) (= x 1))
(define f (factorize n))
(define exponents ((inst map N (List N N)) second f))
(define exponents ((inst map Natural (List Natural Natural)) second f))
(cond
[(every one? exponents)
(define primes ((inst map N (List N N)) first f))
(define primes ((inst map Natural (List Natural Natural)) first f))
(if (even? (length primes))
1 -1)]
[else 0]))
(: divisor-sum : (case-> (N -> N) (N N -> N)))
(: divisor-sum : (case-> (Natural -> Natural) (Natural Natural -> Natural)))
(define divisor-sum
; returns the sum of the kth power of all divisors of n
(let ()
(case-lambda
[(n) (divisor-sum n 1)]
[(n k) (let* ([f (factorize n)]
[ps ((inst map N (List N N)) first f)]
[es ((inst map N (List N N)) second f)])
(: divisor-sum0 : Any N -> N)
[ps ((inst map Natural (List Natural Natural)) first f)]
[es ((inst map Natural (List Natural Natural)) second f)])
(: divisor-sum0 : Any Natural -> Natural)
(define (divisor-sum0 p e) (+ e 1))
(: divisor-sum1 : N N -> N)
(: divisor-sum1 : Natural Natural -> Natural)
(define (divisor-sum1 p e)
(let: loop : N
([sum : N 1]
[n : N 0]
[p-to-n : N 1])
(let: loop : Natural
([sum : Natural 1]
[n : Natural 0]
[p-to-n : Natural 1])
(cond [(= n e) sum]
[else (let ([t (* p p-to-n)])
(loop (+ t sum) (+ n 1) t))])))
(: divisor-sumk : N N -> N)
(: divisor-sumk : Natural Natural -> Natural)
(define (divisor-sumk p e)
(let ([p-to-k (expt p k)])
(let: loop : N
([sum : N 1]
[n : N 0]
[p-to-kn : N 1])
(let: loop : Natural
([sum : Natural 1]
[n : Natural 0]
[p-to-kn : Natural 1])
(cond [(= n e) sum]
[else (let ([t (* p-to-k p-to-kn)])
(loop (+ t sum) (+ n 1) t))]))))
@ -691,7 +693,7 @@
ps es))
natural?))])))
(: mangoldt-lambda : Z -> Real)
(: mangoldt-lambda : Integer -> Real)
(define (mangoldt-lambda n)
(cond
[(<= n 0) (raise-argument-error 'mangoldt-lambda "Natural" n)]

27
collects/math/private/number-theory/primitive-roots.rkt
View File

@ -12,8 +12,8 @@
prime-power))
(provide unit-group
order
orders
unit-group-order
unit-group-orders
exists-primitive-root?
primitive-root?
primitive-root
@ -34,11 +34,12 @@
[else (filter (λ: ([m : Natural]) (coprime? m n))
(build-list (- n 1) add1))]))
(: order : Integer Integer -> Positive-Integer)
(define (order g n)
(cond [(g . <= . 0) (raise-argument-error 'order "Positive-Integer" 0 g n)]
[(n . <= . 0) (raise-argument-error 'order "Positive-Integer" 1 g n)]
[(not (coprime? g n)) (error 'order "expected coprime arguments; given ~e and ~e" g n)]
(: unit-group-order : Integer Integer -> Positive-Integer)
(define (unit-group-order g n)
(cond [(g . <= . 0) (raise-argument-error 'unit-group-order "Positive-Integer" 0 g n)]
[(n . <= . 0) (raise-argument-error 'unit-group-order "Positive-Integer" 1 g n)]
[(not (coprime? g n))
(error 'unit-group-order "expected coprime arguments; given ~e and ~e" g n)]
[else
(with-modulus n
(let: loop : Positive-Integer ([k : Positive-Integer 1]
@ -46,10 +47,10 @@
(cond [(mod= a 1) k]
[else (loop (+ k 1) (mod* a g))])))]))
(: orders : Integer -> (Listof Positive-Integer))
(define (orders n)
(cond [(n . <= . 0) (raise-argument-error 'orders "Positive-Integer" n)]
[else (map (λ: ([m : Positive-Integer]) (order m n))
(: unit-group-orders : Integer -> (Listof Positive-Integer))
(define (unit-group-orders n)
(cond [(n . <= . 0) (raise-argument-error 'unit-group-orders "Positive-Integer" n)]
[else (map (λ: ([m : Positive-Integer]) (unit-group-order m n))
(unit-group n))]))
; DEFINITION (Primitive Root)
@ -60,7 +61,7 @@
(define (primitive-root? g n)
(if (not (coprime? g n))
(error 'primitive-root? "expected coprime arguments; given ~e and ~e" g n)
(= (order g n) (phi n))))
(= (unit-group-order g n) (phi n))))
; THEOREM (Existence of primitive roots)
; Un is cyclic (i.e. have a primitive root)
@ -149,7 +150,7 @@
(define p (if pp (first pp) (error 'primitive-root "internal error")))
(define gg (primitive-root p))
(define g (or gg (error 'primitive-root "internal error")))
(if (= (order g (* p p)) (totient (* p p)))
(if (= (unit-group-order g (* p p)) (totient (* p p)))
g
(modulo (+ g p) n))]
; U_2p^e , p odd

27
collects/math/private/number-theory/types.rkt
View File

@ -1,26 +1,6 @@
#lang typed/racket
(provide N N+ Z Q Ns Zs Base-Exponent Factorization Prime
cast
natural? naturals? Integer? integers? exact-zero?)
;;;
;;; Types and Predicates
;;;
(define-type N Natural)
(define-type N+ Exact-Positive-Integer)
(define-type Z Integer)
(define-type Q Exact-Rational)
(define-type Ns (Listof N))
(define-type Zs (Listof Z))
(define-type Base-Exponent (List N N))
(define-type BE Base-Exponent)
(define-type Factorization (List Base-Exponent))
(define-type Prime N) ; non-checked (for documentation purposes)
(provide cast natural? exact-zero?)
(define-syntax (cast stx) (syntax-case stx () [(_ . more) #'(assert . more)]))
; Note: (cast val predicate) is used in the code, where
@ -28,8 +8,5 @@
; can prove it. Replace assert with a "proper" cast when
; it appears in Typed Racket.
(define-predicate natural? N) ; Note: 0 is natural
(define-predicate naturals? Ns)
(define-predicate Integer? Z)
(define-predicate integers? Zs)
(define natural? exact-nonnegative-integer?)
(define-predicate exact-zero? Zero)

36
collects/math/scribblings/math-bigfloat.scrbl
View File

@ -26,12 +26,10 @@ a C library that provides
The arbitrary-precision floating-point numbers MPFR provides and operates on are
represented by the type @racket[Bigfloat].
MPFR is free, license-compatible with commercial software, and
@hyperlink["http://www.mpfr.org/ports.html"]{easy to install} on all
major platforms. Even so, not all systems have it, so @racketmodname[math/bigfloat] links
against MPFR lazily.
If you know MPFR is installed but can't use @racketmodname[math/bigfloat], see @secref{loading}.
MPFR is free and license-compatible with commercial software. It is distributed with Racket
for Windows and Mac OS X, is installed on most Linux systems, and is
@hyperlink["http://www.mpfr.org/ports.html"]{easy to install} on all major platforms.
See @secref{loading} for details.
@section[#:tag "quick"]{Quick Start}
@ -711,23 +709,19 @@ Canonicalizing bigfloats won't change answers computed from them.
@section[#:tag "loading"]{Loading the MPFR Library}
The @racketmodname[math/bigfloat] library depends directly on one library: MPFR, which on the filesystem is
named something like @filepath{libmpfr.so}, but with a platform-specific suffix.
MPFR itself depends on @hyperlink["http://gmplib.org/"]{GMP}, the
GNU Multiple Precision Arithmetic Library, which on the filesystem is named @filepath{libgmp.so}.
If @racket[(mpfr-available?)] is @racket[#t], both of these libraries have been loaded
by the Racket runtime, and everything in @racketmodname[math/bigfloat] should work.
(If something doesn't, it's likely an error in the foreign function interface. In that
case, please email the author of this library or submit a bug report.)
If @racket[(mpfr-available?)] is @racket[#f], one of these libraries can't be loaded. The
most common reason is that one or both needs to be installed. Another possibility is that
both are installed, but one is not in the operating system's search path. In this case,
see @racket[set-mpfr-lib-dirs!].
The @racketmodname[math/bigfloat] library depends directly on external one library:
MPFR, which on the filesystem is named something like @filepath{libmpfr.so}, but with a
platform-specific suffix. MPFR itself depends on @hyperlink["http://gmplib.org/"]{GMP}, the
GNU Multiple Precision Arithmetic Library, which on the filesystem is named something like
@filepath{libgmp.so}.
@defproc[(mpfr-available?) Boolean]{
Returns @racket[#t] when both GMP and MPFR are available; @racket[#f] otherwise.
If @racket[(mpfr-available?)] is @racket[#t], both MPFR and GMP have been loaded
by the Racket runtime, and everything in @racketmodname[math/bigfloat] should work.
If @racket[(mpfr-available?)] is @racket[#t] and something doesn't work, it's likely an error
in the foreign function interface. In that case, please email the author of this library or
submit a bug report.
}
@(close-eval untyped-eval)

5
collects/math/scribblings/math-flonum.scrbl
View File

@ -230,12 +230,11 @@ For example, say we want the probability density of the standard normal distribu
(the bell curve) at 50 standard deviations from zero:
@interaction[#:eval untyped-eval
(require math/distributions)
(define d (normal-dist))
((dist-pdf d) 50.0)]
(pdf (normal-dist) 50.0)]
Mathematically, the density is nonzero everywhere, but the density at 50 is less than
@racket[+min.0]. However, its density in log space, or its log-density, is representable:
@interaction[#:eval untyped-eval
((dist-pdf d) 50.0 #t)]
(pdf (normal-dist) 50.0 #t)]
While this example may seem contrived, it is very common, when computing the density
of a @italic{vector} of data, for the product of the densities to be too small to represent directly.

49
collects/math/scribblings/math-number-theory.scrbl
View File

@ -586,13 +586,13 @@ Otherwise 0 is returned.
@section[#:tag "number-sequences"]{Number Sequences}
@margin-note{Wikipedia: @hyperlink["http://en.wikipedia.org/wiki/Bernoulli_number"]{Bernoulli Number}}
@defproc[(bernoulli [n Integer]) Exact-Rational]{
@defproc[(bernoulli-number [n Integer]) Exact-Rational]{
Returns the @racket[n]th Bernoulli number; @racket[n] must be nonnegative.
@interaction[#:eval untyped-eval
(map bernoulli (range 9))]
(map bernoulli-number (range 9))]
Note that these are the @italic{first} Bernoulli numbers, since @racket[(bernoulli 1) = -1/2].
Note that these are the @italic{first} Bernoulli numbers, since @racket[(bernoulli-number 1) = -1/2].
}
@margin-note{MathWorld: @hyperlink["http://mathworld.wolfram.com/EulerianNumber.html"]{Eulerian Number}}
@ -612,26 +612,41 @@ Otherwise 0 is returned.
(map fibonacci (range 10))]
}
@defproc[(fibonacci/mod [n Integer] [m Integer]) Natural]{
@defproc[(make-fibonacci [a Integer] [b Integer]) (Integer -> Integer)]{
Returns a function representing @italic{a} Fibonacci sequence with the first two numbers
@racket[a] and @racket[b]. The @racket[fibonacci] function is defined as
@racket[(make-fibonacci 0 1)].
@margin-note{Wikipedia: @hyperlink["http://wikipedia.org/wiki/Lucas_number"]{Lucas Number}}
The Lucas numbers are defined as a Fibonacci sequence starting with 2 and 1:
@interaction[#:eval untyped-eval
(map (make-fibonacci 2 1) (range 10))]
}
@defproc[(modular-fibonacci [n Integer] [m Integer]) Natural]{
Returns the @racket[n]th Fibonacci number modulo @racket[m]; @racket[n] must be nonnegative
and @racket[m] must be positive.
The ten first Fibonacci numbers modulo 5.
@interaction[#:eval untyped-eval
(map (λ (n) (fibonacci/mod n 5)) (range 10))]
(map (λ (n) (modular-fibonacci n 5)) (range 10))]
}
@defproc[(make-modular-fibonacci [a Integer] [b Integer]) (Integer Integer -> Integer)]{
Like @racket[make-fibonacci], but makes a modular Fibonacci sequence.
}
@margin-note{Wikipedia: @hyperlink["http://en.wikipedia.org/wiki/Farey_sequence"]{Farey Sequence}}
@defproc[(farey [n Integer]) (Listof Exact-Rational)]{
@defproc[(farey-sequence [n Integer]) (Listof Exact-Rational)]{
Returns a list of the numbers in the @racket[n]th Farey sequence; @racket[n] must be positive.
The @racket[n]th Farey sequence is the sequence of all
completely reduced rational numbers from 0 to 1 which denominators
are less than or equal to @racket[n].
@interaction[#:eval untyped-eval
(farey 1)
(farey 2)
(farey 3)]
(farey-sequence 1)
(farey-sequence 2)
(farey-sequence 3)]
}
@margin-note{MathWorld: @hyperlink["http://mathworld.wolfram.com/TangentNumber.html"]{Tangent Number}}
@ -806,21 +821,21 @@ modulus @racket[n] must be positive.
(unit-group 6)]
}
@defproc[(order [x Integer] [n Integer]) Positive-Integer]{
@defproc[(unit-group-order [x Integer] [n Integer]) Positive-Integer]{
Returns the order of @racket[x] in the group @math-style{Un}; both arguments must be positive.
If @racket[x] and @racket[n] are not coprime, @racket[(order x n)] raises an error.
If @racket[x] and @racket[n] are not coprime, @racket[(unit-group-order x n)] raises an error.
@interaction[#:eval untyped-eval
(order 2 5)
(order 2 6)]
(unit-group-order 2 5)
(unit-group-order 2 6)]
}
@defproc[(orders [n Integer]) (Listf Positive-Integer)]{
Returns a list @racket[(list (order x0 n) (order x1 n) ...)] where
@defproc[(unit-group-orders [n Integer]) (Listf Positive-Integer)]{
Returns a list @racket[(list (unit-group-order x0 n) (unit-group-order x1 n) ...)] where
@racket[x0], @racket[x1], ... are the elements of @math-style{Un}.
The modulus @racket[n] must be positive.
@interaction[#:eval untyped-eval
(orders 5)
(map (curryr order 5) (unit-group 5))]
(unit-group-orders 5)
(map (curryr unit-group-order 5) (unit-group 5))]
}
@defproc[(primitive-root? [x Integer] [n Integer]) Boolean]{

16
collects/math/tests/number-theory-tests.rkt
View File

@ -19,13 +19,13 @@
; "primitive-roots.rkt"
(check-equal? (unit-group 20) '(1 3 7 9 11 13 17 19)) ; 19 !!!!
(check-equal? (order 19 20) 2)
(check-equal? (order 3 20) 4)
(check-equal? (orders 20) '(1 4 4 2 2 4 4 2)) ; (order 3 20)=4, ...
(check-equal? (unit-group-order 19 20) 2)
(check-equal? (unit-group-order 3 20) 4)
(check-equal? (unit-group-orders 20) '(1 4 4 2 2 4 4 2)) ; (unit-group-order 3 20)=4, ...
(check-true (andmap exists-primitive-root? '(1 2 4 3 9 6 18)))
(check-false (ormap exists-primitive-root? '(8 16 12)))
(check-equal? (primitive-root 20) #f)
(check-equal? (primitive-root 10) 7) ; (length (unit-group 10)) = (order 7 10)
(check-equal? (primitive-root 10) 7) ; (length (unit-group 10)) = (unit-group-order 7 10)
(check-true (primitive-root? 7 10))
(check-false (primitive-root? 7 20))
(check-equal? (primitive-roots 10) '(3 7))
@ -33,7 +33,7 @@
(define (find-and-check-root n)
(define r (primitive-root n))
(cond [(not r) #t]
[else (= (length (unit-group n)) (order r n))]))
[else (= (length (unit-group n)) (unit-group-order r n))]))
(check-true (andmap find-and-check-root '(1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 78125)))
;"polygonal.rkt"
@ -50,7 +50,7 @@
(check-true (andmap octagonal-number? '(0 1 8 21 40 65)))
; "farey.rkt"
(check-equal? (farey 5) '(0 1/5 1/4 1/3 2/5 1/2 3/5 2/3 3/4 4/5 1))
(check-equal? (farey-sequence 5) '(0 1/5 1/4 1/3 2/5 1/2 3/5 2/3 3/4 4/5 1))
(check-equal? (mediant 1/1 1/2) 2/3)
; "fibonacci.rkt"
@ -59,7 +59,7 @@
(for*: ([a (in-range -5 6)]
[b (in-range -5 6)]
[mod (in-range 1 8)])
(check-equal? (build-list 20 (λ: ([n : Integer]) ((make-fibonacci/mod a b) n mod)))
(check-equal? (build-list 20 (λ: ([n : Integer]) ((make-modular-fibonacci a b) n mod)))
(build-list 20 (λ: ([n : Integer]) (modulo ((make-fibonacci a b) n) mod)))))
; "partitions.rkt"
@ -68,7 +68,7 @@
; "bernoulli.rkt"
(check-equal? (map bernoulli '(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18))
(check-equal? (map bernoulli-number '(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18))
'(1 -1/2 1/6 0 -1/30 0 1/42 0 -1/30 0 5/66 0 -691/2730 0 7/6 0 -3617/510 0 43867/798))
; "tangent-number.rkt"