Average loop length: Difference between revisions
added Easylang
m (→{{header|REXX}}: elided an extra blank line.) |
(added Easylang) |
||
| (65 intermediate revisions by 39 users not shown) | |||
Line 2:
Let <code>f</code> be a uniformly-randomly chosen mapping from the numbers 1..N to the numbers 1..N (note: not necessarily a permutation of 1..N; the mapping could produce a number in more than one way or not at all). At some point, the sequence <code>1, f(1), f(f(1))...</code> will contain a <em>repetition</em>, a number that occurring for the second time in the sequence.
;Task:
Write a program or a script that estimates, for each <code>N</code>, the average length until the first such repetition.
Also calculate this expected length using an analytical formula, and optionally compare the simulated result with the theoretical one.
This problem comes from the end of Donald Knuth's [http://www.youtube.com/watch?v=cI6tt9QfRdo Christmas tree lecture 2011].
Line 32 ⟶ 35:
19 5.1312 5.1522 ( 0.41%)
20 5.2699 5.2936 ( 0.45%)</pre>
<br>
=={{header|11l}}==
{{trans|Python}}
<syntaxhighlight lang="11l">F ffactorial(n)
V result = 1.0
L(i) 2..n
result *= i
R result
V MAX_N = 20
V TIMES = 1000000
F analytical(n)
R sum((1..n).map(i -> ffactorial(@n) / pow(Float(@n), Float(i)) / ffactorial(@n - i)))
F test(n, times)
V count = 0
L(i) 0 .< times
V (x, bits) = (1, 0)
L (bits [&] x) == 0
count++
bits [|]= x
x = 1 << random:(n)
R Float(count) / times
print(" n avg exp. diff\n-------------------------------")
L(n) 1 .. MAX_N
V avg = test(n, TIMES)
V theory = analytical(n)
V diff = (avg / theory - 1) * 100
print(‘#2 #3.4 #3.4 #2.3%’.format(n, avg, theory, diff))</syntaxhighlight>
{{out}}
<pre>
n avg exp. diff
-------------------------------
1 1.0000 1.0000 0.000%
2 1.5003 1.5000 0.022%
3 1.8897 1.8889 0.044%
4 2.2170 2.2187 -0.080%
5 2.5099 2.5104 -0.022%
6 2.7736 2.7747 -0.040%
7 3.0182 3.0181 0.001%
8 3.2438 3.2450 -0.037%
9 3.4589 3.4583 0.018%
10 3.6605 3.6602 0.008%
11 3.8517 3.8524 -0.017%
12 4.0373 4.0361 0.032%
13 4.2159 4.2123 0.085%
14 4.3828 4.3820 0.017%
15 4.5465 4.5458 0.016%
16 4.7048 4.7043 0.013%
17 4.8585 4.8579 0.012%
18 5.0042 5.0071 -0.057%
19 5.1465 5.1522 -0.110%
20 5.2907 5.2936 -0.054%
</pre>
=={{header|Ada}}==
<
with Ada.Numerics.Generic_Elementary_Functions;
with Ada.Numerics.Discrete_Random;
Line 86 ⟶ 148:
Put(err, Fore=>3, Aft=>3, exp=>0); New_line;
end loop;
end Avglen;</
{{out}}
<pre>
Line 113 ⟶ 175:
=={{header|BBC BASIC}}==
<
MAX_N = 20
TIMES = 1000000
Line 145 ⟶ 207:
DEF FNfactorial(n)
IF n=1 OR n=0 THEN =1 ELSE = n * FNfactorial(n-1)</
{{out}}
<pre>
Line 171 ⟶ 233:
=={{header|C}}==
<syntaxhighlight lang="c">#include <stdio.h>
#include <stdlib.h>
#include <math.h>
Line 226 ⟶ 289:
}
return 0;
}</
{{out}}
<pre>
Line 251 ⟶ 314:
19 5.1479 5.1522 -0.084%
20 5.2957 5.2936 0.040%
</pre>
=={{header|C sharp|C#}}==
{{trans|Java}}
<syntaxhighlight lang="csharp">public class AverageLoopLength {
private static int N = 100000;
private static double analytical(int n) {
double[] factorial = new double[n + 1];
double[] powers = new double[n + 1];
powers[0] = 1.0;
factorial[0] = 1.0;
for (int i = 1; i <= n; i++) {
factorial[i] = factorial[i - 1] * i;
powers[i] = powers[i - 1] * n;
}
double sum = 0;
for (int i = 1; i <= n; i++) {
sum += factorial[n] / factorial[n - i] / powers[i];
}
return sum;
}
private static double average(int n) {
Random rnd = new Random();
double sum = 0.0;
for (int a = 0; a < N; a++) {
int[] random = new int[n];
for (int i = 0; i < n; i++) {
random[i] = rnd.Next(n);
}
var seen = new HashSet<double>(n);
int current = 0;
int length = 0;
while (seen.Add(current)) {
length++;
current = random[current];
}
sum += length;
}
return sum / N;
}
public static void Main(string[] args) {
Console.WriteLine(" N average analytical (error)");
Console.WriteLine("=== ========= ============ =========");
for (int i = 1; i <= 20; i++) {
var average = AverageLoopLength.average(i);
var analytical = AverageLoopLength.analytical(i);
Console.WriteLine("{0,3} {1,10:N4} {2,13:N4} {3,8:N2}%", i, average, analytical, (analytical - average) / analytical * 100);
}
}
}
</syntaxhighlight>
{{out}}
<pre>
N average analytical (error)
=== ========= ============ =========
1 1.0000 1.0000 0.00%
2 1.4999 1.5000 0.01%
3 1.8860 1.8889 0.15%
4 2.2235 2.2188 -0.22%
5 2.5115 2.5104 -0.04%
6 2.7793 2.7747 -0.17%
7 3.0149 3.0181 0.11%
8 3.2457 3.2450 -0.02%
9 3.4559 3.4583 0.07%
10 3.6558 3.6602 0.12%
11 3.8428 3.8524 0.25%
12 4.0270 4.0361 0.22%
13 4.2111 4.2123 0.03%
14 4.3766 4.3820 0.12%
15 4.5535 4.5458 -0.17%
16 4.6989 4.7043 0.11%
17 4.8590 4.8579 -0.02%
18 4.9972 5.0071 0.20%
19 5.1542 5.1522 -0.04%
20 5.3024 5.2936 -0.17%
</pre>
=={{header|C++}}==
Partial translation of C using stl and std.
<syntaxhighlight lang="cpp">#include <random>
#include <random>
#include <vector>
#include <iostream>
#define MAX_N 20
#define TIMES 1000000
/**
* Used to generate a uniform random distribution
*/
static std::random_device rd; //Will be used to obtain a seed for the random number engine
static std::mt19937 gen(rd()); //Standard mersenne_twister_engine seeded with rd()
static std::uniform_int_distribution<> dis;
int randint(int n) {
int r, rmax = RAND_MAX / n * n;
dis=std::uniform_int_distribution<int>(0,rmax) ;
r = dis(gen);
return r / (RAND_MAX / n);
}
unsigned long long factorial(size_t n) {
//Factorial using dynamic programming to memoize the values.
static std::vector<unsigned long long>factorials{1,1,2};
for (;factorials.size() <= n;)
factorials.push_back(((unsigned long long) factorials.back())*factorials.size());
return factorials[n];
}
long double expected(size_t n) {
long double sum = 0;
for (size_t i = 1; i <= n; i++)
sum += factorial(n) / pow(n, i) / factorial(n - i);
return sum;
}
int test(int n, int times) {
int i, count = 0;
for (i = 0; i < times; i++) {
unsigned int x = 1, bits = 0;
while (!(bits & x)) {
count++;
bits |= x;
x = static_cast<unsigned int>(1 << randint(n));
}
}
return count;
}
int main() {
puts(" n\tavg\texp.\tdiff\n-------------------------------");
int n;
for (n = 1; n <= MAX_N; n++) {
int cnt = test(n, TIMES);
long double avg = (double)cnt / TIMES;
long double theory = expected(static_cast<size_t>(n));
long double diff = (avg / theory - 1) * 100;
printf("%2d %8.4f %8.4f %6.3f%%\n", n, static_cast<double>(avg), static_cast<double>(theory), static_cast<double>(diff));
}
return 0;
}
</syntaxhighlight>
{{out}}
<pre>
n avg exp. diff
-------------------------------
1 1.0000 1.0000 0.002%
2 1.4999 1.5000 -0.006%
3 1.8897 1.8889 0.042%
4 2.2177 2.2188 -0.046%
5 2.5109 2.5104 0.018%
6 2.7768 2.7747 0.077%
7 3.0187 3.0181 0.019%
8 3.2448 3.2450 -0.008%
9 3.4600 3.4583 0.049%
10 3.6619 3.6602 0.046%
11 3.8526 3.8524 0.006%
12 4.0391 4.0361 0.076%
13 4.2129 4.2123 0.012%
14 4.3858 4.3820 0.087%
15 4.5469 4.5458 0.023%
16 4.7045 4.7043 0.006%
17 4.8587 4.8579 0.016%
18 5.0071 5.0071 0.001%
19 5.1529 5.1522 0.013%
20 5.2931 5.2936 -0.010%
</pre>
=={{header|Clojure}}==
{{trans|Python}}
<syntaxhighlight lang="lisp">(ns cyclelengths
(:gen-class))
(defn factorial [n]
" n! "
(apply *' (range 1 (inc n)))) ; Use *' (vs. *) to allow arbitrary length arithmetic
(defn pow [n i]
" n^i"
(apply *' (repeat i n)))
(defn analytical [n]
" Analytical Computation "
(->>(range 1 (inc n))
(map #(/ (factorial n) (pow n %) (factorial (- n %)))) ;calc n %))
(reduce + 0)))
;; Number of random times to test each n
(def TIMES 1000000)
(defn single-test-cycle-length [n]
" Single random test of cycle length "
(loop [count 0
bits 0
x 1]
(if (zero? (bit-and x bits))
(recur (inc count) (bit-or bits x) (bit-shift-left 1 (rand-int n)))
count)))
(defn avg-cycle-length [n times]
" Average results of single tests of cycle lengths "
(/
(reduce +
(for [i (range times)]
(single-test-cycle-length n)))
times))
;; Show Results
(println "\tAvg\t\tExp\t\tDiff")
(doseq [q (range 1 21)
:let [anal (double (analytical q))
avg (double (avg-cycle-length q TIMES))
diff (Math/abs (* 100 (- 1 (/ avg anal))))]]
(println (format "%3d\t%.4f\t%.4f\t%.2f%%" q avg anal diff)))
</syntaxhighlight>
{{Output}}
<pre>
Avg Exp Diff
1 1.0000 1.0000 0.00%
2 1.4995 1.5000 0.03%
3 1.8899 1.8889 0.05%
4 2.2178 2.2188 0.04%
5 2.5118 2.5104 0.06%
6 2.7773 2.7747 0.09%
7 3.0177 3.0181 0.02%
8 3.2448 3.2450 0.01%
9 3.4587 3.4583 0.01%
10 3.6594 3.6602 0.02%
11 3.8553 3.8524 0.08%
12 4.0335 4.0361 0.06%
13 4.2113 4.2123 0.03%
14 4.3823 4.3820 0.01%
15 4.5491 4.5458 0.07%
16 4.7035 4.7043 0.02%
17 4.8580 4.8579 0.00%
18 5.0050 5.0071 0.04%
19 5.1543 5.1522 0.04%
20 5.2956 5.2936 0.04%
</pre>
=={{header|D}}==
{{trans|
<
real analytical(in int n) pure nothrow @safe /*@nogc*/ {
Line 295 ⟶ 603:
n, average, an, errorS);
}
}</
{{out}}
<pre> n average analytical (error)
Line 339 ⟶ 647:
39 7.51864 7.50959 ( 0.1204%)
40 7.60255 7.60911 ( 0.0863%)</pre>
=={{header|Delphi}}==
{{libheader| System.SysUtils}}
{{libheader| System.Math}}
{{Trans|C}}
<syntaxhighlight lang="delphi">
program Average_loop_length;
{$APPTYPE CONSOLE}
uses
System.SysUtils,
System.Math;
const
MAX_N = 20;
TIMES = 1000000;
function Factorial(const n: Double): Double;
begin
Result := 1;
if n > 1 then
Result := n * Factorial(n - 1);
end;
function Expected(const n: Integer): Double;
var
i: Integer;
begin
Result := 0;
for i := 1 to n do
Result := Result + (factorial(n) / Power(n, i) / factorial(n - i));
end;
function Test(const n, times: Integer): integer;
var
i, x, bits: Integer;
begin
Result := 0;
for i := 0 to times - 1 do
begin
x := 1;
bits := 0;
while ((bits and x) = 0) do
begin
inc(Result);
bits := bits or x;
x := 1 shl random(n);
end;
end;
end;
var
n, cnt: Integer;
avg, theory, diff: Double;
begin
Randomize;
Writeln(#10' tavg'^I'exp.'^I'diff'#10'-------------------------------');
for n := 1 to MAX_N do
begin
cnt := test(n, times);
avg := cnt / times;
theory := expected(n);
diff := (avg / theory - 1) * 100;
writeln(format('%2d %8.4f %8.4f %6.3f%%', [n, avg, theory, diff]));
end;
readln;
end.
</syntaxhighlight>
{{out}}
<pre>
tavg exp. diff
-------------------------------
1 1,0000 1,0000 0,000%
2 1,4985 1,5000 -0,101%
3 1,8896 1,8889 0,037%
4 2,2195 2,2188 0,035%
5 2,5103 2,5104 -0,003%
6 2,7746 2,7747 -0,005%
7 3,0176 3,0181 -0,017%
8 3,2458 3,2450 0,023%
9 3,4572 3,4583 -0,032%
10 3,6623 3,6602 0,057%
11 3,8494 3,8524 -0,078%
12 4,0373 4,0361 0,029%
13 4,2114 4,2123 -0,023%
14 4,3834 4,3820 0,032%
15 4,5449 4,5458 -0,020%
16 4,7030 4,7043 -0,027%
17 4,8574 4,8579 -0,009%
18 5,0063 5,0071 -0,014%
19 5,1506 5,1522 -0,030%
20 5,2960 5,2936 0,046%
</pre>
=={{header|EasyLang}}==
{{trans|Lua}}
<syntaxhighlight>
func average n reps .
for r to reps
f[] = [ ]
for i to n
f[] &= randint n
.
seen[] = [ ]
len seen[] n
x = 1
while seen[x] = 0
seen[x] = 1
x = f[x]
count += 1
.
.
return count / reps
.
func analytical n .
s = 1
t = 1
for i = n - 1 downto 1
t = t * i / n
s += t
.
return s
.
print " N average analytical (error)"
print "=== ======= ========== ======="
for n to 20
avg = average n 1e6
ana = analytical n
err = (avg - ana) / ana * 100
numfmt 0 2
write n
numfmt 4 9
print avg & ana & err & "%"
.
</syntaxhighlight>
=={{header|EchoLisp}}==
<
(lib 'math) ;; Σ aka (sigma f(n) nfrom nto)
Line 371 ⟶ 819:
(define fa (f-anal N))
(printf "%3d %10d %10d %10.2d %%" N fc fa (// (abs (- fa fc)) fc 0.01))))
</syntaxhighlight>
{{out}}
<pre>
Line 396 ⟶ 844:
20 5.30243 5.29358 0.17 %
</pre>
=={{header|Elixir}}==
{{trans|Ruby}}
{{works with|Elixir|1.1+}}
<
def factorial(0), do: 1
def factorial(n), do: Enum.reduce(1..n, 1, &(&1 * &2))
Line 426 ⟶ 873:
runs = 1_000_000
RC.task(runs)</
{{out}}
Line 452 ⟶ 899:
19 5.1533 5.1522 ( 0.02%)
20 5.2951 5.2936 ( 0.03%)
</pre>
=={{header|F_Sharp|F#}}==
{{trans|Scala}}
<p>But uses the Gamma function instead of factorials.</p>
<syntaxhighlight lang="fsharp">open System
let gamma z =
let lanczosCoefficients = [76.18009172947146;-86.50532032941677;24.01409824083091;-1.231739572450155;0.1208650973866179e-2;-0.5395239384953e-5]
let rec sumCoefficients acc i coefficients =
match coefficients with
| [] -> acc
| h::t -> sumCoefficients (acc + (h/i)) (i+1.0) t
let gamma = 5.0
let x = z - 1.0
Math.Pow(x + gamma + 0.5, x + 0.5) * Math.Exp( -(x + gamma + 0.5) ) * Math.Sqrt( 2.0 * Math.PI ) * sumCoefficients 1.000000000190015 (x + 1.0) lanczosCoefficients
let factorial n = gamma ((float n) + 1.)
let expected n =
seq {for i in 1 .. n do yield (factorial n) / System.Math.Pow((float n), (float i)) / (factorial (n - i)) }
|> Seq.sum
let r = System.Random()
let trial n =
let count = ref 0
let x = ref 1
let bits = ref 0
while (!bits &&& !x) = 0 do
count := !count + 1
bits := !bits ||| !x
x := 1 <<< r.Next(n)
!count
let tested n times = (float (Seq.sum (seq { for i in 1 .. times do yield (trial n) }))) / (float times)
let results = seq {
for n in 1 .. 20 do
let avg = tested n 1000000
let theory = expected n
yield n, avg, theory
}
[<EntryPoint>]
let main argv =
printfn " N average analytical (error)"
printfn "------------------------------------"
results
|> Seq.iter (fun (n, avg, theory) ->
printfn "%2i %2.6f %2.6f %+2.3f%%" n avg theory ((avg / theory - 1.) * 100.))
0
</syntaxhighlight>
{{out}}
<pre> N average analytical (error)
------------------------------------
1 1.000000 1.000000 +0.000%
2 1.498934 1.500000 -0.071%
3 1.889318 1.888889 +0.023%
4 2.219397 2.218750 +0.029%
5 2.510618 2.510400 +0.009%
6 2.771914 2.774691 -0.100%
7 3.014726 3.018139 -0.113%
8 3.245022 3.245018 +0.000%
9 3.457096 3.458316 -0.035%
10 3.660337 3.660216 +0.003%
11 3.849770 3.852372 -0.068%
12 4.038977 4.036074 +0.072%
13 4.213248 4.212348 +0.021%
14 4.380451 4.382029 -0.036%
15 4.541868 4.545807 -0.087%
16 4.704117 4.704258 -0.003%
17 4.858934 4.857871 +0.022%
18 5.004236 5.007063 -0.056%
19 5.154166 5.152196 +0.038%
20 5.298119 5.293585 +0.086%</pre>
=={{header|Factor}}==
The <code>loop-length</code> word is more or less a translation of the inner loop of C's <code>test</code> function.
{{works with|Factor|0.99 2020-01-23}}
<syntaxhighlight lang="factor">USING: formatting fry io kernel locals math math.factorials
math.functions math.ranges random sequences ;
: (analytical) ( m n -- x )
[ drop factorial ] [ ^ /f ] [ - factorial / ] 2tri ;
: analytical ( n -- x )
dup [1,b] [ (analytical) ] with map-sum ;
: loop-length ( n -- x )
[ 0 0 1 [ 2dup bitand zero? ] ] dip
'[ [ 1 + ] 2dip bitor 1 _ random shift ] while 2drop ;
:: average-loop-length ( n #tests -- x )
0 #tests [ n loop-length + ] times #tests / ;
: stats ( n -- avg exp )
[ 1,000,000 average-loop-length ] [ analytical ] bi ;
: .line ( n -- )
dup stats 2dup / 1 - 100 *
"%2d %8.4f %8.4f %6.3f%%\n" printf ;
" n\tavg\texp.\tdiff\n-------------------------------" print
20 [1,b] [ .line ] each</syntaxhighlight>
{{out}}
<pre>
n avg exp. diff
-------------------------------
1 1.0000 1.0000 0.000%
2 1.4993 1.5000 -0.044%
3 1.8877 1.8889 -0.064%
4 2.2193 2.2188 0.023%
5 2.5099 2.5104 -0.021%
6 2.7728 2.7747 -0.068%
7 3.0165 3.0181 -0.056%
8 3.2442 3.2450 -0.026%
9 3.4574 3.4583 -0.027%
10 3.6622 3.6602 0.054%
11 3.8537 3.8524 0.033%
12 4.0365 4.0361 0.010%
13 4.2094 4.2123 -0.070%
14 4.3819 4.3820 -0.004%
15 4.5469 4.5458 0.023%
16 4.7028 4.7043 -0.031%
17 4.8571 4.8579 -0.016%
18 5.0049 5.0071 -0.043%
19 5.1519 5.1522 -0.005%
20 5.2927 5.2936 -0.017%
</pre>
=={{header|FreeBASIC}}==
<syntaxhighlight lang="freebasic">Const max_N = 20, max_ciclos = 1000000
Function Factorial(Byval N As Integer) As Double
Dim As Double d: d = 1
If N = 0 Then Factorial = 1: Exit Function
While (N > 1)
d *= N
N -= 1
Wend
Factorial = d
End Function
Function Analytical(N As Integer) As Double
Dim As Double i, sum = 0
For i = 1 To N
sum += Factorial(N) / N^i / Factorial(N-i)
Next i
Return sum
End Function
Function Average(N As Integer, ciclos As Double) As Double
Dim As Integer i, x, bits, sum = 0
For i = 0 To ciclos - 1
x = 1 : bits = 0
While (bits And x) = 0
sum += 1
bits Or= x
x = 1 Shl (Rnd * (N - 1))
Wend
Next i
Return sum / ciclos
End Function
Randomize Timer
Print " N promedio analitico (error)"
Print "--- ---------- ----------- ----------"
For N As Integer = 1 To max_N
Dim As Double avg = Average(N, max_ciclos)
Dim As Double ana = Analytical(N)
Dim As Double diff = abs(avg-ana) / ana * 100
Print Using " ## #####.###0 #####.###0 ###.#0%"; N; avg; ana; diff
Next N
Sleep
</syntaxhighlight>
=={{header|FutureBasic}}==
<syntaxhighlight lang="futurebasic">
_nmax = 20
_times = 1000000
local fn Average( n as long, times as long ) as double
long i, x
double b, c = 0
for i = 0 to times
x = 1 : b = 0
while ( b and x ) == 0
c++
b = b || x
x = 1 << ( rnd(n) - 1 )
wend
next
end fn = c / times
local fn Analyltic( n as long ) as double
double nn = (double)n
double term = 1.0
double sum = 1.0
long i
for i = nn - 1 to i >= 1 step -1
term = term * i / nn
sum = sum + term
next
end fn = sum
local fn DoIt
long n
double average, theory, difference
window 1
printf @"\nSamples tested: %ld\n", _times
print " N Average Analytical (error)"
print "=== ========= ============ ========="
for n = 1 to _nmax
average = fn Average( n, _times )
theory = fn Analyltic( n )
difference = ( average / theory - 1) * 100
printf @"%3d %9.4f %9.4f %10.4f%%", n, average, theory, difference
next
end fn
randomize
fn DoIt
HandleEvents
</syntaxhighlight>
{{output}}
<pre>
Number of tests performed: 1000000
N Average Analytical (error)
=== ========= ============ =========
1 1.0000 1.0000 0.0001%
2 1.4999 1.5000 -0.0070%
3 1.8877 1.8889 -0.0630%
4 2.2187 2.2188 -0.0011%
5 2.5102 2.5104 -0.0065%
6 2.7735 2.7747 -0.0438%
7 3.0173 3.0181 -0.0273%
8 3.2478 3.2450 0.0854%
9 3.4569 3.4583 -0.0424%
10 3.6604 3.6602 0.0060%
11 3.8506 3.8524 -0.0449%
12 4.0361 4.0361 0.0003%
13 4.2148 4.2123 0.0582%
14 4.3845 4.3820 0.0559%
15 4.5454 4.5458 -0.0079%
16 4.7056 4.7043 0.0288%
17 4.8558 4.8579 -0.0428%
18 5.0143 5.0071 0.1450%
19 5.1509 5.1522 -0.0254%
20 5.2985 5.2936 0.0929%
</pre>
=={{header|Go}}==
<
import (
Line 498 ⟶ 1,204:
}
return sum
}</
{{out}}
<pre>
Line 526 ⟶ 1,232:
=={{header|Haskell}}==
<
import qualified Data.Set as S
import Text.Printf
Line 579 ⟶ 1,285:
range = [1..20] :: [Integer]
_ <- test samples range $ mkStdGen 0
return ()</
<pre> N average analytical (error)
=== ========= ============ =========
Line 609 ⟶ 1,315:
We can implement f as {&LIST where LIST is an arbitrary list of N numbers, each picked independently from the range 0..(N-1). We can incrementally build the described sequence using (, f@{:) - here we extend the sequence by applying f to the last element of the sequence. Since we are only concerned with the sequence up to the point of the first repeat, we can select the unique values giving us (~.@, f@{:). This routine stops changing when we reach the desired length, so we can repeatedly apply it forever. For example:
<
0
(~.@, {&0 0@{:)^:_] 1
1 0</
Once we have the sequence, we can count how many elements are in it.
<
2</
Meanwhile, we can also generate all possible values of 1..N by counting out N^N values and breaking out the result as a base N list of digits.
<
0 0
0 1
1 0
1 1</
All that's left is to count the lengths of all possible sequences for all possible distinct instances of f and average the results:
<
1
(+/ % #)@,@((#.inv i.@^~) ([: # (] ~.@, {:@] { [)^:_)"1 0/ i.)2
Line 634 ⟶ 1,340:
2.5104
(+/ % #)@,@((#.inv i.@^~) ([: # (] ~.@, {:@] { [)^:_)"1 0/ i.)6
2.77469</
Meanwhile the analytic solution (derived by reading the Ada implementation) looks like this:
<
ana"0]1 2 3 4 5 6
1 1.5 1.88889 2.21875 2.5104 2.77469</
To get our simulation, we can take the exact approach and replace the part that generates all possible values for f with a random mechanism. Since the task does not specify how long to run the simulation, and to make this change easy, we'll use N*1e4 tests.
<
sim"0]1 2 3 4 5 6
1 1.5034 1.8825 2.22447 2.51298 2.76898</
The simulation approach is noticeably slower than the analytic approach, while being less accurate.
Finally, we can generate our desired results:
<
+--+-------+--------+-----------+
|N |average|analytic|error |
Line 670 ⟶ 1,376:
|19| 5.1522|5.14785 |0.000843052|
|20|5.29358|5.28587 | 0.00145829|
+--+-------+--------+-----------+</
Here, error is the difference between the two values divided by the analytic value.
Line 677 ⟶ 1,383:
This uses a 0-based index (0, 1, ..., n-1) as opposed to the 1-based index (1, 2, ..., n) specified in the question, because it fits better with the native structure of Java.
<
import java.util.Random;
import java.util.Set;
public class AverageLoopLength {
private static final int N = 100000;
//analytical(n) = sum_(i=1)^n (n!/(n-i)!/n**i)
private static double analytical(int n) {
double[] powers = new double[n + 1];
powers[0] = 1.0;
factorial[0] = 1.0;
for (int i = 1; i <= n; i++) {
factorial[i] = factorial[i - 1] * i;
powers[i] = powers[i - 1] * n;
}
double sum = 0;
//memoized factorial and powers
for (int i = 1; i <= n; i++) {
sum += factorial[n] / factorial[n - i] / powers[i];
}
return sum;
}
private static double average(int n) {
Random rnd = new Random();
double sum = 0.0;
for (int a = 0; a < N; a++) {
for (int i = 0; i < n; i++) {
random[i] = rnd.nextInt(n);
}
Set<Integer> seen = new HashSet<>(n);
int length = 0;
while (seen.add(current)) {
length++;
current = random[current];
}
sum += length;
}
return sum / N;
}
public static void main(String[] args) {
System.out.println(" N average analytical (error)");
System.out.println("=== ========= ============ =========");
for (int i = 1; i <= 20; i++) {
double avg = average(i);
double ana = analytical(i);
System.out.println(String.format("%3d %9.4f %12.4f (%6.2f%%)", i, avg, ana, ((ana - avg) / ana * 100)));
}
}
}</syntaxhighlight>
=={{header|Julia}}==
{{trans|Python}}
<syntaxhighlight lang="julia">using Printf
analytical(n::Integer) = sum(factorial(n) / big(n) ^ i / factorial(n - i) for i = 1:n)
function test(n::Integer, times::Integer = 1000000)
c = 0
for i = range(0, times)
x, bits = 1, 0
while (bits & x) == 0
c += 1
bits |= x
x = 1 << rand(0:(n - 1))
end
end
return c / times
end
function main(n::Integer)
println(" n\tavg\texp.\tdiff\n-------------------------------")
for n in 1:n
avg = test(n)
theory = analytical(n)
diff = (avg / theory - 1) * 100
@printf(STDOUT, "%2d %8.4f %8.4f %6.3f%%\n", n, avg, theory, diff)
end
end
main(20)
</syntaxhighlight>
{{out}}
<pre>
n avg exp. diff
-------------------------------
1 1.0000 1.0000 0.000%
2 1.4998 1.5000 -0.015%
3 1.8895 1.8889 0.034%
4 2.2171 2.2188 -0.075%
5 2.5082 2.5104 -0.088%
6 2.7729 2.7747 -0.063%
7 3.0171 3.0181 -0.033%
8 3.2439 3.2450 -0.034%
9 3.4578 3.4583 -0.016%
10 3.6616 3.6602 0.038%
11 3.8525 3.8524 0.004%
12 4.0353 4.0361 -0.020%
13 4.2126 4.2123 0.006%
14 4.3835 4.3820 0.034%
15 4.5428 4.5458 -0.067%
16 4.7027 4.7043 -0.033%
17 4.8560 4.8579 -0.039%
18 5.0054 5.0071 -0.033%
19 5.1492 5.1522 -0.058%
20 5.2896 5.2936 -0.076%
</pre>
=={{header|Kotlin}}==
{{trans|Go}}
<syntaxhighlight lang="scala">const val NMAX = 20
const val TESTS = 1000000
val rand = java.util.Random()
fun avg(n: Int): Double {
var sum = 0
for (t in 0 until TESTS) {
val v = BooleanArray(NMAX)
var x = 0
while (!v[x]) {
v[x] = true
sum++
x = rand.nextInt(n)
}
}
return sum.toDouble() / TESTS
}
fun ana(n: Int): Double {
val nn = n.toDouble()
var term = 1.0
var sum = 1.0
for (i in n - 1 downTo 1) {
term *= i / nn
sum += term
}
return sum
}
fun main(args: Array<String>) {
println(" N average analytical (error)")
println("=== ========= ============ =========")
for (n in 1..NMAX) {
val a = avg(n)
val b = ana(n)
println(String.format("%3d %6.4f %10.4f (%4.2f%%)", n, a, b, Math.abs(a - b) / b * 100.0))
}
}</syntaxhighlight>
Sample output:
{{out}}
<pre>
N average analytical (error)
=== ========= ============ =========
1 1.0000 1.0000 (0.00%)
2 1.5004 1.5000 (0.03%)
3 1.8890 1.8889 (0.00%)
4 2.2179 2.2188 (0.04%)
5 2.5108 2.5104 (0.02%)
6 2.7738 2.7747 (0.03%)
7 3.0178 3.0181 (0.01%)
8 3.2482 3.2450 (0.10%)
9 3.4572 3.4583 (0.03%)
10 3.6608 3.6602 (0.02%)
11 3.8545 3.8524 (0.06%)
12 4.0378 4.0361 (0.04%)
13 4.2131 4.2123 (0.02%)
14 4.3795 4.3820 (0.06%)
15 4.5481 4.5458 (0.05%)
16 4.7044 4.7043 (0.00%)
17 4.8610 4.8579 (0.06%)
18 5.0027 5.0071 (0.09%)
19 5.1498 5.1522 (0.05%)
20 5.2941 5.2936 (0.01%)
</pre>
=={{header|Liberty BASIC}}==
{{trans|BBC BASIC}}
<syntaxhighlight lang="lb">
MAXN = 20
TIMES = 10000'00
Line 769 ⟶ 1,604:
IF n=1 OR n=0 THEN FNfactorial=1 ELSE FNfactorial= n * FNfactorial(n-1)
end function
</syntaxhighlight>
{{out}}
Line 794 ⟶ 1,629:
20 5.2792 5.29358459 -0.2717%
</pre>
=={{header|Lua}}==
<syntaxhighlight lang="lua">function average(n, reps)
local count = 0
for r = 1, reps do
local f = {}
for i = 1, n do f[i] = math.random(n) end
local seen, x = {}, 1
while not seen[x] do
seen[x], x, count = true, f[x], count+1
end
end
return count / reps
end
function analytical(n)
local s, t = 1, 1
for i = n-1, 1, -1 do t=t*i/n s=s+t end
return s
end
print(" N average analytical (error)")
print("=== ========= ============ =========")
for n = 1, 20 do
local avg, ana = average(n, 1e6), analytical(n)
local err = (avg-ana) / ana * 100
print(string.format("%3d %9.4f %12.4f (%6.3f%%)", n, avg, ana, err))
end</syntaxhighlight>
{{out}}
<pre> N average analytical (error)
=== ========= ============ =========
1 1.0000 1.0000 ( 0.000%)
2 1.5002 1.5000 ( 0.014%)
3 1.8896 1.8889 ( 0.037%)
4 2.2176 2.2188 (-0.054%)
5 2.5094 2.5104 (-0.038%)
6 2.7732 2.7747 (-0.054%)
7 3.0186 3.0181 ( 0.016%)
8 3.2440 3.2450 (-0.031%)
9 3.4554 3.4583 (-0.085%)
10 3.6625 3.6602 ( 0.063%)
11 3.8534 3.8524 ( 0.026%)
12 4.0354 4.0361 (-0.016%)
13 4.2111 4.2123 (-0.031%)
14 4.3839 4.3820 ( 0.043%)
15 4.5453 4.5458 (-0.012%)
16 4.7054 4.7043 ( 0.024%)
17 4.8596 4.8579 ( 0.035%)
18 5.0099 5.0071 ( 0.056%)
19 5.1553 5.1522 ( 0.060%)
20 5.2901 5.2936 (-0.066%)</pre>
=={{header|Mathematica}} / {{header|Wolfram Language}}==
<
Table[{n, #[[1]], #[[2]],
Row[{Round[10000 Abs[#[[1]] - #[[2]]]/#[[2]]]/100., "%"}]} &@
Line 804 ⟶ 1,690:
RandomInteger[{1, n}, n]] &] &, 10000]],
Sum[n! n^(n - k - 1)/(n - k)!, {k, n}]/n^(n - 1)}, 5], {n, 1,
20}], {"N", "average", "analytical", "error"}]</
{{Out}}
<pre>N average analytical error
Line 830 ⟶ 1,716:
=={{header|Nim}}==
{{trans|C}}
<
randomize()
const
maxN = 20
times = 1_000_000
proc factorial(n: int): float =
result = 1
for i in 1 .. n:
result *= i.float
proc expected(n: int): float =
for i in 1 .. n:
result += factorial(n) / pow(n.float, i.float) / factorial(n - i)
proc test(n, times: int): int =
for i in 1 .. times:
var
Line 854 ⟶ 1,740:
inc result
bits = bits or x
x = 1 shl
echo " n\tavg\texp.\tdiff"
echo "-------------------------------"
Line 863 ⟶ 1,749:
let theory = expected(n)
let diff = (avg / theory - 1) * 100
{{out}}
<pre> n avg exp. diff
Line 887 ⟶ 1,773:
19 5.1554 5.1522 0.061%
20 5.2915 5.2936 -0.040%</pre>
=={{header|Oberon-2}}==
<syntaxhighlight lang="oberon2">
MODULE AvgLoopLen;
(* Oxford Oberon-2 *)
IMPORT Random, Out;
PROCEDURE Fac(n: INTEGER; f: REAL): REAL;
BEGIN
IF n = 0 THEN
RETURN f
ELSE
RETURN Fac(n - 1,n*f)
END
END Fac;
PROCEDURE Power(n,i: INTEGER): REAL;
VAR
p: REAL;
BEGIN
p := 1.0;
WHILE i > 0 DO p := p * n; DEC(i) END;
RETURN p
END Power;
PROCEDURE Abs(x: REAL): REAL;
BEGIN
IF x < 0 THEN RETURN -x ELSE RETURN x END
END Abs;
PROCEDURE Analytical(n: INTEGER): REAL;
VAR
i: INTEGER;
res: REAL;
BEGIN
res := 0.0;
FOR i := 1 TO n DO
res := res + (Fac(n,1.0) / Power(n,i) / Fac(n - i,1.0));
END;
RETURN res
END Analytical;
PROCEDURE Averages(n: INTEGER): REAL;
CONST
times = 100000;
VAR
rnds: SET;
r,count,i: INTEGER;
BEGIN
count := 0; i := 0;
WHILE i < times DO
rnds := {};
LOOP
r := Random.Roll(n);
IF r IN rnds THEN EXIT ELSE INCL(rnds,r); INC(count) END
END;
INC(i)
END;
RETURN count / times
END Averages;
VAR
i: INTEGER;
av,an,df: REAL;
BEGIN
Random.Randomize;
Out.String(" Averages Analytical Diff% ");Out.Ln;
FOR i := 1 TO 20 DO
Out.Int(i,3); Out.String(": ");
av := Averages(i);an := Analytical(i);df := Abs(av - an) / an * 100.0;
Out.Fixed(av,10,4);Out.Fixed(an,11,4);Out.Fixed(df,10,4);Out.Ln
END
END AvgLoopLen.
</syntaxhighlight>
{{Out}}
<pre>
Averages Analytical Diff%
1: 1.0000 1.0000 0.0000
2: 1.5015 1.5000 0.0993
3: 1.8868 1.8889 0.1085
4: 2.2187 2.2188 0.0005
5: 2.5119 2.5104 0.0578
6: 2.7785 2.7747 0.1366
7: 3.0184 3.0181 0.0090
8: 3.2435 3.2450 0.0471
9: 3.4585 3.4583 0.0056
10: 3.6549 3.6602 0.1463
11: 3.8559 3.8524 0.0918
12: 4.0452 4.0361 0.2264
13: 4.2097 4.2123 0.0628
14: 4.3740 4.3820 0.1830
15: 4.5583 4.5458 0.2739
16: 4.7001 4.7043 0.0882
17: 4.8654 4.8579 0.1556
18: 5.0157 5.0071 0.1731
19: 5.1515 5.1522 0.0135
20: 5.2930 5.2936 0.0105
</pre>
=={{header|PARI/GP}}==
{{trans|C}}
<
test(n, times)={
my(ct);
Line 903 ⟶ 1,888:
my(cnt=test(n, TIMES),avg=cnt/TIMES,ex=expected(n),diff=(avg/ex-1)*100.);
print(n"\t"avg*1."\t"ex*1."\t"diff);
)}</
{{out}}
<pre>1 1.0000 1.0000 0.E-7
Line 926 ⟶ 1,911:
20 5.2970 5.2936 0.065143</pre>
=={{header|Perl
<syntaxhighlight lang="perl">use List::Util qw(sum reduce);
sub find_loop {
my($n) = @_;
my($r,@seen);
while () { $seen[$r] = $seen[($r = int(1+rand $n))] ? return sum @seen : 1 }
}
print " N empiric theoric (error)\n";
print "=== ========= ============ =========\n";
my $MAX = 20;
my $TRIALS = 1000;
for my $n (1 .. $MAX) {
my $empiric = ( sum map { find_loop($n) } 1..$TRIALS ) / $TRIALS;
my $theoric = sum map { (reduce { $a*$b } $_**2, ($n-$_+1)..$n ) / $n ** ($_+1) } 1..$n;
printf "%3d %9.4f %12.4f (%5.2f%%)\n",
$n, $empiric, $theoric, 100 * ($empiric - $theoric) / $theoric;
}</syntaxhighlight>
{{out}}
<pre> N empiric
=== ========= ============ =========
1 1.0000 1.0000 ( 0.00%)
2 1.4950
3 1.9190
4 2.2400
5 2.5120
6 2.7500 2.7747 (-0.89%)
7 3.0360 3.0181 ( 0.59%)
8 3.2600 3.2450 ( 0.46%)
9 3.4440 3.4583 (-0.41%)
10 3.6670 3.6602 ( 0.19%)
11
12
13 4.2160 4.2123 ( 0.09%)
14 4.4420
15 4.5600
16 4.7940 4.7043 ( 1.91%)
17 4.7830 4.8579 (-1.54%)
18 4.9140 5.0071 (-1.86%)
19 5.2490 5.1522 ( 1.88%)
20 5.2930 5.2936 (-0.01%)</pre>
=={{header|Phix}}==
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">constant</span> <span style="color: #000000;">MAX</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">20</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">ITER</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1000000</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">expected</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">atom</span> <span style="color: #7060A8;">sum</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span>
<span style="color: #7060A8;">sum</span> <span style="color: #0000FF;">+=</span> <span style="color: #7060A8;">factorial</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">i</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #7060A8;">factorial</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">i</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #7060A8;">sum</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">test</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">count</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">bits</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">ITER</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span>
<span style="color: #000000;">bits</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #008080;">while</span> <span style="color: #008080;">not</span> <span style="color: #7060A8;">and_bits</span><span style="color: #0000FF;">(</span><span style="color: #000000;">bits</span><span style="color: #0000FF;">,</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">count</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
<span style="color: #000000;">bits</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">or_bits</span><span style="color: #0000FF;">(</span><span style="color: #000000;">bits</span><span style="color: #0000FF;">,</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">rand</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">count</span><span style="color: #0000FF;">/</span><span style="color: #000000;">ITER</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">av</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">ex</span>
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" n avg. exp. (error%)\n"</span><span style="color: #0000FF;">);</span>
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"== ====== ====== ========\n"</span><span style="color: #0000FF;">);</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">MAX</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">av</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">ex</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">expected</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%2d %8.4f %8.4f (%5.3f%%)\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">av</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ex</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">-</span><span style="color: #000000;">av</span><span style="color: #0000FF;">/</span><span style="color: #000000;">ex</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">100</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
Line 1,029 ⟶ 2,018:
20 5.2955 5.2936 (0.037%)
</pre>
=={{header|Phixmonti}}==
{{trans|Phix}}
<syntaxhighlight lang="phixmonti">include ..\Utilitys.pmt
20 var MAX
100000 var ITER
def factorial 1 swap for * endfor enddef
def expected /# n -- n #/
>ps
0
tps for var i
tps factorial tps i power / tps i - factorial / +
endfor
ps> drop
enddef
def condition over over bitand not enddef
def test /# n -- n #/
0 >ps
ITER for var i
0 1
condition while
ps> 1 + >ps
bitor
over rand * 1 + int 1 - 2 swap power
condition endwhile
drop drop
endfor
drop ps> ITER /
enddef
def printAll len for get print 9 tochar print endfor enddef
( "n" "avg." "exp." "(error%)" ) printAll drop nl
( "==" "======" "======" "========" ) printAll drop nl
MAX for var n
n test
n expected
n rot rot over over / 1 swap - abs 100 * 4 tolist printAll drop nl
endfor</syntaxhighlight>
{{out}}
<pre>n avg. exp. (error%)
== ====== ====== ========
1 1 1 0
2 1.50119 1.5 0.0793333
3 1.89076 1.88889 0.0990588
4 2.22164 2.21875 0.130254
5 2.50989 2.5104 0.0203155
6 2.78108 2.77469 0.230247
7 3.02431 3.01814 0.204474
8 3.24594 3.24502 0.0284126
9 3.46167 3.45832 0.096991
10 3.66691 3.66022 0.182894
11 3.84558 3.85237 0.176308
12 4.03174 4.03607 0.107374
13 4.21113 4.21235 0.0289129
14 4.37294 4.38203 0.207425
15 4.54199 4.54581 0.0839738
16 4.69651 4.70426 0.164707
17 4.8463 4.85787 0.238187
18 5.01786 5.00706 0.215633
19 5.15783 5.1522 0.109348
20 5.29575 5.29358 0.0409064
=== Press any key to exit ===</pre>
=={{header|PicoLisp}}==
{{trans|Python}}
<
(seed (in "/dev/urandom" (rd 8)))
Line 1,073 ⟶ 2,132:
2 ) ) ) ) )
(bye)</
=={{header|PowerShell}}==
{{works with|PowerShell|2}}
<syntaxhighlight lang="powershell">
function Get-AnalyticalLoopAverage ( [int]$N )
{
# Expected loop average = sum from i = 1 to N of N! / (N-i)! / N^(N-i+1)
# Equivalently, Expected loop average = sum from i = 1 to N of F(i)
# where F(N) = 1, and F(i) = F(i+1)*i/N
$LoopAverage = $Fi = 1
If ( $N -eq 1 ) { return $LoopAverage }
ForEach ( $i in ($N-1)..1 )
{
$Fi *= $i / $N
$LoopAverage += $Fi
}
return $LoopAverage
}
function Get-ExperimentalLoopAverage ( [int]$N, [int]$Tests = 100000 )
{
If ( $N -eq 1 ) { return 1 }
# Using 0 through N-1 instead of 1 through N for speed and simplicity
$NMO = $N - 1
# Create array to hold mapping function
$F = New-Object int[] ( $N )
$Count = 0
$Random = New-Object System.Random
ForEach ( $Test in 1..$Tests )
{
# Map each number to a random number
ForEach ( $i in 0..$NMO )
{
$F[$i] = $Random.Next( $N )
}
# For each number...
ForEach ( $i in 0..$NMO )
{
# Add the number to the list
$List = @()
$Count++
$List += $X = $i
# If loop does not yet exist in list...
While ( $F[$X] -notin $List )
{
# Go to the next mapped number and add it to the list
$Count++
$List += $X = $F[$X]
}
}
}
$LoopAvereage = $Count / $N / $Tests
return $LoopAvereage
}
</syntaxhighlight>
Note: The use of the [pscustomobject] type accelerator to simplify making the test result table look pretty requires PowerShell 3.0.
<syntaxhighlight lang="powershell">
# Display results for N = 1 through 20
ForEach ( $N in 1..20 )
{
$AnalyticalAverage = Get-AnalyticalLoopAverage $N
$ExperimentalAverage = Get-ExperimentalLoopAverage $N
[pscustomobject] @{
N = $N.ToString().PadLeft( 2, ' ' )
Analytical = $AnalyticalAverage.ToString( '0.00000000' )
Experimental = $ExperimentalAverage.ToString( '0.00000000' )
'Error (%)' = ( [math]::Abs( $AnalyticalAverage - $ExperimentalAverage ) / $AnalyticalAverage * 100 ).ToString( '0.00000000' )
}
}
</syntaxhighlight>
{{out}}
<pre>
N Analytical Experimental Error (%)
- ---------- ------------ ---------
1 1.00000000 1.00000000 0.00000000
2 1.50000000 1.49985500 0.00966667
3 1.88888889 1.88713000 0.09311765
4 2.21875000 2.22103500 0.10298592
5 2.51040000 2.51069200 0.01163161
6 2.77469136 2.77264833 0.07363070
7 3.01813870 3.01547143 0.08837474
8 3.24501801 3.25003875 0.15472163
9 3.45831574 3.45067667 0.22089013
10 3.66021568 3.65659000 0.09905646
11 3.85237205 3.85669273 0.11215626
12 4.03607368 4.03813500 0.05107253
13 4.21234791 4.20946231 0.06850349
14 4.38202942 4.38458786 0.05838465
15 4.54580729 4.54466400 0.02515032
16 4.70425825 4.70146375 0.05940356
17 4.85787082 4.86807647 0.21008483
18 5.00706310 5.01939278 0.24624572
19 5.15219620 5.15179263 0.00783296
20 5.29358459 5.29214950 0.02710991
</pre>
=={{header|Python}}==
{{trans|C}}
<
from math import factorial
from random import randrange
Line 1,103 ⟶ 2,266:
theory = analytical(n)
diff = (avg / theory - 1) * 100
print("%2d %8.4f %8.4f %6.3f%%" % (n, avg, theory, diff))</
{{out}}
<pre> n avg exp. diff
Line 1,127 ⟶ 2,290:
19 5.1534 5.1522 0.024%
20 5.2927 5.2936 -0.017%</pre>
=={{header|Quackery}}==
<syntaxhighlight lang="Quackery"> [ $ "bigrat.qky" loadfile ] now!
[ tuck space swap of
join
swap split drop echo$ ] is lecho$ ( $ n --> )
[ 1 swap times [ i 1+ * ] ] is ! ( n --> n )
[ 0 n->v rot
dup temp put
times
[ temp share ! n->v
temp share i 1+ - ! n->v
v/
temp share i 1+ ** n->v
v/ v+ ]
temp release ] is expected ( n --> n/d )
[ -1 temp put
0
[ 1 temp tally
over random bit
2dup & not while
| again ]
2drop drop
temp take ] is trial ( n --> n )
[ tuck 0 swap
times
[ over trial + ]
nip swap reduce ] is trials ( n n --> n/d )
[ say " n average expected difference"
cr
say "-- ------- -------- ----------"
cr
20 times
[ i^ 1+ dup 10 < if sp echo
2 times sp
i^ 1+ 1000000 trials
2dup 7 point$ 10 lecho$
i^ 1+ expected
2dup 7 point$ 11 lecho$
v/ 1 n->v v- 100 1 v* vabs
7 point$ echo$ say "%" cr ] ] is task ( --> )</syntaxhighlight>
{{out}}
<pre> n average expected difference
-- ------- -------- ----------
1 1 1 0%
2 1.499195 1.5 0.0536667%
3 1.88936 1.8888889 0.0249412%
4 2.220728 2.21875 0.0891493%
5 2.508183 2.5104 0.0883126%
6 2.773072 2.7746914 0.0583617%
7 3.019331 3.0181387 0.0395045%
8 3.243534 3.245018 0.0457318%
9 3.45625 3.4583157 0.0597327%
10 3.658848 3.6602157 0.0373661%
11 3.850874 3.8523721 0.0388865%
12 4.032375 4.0360737 0.0916404%
13 4.212238 4.2123479 0.0026093%
14 4.383076 4.3820294 0.0238834%
15 4.544029 4.5458073 0.0391192%
16 4.706797 4.7042582 0.0539671%
17 4.856011 4.8578708 0.0382847%
18 5.004107 5.0070631 0.0590386%
19 5.152561 5.1521962 0.0070805%
20 5.288056 5.2935846 0.1044394%
</pre>
=={{header|R}}==
<syntaxhighlight lang="r">
expected <- function(size) {
result <- 0
for (i in 1:size) {
result <- result + factorial(size) / size^i / factorial(size -i)
}
result
}
knuth <- function(size) {
v <- sample(1:size, size, replace = TRUE)
visit <- vector('logical',size)
place <- 1
visit[[1]] <- TRUE
steps <- 0
repeat {
place <- v[[place]]
steps <- steps + 1
if (visit[[place]]) break
visit[[place]] <- TRUE
}
steps
}
cat(" N average analytical (error)\n")
cat("=== ========= ============ ==========\n")
for (num in 1:20) {
average <- mean(replicate(1e6, knuth(num)))
analytical <- expected(num)
error <- abs(average/analytical-1)*100
cat(sprintf("%3d%11.4f%14.4f ( %4.4f%%)\n", num, round(average,4), round(analytical,4), round(error,2)))
}
</syntaxhighlight>
{{out}}
<pre>
N average analytical (error)
=== ========= ============ ==========
1 1.0000 1.0000 ( 0.0000%)
2 1.5002 1.5000 ( 0.0100%)
3 1.8892 1.8889 ( 0.0100%)
4 2.2190 2.2188 ( 0.0100%)
5 2.5108 2.5104 ( 0.0200%)
6 2.7751 2.7747 ( 0.0200%)
7 3.0177 3.0181 ( 0.0100%)
8 3.2472 3.2450 ( 0.0700%)
9 3.4582 3.4583 ( 0.0000%)
10 3.6600 3.6602 ( 0.0100%)
11 3.8530 3.8524 ( 0.0200%)
12 4.0366 4.0361 ( 0.0100%)
13 4.2085 4.2123 ( 0.0900%)
14 4.3814 4.3820 ( 0.0100%)
15 4.5446 4.5458 ( 0.0300%)
16 4.7063 4.7043 ( 0.0400%)
17 4.8555 4.8579 ( 0.0500%)
18 5.0099 5.0071 ( 0.0600%)
19 5.1567 5.1522 ( 0.0900%)
20 5.2940 5.2936 ( 0.0100%)
</pre>
=={{header|Racket}}==
<
#lang racket
(require (only-in math factorial))
Line 1,157 ⟶ 2,458:
(test-table 20 10000)
</syntaxhighlight>
{{out}}
Line 1,184 ⟶ 2,485:
20 5.3316 5.2936 0.7181%
</pre>
=={{header|Raku}}==
(formerly Perl 6)
<syntaxhighlight lang="raku" line>constant MAX_N = 20;
constant TRIALS = 100;
for 1 .. MAX_N -> $N {
my $empiric = TRIALS R/ [+] find-loop(random-mapping $N).elems xx TRIALS;
my $theoric = [+]
map -> $k { $N ** ($k + 1) R/ [×] flat $k**2, $N - $k + 1 .. $N }, 1 .. $N;
FIRST say " N empiric theoric (error)";
FIRST say "=== ========= ============ =========";
printf "%3d %9.4f %12.4f (%4.2f%%)\n",
$N, $empiric, $theoric, 100 × abs($theoric - $empiric) / $theoric;
}
sub random-mapping { hash .list Z=> .roll($_) given ^$^size }
sub find-loop { 0, | %^mapping{*} ...^ { (%){$_}++ } }</syntaxhighlight>
{{out|Example}}
<pre> N empiric theoric (error)
=== ========= ============ =========
1 1.0000 1.0000 (0.00%)
2 1.5600 1.5000 (4.00%)
3 1.7800 1.8889 (5.76%)
4 2.1800 2.2188 (1.75%)
5 2.6200 2.5104 (4.37%)
6 2.8300 2.7747 (1.99%)
7 3.1200 3.0181 (3.37%)
8 3.1400 3.2450 (3.24%)
9 3.4500 3.4583 (0.24%)
10 3.6700 3.6602 (0.27%)
11 3.8300 3.8524 (0.58%)
12 4.3600 4.0361 (8.03%)
13 3.9000 4.2123 (7.42%)
14 4.4900 4.3820 (2.46%)
15 4.9500 4.5458 (8.89%)
16 4.9800 4.7043 (5.86%)
17 4.9100 4.8579 (1.07%)
18 4.9700 5.0071 (0.74%)
19 5.1000 5.1522 (1.01%)
20 5.2300 5.2936 (1.20%)</pre>
=={{header|REXX}}==
Line 1,190 ⟶ 2,534:
Also note that the <big>'''!'''</big> (factorial function) uses memoization for optimization.
<
parse arg runs tests seed . /*obtain optional arguments from the CL*/
if runs =='' | runs =="," then runs = 40 /*Not specified? Then use the default.*/
if tests =='' | tests =="," then tests= 1000000 /* " " " " " " */
if datatype(seed, 'W')
!.=0;
numeric digits 100000 /*be able to calculate 25k! if need be.*/
numeric digits max(9, length( !(runs) ) ) /*set the NUMERIC DIGITS for !(runs). */
say right( runs, 24)
say right( tests, 24)
say right( digits(), 24)
say
say " N average exact % error " /* ◄─── title, header ►────────┐ */
hdr=" ═══ ═════════ ═════════ ═════════"; pad=left('',3) /* ◄────────┘ */
say hdr
do #=1 for runs; av=fmtD( exact(#
xa=fmtD( exper(#
say right(#,9) pad xa pad av pad fmtD( abs(xa-av) * 100 / av)
end /*#*/
say hdr /*display the final header (some bars).*/
Line 1,213 ⟶ 2,557:
/*──────────────────────────────────────────────────────────────────────────────────────*/
!: procedure expose !.; parse arg z; if !.z\==0 then return !.z
!=1; do j=2 for z -1; !=!*j; !.j=!; end;
/*──────────────────────────────────────────────────────────────────────────────────────*/
exact: parse arg x; s=0; do j=1 for x; s=s + !(x) / !(x-j) / x**j; end;
/*──────────────────────────────────────────────────────────────────────────────────────*/
exper: parse arg n; k=0; do tests; $.=0
do n; r=random(1, n); if $.r then leave
$.r=1; k=k +
end /*n*/
end /*tests*/
return k/tests
/*──────────────────────────────────────────────────────────────────────────────────────*/
fmtD: parse arg y,d; d=word(d 4, 1); y=format(y, , d); parse var y w '.' f
if f=0 then return w || left('', d +1); return y</
<pre style="height:90ex">
40 runs
Line 1,275 ⟶ 2,619:
40 7.6151 7.6091 0.0789
═══ ═════════ ═════════ ═════════
</pre>
=={{header|RPL}}==
This task is an opportunity to showcase several useful instructions - <code>CON, SEQ</code> and <code>∑</code> - which avoid to use a <code>FOR..NEXT</code> loop, to create a constant array, to generate a list of random numbers and to calculate a finite sum.
{{works with|HP|48G}}
« 0 → n times count
« 1 times '''FOR''' j
n { } + 0 CON
« n RAND * CEIL » 'z' 1 n 1 SEQ
1
'''WHILE''' 3 PICK OVER GET NOT '''REPEAT'''
ROT OVER 1 PUT ROT ROT
OVER SWAP GET
'count' 1 STO+
'''END''' 3 DROPN
'''NEXT'''
count times /
» » '<span style="color:blue">XPRMT</span>' STO
« { } DUP
1 20 '''FOR''' n
SWAP n 1000 <span style="color:blue">XPRMT</span> +
SWAP 'j' 1 n 'n!/n^j/(n-j)!' ∑ +
'''NEXT'''
DUP2 SWAP %CH ABS 2 FIX "%" ADD STD
» '<span style="color:blue">TASK</span>' STO
{{out}}
<pre>
3: { 1 1.497 1.882 2.2 2.468 2.823 3 3.258 3.475 3.647 3.854 4.003 4.234 4.46 4.589 4.767 4.852 4.929 5.154 5.251 }
2: { 1 1.5 1.88888888889 2.21875 2.5104 2.77469135802 3.0181387007 3.24501800538 3.45831574488 3.66021568 3.85237205073 4.03607367511 4.21234791298 4.38202942438 4.54580728514 4.70425824709 4.85787082081 5.00706309901 5.15219620097 5.29358458601 }
1:{ "0.00%" "0.20%" "0.36%" "0.85%" "1.69%" "1.74%" "0.60%" "0.40%" "0.48%" "0.36%" "0.04%" "0.82%" "0.51%" "1.78%" "0.95%" "1.33%" "0.12%" "1.56%" "0.04%" "0.80%" }
</pre>
=={{header|Ruby}}==
Ruby does not have a factorial method, not even in it's math library.
<
def factorial
self == 0 ? 1 : (1..self).inject(:*)
Line 1,303 ⟶ 2,678:
analytical = (1..n).inject(0){|sum, i| sum += (n.factorial / (n**i).to_f / (n-i).factorial)}
puts "%3d %8.4f %8.4f (%8.4f%%)" % [n, avg, analytical, (avg/analytical - 1)*100]
end</
{{out}}
<pre>
Line 1,331 ⟶ 2,706:
=={{header|Rust}}==
{{libheader|rand}}
<syntaxhighlight lang="rust">extern crate rand;
use rand::{ThreadRng, thread_rng};
Line 1,407 ⟶ 2,783:
</syntaxhighlight>
{{out}}
Using default arguments:
Line 1,437 ⟶ 2,813:
=={{header|Scala}}==
<syntaxhighlight lang="scala">
import scala.util.Random
object AverageLoopLength extends App {
val factorial:
val results = for (n <- 1 to 20;
avg = tested(n, 1000000);
theory = expected(n)
) yield (n, avg, theory, (avg / theory - 1) * 100)
def expected(n: Int): Double = (for (i <- 1 to n) yield factorial(n) / Math.pow(n, i) / factorial(n - i)).sum
def
def trial(n: Int): Double = {
var count = 0
var x = 1
Line 1,458 ⟶ 2,840:
count
}
Line 1,470 ⟶ 2,845:
println("------------------------------------")
results foreach { n => {
}
}
</syntaxhighlight>
{{out}}
<pre>
Line 1,500 ⟶ 2,875:
19 5.148767 5.152196 -0.067%
20 5.292875 5.293585 -0.013%
</pre>
=={{header|Scheme}}==
<syntaxhighlight lang="scheme">
(import (scheme base)
(scheme write)
(srfi 1 lists)
(only (srfi 13 strings) string-pad-right)
(srfi 27 random-bits))
(define (analytical-function n)
(define (factorial n)
(fold * 1 (iota n 1)))
;
(fold (lambda (i sum)
(+ sum
(/ (factorial n) (expt n i) (factorial (- n i)))))
0
(iota n 1)))
(define (simulation n runs)
(define (single-simulation)
(random-source-randomize! default-random-source)
(let ((vec (make-vector n #f)))
(let loop ((count 0)
(num (random-integer n)))
(if (vector-ref vec num)
count
(begin (vector-set! vec num #t)
(loop (+ 1 count)
(random-integer n)))))))
;;
(let loop ((total 0)
(run runs))
(if (zero? run)
(/ total runs)
(loop (+ total (single-simulation))
(- run 1)))))
(display " N average formula (error) \n")
(display "=== ========= ========= =========\n")
(for-each
(lambda (n)
(let ((simulation (inexact (simulation n 10000)))
(formula (inexact (analytical-function n))))
(display
(string-append
" "
(string-pad-right (number->string n) 3)
" "
(string-pad-right (number->string simulation) 6)
" "
(string-pad-right (number->string formula) 6)
" ("
(string-pad-right
(number->string (* 100 (/ (- simulation formula) formula)))
5)
"%)"))
(newline)))
(iota 20 1))
</syntaxhighlight>
{{out}}
<pre>
N average formula (error)
=== ========= ========= =========
1 1.0 1.0 (0.0 %)
2 1.5018 1.5 (0.120%)
3 1.8863 1.8888 (-0.13%)
4 2.2154 2.2187 (-0.15%)
5 2.5082 2.5104 (-0.08%)
6 2.7613 2.7746 (-0.48%)
7 3.036 3.0181 (0.591%)
8 3.2656 3.2450 (0.634%)
9 3.455 3.4583 (-0.09%)
10 3.682 3.6602 (0.595%)
11 3.8233 3.8523 (-0.75%)
12 4.0409 4.0360 (0.119%)
13 4.2471 4.2123 (0.825%)
14 4.3577 4.3820 (-0.55%)
15 4.5351 4.5458 (-0.23%)
16 4.7181 4.7042 (0.294%)
17 4.8877 4.8578 (0.614%)
18 5.0239 5.0070 (0.336%)
19 5.1216 5.1521 (-0.59%)
20 5.2717 5.2935 (-0.41%)
</pre>
=={{header|Seed7}}==
<
include "float.s7i";
Line 1,566 ⟶ 3,028:
analytical digits 4 lpad 7 <& err digits 3 lpad 7);
end for;
end func;</
{{out}}
Line 1,592 ⟶ 3,054:
20 5.2907 5.2936 0.054
</pre>
=={{header|Sidef}}==
{{trans|Perl}}
<syntaxhighlight lang="ruby">func find_loop(n) {
var seen = Hash()
loop {
with (irand(1, n)) { |r|
seen.has(r) ? (return seen.len) : (seen{r} = true)
}
}
}
print " N empiric theoric (error)\n";
print "=== ========= ============ =========\n";
define MAX = 20
define TRIALS = 1000
for n in (1..MAX) {
var empiric = (1..TRIALS -> sum { find_loop(n) } / TRIALS)
var theoric = (1..n -> sum {|k| prod(n - k + 1 .. n) * k**2 / n**(k+1) })
printf("%3d %9.4f %12.4f (%5.2f%%)\n",
n, empiric, theoric, 100*(empiric-theoric)/theoric)
}</syntaxhighlight>
{{out}}
<pre>
N empiric theoric (error)
=== ========= ============ =========
1 1.0000 1.0000 ( 0.00%)
2 1.4990 1.5000 (-0.07%)
3 1.8560 1.8889 (-1.74%)
4 2.1730 2.2188 (-2.06%)
5 2.5440 2.5104 ( 1.34%)
6 2.7490 2.7747 (-0.93%)
7 3.0390 3.0181 ( 0.69%)
8 3.1820 3.2450 (-1.94%)
9 3.4520 3.4583 (-0.18%)
10 3.6580 3.6602 (-0.06%)
11 3.9650 3.8524 ( 2.92%)
12 3.9920 4.0361 (-1.09%)
13 4.1270 4.2123 (-2.03%)
14 4.3710 4.3820 (-0.25%)
15 4.5520 4.5458 ( 0.14%)
16 4.7120 4.7043 ( 0.16%)
17 4.9400 4.8579 ( 1.69%)
18 5.0070 5.0071 (-0.00%)
19 5.2370 5.1522 ( 1.65%)
20 5.2940 5.2936 ( 0.01%)
</pre>
=={{header|Simula}}==
<syntaxhighlight lang="simula">BEGIN
REAL PROCEDURE FACTORIAL(N); INTEGER N;
BEGIN
REAL RESULT;
INTEGER I;
RESULT := 1.0;
FOR I := 2 STEP 1 UNTIL N DO
RESULT := RESULT * I;
FACTORIAL := RESULT;
END FACTORIAL;
REAL PROCEDURE ANALYTICAL (N); INTEGER N;
BEGIN
REAL SUM, RN;
INTEGER I;
RN := N;
FOR I := 1 STEP 1 UNTIL N DO
BEGIN
SUM := SUM + FACTORIAL(N) / FACTORIAL(N - I) / RN ** I;
END;
ANALYTICAL := SUM;
END ANALYTICAL;
REAL PROCEDURE EXPERIMENTAL(N); INTEGER N;
BEGIN
INTEGER NUM;
INTEGER COUNT;
INTEGER RUN;
FOR RUN := 1 STEP 1 UNTIL TESTS DO
BEGIN
BOOLEAN ARRAY BITS(1:N);
INTEGER I;
FOR I := 1 STEP 1 UNTIL N DO
BEGIN
NUM := RANDINT(1,N,SEED);
IF BITS(NUM) THEN GOTO L;
BITS(NUM) := TRUE;
COUNT := COUNT + 1;
END FOR I;
L:
END FOR RUN;
EXPERIMENTAL := COUNT / TESTS;
END EXPERIMENTAL;
INTEGER SEED, TESTS;
SEED := ININT;
TESTS := 1000000;
BEGIN
REAL A, E, ERR;
INTEGER I;
OUTTEXT(" N AVG CALC %DIFF"); OUTIMAGE;
FOR I := 1 STEP 1 UNTIL 20 DO
BEGIN
A := ANALYTICAL(I);
E := EXPERIMENTAL(I);
ERR := (ABS(E-A)/A)*100.0;
OUTINT(I, 2);
OUTFIX(E, 4, 7);
OUTFIX(A, 4, 10);
OUTFIX(ERR, 4, 10);
OUTIMAGE;
END FOR I;
END;
END</syntaxhighlight>
{{in}}
<pre>678</pre>
{{out}}
<pre> N AVG CALC %DIFF
1 1.0000 1.0000 0.0000
2 1.4999 1.5000 0.0075
3 1.8890 1.8889 0.0072
4 2.2182 2.2188 0.0243
5 2.5105 2.5104 0.0027
6 2.7746 2.7747 0.0025
7 3.0164 3.0181 0.0590
8 3.2447 3.2450 0.0110
9 3.4567 3.4583 0.0453
10 3.6622 3.6602 0.0539
11 3.8503 3.8524 0.0546
12 4.0373 4.0361 0.0300
13 4.2105 4.2123 0.0445
14 4.3819 4.3820 0.0027
15 4.5475 4.5458 0.0376
16 4.7056 4.7043 0.0295
17 4.8559 4.8579 0.0396
18 5.0105 5.0071 0.0694
19 5.1541 5.1522 0.0376
20 5.2961 5.2936 0.0467</pre>
=={{header|Tcl}}==
<
proc range {a b} {
for {set result {}} {$a <= $b} {incr a} {lappend result $a}
Line 1,633 ⟶ 3,236:
set e [Experimental $n 100000]
puts [format "%3d %9.4f %12.4f (%6.2f%%)" $n $e $a [expr {abs($e-$a)/$a*100.0}]]
}</
{{out}}
<pre>
Line 1,659 ⟶ 3,262:
20 5.2992 5.2936 ( 0.11%)
</pre>
=={{header|uBasic/4tH}}==
{{trans|BBC BASIC}}
This is about the limit of what you can do with uBasic/4tH. Since it is an integer BASIC, it uses what has become known in the Forth community as "Brodie math". The last 14 bits of an 64-bit integer are used to represent the fraction, so basically it is a form of "fixed point math". This, of course, leads inevitably to rounding errors. After step 14 the number is too large to fit in a 64-bit integer - so at that point it simply breaks down. Performance is another issue.
<syntaxhighlight lang="uBasic/4tH">M = 14
T = 100000
If Info("wordsize") < 64 Then Print "This program requires a 64-bit uBasic" : End
Print "N\tavg\tcalc\t%diff"
For n = 1 To M
a = FUNC(_Test(n, T))
h = FUNC(_Analytical(n))
d = FUNC(_Fmul(FUNC(_Fdiv(a, h)) - FUNC(_Ntof(1)), FUNC(_Ntof(100))))
Print n; "\t";
Proc _Fprint (a) : Print "\t";
Proc _Fprint (h) : Print "\t";
Proc _Fprint (d) : Print "%"
Next
End
_Analytical
Param (1)
Local (3)
c@ = 0
For b@ = 1 To a@
d@ = FUNC(_Fdiv(FUNC(_Factorial(a@)), a@^b@))
c@ = c@ + FUNC(_Fdiv (d@, FUNC(_Ntof(FUNC(_Factorial(a@-b@))))))
Next
Return (c@)
_Test
Param (2)
Local (4)
e@ = 0
For c@ = 1 To b@
f@ = 1 : d@ = 0
Do While AND(d@, f@) = 0
e@ = e@ + 1
d@ = OR(d@, f@)
f@ = SHL(1, Rnd(a@))
Loop
Next
Return (FUNC(_Fdiv(e@, b@)))
_Factorial
Param(1)
If (a@ = 1) + (a@ = 0) Then Return (1)
Return (a@ * FUNC(_Factorial(a@-1)))
_Fmul Param (2) : Return ((a@*b@)/16384)
_Fdiv Param (2) : Return ((a@*16384)/b@)
_Ftoi Param (1) : Return ((10000*a@)/16384)
_Itof Param (1) : Return ((16384*a@)/10000)
_Ntof Param (1) : Return (16384*a@)
_Fprint Param (1) : a@ = FUNC(_Ftoi(a@)) : Print Using "+?.####";a@; : Return</syntaxhighlight>
{{Out}}
<pre>N avg calc %diff
1 1.0000 1.0000 0.0000%
2 1.4985 1.5000 -0.0976%
3 1.8869 1.8887 -0.1037%
4 2.2192 2.2187 0.0183%
5 2.5130 2.5103 0.1037%
6 2.7761 2.7745 0.0549%
7 3.0264 3.0180 0.2746%
8 3.2504 3.2449 0.1647%
9 3.4528 3.4581 -0.1525%
10 3.6651 3.6599 0.1403%
11 3.8543 3.8521 0.0549%
12 4.0364 4.0357 0.0122%
13 4.2153 4.2119 0.0793%
14 4.3866 4.3815 0.1098%
0 OK, 0:392</pre>
=={{header|Unicon}}==
{{trans|C}}
<
$define MAX_N 20
Line 1,700 ⟶ 3,383:
}
return 0
end</
{{out}}
<pre> n avg exp. diff
Line 1,724 ⟶ 3,407:
19 5.1533 5.1522 0.021%
20 5.2893 5.2936 -0.081%</pre>
=={{header|VBA}}==
{{trans|Phix}}
<syntaxhighlight lang="vb">Const MAX = 20
Const ITER = 1000000
Function expected(n As Long) As Double
Dim sum As Double
For i = 1 To n
sum = sum + WorksheetFunction.Fact(n) / n ^ i / WorksheetFunction.Fact(n - i)
Next i
expected = sum
End Function
Function test(n As Long) As Double
Dim count As Long
Dim x As Long, bits As Long
For i = 1 To ITER
x = 1
bits = 0
Do While Not bits And x
count = count + 1
bits = bits Or x
x = 2 ^ (Int(n * Rnd()))
Loop
Next i
test = count / ITER
End Function
Public Sub main()
Dim n As Long
Debug.Print " n avg. exp. (error%)"
Debug.Print "== ====== ====== ========"
For n = 1 To MAX
av = test(n)
ex = expected(n)
Debug.Print Format(n, "@@"); " "; Format(av, "0.0000"); " ";
Debug.Print Format(ex, "0.0000"); " ("; Format(Abs(1 - av / ex), "0.000%"); ")"
Next n
End Sub</syntaxhighlight>{{out}}
<pre> n avg. exp. (error%)
== ====== ====== ========
1 1,0000 1,0000 (0,000%)
2 1,4994 1,5000 (0,041%)
3 1,8893 1,8889 (0,023%)
4 2,2187 2,2188 (0,001%)
5 2,5107 2,5104 (0,010%)
6 2,7769 2,7747 (0,080%)
7 3,0162 3,0181 (0,064%)
8 3,2472 3,2450 (0,066%)
9 3,4603 3,4583 (0,056%)
10 3,6577 3,6602 (0,070%)
11 3,8527 3,8524 (0,010%)
12 4,0361 4,0361 (0,001%)
13 4,2121 4,2123 (0,005%)
14 4,3825 4,3820 (0,010%)
15 4,5466 4,5458 (0,016%)
16 4,7023 4,7043 (0,041%)
17 4,8567 4,8579 (0,025%)
18 5,0031 5,0071 (0,079%)
19 5,1530 5,1522 (0,016%)
20 5,2958 5,2936 (0,041%)</pre>
=={{header|VBScript}}==
Ported from the VBA version. I added some precalculations to speed it up
<syntaxhighlight lang="vb">
Const MAX = 20
Const ITER = 100000
Function expected(n)
Dim sum
ni=n
For i = 1 To n
sum = sum + fact(n) / ni / fact(n-i)
ni=ni*n
Next
expected = sum
End Function
Function test(n )
Dim coun,x,bits
For i = 1 To ITER
x = 1
bits = 0
Do While Not bits And x
count = count + 1
bits = bits Or x
x =shift(Int(n * Rnd()))
Loop
Next
test = count / ITER
End Function
'VBScript formats numbers but does'nt align them!
function rf(v,n,s) rf=right(string(n,s)& v,n):end function
'some precalculations to speed things up...
dim fact(20),shift(20)
fact(0)=1:shift(0)=1
for i=1 to 20
fact(i)=i*fact(i-1)
shift(i)=2*shift(i-1)
next
Dim n
Wscript.echo "For " & ITER &" iterations"
Wscript.Echo " n avg. exp. (error%)"
Wscript.Echo "== ====== ====== =========="
For n = 1 To MAX
av = test(n)
ex = expected(n)
Wscript.Echo rf(n,2," ")& " "& rf(formatnumber(av, 4),7," ") & " "& _
rf(formatnumber(ex,4),6," ")& " ("& rf(Formatpercent(1 - av / ex,4),8," ") & ")"
Next
</syntaxhighlight>
Output
<pre>
For 100000 iterations
n avg. exp. (error%)
== ====== ====== ==========
1 1.0000 1.0000 ( 0.0000%)
2 1.4982 1.5000 ( 0.0012%)
3 1.8909 1.8889 (-0.0010%)
4 2.2190 2.2188 (-0.0001%)
5 2.5102 2.5104 ( 0.0001%)
6 2.7789 2.7747 (-0.0015%)
7 3.0230 3.0181 (-0.0016%)
8 3.2449 3.2450 ( 0.0000%)
9 3.4543 3.4583 ( 0.0012%)
10 3.6714 3.6602 (-0.0031%)
11 3.8559 3.8524 (-0.0009%)
12 4.0345 4.0361 ( 0.0004%)
13 4.2141 4.2123 (-0.0004%)
14 4.3762 4.3820 ( 0.0013%)
15 4.5510 4.5458 (-0.0011%)
16 4.6979 4.7043 ( 0.0014%)
17 4.8628 4.8579 (-0.0010%)
18 5.0081 5.0071 (-0.0002%)
19 5.1518 5.1522 ( 0.0001%)
20 5.2906 5.2936 ( 0.0006%)
</pre>
=={{header|V (Vlang)}}==
{{trans|Go}}
<syntaxhighlight lang="v (vlang)">import rand
import math
const nmax = 20
fn main() {
println(" N average analytical (error)")
println("=== ========= ============ =========")
for n := 1; n <= nmax; n++ {
a := avg(n)
b := ana(n)
println("${n:3} ${a:9.4f} ${b:12.4f} (${math.abs(a-b)/b*100:6.2f}%)" )
}
}
fn avg(n int) f64 {
tests := int(1e4)
mut sum := 0
for _ in 0..tests {
mut v := [nmax]bool{}
for x := 0; !v[x]; {
v[x] = true
sum++
x = rand.intn(n) or {0}
}
}
return f64(sum) / tests
}
fn ana(n int) f64 {
nn := f64(n)
mut term := 1.0
mut sum := 1.0
for i := nn - 1; i >= 1; i-- {
term *= i / nn
sum += term
}
return sum
}</syntaxhighlight>
{{out}}
Sample output:
<pre>
N average analytical (error)
=== ========= ============ =========
1 1.0000 1.0000 ( 0.00%)
2 1.4967 1.5000 ( 0.22%)
3 1.8970 1.8889 ( 0.43%)
4 2.2151 2.2188 ( 0.16%)
5 2.5044 2.5104 ( 0.24%)
6 2.7884 2.7747 ( 0.49%)
7 3.0356 3.0181 ( 0.58%)
8 3.2468 3.2450 ( 0.05%)
9 3.4692 3.4583 ( 0.31%)
10 3.6538 3.6602 ( 0.18%)
11 3.8325 3.8524 ( 0.52%)
12 4.0674 4.0361 ( 0.78%)
13 4.2199 4.2123 ( 0.18%)
14 4.3808 4.3820 ( 0.03%)
15 4.5397 4.5458 ( 0.13%)
16 4.6880 4.7043 ( 0.35%)
17 4.8554 4.8579 ( 0.05%)
18 5.0311 5.0071 ( 0.48%)
19 5.1577 5.1522 ( 0.11%)
20 5.2995 5.2936 ( 0.11%)
</pre>
=={{header|Wren}}==
{{trans|Go}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="wren">import "random" for Random
import "./fmt" for Fmt
var nmax = 20
var rand = Random.new()
var avg = Fn.new { |n|
var tests = 1e4
var sum = 0
for (t in 0...tests) {
var v = List.filled(nmax, false)
var x = 0
while (!v[x]) {
v[x] = true
sum = sum + 1
x = rand.int(n)
}
}
return sum/tests
}
var ana = Fn.new { |n|
if (n < 2) return 1
var term = 1
var sum = 1
for (i in n-1..1) {
term = term * i / n
sum = sum + term
}
return sum
}
System.print(" N average analytical (error)")
System.print("=== ========= ============ =========")
for (n in 1..nmax) {
var a = avg.call(n)
var b = ana.call(n)
var e = (a - b).abs/ b * 100
Fmt.print("$3d $9.4f $12.4f ($6.2f\%)", n, a, b, e)
}</syntaxhighlight>
{{out}}
Sample output:
<pre>
N average analytical (error)
=== ========= ============ =========
1 1.0000 1.0000 ( 0.00%)
2 1.4967 1.5000 ( 0.22%)
3 1.8970 1.8889 ( 0.43%)
4 2.2151 2.2188 ( 0.16%)
5 2.5044 2.5104 ( 0.24%)
6 2.7884 2.7747 ( 0.49%)
7 3.0356 3.0181 ( 0.58%)
8 3.2468 3.2450 ( 0.05%)
9 3.4692 3.4583 ( 0.31%)
10 3.6538 3.6602 ( 0.18%)
11 3.8325 3.8524 ( 0.52%)
12 4.0674 4.0361 ( 0.78%)
13 4.2199 4.2123 ( 0.18%)
14 4.3808 4.3820 ( 0.03%)
15 4.5397 4.5458 ( 0.13%)
16 4.6880 4.7043 ( 0.35%)
17 4.8554 4.8579 ( 0.05%)
18 5.0311 5.0071 ( 0.48%)
19 5.1577 5.1522 ( 0.11%)
20 5.2995 5.2936 ( 0.11%)
</pre>
=={{header|zkl}}==
<
(" N average analytical (error)").println();
Line 1,758 ⟶ 3,722:
n=n.toFloat();
(1).reduce(n,'wrap(sum,i){ sum+fact(n)/n.pow(i)/fact(n-i) },0.0);
}</
{{out}}
<pre>
| |||