-
-
Notifications
You must be signed in to change notification settings - Fork 2
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Added the Array.pairmap{} alias for Array.reduce_pairs{}
- new file: scripts/Expensive/numerical_integration.sf - new file: scripts/RosettaCode/ordered_partitions.sf
- Loading branch information
trizen
committed
Dec 16, 2015
1 parent
93af122
commit ce2fc7b
Showing
4 changed files
with
97 additions
and
1 deletion.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,64 @@ | ||
#!/usr/bin/ruby | ||
|
||
# | ||
## http://rosettacode.org/wiki/Numerical_integration/Gauss-Legendre_Quadrature | ||
# | ||
|
||
func legendre_pair((1), x) { (x, 1) } | ||
func legendre_pair( n, x) { | ||
var (m1, m2) = legendre_pair(n - 1, x); | ||
var u = (1 - 1/n); | ||
((1 + u)*x*m1 - u*m2, m1); | ||
} | ||
|
||
func legendre((0), _) { 1 } | ||
func legendre( n, x) { [legendre_pair(n, x)][0] } | ||
|
||
func legendre_prime({ .is_zero }, _) { 0 } | ||
func legendre_prime({ .is_one }, _) { 1 } | ||
|
||
func legendre_prime(n, x) { | ||
var (m0, m1) = legendre_pair(n, x); | ||
(m1 - x*m0) * n / (1 - x**2); | ||
} | ||
|
||
func approximate_legendre_root(n, k) { | ||
# Approximation due to Francesco Tricomi | ||
var t = ((4*k - 1) / (4*n + 2)); | ||
(1 - ((n - 1)/(8 * n**3))) * (Math.pi * t -> cos); | ||
} | ||
|
||
func newton_raphson(f, f_prime, r, eps = 2e-16) { | ||
while (var dr = (-f(r) / f_prime(r)) -> abs >= eps) { | ||
r += dr; | ||
} | ||
return r; | ||
} | ||
|
||
func legendre_root(n, k) { | ||
newton_raphson(legendre.method(:call, n), legendre_prime.method(:call, n), | ||
approximate_legendre_root(n, k)); | ||
} | ||
|
||
func weight(n, r) { 2 / ((1 - r**2) * legendre_prime(n, r)**2) } | ||
|
||
func nodes(n) { | ||
gather { | ||
take(Pair(0, weight(n, 0))) if n.is_odd; | ||
(n/2 -> int).times { |i| | ||
var r = legendre_root(n, i); | ||
var w = weight(n, r); | ||
take(Pair(r, w), Pair(-r, w)); | ||
} | ||
} | ||
} | ||
|
||
func quadrature(n, f, a, b, nds = nodes(n)) { | ||
func scale(x) { (x*(b - a) + a + b) / 2 } | ||
(b - a) / 2 * nds.map{ .second * f(scale(.first)) }.sum | ||
} | ||
|
||
[@(5..10), 20].each { |i| | ||
printf("Gauss-Legendre %2d-point quadrature ∫₋₃⁺³ exp(x) dx ≈ %.15f\n", | ||
i, quadrature(i, {.exp}, -3, +3)) | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,28 @@ | ||
#!/usr/bin/ruby | ||
|
||
# | ||
## http://rosettacode.org/wiki/Ordered_Partitions#Sidef | ||
# | ||
|
||
func part(_, {.is_empty}) { [[]] } | ||
func partitions({.is_empty}) { [[]] } | ||
|
||
func part(s, args) { | ||
var res = []; | ||
var combs = s.combinations(args[0]); | ||
combs << [] if combs.is_empty; | ||
combs.each { |c| | ||
part(s - c, args.ft(1)).each{|r| res << ([c] + r)} | ||
} | ||
res | ||
} | ||
|
||
func partitions(args) { | ||
part(1..args.sum, args) | ||
} | ||
|
||
[[],[0,0,0],[1,1,1],[2,0,2]].each { |test_case| | ||
say "partitions #{test_case.dump}:" | ||
partitions(test_case).each{|part| say part.dump } | ||
print "\n" | ||
} |