Calculate the Mean mean of two numbersAGM Series Hole 1: Calculate the Arithmetic–Geometric MeanSum the...
How to type dʒ symbol (IPA) on Mac?
What do you call a Matrix-like slowdown and camera movement effect?
How can I fix this gap between bookcases I made?
Non-Jewish family in an Orthodox Jewish Wedding
Is Social Media Science Fiction?
How is this relation reflexive?
How does one intimidate enemies without having the capacity for violence?
Should I join an office cleaning event for free?
How can bays and straits be determined in a procedurally generated map?
Why is an old chain unsafe?
What does "enim et" mean?
Infinite past with a beginning?
Are white and non-white police officers equally likely to kill black suspects?
My colleague's body is amazing
Email Account under attack (really) - anything I can do?
Why do we use polarized capacitor?
The magic money tree problem
Can a German sentence have two subjects?
Shell script can be run only with sh command
Why doesn't Newton's third law mean a person bounces back to where they started when they hit the ground?
New order #4: World
Why has Russell's definition of numbers using equivalence classes been finally abandoned? ( If it has actually been abandoned).
Circuitry of TV splitters
What makes Graph invariants so useful/important?
Calculate the Mean mean of two numbers
AGM Series Hole 1: Calculate the Arithmetic–Geometric MeanSum the means of the two integersConjugation in Real LifeRandom number generation, without built-in random number generation librariesMusic: what's in this chord?Rosetta Stone Challenge: What's Average Anyways?7-segment differencesPretty Smooth MovesFind the Fibonacci KernelNot just a four fours puzzleStandardize the Samples (Compute the z-Score)Eight coins for the fair kingSquare Root Distance from Integers
$begingroup$
disclaimer: the Mean mean is made up by me
Define the arithmetic mean of $n$ numbers as
$$M_1(x_1,...,x_n)=frac{x_1+x_2+...+x_n}{n}$$
Define the geometric mean of $n$ numbers as
$$M_0(x_1,...,x_n)=root{n}of{x_1x_2...x_n}$$
Define the harmonic mean of $n$ numbers as
$$M_{-1}(x_1,...,x_n)=frac{n}{frac{1}{x_2} + frac{1}{x_2} + ... + frac{1}{x_n}}$$
Define the quadratic mean of $n$ numbers as
$$M_2(x_1,...,x_n)=rootof{frac{x_1^2+x_2^2+...+x_n^2}{n}}$$
The Mean mean ($M_M$) is defined as follows:
Define four sequences ($a_k, b_k, c_k, d_k$) as $$a_0=M_1(x_1,...,x_n),\b_0=M_0(x_1,...,x_n),\c_0=M_{-1}(x_1,...,x_n),\d_0=M_2(x_1,...,x_n),\
a_{k+1}=M_1(a_k,b_k,c_k,d_k),\b_{k+1}=M_0(a_k,b_k,c_k,d_k),\c_{k+1}=M_{-1}(a_k,b_k,c_k,d_k),\d_{k+1}=M_2(a_k,b_k,c_k,d_k)$$
All four sequences converge to the same number, $M_M(x_1,x_2,...,x_n)$.
Example
The Mean mean of 1 and 2 is calculated as follows: start with $$a_0 = (1+2)/2 = 1.5, b_0 = rootof{1 * 2} = rootof2 approx 1.4142,\ c_0 = frac2{frac1{1}+frac1{2}} = frac4{3} approx 1.3333, d_0 = rootof{frac{1^2+2^2}2} = rootof{frac5{2}} approx 1.5811.$$
Then
$$a_1 = frac{1.5+1.4142+1.3333+1.5811}4 approx 1.4571,\ b_1 = root^4of{1.5*1.4142*1.3333*1.5811} approx 1.4542,\ c_1 = frac4{frac1{1.5}+frac1{1.4142}+frac1{1.3333}+frac1{1.5811}} approx 1.4512,\ d_1 = rootof{frac{1.5^2+1.4142^2+1.3333^2+1.5811^2}4} approx 1.4601.$$
The further calculation of the sequences should be clear. It can be seen that they converge to the same number, approximately $1.45568889$.
Challenge
Given two positive real numbers, $a$ and $b$ ($a<b$), calculate their Mean mean $M_M(a,b)$.
Test cases
1 1 => 1
1 2 => 1.45568889
100 200 => 145.568889
2.71 3.14 => 2.92103713
0.57 1.78 => 1.0848205
1.61 2.41 => 1.98965438
0.01 100 => 6.7483058
Notes
- Your program is valid if the difference between its output and the correct output is not greater than 1/100000 of the absolute value of the difference between input numbers.
- The output should be a single number.
This is code-golf, so the shortest code wins!
code-golf math
$endgroup$
add a comment |
$begingroup$
disclaimer: the Mean mean is made up by me
Define the arithmetic mean of $n$ numbers as
$$M_1(x_1,...,x_n)=frac{x_1+x_2+...+x_n}{n}$$
Define the geometric mean of $n$ numbers as
$$M_0(x_1,...,x_n)=root{n}of{x_1x_2...x_n}$$
Define the harmonic mean of $n$ numbers as
$$M_{-1}(x_1,...,x_n)=frac{n}{frac{1}{x_2} + frac{1}{x_2} + ... + frac{1}{x_n}}$$
Define the quadratic mean of $n$ numbers as
$$M_2(x_1,...,x_n)=rootof{frac{x_1^2+x_2^2+...+x_n^2}{n}}$$
The Mean mean ($M_M$) is defined as follows:
Define four sequences ($a_k, b_k, c_k, d_k$) as $$a_0=M_1(x_1,...,x_n),\b_0=M_0(x_1,...,x_n),\c_0=M_{-1}(x_1,...,x_n),\d_0=M_2(x_1,...,x_n),\
a_{k+1}=M_1(a_k,b_k,c_k,d_k),\b_{k+1}=M_0(a_k,b_k,c_k,d_k),\c_{k+1}=M_{-1}(a_k,b_k,c_k,d_k),\d_{k+1}=M_2(a_k,b_k,c_k,d_k)$$
All four sequences converge to the same number, $M_M(x_1,x_2,...,x_n)$.
Example
The Mean mean of 1 and 2 is calculated as follows: start with $$a_0 = (1+2)/2 = 1.5, b_0 = rootof{1 * 2} = rootof2 approx 1.4142,\ c_0 = frac2{frac1{1}+frac1{2}} = frac4{3} approx 1.3333, d_0 = rootof{frac{1^2+2^2}2} = rootof{frac5{2}} approx 1.5811.$$
Then
$$a_1 = frac{1.5+1.4142+1.3333+1.5811}4 approx 1.4571,\ b_1 = root^4of{1.5*1.4142*1.3333*1.5811} approx 1.4542,\ c_1 = frac4{frac1{1.5}+frac1{1.4142}+frac1{1.3333}+frac1{1.5811}} approx 1.4512,\ d_1 = rootof{frac{1.5^2+1.4142^2+1.3333^2+1.5811^2}4} approx 1.4601.$$
The further calculation of the sequences should be clear. It can be seen that they converge to the same number, approximately $1.45568889$.
Challenge
Given two positive real numbers, $a$ and $b$ ($a<b$), calculate their Mean mean $M_M(a,b)$.
Test cases
1 1 => 1
1 2 => 1.45568889
100 200 => 145.568889
2.71 3.14 => 2.92103713
0.57 1.78 => 1.0848205
1.61 2.41 => 1.98965438
0.01 100 => 6.7483058
Notes
- Your program is valid if the difference between its output and the correct output is not greater than 1/100000 of the absolute value of the difference between input numbers.
- The output should be a single number.
This is code-golf, so the shortest code wins!
code-golf math
$endgroup$
$begingroup$
Sandbox link
$endgroup$
– someone
Mar 31 at 7:41
$begingroup$
Somewhat related: Calculate the Arithmetic–Geometric Mean
$endgroup$
– user202729
Mar 31 at 13:38
8
$begingroup$
How precise are we supposed to be?
$endgroup$
– Embodiment of Ignorance
Mar 31 at 17:52
2
$begingroup$
Closely related
$endgroup$
– Giuseppe
Apr 1 at 13:18
1
$begingroup$
Can we assume the first input is always smaller than the second, as in all your test cases? (If not I'll rollback my Java answer.)
$endgroup$
– Kevin Cruijssen
Apr 2 at 7:59
add a comment |
$begingroup$
disclaimer: the Mean mean is made up by me
Define the arithmetic mean of $n$ numbers as
$$M_1(x_1,...,x_n)=frac{x_1+x_2+...+x_n}{n}$$
Define the geometric mean of $n$ numbers as
$$M_0(x_1,...,x_n)=root{n}of{x_1x_2...x_n}$$
Define the harmonic mean of $n$ numbers as
$$M_{-1}(x_1,...,x_n)=frac{n}{frac{1}{x_2} + frac{1}{x_2} + ... + frac{1}{x_n}}$$
Define the quadratic mean of $n$ numbers as
$$M_2(x_1,...,x_n)=rootof{frac{x_1^2+x_2^2+...+x_n^2}{n}}$$
The Mean mean ($M_M$) is defined as follows:
Define four sequences ($a_k, b_k, c_k, d_k$) as $$a_0=M_1(x_1,...,x_n),\b_0=M_0(x_1,...,x_n),\c_0=M_{-1}(x_1,...,x_n),\d_0=M_2(x_1,...,x_n),\
a_{k+1}=M_1(a_k,b_k,c_k,d_k),\b_{k+1}=M_0(a_k,b_k,c_k,d_k),\c_{k+1}=M_{-1}(a_k,b_k,c_k,d_k),\d_{k+1}=M_2(a_k,b_k,c_k,d_k)$$
All four sequences converge to the same number, $M_M(x_1,x_2,...,x_n)$.
Example
The Mean mean of 1 and 2 is calculated as follows: start with $$a_0 = (1+2)/2 = 1.5, b_0 = rootof{1 * 2} = rootof2 approx 1.4142,\ c_0 = frac2{frac1{1}+frac1{2}} = frac4{3} approx 1.3333, d_0 = rootof{frac{1^2+2^2}2} = rootof{frac5{2}} approx 1.5811.$$
Then
$$a_1 = frac{1.5+1.4142+1.3333+1.5811}4 approx 1.4571,\ b_1 = root^4of{1.5*1.4142*1.3333*1.5811} approx 1.4542,\ c_1 = frac4{frac1{1.5}+frac1{1.4142}+frac1{1.3333}+frac1{1.5811}} approx 1.4512,\ d_1 = rootof{frac{1.5^2+1.4142^2+1.3333^2+1.5811^2}4} approx 1.4601.$$
The further calculation of the sequences should be clear. It can be seen that they converge to the same number, approximately $1.45568889$.
Challenge
Given two positive real numbers, $a$ and $b$ ($a<b$), calculate their Mean mean $M_M(a,b)$.
Test cases
1 1 => 1
1 2 => 1.45568889
100 200 => 145.568889
2.71 3.14 => 2.92103713
0.57 1.78 => 1.0848205
1.61 2.41 => 1.98965438
0.01 100 => 6.7483058
Notes
- Your program is valid if the difference between its output and the correct output is not greater than 1/100000 of the absolute value of the difference between input numbers.
- The output should be a single number.
This is code-golf, so the shortest code wins!
code-golf math
$endgroup$
disclaimer: the Mean mean is made up by me
Define the arithmetic mean of $n$ numbers as
$$M_1(x_1,...,x_n)=frac{x_1+x_2+...+x_n}{n}$$
Define the geometric mean of $n$ numbers as
$$M_0(x_1,...,x_n)=root{n}of{x_1x_2...x_n}$$
Define the harmonic mean of $n$ numbers as
$$M_{-1}(x_1,...,x_n)=frac{n}{frac{1}{x_2} + frac{1}{x_2} + ... + frac{1}{x_n}}$$
Define the quadratic mean of $n$ numbers as
$$M_2(x_1,...,x_n)=rootof{frac{x_1^2+x_2^2+...+x_n^2}{n}}$$
The Mean mean ($M_M$) is defined as follows:
Define four sequences ($a_k, b_k, c_k, d_k$) as $$a_0=M_1(x_1,...,x_n),\b_0=M_0(x_1,...,x_n),\c_0=M_{-1}(x_1,...,x_n),\d_0=M_2(x_1,...,x_n),\
a_{k+1}=M_1(a_k,b_k,c_k,d_k),\b_{k+1}=M_0(a_k,b_k,c_k,d_k),\c_{k+1}=M_{-1}(a_k,b_k,c_k,d_k),\d_{k+1}=M_2(a_k,b_k,c_k,d_k)$$
All four sequences converge to the same number, $M_M(x_1,x_2,...,x_n)$.
Example
The Mean mean of 1 and 2 is calculated as follows: start with $$a_0 = (1+2)/2 = 1.5, b_0 = rootof{1 * 2} = rootof2 approx 1.4142,\ c_0 = frac2{frac1{1}+frac1{2}} = frac4{3} approx 1.3333, d_0 = rootof{frac{1^2+2^2}2} = rootof{frac5{2}} approx 1.5811.$$
Then
$$a_1 = frac{1.5+1.4142+1.3333+1.5811}4 approx 1.4571,\ b_1 = root^4of{1.5*1.4142*1.3333*1.5811} approx 1.4542,\ c_1 = frac4{frac1{1.5}+frac1{1.4142}+frac1{1.3333}+frac1{1.5811}} approx 1.4512,\ d_1 = rootof{frac{1.5^2+1.4142^2+1.3333^2+1.5811^2}4} approx 1.4601.$$
The further calculation of the sequences should be clear. It can be seen that they converge to the same number, approximately $1.45568889$.
Challenge
Given two positive real numbers, $a$ and $b$ ($a<b$), calculate their Mean mean $M_M(a,b)$.
Test cases
1 1 => 1
1 2 => 1.45568889
100 200 => 145.568889
2.71 3.14 => 2.92103713
0.57 1.78 => 1.0848205
1.61 2.41 => 1.98965438
0.01 100 => 6.7483058
Notes
- Your program is valid if the difference between its output and the correct output is not greater than 1/100000 of the absolute value of the difference between input numbers.
- The output should be a single number.
This is code-golf, so the shortest code wins!
code-golf math
code-golf math
edited Apr 2 at 8:12
someone
asked Mar 31 at 7:39
someonesomeone
732720
732720
$begingroup$
Sandbox link
$endgroup$
– someone
Mar 31 at 7:41
$begingroup$
Somewhat related: Calculate the Arithmetic–Geometric Mean
$endgroup$
– user202729
Mar 31 at 13:38
8
$begingroup$
How precise are we supposed to be?
$endgroup$
– Embodiment of Ignorance
Mar 31 at 17:52
2
$begingroup$
Closely related
$endgroup$
– Giuseppe
Apr 1 at 13:18
1
$begingroup$
Can we assume the first input is always smaller than the second, as in all your test cases? (If not I'll rollback my Java answer.)
$endgroup$
– Kevin Cruijssen
Apr 2 at 7:59
add a comment |
$begingroup$
Sandbox link
$endgroup$
– someone
Mar 31 at 7:41
$begingroup$
Somewhat related: Calculate the Arithmetic–Geometric Mean
$endgroup$
– user202729
Mar 31 at 13:38
8
$begingroup$
How precise are we supposed to be?
$endgroup$
– Embodiment of Ignorance
Mar 31 at 17:52
2
$begingroup$
Closely related
$endgroup$
– Giuseppe
Apr 1 at 13:18
1
$begingroup$
Can we assume the first input is always smaller than the second, as in all your test cases? (If not I'll rollback my Java answer.)
$endgroup$
– Kevin Cruijssen
Apr 2 at 7:59
$begingroup$
Sandbox link
$endgroup$
– someone
Mar 31 at 7:41
$begingroup$
Sandbox link
$endgroup$
– someone
Mar 31 at 7:41
$begingroup$
Somewhat related: Calculate the Arithmetic–Geometric Mean
$endgroup$
– user202729
Mar 31 at 13:38
$begingroup$
Somewhat related: Calculate the Arithmetic–Geometric Mean
$endgroup$
– user202729
Mar 31 at 13:38
8
8
$begingroup$
How precise are we supposed to be?
$endgroup$
– Embodiment of Ignorance
Mar 31 at 17:52
$begingroup$
How precise are we supposed to be?
$endgroup$
– Embodiment of Ignorance
Mar 31 at 17:52
2
2
$begingroup$
Closely related
$endgroup$
– Giuseppe
Apr 1 at 13:18
$begingroup$
Closely related
$endgroup$
– Giuseppe
Apr 1 at 13:18
1
1
$begingroup$
Can we assume the first input is always smaller than the second, as in all your test cases? (If not I'll rollback my Java answer.)
$endgroup$
– Kevin Cruijssen
Apr 2 at 7:59
$begingroup$
Can we assume the first input is always smaller than the second, as in all your test cases? (If not I'll rollback my Java answer.)
$endgroup$
– Kevin Cruijssen
Apr 2 at 7:59
add a comment |
19 Answers
19
active
oldest
votes
$begingroup$
Wolfram Language (Mathematica), 52 bytes
#//.x_:>N@{M@x,E^M@Log@x,1/M[1/x],M[x^2]^.5}&
M=Mean
Try it online!
In my first approach I used these builtinsMean
GeometricMean
HarmonicMean
and RootMeanSquare
Here are some substitutions for saving bytes
HarmonicMean
--> 1/Mean[1/x]
by @Robin Ryder (3 bytes saved)GeometricMean
--> E^Mean@Log@x
by @A. Rex (2 bytes saved)RootMeanSquare
--> Mean[x^2]^.5
by @A. Rex (4 bytes saved)
finally we can assign Mean
to M
(as proposed by @ovs) and save 5 more bytes
$endgroup$
$begingroup$
Save 2 bytes by recoding GeometricMean
$endgroup$
– Robin Ryder
Apr 1 at 9:58
$begingroup$
@RobinRyder I believe you mean Harmonic.. nice!
$endgroup$
– J42161217
Apr 1 at 10:04
1
$begingroup$
Save 8 more bytes:#//.x_:>N@{Mean@x,E^Mean@Log@x,1/Mean[1/x],Mean[x^2]^.5}&
$endgroup$
– A. Rex
Apr 1 at 10:41
$begingroup$
@ovs edited.....
$endgroup$
– J42161217
Apr 1 at 11:27
add a comment |
$begingroup$
R, 70 69 67 bytes
x=scan();`?`=mean;while(x-?x)x=c((?x^2)^.5,?x,2^?log2(x),1/?1/x);?x
Try it online!
-1 byte with better conditioning.
-2 bytes by switching to base 2.
Like some other answers, this uses the expression of the geometric mean as an arithmetic mean on the log scale (here in base 2):
$$M_0(x_1,ldots, x_n) = 2^{M_1(log_2 x_1,ldots,log_2 x_n)}.$$
It also uses the fact that $forall k, d_kgeq a_k geq b_kgeq c_k$, i.e. $d_k=max(a_k, b_k, c_k, d_k)$. The condition $a_k=b_k=c_k=d_k$ is therefore equivalent to $d_k = M_1(a_k, b_k, c_k, d_k)$, which is what I use in the while loop; this is achieved by abusing the syntax of while
which only considers the first element when the condition is a vector, hence the order in which the means are stored. (Note that we could also use $c_k$ instead since it is the minimum of the four, but we could not use $a_k$ or $b_k$ in the condition.)
When we exit the while loop, x
is a constant vector. The final ?x
computes its mean to reduce it to a scalar.
$endgroup$
1
$begingroup$
Shouldn't it be $ln x_n$ instead of $log x_n$ ?
$endgroup$
– Tau
Mar 31 at 23:40
$begingroup$
@Tau Yes, I was denoting natural logarithm by $ log$, which is the default in R. Anyway, I have now changed it to base 2 logarithm for -2 bytes.
$endgroup$
– Robin Ryder
Apr 1 at 7:14
add a comment |
$begingroup$
J, 34 bytes
(31 as an expression without the assignment to variable f
)
f=:1{(^.z,%z,*:z,[z=:(+/%#)&.:)^:_
Try it online!
For functions a
and b
, a &.: b
("a under b" (related challenge)) is equivalent to (b inv) a b
-- apply b, then a, then inverse of b. In this case, geometric/harmonic/quadratic mean is the arithmetic mean "under" logarithm, inversion, and square respectively.
$endgroup$
add a comment |
$begingroup$
TI-BASIC, 42 35 34 bytes
-1 byte thanks to @SolomonUcko
While max(ΔList(Ans:{mean(Ans),√(mean(Ans²)),mean(Ans^-1)^-1,e^(mean(ln(Ans:End:Ans(1
Input is a list of two integers in Ans
.
Output is stored in Ans
and is automatically printed out when the program completes.
Formulas used for geometric, harmonic, and quadratic means are based off of user202729's explanation.
Example:
{1,2
{1 2}
prgmCDGFB
1.455688891
{100,200
{100 200}
prgmCDGFB
145.5688891
Explanation:
(Newlines have been added for clarification. They do NOT appear in the code.)
While max(ΔList(Ans ;loop until all elements of the current list are equal
; the maximum of the change in each element will be 0
{ ;create a list containing...
mean(Ans), ; the arithmetic mean
√(mean(Ans²)), ; the quadratic mean
mean(Ans^-1)^-1, ; the harmonic mean
e^(mean(ln(Ans ; and the geometric mean
End
Ans(1 ;keep the first element in "Ans" and implicitly print it
Notes:
TI-BASIC is a tokenized language. Character count does not equal byte count.
e^(
is this one-byte token.
^-1
is used for this one-byte token.
I opted for writing ^-1
instead because the token looks like ֿ¹
when in a code block.
√(
is this one-byte token.
ΔList(
is this two-byte token.
$endgroup$
$begingroup$
I think you can save a parenthesis by putting the geometric mean last.
$endgroup$
– Solomon Ucko
Apr 1 at 11:06
$begingroup$
@SolomonUcko ah, thanks for noticing! Didn't consider that before.
$endgroup$
– Tau
Apr 1 at 13:01
add a comment |
$begingroup$
Java 10, 234 229 214 211 215 206 203 196 180 177 bytes
a->{for(;a[1]-a[0]>4e-9;){double l=a.length,A[]={0,0,0,1};for(var d:a){A[2]+=d/l;A[3]*=Math.pow(d,1/l);A[0]+=1/d;A[1]+=d*d;}A[0]=l/A[0];A[1]=Math.sqrt(A[1]/l);a=A;}return a[0];}
-5 bytes thanks to @PeterCordes.
-15 more bytes thanks to @PeterCordes, inspired by @RobinRyder's R answer.
+4 bytes because I assumed the inputs are pre-ordered.
-27 bytes thanks to @OlivierGrégoire.
Try it online.
Explanation:
a->{ // Method with double-array parameter and double return-type
for(;a[1]-a[0] // Loop as long as the difference between the 2nd and 1st items
>4e-9;){ // is larger than 0.000000004:
double l=a.length, // Set `l` to the amount of values in the array `a`
A[]={0,0,0,1}; // Create an array `A`, filled with the values [0,0,0,1]
for(var d:a){ // Inner loop over the values of `a`:
A[2]+=d/l; // Calculate the sum divided by the length in the third spot
A[3]*=Math.pow(d,1/l);// The product of the power of 1/length in the fourth spot
A[0]+=1/d; // The sum of 1/value in the first spot
A[1]+=d*d; // And the sum of squares in the second spot
} // After the inner loop:
// (the third spot of the array now holds the Arithmetic Mean)
// (the fourth spot of the array now holds the Geometric Mean)
A[0]=l/A[0]; // Divide the length by the first spot
// (this first spot of the array now holds the Harmonic Mean)
A[1]=Math.sqrt(A[1]/l); // Take the square of the second spot divided by the length
// (this second spot of the array now holds the Quadratic Mean)
a=A; // And then replace input `a` with array `A`
} // After the outer loop when all values are approximately equal:
return a[0];} // Return the value in the first spot as result
$endgroup$
$begingroup$
In C you couldf+=Math.abs(d-D)<1e-9;
and get implicit conversion from a boolean compare result to a 0 / 1 integer and thendouble
. Does Java have any compact syntax for that? Or is it possible to dof+=Math.abs(d-D)
and then check that the sum of absolute differences is small enough?
$endgroup$
– Peter Cordes
Apr 2 at 2:28
1
$begingroup$
Yup, for your test cases,f>1e-8
works as a loop condition: 229 bytes.a->{for(double f=1,D,A[],l;f>1e-8;a=A){D=a[0];A=new double[]{f=0,1,0,0};for(var d:a){f+=Math.abs(d-D);A[0]+=d;A[1]*=d;A[2]+=1/d;A[3]+=d*d;}A[0]/=l=a.length;A[1]=Math.pow(A[1],1/l);A[2]=l/A[2];A[3]=Math.sqrt(A[3]/l);}return a[0];}
. With1e-9
, it runs slower (about twice the CPU time), having to do more iterations to get essentially 4 *d-D
down that small. With1e-7
, it's about the same speed as 1e-8. With1e-6
, some of the trailing digits differ for one case.
$endgroup$
– Peter Cordes
Apr 2 at 2:36
1
$begingroup$
@RobinRyder's answer points out that the quadratic mean is always the largest, and harmonic always the smallest, so perhaps you can ditchf
entirely and only checka[3]-a[2]<4e-9
.
$endgroup$
– Peter Cordes
Apr 2 at 7:32
1
$begingroup$
@PeterCordesl==2||
you mean (golfed tol<3|
). But yes, good point; I've added it. :)
$endgroup$
– Kevin Cruijssen
Apr 2 at 8:03
2
$begingroup$
180 bytes by aggregating aggregable reducers.
$endgroup$
– Olivier Grégoire
Apr 2 at 8:45
|
show 12 more comments
$begingroup$
Charcoal, 40 bytes
W‹⌊θ⌈θ≔⟦∕ΣθLθXΠθ∕¹Lθ∕LθΣ∕¹θ₂∕ΣXθ²Lθ⟧θI⊟θ
Try it online! Link is to verbose version of code. Takes input as an array of numbers. Explanation:
W‹⌊θ⌈θ
Repeat while the array contains different values...
≔⟦....⟧θ
... replace the array with a list of values:
∕ΣθLθ
... the mean...
XΠθ∕¹Lθ
... the geometric mean...
∕LθΣ∕¹θ
... the harmonic mean...
₂∕ΣXθ²Lθ
... and the root mean square.
I⊟θ
Cast an element of the array to string and implicitly print it.
$endgroup$
add a comment |
$begingroup$
Jelly, 24 bytes
;µ*Ɱ-;ؽ¤©Æm*İɼ;P½½ƊµÐLḢ
Try it online!
$endgroup$
add a comment |
$begingroup$
PowerShell, 182 180 183 bytes
$f={$a=$c=$d=$n=0
$b=1
$m=[math]
$args|%{$n++
$a+=$_
$b*=$_
$c+=1/$_
$d+=+$_*$_}
$p=($a/$n),$m::pow($b,1/$n),($n/$c),$m::sqrt($d/$n)|%{$m::round($_,9)}
if($p-ne$p[0]){$p=&$f @p}$p[0]}
Try it online!
$endgroup$
$begingroup$
171 bytes
$endgroup$
– mazzy
Apr 1 at 12:10
$begingroup$
@mazzy, impressively :)
$endgroup$
– Andrei Odegov
Apr 1 at 14:14
add a comment |
$begingroup$
Jelly, 25 24 bytes
Wẋ4¹ÆlÆeƭ²½ƭİ4ƭÆm$€⁺µÐLḢ
Try it online!
Explanation
µÐL | Repeat until unchanged:
W | Wrap as a list
ẋ4 | Copy list 4 times
⁺ | Do twice:
$€ | For each copy of the list:
4ƭ | One of these functions, cycling between them:
¹ | Identity
ÆlÆeƭ | Alternate between log and exp
²½ƭ | Alternate between square and square root
İ | Reciprocal
Æm | Then take the mean
Ḣ| Finally take the first item
$endgroup$
$begingroup$
I am fairly bad at Jelly, but could something similar toP*İL
work for the geometric mean?
$endgroup$
– someone
Mar 31 at 15:08
$begingroup$
@someone it would need to beP*Lİ$
so wouldn’t save bytes. It would mean I could bringÆm
back down a line without costing bytes, but I quite like the fact each one currently has an arithmetic mean at its core.
$endgroup$
– Nick Kennedy
Mar 31 at 15:19
add a comment |
$begingroup$
05AB1E (legacy), 26 24 bytes
Δ©ÅA®PYzmY®zO/®nÅAt)4V}н
Try it online or see the steps of all test cases.
Explanation:
Δ # Loop until the list no longer changes:
© # Store the current list in the register (without popping)
# (which is the implicit input-list in the first iteration)
# Arithmetic mean:
ÅA # Builtin to calculate the arithmetic mean of this list
# Geometric mean:
®P # Take the product of the list from the register
Yz # And 1/`Y` (NOTE: `Y` is 2 by default)
m # Then take the result of the product to the power of 1/`Y`
# (NOTE: a^(1/4) is the same as ⁴√a)
# Harmonic mean:
®z # Get 1/x for each value x in the list from the register
O # Take the sum of that
Y / # And then divide `Y` by that sum
# Quadratic mean:
®n # Take 2 to the power of each number x in the list from the register
ÅA # Calculate the arithmetic mean of this list
t # And then take the square-root of that
) # Wrap all four results into a list
4V # And set variable `Y` now to 4 instead of 2
}н # After the list no longer changes: pop and push its first value
# (which is output implicitly as result)
$endgroup$
add a comment |
$begingroup$
C#, 173 bytes
double m(int n,params double[]a)=>(n<1?a[0]:m(n-1,a.Sum()/a.Length,Math.Pow(a.Aggregate((t,x)=>t*x),1.0/a.Length),a.Length/a.Sum(x=>1/x),Math.Sqrt(a.Sum(x=>x*x)/a.Length)));
Try it online!
$endgroup$
2
$begingroup$
This seems to really on a variable that must be passed. Also, you have to includeusing System
andusing System.Linq
in your byte count, since they are required for the program to run. You can change your compiler to the C# Visual Interactive Compiler, which doesn't need those imports. Also,1.0
->1d
$endgroup$
– Embodiment of Ignorance
Mar 31 at 20:59
add a comment |
$begingroup$
Clean, 124 bytes
import StdEnv
f=avg o limit o iteratel=let n=toReal(length l)in[avg l,prod l^(1.0/n),n/sum[1.0/x\x<-l],avg[x*x\x<-l]^0.5]
Try it online!
Performs the operation until the result stops changing.
Hurray for limited precision floating point!
$endgroup$
add a comment |
$begingroup$
Pyth, 32 bytes
h.Wt{H[.OZ@*FZJlZcJscL1Z@.O^R2Z2
Try it online here, or verify all the test cases (bar two, see note below) at once here. Accepts input as a list.
There seem to be some issues with rounding, as certain inputs don't converge correctly when they otherwise should. In particular, test case 0.01 100
gets stuck at values [6.748305820749738, 6.748305820749738, 6.748305820749739, 6.748305820749738]
, and test case 1.61 2.41
gets stuck at [1.9896543776640825, 1.9896543776640825, 1.9896543776640827, 1.9896543776640825]
- note in both cases that the 3rd mean (harmonic mean) differs from the others.
I'm not sure if this problem invalidates my entry, but I'm posting it anyway as it should work. If this isn't acceptable, it can be fixed by instering .RRT
before the [
, to round each of the means to 10 decimal places, as seen in this test suite.
h.Wt{H[.OZ@*FZJlZcJscL1Z@.O^R2Z2)Q Implicit: Q=eval(input())
Trailing )Q inferred
.W Q Funcitonal while: While condition is true, call inner. Starting value Q
t{H Condition function: current input H
{H Deduplicate H
t Discard first value
Empty list is falsey, so while is terminated when means converge
[.OZ@*FZJlZcJscL1Z@.O^R2Z2) Inner function: current input Z
JlZ Take length of Z, store in J
.OZ (1) Arithmetic mean of Z
*FZ Product of Z
@ J (2) Jth root of the above
L Z Map each element of Z...
c 1 ... to its reciprocal
s Sum the above
cJ (3) J / the above
R Z Map each element of Z...
^ 2 ... to its square
.O Arithmetic mean of the above
@ 2 (4) Square root of the above
[ ) Wrap results (1), (2), (3), and (4) in a list
This is used as the input for the next iteration of the loop
h Take the first element of the result, implicit print
$endgroup$
add a comment |
$begingroup$
Japt v2.0a0 -g
, 42 38 bytes
â ÊÉ?ß[Ux²÷(V=UÊ)¬Ux÷V U×qV V÷Ux!÷1]:U
There's got to be a shorter way... This is a monstrosity! Saved 4 bytes thanks to @Shaggy!
Try it
$endgroup$
$begingroup$
38 bytes. But I agree, there must be a shorter way!
$endgroup$
– Shaggy
Apr 1 at 17:12
add a comment |
$begingroup$
C# (Visual C# Interactive Compiler), 177 bytes
double f(double[]g)=>g.All(c=>Math.Abs(c-g[0])<1e-9)?g[0]:f(new[]{g.Sum()/(z=g.Length),Math.Pow(g.Aggregate((a,b)=>a*b),1d/z),z/g.Sum(x=>1/x),Math.Sqrt(g.Sum(x=>x*x)/z)});int z;
Thanks to @KevinCruijjsen for pointing out that using floating point precision was causing problems! Would be 163 bytes if doubles were perfectly precise.
Try it online!
$endgroup$
$begingroup$
The last two test cases give aStackOverflowException
due to floating point precision. Instead ofc==g[0]
you could do something likeMath.Abs(c-g[0])<1e-9
. Try it online.
$endgroup$
– Kevin Cruijssen
Apr 1 at 9:01
$begingroup$
@KevinCruijssen Thanks, it's such a pain dealing with floating point numbers
$endgroup$
– Embodiment of Ignorance
Apr 2 at 15:32
add a comment |
$begingroup$
x86 machine code (SIMD 4x float using 128-bit SSE1&AVX) 94 bytes
x86 machine code (SIMD 4x double using 256-bit AVX) 123 bytes
float
passes the test cases in the question, but with a loop-exit threshold small enough to make that happen, it's easy for it to get stuck in an infinite loop with random inputs.
SSE1 packed-single-precision instructions are 3 bytes long, but SSE2 and simple AVX instructions are 4 bytes long. (Scalar-single instructions like sqrtss
are also 4 bytes long, which is why I use sqrtps
even though I only care about the low element. It's not even slower than sqrtss on modern hardware). I used AVX for non-destructive destination to save 2 bytes vs. movaps+op.
In the double version we can still do a couple movlhps
to copy 64-bit chunks (because often we only care about the low element of a horizontal sum). Horizontal sum of a 256-bit SIMD vector also requires an extra vextractf128
to get the high half, vs. the slow but small 2x haddps
strategy for float. The double
version also needs 2x 8-byte constants, instead of 2x 4-byte. Overall it comes out at close to 4/3 the size of the float
version.
mean(a,b) = mean(a,a,b,b)
for all 4 of these means, so we can simply duplicate the input up to 4 elements and never have to implement length=2. Thus we can hardcode geometric mean as 4th-root = sqrt(sqrt), for example. And we only need one FP constant, 4.0
.
We have a single SIMD vector of all 4 [a_i, b_i, c_i, d_i]
. From that, we calculate the 4 means as scalars in separate registers, and shuffle them back together for the next iteration. (Horizontal operations on SIMD vectors are inconvenient, but we need to do the same thing for all 4 elements in enough cases that it balances out. I started on an x87 version of this, but it was getting very long and not fun.)
The loop-exit condition of }while(quadratic - harmonic > 4e-5)
(or a smaller constant for double
) is based on @RobinRyder's R answer, and Kevin Cruijssen's Java answer: quadratic mean is always the largest magnitude, and harmonic mean is always the smallest (ignoring rounding errors). So we can check the delta between those two to detect convergence. We return the arithmetic mean as the scalar result. It's usually between those two and is probably the least susceptible to rounding errors.
Float version: callable as float meanmean_float_avx(__m128);
with the arg and return value in xmm0. (So x86-64 System V, or Windows x64 vectorcall, but not x64 fastcall.) Or declare the return-type as __m128
so you can get at the quadratic and harmonic mean for testing.
Letting this take 2 separate float
args in xmm0 and xmm1 would cost 1 extra byte: we'd need a shufps
with an imm8 (instead of just unpcklps xmm0,xmm0
) to shuffle together and duplicate 2 inputs.
40 address align 32
41 code bytes global meanmean_float_avx
42 meanmean_float_avx:
43 00000000 B9[52000000] mov ecx, .arith_mean ; allows 2-byte call reg, and a base for loading constants
44 00000005 C4E2791861FC vbroadcastss xmm4, [rcx-4] ; float 4.0
45
46 ;; mean(a,b) = mean(a,b,a,b) for all 4 types of mean
47 ;; so we only ever have to do the length=4 case
48 0000000B 0F14C0 unpcklps xmm0,xmm0 ; [b,a] => [b,b,a,a]
49
50 ; do{ ... } while(quadratic - harmonic > threshold);
51 .loop:
52 ;;; XMM3 = geometric mean: not based on addition. (Transform to log would be hard. AVX512ER has exp with 23-bit accuracy, but not log. vgetexp = floor(lofg2(x)), so that's no good.)
53 ;; sqrt once *first*, making magnitudes closer to 1.0 to reduce rounding error. Numbers are all positive so this is safe.
54 ;; both sqrts first was better behaved, I think.
55 0000000E 0F51D8 sqrtps xmm3, xmm0 ; xmm3 = 4th root(x)
56 00000011 F30F16EB movshdup xmm5, xmm3 ; bring odd elements down to even
57 00000015 0F59EB mulps xmm5, xmm3
58 00000018 0F12DD movhlps xmm3, xmm5 ; high half -> low
59 0000001B 0F59DD mulps xmm3, xmm5 ; xmm3[0] = hproduct(sqrt(xmm))
60 ; sqrtps xmm3, xmm3 ; sqrt(hprod(sqrt)) = 4th root(hprod)
61 ; common final step done after interleaving with quadratic mean
62
63 ;;; XMM2 = quadratic mean = max of the means
64 0000001E C5F859E8 vmulps xmm5, xmm0,xmm0
65 00000022 FFD1 call rcx ; arith mean of squares
66 00000024 0F14EB unpcklps xmm5, xmm3 ; [quad^2, geo^2, ?, ?]
67 00000027 0F51D5 sqrtps xmm2, xmm5 ; [quad, geo, ?, ?]
68
69 ;;; XMM1 = harmonic mean = min of the means
70 0000002A C5D85EE8 vdivps xmm5, xmm4, xmm0 ; 4/x
71 0000002E FFD1 call rcx ; arithmetic mean (under inversion)
72 00000030 C5D85ECD vdivps xmm1, xmm4, xmm5 ; 4/. (the factor of 4 cancels out)
73
74 ;;; XMM5 = arithmetic mean
75 00000034 0F28E8 movaps xmm5, xmm0
76 00000037 FFD1 call rcx
77
78 00000039 0F14E9 unpcklps xmm5, xmm1 ; [arith, harm, ?,?]
79 0000003C C5D014C2 vunpcklps xmm0, xmm5,xmm2 ; x = [arith, harm, quad, geo]
80
81 00000040 0F5CD1 subps xmm2, xmm1 ; largest - smallest mean: guaranteed non-negative
82 00000043 0F2E51F8 ucomiss xmm2, [rcx-8] ; quad-harm > convergence_threshold
83 00000047 73C5 jae .loop
84
85 ; return with the arithmetic mean in the low element of xmm0 = scalar return value
86 00000049 C3 ret
87
88 ;;; "constant pool" between the main function and the helper, like ARM literal pools
89 0000004A ACC52738 .fpconst_threshold: dd 4e-5 ; 4.3e-5 is the highest we can go and still pass the main test cases
90 0000004E 00008040 .fpconst_4: dd 4.0
91 .arith_mean: ; returns XMM5 = hsum(xmm5)/4.
92 00000052 C5D37CED vhaddps xmm5, xmm5 ; slow but small
93 00000056 C5D37CED vhaddps xmm5, xmm5
94 0000005A 0F5EEC divps xmm5, xmm4 ; divide before/after summing doesn't matter mathematically or numerically; divisor is a power of 2
95 0000005D C3 ret
96 0000005E 5E000000 .size: dd $ - meanmean_float_avx
0x5e = 94 bytes
(NASM listing created with nasm -felf64 mean-mean.asm -l/dev/stdout | cut -b -34,$((34+6))-
. Strip the listing part and recover the source with cut -b 34- > mean-mean.asm
)
SIMD horizontal sum and divide by 4 (i.e. arithmetic mean) is implemented in a separate function that we call
(with a function pointer to amortize the cost of the address). With 4/x
before/after, or x^2
before and sqrt after, we get the harmonic mean and quadratic mean. (It was painful to write these div
instructions instead of multiplying by an exactly-representable 0.25
.)
Geometric mean is implemented separately with multiply and chained sqrt. Or with one sqrt first to reduce exponent magnitude and maybe help numerical precision. log is not available, only floor(log2(x))
via AVX512 vgetexpps/pd
. Exp is sort of available via AVX512ER (Xeon Phi only), but with only 2^-23 precision.
Mixing 128-bit AVX instructions and legacy SSE is not a performance problem. Mixing 256-bit AVX with legacy SSE can be on Haswell, but on Skylake it just potentially creates a potential false dependency for SSE instructions. I think my double
version avoids any unnecessary loop-carried dep chains, and bottlenecks
on div/sqrt latency/throughput.
Double version:
108 global meanmean_double_avx
109 meanmean_double_avx:
110 00000080 B9[E8000000] mov ecx, .arith_mean
111 00000085 C4E27D1961F8 vbroadcastsd ymm4, [rcx-8] ; float 4.0
112
113 ;; mean(a,b) = mean(a,b,a,b) for all 4 types of mean
114 ;; so we only ever have to do the length=4 case
115 0000008B C4E37D18C001 vinsertf128 ymm0, ymm0, xmm0, 1 ; [b,a] => [b,a,b,a]
116
117 .loop:
118 ;;; XMM3 = geometric mean: not based on addition.
119 00000091 C5FD51D8 vsqrtpd ymm3, ymm0 ; sqrt first to get magnitude closer to 1.0 for better(?) numerical precision
120 00000095 C4E37D19DD01 vextractf128 xmm5, ymm3, 1 ; extract high lane
121 0000009B C5D159EB vmulpd xmm5, xmm3
122 0000009F 0F12DD movhlps xmm3, xmm5 ; extract high half
123 000000A2 F20F59DD mulsd xmm3, xmm5 ; xmm3 = hproduct(sqrt(xmm0))
124 ; sqrtsd xmm3, xmm3 ; xmm3 = 4th root = geomean(xmm0) ;deferred until quadratic
125
126 ;;; XMM2 = quadratic mean = max of the means
127 000000A6 C5FD59E8 vmulpd ymm5, ymm0,ymm0
128 000000AA FFD1 call rcx ; arith mean of squares
129 000000AC 0F16EB movlhps xmm5, xmm3 ; [quad^2, geo^2]
130 000000AF 660F51D5 sqrtpd xmm2, xmm5 ; [quad , geo]
131
132 ;;; XMM1 = harmonic mean = min of the means
133 000000B3 C5DD5EE8 vdivpd ymm5, ymm4, ymm0 ; 4/x
134 000000B7 FFD1 call rcx ; arithmetic mean under inversion
135 000000B9 C5DB5ECD vdivsd xmm1, xmm4, xmm5 ; 4/. (the factor of 4 cancels out)
136
137 ;;; XMM5 = arithmetic mean
138 000000BD C5FC28E8 vmovaps ymm5, ymm0
139 000000C1 FFD1 call rcx
140
141 000000C3 0F16E9 movlhps xmm5, xmm1 ; [arith, harm]
142 000000C6 C4E35518C201 vinsertf128 ymm0, ymm5, xmm2, 1 ; x = [arith, harm, quad, geo]
143
144 000000CC C5EB5CD1 vsubsd xmm2, xmm1 ; largest - smallest mean: guaranteed non-negative
145 000000D0 660F2E51F0 ucomisd xmm2, [rcx-16] ; quad - harm > threshold
146 000000D5 77BA ja .loop
147
148 ; vzeroupper ; not needed for correctness, only performance
149 ; return with the arithmetic mean in the low element of xmm0 = scalar return value
150 000000D7 C3 ret
151
152 ; "literal pool" between the function
153 000000D8 95D626E80B2E113E .fpconst_threshold: dq 1e-9
154 000000E0 0000000000001040 .fpconst_4: dq 4.0 ; TODO: golf these zeros? vpbroadcastb and convert?
155 .arith_mean: ; returns YMM5 = hsum(ymm5)/4.
156 000000E8 C4E37D19EF01 vextractf128 xmm7, ymm5, 1
157 000000EE C5D158EF vaddpd xmm5, xmm7
158 000000F2 C5D17CED vhaddpd xmm5, xmm5 ; slow but small
159 000000F6 C5D35EEC vdivsd xmm5, xmm4 ; only low element matters
160 000000FA C3 ret
161 000000FB 7B000000 .size: dd $ - meanmean_double_avx
0x7b = 123 bytes
C test harness
#include <immintrin.h>
#include <stdio.h>
#include <math.h>
static const struct ab_avg {
double a,b;
double mean;
} testcases[] = {
{1, 1, 1},
{1, 2, 1.45568889},
{100, 200, 145.568889},
{2.71, 3.14, 2.92103713},
{0.57, 1.78, 1.0848205},
{1.61, 2.41, 1.98965438},
{0.01, 100, 6.7483058},
};
// see asm comments for order of arith, harm, quad, geo
__m128 meanmean_float_avx(__m128); // or float ...
__m256d meanmean_double_avx(__m128d); // or double ...
int main(void) {
int len = sizeof(testcases) / sizeof(testcases[0]);
for(int i=0 ; i<len ; i++) {
const struct ab_avg *p = &testcases[i];
#if 1
__m128 arg = _mm_set_ps(0,0, p->b, p->a);
double res = meanmean_float_avx(arg)[0];
#else
__m128d arg = _mm_loadu_pd(&p->a);
double res = meanmean_double_avx(arg)[0];
#endif
double allowed_diff = (p->b - p->a) / 100000.0;
double delta = fabs(p->mean - res);
if (delta > 1e-3 || delta > allowed_diff) {
printf("%f %f => %.9f but we got %.9f. delta = %g allowed=%gn",
p->a, p->b, p->mean, res, p->mean - res, allowed_diff);
}
}
while(1) {
double a = drand48(), b = drand48(); // range= [0..1)
if (a>b) {
double tmp=a;
a=b;
b=tmp; // sorted
}
// a *= 0.00000001;
// b *= 123156;
// a += 1<<11; b += (1<<12)+1; // float version gets stuck inflooping on 2048.04, 4097.18 at fpthreshold = 4e-5
// a *= 1<<11 ; b *= 1<<11; // scaling to large magnitude makes sum of squares loses more precision
//a += 1<<11; b+= 1<<11; // adding to large magnitude is hard for everything, catastrophic cancellation
#if 1
printf("testing float %g, %gn", a, b);
__m128 arg = _mm_set_ps(0,0, b, a);
__m128 res = meanmean_float_avx(arg);
double quad = res[2], harm = res[1]; // same order as double... for now
#else
printf("testing double %g, %gn", a, b);
__m128d arg = _mm_set_pd(b, a);
__m256d res = meanmean_double_avx(arg);
double quad = res[2], harm = res[1];
#endif
double delta = fabs(quad - harm);
double allowed_diff = (b - a) / 100000.0; // calculated in double even for the float case.
// TODO: use the double res as a reference for float res
// instead of just checking quadratic vs. harmonic mean
if (delta > 1e-3 || delta > allowed_diff) {
printf("%g %g we got q=%g, h=%g, a=%g. delta = %g, allowed=%gn",
a, b, quad, harm, res[0], quad-harm, allowed_diff);
}
}
}
Build with:
nasm -felf64 mean-mean.asm &&
gcc -no-pie -fno-pie -g -O2 -march=native mean-mean.c mean-mean.o
Obviously you need a CPU with AVX support, or an emulator like Intel SDE. To compile on a host without native AVX support, use -march=sandybridge
or -mavx
Run: passes the hard-coded test cases, but for the float version, random test cases often fail the (b-a)/10000
threshold set in the question.
$ ./a.out
(note: empty output before the first "testing float" means clean pass on the constant test cases)
testing float 3.90799e-14, 0.000985395
3.90799e-14 0.000985395 we got q=3.20062e-10, h=3.58723e-05, a=2.50934e-05. delta = -3.5872e-05, allowed=9.85395e-09
testing float 0.041631, 0.176643
testing float 0.0913306, 0.364602
testing float 0.0922976, 0.487217
testing float 0.454433, 0.52675
0.454433 0.52675 we got q=0.48992, h=0.489927, a=0.489925. delta = -6.79493e-06, allowed=7.23169e-07
testing float 0.233178, 0.831292
testing float 0.56806, 0.931731
testing float 0.0508319, 0.556094
testing float 0.0189148, 0.767051
0.0189148 0.767051 we got q=0.210471, h=0.210484, a=0.21048. delta = -1.37389e-05, allowed=7.48136e-06
testing float 0.25236, 0.298197
0.25236 0.298197 we got q=0.274796, h=0.274803, a=0.274801. delta = -6.19888e-06, allowed=4.58374e-07
testing float 0.531557, 0.875981
testing float 0.515431, 0.920261
testing float 0.18842, 0.810429
testing float 0.570614, 0.886314
testing float 0.0767746, 0.815274
testing float 0.118352, 0.984891
0.118352 0.984891 we got q=0.427845, h=0.427872, a=0.427863. delta = -2.66135e-05, allowed=8.66539e-06
testing float 0.784484, 0.893906
0.784484 0.893906 we got q=0.838297, h=0.838304, a=0.838302. delta = -7.09295e-06, allowed=1.09422e-06
FP errors are enough that quad-harm comes out less than zero for some inputs.
Or with a += 1<<11; b += (1<<12)+1;
uncommented:
testing float 2048, 4097
testing float 2048.04, 4097.18
^C (stuck in an infinite loop).
None of these problems happen with double
. Comment out the printf
before each test to see that the output is empty (nothing from the if(delta too high)
block).
TODO: use the double
version as a reference for the float
version, instead of just looking at how they converging with quad-harm.
$endgroup$
add a comment |
$begingroup$
Javascript - 186 bytes
Takes input as an array of numbers. Uses the mean transformations in J42161217's answer to shorten the code.
Try It Online
f=(v,l=[m=(w,k=0)=>w.map(x=>k+=x)&&k/w.length,w=>1/m(w.map(x=>1/x)),w=>Math.E**m(w.map(x=>Math.log(x))),w=>m(w.map(x=>x**2))**.5].map(x=>x(v)).sort((a,b)=>a-b))=>l[3]-l[0]>1e-5?f(l):l[0]
Explanation
f = (
v,
l=[
m=(w,k=0)=>w.map(x=>k+=x)&&k/w.length, // m = w => arithmetic mean of values in w
w=>1/m(w.map(x=>1/x)), // w => harmonic mean of values in w
w=>Math.E**m(w.map(x=>Math.log(x))), // w => geometric mean of values in w
w=>m(w.map(x=>x**2))**.5 // w => quadratic mean of values in w
].map(x=>x(v)) // get array of each mean using input v, stored in l
.sort((a,b)=>a-b) // sort the outputs
) =>
l[3] - l[0] > 1e-5 ? // is the difference between the largest
// and smallest means > 1/100000?
f(l) : // if yes, get the mean mean of the means
l[0] // if no, arbitrarily return the smallest value
// as close enough
$endgroup$
$begingroup$
I thought I was going to be clever and implement the relationship with logarithms but it looks like you and J42161217 got there first!
$endgroup$
– Pureferret
Apr 3 at 11:08
$begingroup$
@Pureferret I take no credit for that,I blatantly stole it :D
$endgroup$
– asgallant
Apr 3 at 15:58
$begingroup$
you wrote it in JavaScript though!
$endgroup$
– Pureferret
Apr 3 at 15:59
1
$begingroup$
That was the easy part. Golfing it was hard.
$endgroup$
– asgallant
Apr 3 at 16:13
1
$begingroup$
The TIL was not configured correctly. I added a TIL link to the answer.
$endgroup$
– asgallant
Apr 4 at 16:27
|
show 2 more comments
$begingroup$
Perl 5, 92 72 bytes
sub f{@_=map{$m=$_;sum(map$_**$m/@_,@_)**(1/$m)}1,1e-9,-1,2for 1..9;pop}
Try it online!
...using some dirty tricks.
$endgroup$
add a comment |
$begingroup$
SNOBOL4 (CSNOBOL4), 296 bytes
X =INPUT
Y =INPUT
A =(X + Y) / 2.
P =X * Y
G =P ^ 0.5
H =P / A
Q =(2 * (A ^ 2) - P) ^ 0.5
O OUTPUT =EQ(Q,A) Q :S(END)
M =A
N =G
O =H
P =Q
A =(M + N + O + P) / 4
G =(M * N * O * P) ^ 0.25
H =4 / (1 / M + 1 / N + 1 / O + 1 / P)
Q =((M ^ 2 + N ^ 2 + O ^ 2 + P ^ 2) / 4) ^ 0.5 :(O)
END
Try it online!
Straightforward implementation. Uses a trick from my answer to a related question to golf a bit more.
$endgroup$
add a comment |
Your Answer
StackExchange.ifUsing("editor", function () {
return StackExchange.using("mathjaxEditing", function () {
StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix) {
StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["\$", "\$"]]);
});
});
}, "mathjax-editing");
StackExchange.ifUsing("editor", function () {
StackExchange.using("externalEditor", function () {
StackExchange.using("snippets", function () {
StackExchange.snippets.init();
});
});
}, "code-snippets");
StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "200"
};
initTagRenderer("".split(" "), "".split(" "), channelOptions);
StackExchange.using("externalEditor", function() {
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled) {
StackExchange.using("snippets", function() {
createEditor();
});
}
else {
createEditor();
}
});
function createEditor() {
StackExchange.prepareEditor({
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: false,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: null,
bindNavPrevention: true,
postfix: "",
imageUploader: {
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
},
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});
}
});
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fcodegolf.stackexchange.com%2fquestions%2f182460%2fcalculate-the-mean-mean-of-two-numbers%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
19 Answers
19
active
oldest
votes
19 Answers
19
active
oldest
votes
active
oldest
votes
active
oldest
votes
$begingroup$
Wolfram Language (Mathematica), 52 bytes
#//.x_:>N@{M@x,E^M@Log@x,1/M[1/x],M[x^2]^.5}&
M=Mean
Try it online!
In my first approach I used these builtinsMean
GeometricMean
HarmonicMean
and RootMeanSquare
Here are some substitutions for saving bytes
HarmonicMean
--> 1/Mean[1/x]
by @Robin Ryder (3 bytes saved)GeometricMean
--> E^Mean@Log@x
by @A. Rex (2 bytes saved)RootMeanSquare
--> Mean[x^2]^.5
by @A. Rex (4 bytes saved)
finally we can assign Mean
to M
(as proposed by @ovs) and save 5 more bytes
$endgroup$
$begingroup$
Save 2 bytes by recoding GeometricMean
$endgroup$
– Robin Ryder
Apr 1 at 9:58
$begingroup$
@RobinRyder I believe you mean Harmonic.. nice!
$endgroup$
– J42161217
Apr 1 at 10:04
1
$begingroup$
Save 8 more bytes:#//.x_:>N@{Mean@x,E^Mean@Log@x,1/Mean[1/x],Mean[x^2]^.5}&
$endgroup$
– A. Rex
Apr 1 at 10:41
$begingroup$
@ovs edited.....
$endgroup$
– J42161217
Apr 1 at 11:27
add a comment |
$begingroup$
Wolfram Language (Mathematica), 52 bytes
#//.x_:>N@{M@x,E^M@Log@x,1/M[1/x],M[x^2]^.5}&
M=Mean
Try it online!
In my first approach I used these builtinsMean
GeometricMean
HarmonicMean
and RootMeanSquare
Here are some substitutions for saving bytes
HarmonicMean
--> 1/Mean[1/x]
by @Robin Ryder (3 bytes saved)GeometricMean
--> E^Mean@Log@x
by @A. Rex (2 bytes saved)RootMeanSquare
--> Mean[x^2]^.5
by @A. Rex (4 bytes saved)
finally we can assign Mean
to M
(as proposed by @ovs) and save 5 more bytes
$endgroup$
$begingroup$
Save 2 bytes by recoding GeometricMean
$endgroup$
– Robin Ryder
Apr 1 at 9:58
$begingroup$
@RobinRyder I believe you mean Harmonic.. nice!
$endgroup$
– J42161217
Apr 1 at 10:04
1
$begingroup$
Save 8 more bytes:#//.x_:>N@{Mean@x,E^Mean@Log@x,1/Mean[1/x],Mean[x^2]^.5}&
$endgroup$
– A. Rex
Apr 1 at 10:41
$begingroup$
@ovs edited.....
$endgroup$
– J42161217
Apr 1 at 11:27
add a comment |
$begingroup$
Wolfram Language (Mathematica), 52 bytes
#//.x_:>N@{M@x,E^M@Log@x,1/M[1/x],M[x^2]^.5}&
M=Mean
Try it online!
In my first approach I used these builtinsMean
GeometricMean
HarmonicMean
and RootMeanSquare
Here are some substitutions for saving bytes
HarmonicMean
--> 1/Mean[1/x]
by @Robin Ryder (3 bytes saved)GeometricMean
--> E^Mean@Log@x
by @A. Rex (2 bytes saved)RootMeanSquare
--> Mean[x^2]^.5
by @A. Rex (4 bytes saved)
finally we can assign Mean
to M
(as proposed by @ovs) and save 5 more bytes
$endgroup$
Wolfram Language (Mathematica), 52 bytes
#//.x_:>N@{M@x,E^M@Log@x,1/M[1/x],M[x^2]^.5}&
M=Mean
Try it online!
In my first approach I used these builtinsMean
GeometricMean
HarmonicMean
and RootMeanSquare
Here are some substitutions for saving bytes
HarmonicMean
--> 1/Mean[1/x]
by @Robin Ryder (3 bytes saved)GeometricMean
--> E^Mean@Log@x
by @A. Rex (2 bytes saved)RootMeanSquare
--> Mean[x^2]^.5
by @A. Rex (4 bytes saved)
finally we can assign Mean
to M
(as proposed by @ovs) and save 5 more bytes
edited Apr 1 at 11:26
answered Mar 31 at 8:39
J42161217J42161217
13.8k21253
13.8k21253
$begingroup$
Save 2 bytes by recoding GeometricMean
$endgroup$
– Robin Ryder
Apr 1 at 9:58
$begingroup$
@RobinRyder I believe you mean Harmonic.. nice!
$endgroup$
– J42161217
Apr 1 at 10:04
1
$begingroup$
Save 8 more bytes:#//.x_:>N@{Mean@x,E^Mean@Log@x,1/Mean[1/x],Mean[x^2]^.5}&
$endgroup$
– A. Rex
Apr 1 at 10:41
$begingroup$
@ovs edited.....
$endgroup$
– J42161217
Apr 1 at 11:27
add a comment |
$begingroup$
Save 2 bytes by recoding GeometricMean
$endgroup$
– Robin Ryder
Apr 1 at 9:58
$begingroup$
@RobinRyder I believe you mean Harmonic.. nice!
$endgroup$
– J42161217
Apr 1 at 10:04
1
$begingroup$
Save 8 more bytes:#//.x_:>N@{Mean@x,E^Mean@Log@x,1/Mean[1/x],Mean[x^2]^.5}&
$endgroup$
– A. Rex
Apr 1 at 10:41
$begingroup$
@ovs edited.....
$endgroup$
– J42161217
Apr 1 at 11:27
$begingroup$
Save 2 bytes by recoding GeometricMean
$endgroup$
– Robin Ryder
Apr 1 at 9:58
$begingroup$
Save 2 bytes by recoding GeometricMean
$endgroup$
– Robin Ryder
Apr 1 at 9:58
$begingroup$
@RobinRyder I believe you mean Harmonic.. nice!
$endgroup$
– J42161217
Apr 1 at 10:04
$begingroup$
@RobinRyder I believe you mean Harmonic.. nice!
$endgroup$
– J42161217
Apr 1 at 10:04
1
1
$begingroup$
Save 8 more bytes:
#//.x_:>N@{Mean@x,E^Mean@Log@x,1/Mean[1/x],Mean[x^2]^.5}&
$endgroup$
– A. Rex
Apr 1 at 10:41
$begingroup$
Save 8 more bytes:
#//.x_:>N@{Mean@x,E^Mean@Log@x,1/Mean[1/x],Mean[x^2]^.5}&
$endgroup$
– A. Rex
Apr 1 at 10:41
$begingroup$
@ovs edited.....
$endgroup$
– J42161217
Apr 1 at 11:27
$begingroup$
@ovs edited.....
$endgroup$
– J42161217
Apr 1 at 11:27
add a comment |
$begingroup$
R, 70 69 67 bytes
x=scan();`?`=mean;while(x-?x)x=c((?x^2)^.5,?x,2^?log2(x),1/?1/x);?x
Try it online!
-1 byte with better conditioning.
-2 bytes by switching to base 2.
Like some other answers, this uses the expression of the geometric mean as an arithmetic mean on the log scale (here in base 2):
$$M_0(x_1,ldots, x_n) = 2^{M_1(log_2 x_1,ldots,log_2 x_n)}.$$
It also uses the fact that $forall k, d_kgeq a_k geq b_kgeq c_k$, i.e. $d_k=max(a_k, b_k, c_k, d_k)$. The condition $a_k=b_k=c_k=d_k$ is therefore equivalent to $d_k = M_1(a_k, b_k, c_k, d_k)$, which is what I use in the while loop; this is achieved by abusing the syntax of while
which only considers the first element when the condition is a vector, hence the order in which the means are stored. (Note that we could also use $c_k$ instead since it is the minimum of the four, but we could not use $a_k$ or $b_k$ in the condition.)
When we exit the while loop, x
is a constant vector. The final ?x
computes its mean to reduce it to a scalar.
$endgroup$
1
$begingroup$
Shouldn't it be $ln x_n$ instead of $log x_n$ ?
$endgroup$
– Tau
Mar 31 at 23:40
$begingroup$
@Tau Yes, I was denoting natural logarithm by $ log$, which is the default in R. Anyway, I have now changed it to base 2 logarithm for -2 bytes.
$endgroup$
– Robin Ryder
Apr 1 at 7:14
add a comment |
$begingroup$
R, 70 69 67 bytes
x=scan();`?`=mean;while(x-?x)x=c((?x^2)^.5,?x,2^?log2(x),1/?1/x);?x
Try it online!
-1 byte with better conditioning.
-2 bytes by switching to base 2.
Like some other answers, this uses the expression of the geometric mean as an arithmetic mean on the log scale (here in base 2):
$$M_0(x_1,ldots, x_n) = 2^{M_1(log_2 x_1,ldots,log_2 x_n)}.$$
It also uses the fact that $forall k, d_kgeq a_k geq b_kgeq c_k$, i.e. $d_k=max(a_k, b_k, c_k, d_k)$. The condition $a_k=b_k=c_k=d_k$ is therefore equivalent to $d_k = M_1(a_k, b_k, c_k, d_k)$, which is what I use in the while loop; this is achieved by abusing the syntax of while
which only considers the first element when the condition is a vector, hence the order in which the means are stored. (Note that we could also use $c_k$ instead since it is the minimum of the four, but we could not use $a_k$ or $b_k$ in the condition.)
When we exit the while loop, x
is a constant vector. The final ?x
computes its mean to reduce it to a scalar.
$endgroup$
1
$begingroup$
Shouldn't it be $ln x_n$ instead of $log x_n$ ?
$endgroup$
– Tau
Mar 31 at 23:40
$begingroup$
@Tau Yes, I was denoting natural logarithm by $ log$, which is the default in R. Anyway, I have now changed it to base 2 logarithm for -2 bytes.
$endgroup$
– Robin Ryder
Apr 1 at 7:14
add a comment |
$begingroup$
R, 70 69 67 bytes
x=scan();`?`=mean;while(x-?x)x=c((?x^2)^.5,?x,2^?log2(x),1/?1/x);?x
Try it online!
-1 byte with better conditioning.
-2 bytes by switching to base 2.
Like some other answers, this uses the expression of the geometric mean as an arithmetic mean on the log scale (here in base 2):
$$M_0(x_1,ldots, x_n) = 2^{M_1(log_2 x_1,ldots,log_2 x_n)}.$$
It also uses the fact that $forall k, d_kgeq a_k geq b_kgeq c_k$, i.e. $d_k=max(a_k, b_k, c_k, d_k)$. The condition $a_k=b_k=c_k=d_k$ is therefore equivalent to $d_k = M_1(a_k, b_k, c_k, d_k)$, which is what I use in the while loop; this is achieved by abusing the syntax of while
which only considers the first element when the condition is a vector, hence the order in which the means are stored. (Note that we could also use $c_k$ instead since it is the minimum of the four, but we could not use $a_k$ or $b_k$ in the condition.)
When we exit the while loop, x
is a constant vector. The final ?x
computes its mean to reduce it to a scalar.
$endgroup$
R, 70 69 67 bytes
x=scan();`?`=mean;while(x-?x)x=c((?x^2)^.5,?x,2^?log2(x),1/?1/x);?x
Try it online!
-1 byte with better conditioning.
-2 bytes by switching to base 2.
Like some other answers, this uses the expression of the geometric mean as an arithmetic mean on the log scale (here in base 2):
$$M_0(x_1,ldots, x_n) = 2^{M_1(log_2 x_1,ldots,log_2 x_n)}.$$
It also uses the fact that $forall k, d_kgeq a_k geq b_kgeq c_k$, i.e. $d_k=max(a_k, b_k, c_k, d_k)$. The condition $a_k=b_k=c_k=d_k$ is therefore equivalent to $d_k = M_1(a_k, b_k, c_k, d_k)$, which is what I use in the while loop; this is achieved by abusing the syntax of while
which only considers the first element when the condition is a vector, hence the order in which the means are stored. (Note that we could also use $c_k$ instead since it is the minimum of the four, but we could not use $a_k$ or $b_k$ in the condition.)
When we exit the while loop, x
is a constant vector. The final ?x
computes its mean to reduce it to a scalar.
edited Apr 2 at 7:55
answered Mar 31 at 19:35
Robin RyderRobin Ryder
5417
5417
1
$begingroup$
Shouldn't it be $ln x_n$ instead of $log x_n$ ?
$endgroup$
– Tau
Mar 31 at 23:40
$begingroup$
@Tau Yes, I was denoting natural logarithm by $ log$, which is the default in R. Anyway, I have now changed it to base 2 logarithm for -2 bytes.
$endgroup$
– Robin Ryder
Apr 1 at 7:14
add a comment |
1
$begingroup$
Shouldn't it be $ln x_n$ instead of $log x_n$ ?
$endgroup$
– Tau
Mar 31 at 23:40
$begingroup$
@Tau Yes, I was denoting natural logarithm by $ log$, which is the default in R. Anyway, I have now changed it to base 2 logarithm for -2 bytes.
$endgroup$
– Robin Ryder
Apr 1 at 7:14
1
1
$begingroup$
Shouldn't it be $ln x_n$ instead of $log x_n$ ?
$endgroup$
– Tau
Mar 31 at 23:40
$begingroup$
Shouldn't it be $ln x_n$ instead of $log x_n$ ?
$endgroup$
– Tau
Mar 31 at 23:40
$begingroup$
@Tau Yes, I was denoting natural logarithm by $ log$, which is the default in R. Anyway, I have now changed it to base 2 logarithm for -2 bytes.
$endgroup$
– Robin Ryder
Apr 1 at 7:14
$begingroup$
@Tau Yes, I was denoting natural logarithm by $ log$, which is the default in R. Anyway, I have now changed it to base 2 logarithm for -2 bytes.
$endgroup$
– Robin Ryder
Apr 1 at 7:14
add a comment |
$begingroup$
J, 34 bytes
(31 as an expression without the assignment to variable f
)
f=:1{(^.z,%z,*:z,[z=:(+/%#)&.:)^:_
Try it online!
For functions a
and b
, a &.: b
("a under b" (related challenge)) is equivalent to (b inv) a b
-- apply b, then a, then inverse of b. In this case, geometric/harmonic/quadratic mean is the arithmetic mean "under" logarithm, inversion, and square respectively.
$endgroup$
add a comment |
$begingroup$
J, 34 bytes
(31 as an expression without the assignment to variable f
)
f=:1{(^.z,%z,*:z,[z=:(+/%#)&.:)^:_
Try it online!
For functions a
and b
, a &.: b
("a under b" (related challenge)) is equivalent to (b inv) a b
-- apply b, then a, then inverse of b. In this case, geometric/harmonic/quadratic mean is the arithmetic mean "under" logarithm, inversion, and square respectively.
$endgroup$
add a comment |
$begingroup$
J, 34 bytes
(31 as an expression without the assignment to variable f
)
f=:1{(^.z,%z,*:z,[z=:(+/%#)&.:)^:_
Try it online!
For functions a
and b
, a &.: b
("a under b" (related challenge)) is equivalent to (b inv) a b
-- apply b, then a, then inverse of b. In this case, geometric/harmonic/quadratic mean is the arithmetic mean "under" logarithm, inversion, and square respectively.
$endgroup$
J, 34 bytes
(31 as an expression without the assignment to variable f
)
f=:1{(^.z,%z,*:z,[z=:(+/%#)&.:)^:_
Try it online!
For functions a
and b
, a &.: b
("a under b" (related challenge)) is equivalent to (b inv) a b
-- apply b, then a, then inverse of b. In this case, geometric/harmonic/quadratic mean is the arithmetic mean "under" logarithm, inversion, and square respectively.
edited Mar 31 at 15:26
answered Mar 31 at 14:00
user202729user202729
14.2k12653
14.2k12653
add a comment |
add a comment |
$begingroup$
TI-BASIC, 42 35 34 bytes
-1 byte thanks to @SolomonUcko
While max(ΔList(Ans:{mean(Ans),√(mean(Ans²)),mean(Ans^-1)^-1,e^(mean(ln(Ans:End:Ans(1
Input is a list of two integers in Ans
.
Output is stored in Ans
and is automatically printed out when the program completes.
Formulas used for geometric, harmonic, and quadratic means are based off of user202729's explanation.
Example:
{1,2
{1 2}
prgmCDGFB
1.455688891
{100,200
{100 200}
prgmCDGFB
145.5688891
Explanation:
(Newlines have been added for clarification. They do NOT appear in the code.)
While max(ΔList(Ans ;loop until all elements of the current list are equal
; the maximum of the change in each element will be 0
{ ;create a list containing...
mean(Ans), ; the arithmetic mean
√(mean(Ans²)), ; the quadratic mean
mean(Ans^-1)^-1, ; the harmonic mean
e^(mean(ln(Ans ; and the geometric mean
End
Ans(1 ;keep the first element in "Ans" and implicitly print it
Notes:
TI-BASIC is a tokenized language. Character count does not equal byte count.
e^(
is this one-byte token.
^-1
is used for this one-byte token.
I opted for writing ^-1
instead because the token looks like ֿ¹
when in a code block.
√(
is this one-byte token.
ΔList(
is this two-byte token.
$endgroup$
$begingroup$
I think you can save a parenthesis by putting the geometric mean last.
$endgroup$
– Solomon Ucko
Apr 1 at 11:06
$begingroup$
@SolomonUcko ah, thanks for noticing! Didn't consider that before.
$endgroup$
– Tau
Apr 1 at 13:01
add a comment |
$begingroup$
TI-BASIC, 42 35 34 bytes
-1 byte thanks to @SolomonUcko
While max(ΔList(Ans:{mean(Ans),√(mean(Ans²)),mean(Ans^-1)^-1,e^(mean(ln(Ans:End:Ans(1
Input is a list of two integers in Ans
.
Output is stored in Ans
and is automatically printed out when the program completes.
Formulas used for geometric, harmonic, and quadratic means are based off of user202729's explanation.
Example:
{1,2
{1 2}
prgmCDGFB
1.455688891
{100,200
{100 200}
prgmCDGFB
145.5688891
Explanation:
(Newlines have been added for clarification. They do NOT appear in the code.)
While max(ΔList(Ans ;loop until all elements of the current list are equal
; the maximum of the change in each element will be 0
{ ;create a list containing...
mean(Ans), ; the arithmetic mean
√(mean(Ans²)), ; the quadratic mean
mean(Ans^-1)^-1, ; the harmonic mean
e^(mean(ln(Ans ; and the geometric mean
End
Ans(1 ;keep the first element in "Ans" and implicitly print it
Notes:
TI-BASIC is a tokenized language. Character count does not equal byte count.
e^(
is this one-byte token.
^-1
is used for this one-byte token.
I opted for writing ^-1
instead because the token looks like ֿ¹
when in a code block.
√(
is this one-byte token.
ΔList(
is this two-byte token.
$endgroup$
$begingroup$
I think you can save a parenthesis by putting the geometric mean last.
$endgroup$
– Solomon Ucko
Apr 1 at 11:06
$begingroup$
@SolomonUcko ah, thanks for noticing! Didn't consider that before.
$endgroup$
– Tau
Apr 1 at 13:01
add a comment |
$begingroup$
TI-BASIC, 42 35 34 bytes
-1 byte thanks to @SolomonUcko
While max(ΔList(Ans:{mean(Ans),√(mean(Ans²)),mean(Ans^-1)^-1,e^(mean(ln(Ans:End:Ans(1
Input is a list of two integers in Ans
.
Output is stored in Ans
and is automatically printed out when the program completes.
Formulas used for geometric, harmonic, and quadratic means are based off of user202729's explanation.
Example:
{1,2
{1 2}
prgmCDGFB
1.455688891
{100,200
{100 200}
prgmCDGFB
145.5688891
Explanation:
(Newlines have been added for clarification. They do NOT appear in the code.)
While max(ΔList(Ans ;loop until all elements of the current list are equal
; the maximum of the change in each element will be 0
{ ;create a list containing...
mean(Ans), ; the arithmetic mean
√(mean(Ans²)), ; the quadratic mean
mean(Ans^-1)^-1, ; the harmonic mean
e^(mean(ln(Ans ; and the geometric mean
End
Ans(1 ;keep the first element in "Ans" and implicitly print it
Notes:
TI-BASIC is a tokenized language. Character count does not equal byte count.
e^(
is this one-byte token.
^-1
is used for this one-byte token.
I opted for writing ^-1
instead because the token looks like ֿ¹
when in a code block.
√(
is this one-byte token.
ΔList(
is this two-byte token.
$endgroup$
TI-BASIC, 42 35 34 bytes
-1 byte thanks to @SolomonUcko
While max(ΔList(Ans:{mean(Ans),√(mean(Ans²)),mean(Ans^-1)^-1,e^(mean(ln(Ans:End:Ans(1
Input is a list of two integers in Ans
.
Output is stored in Ans
and is automatically printed out when the program completes.
Formulas used for geometric, harmonic, and quadratic means are based off of user202729's explanation.
Example:
{1,2
{1 2}
prgmCDGFB
1.455688891
{100,200
{100 200}
prgmCDGFB
145.5688891
Explanation:
(Newlines have been added for clarification. They do NOT appear in the code.)
While max(ΔList(Ans ;loop until all elements of the current list are equal
; the maximum of the change in each element will be 0
{ ;create a list containing...
mean(Ans), ; the arithmetic mean
√(mean(Ans²)), ; the quadratic mean
mean(Ans^-1)^-1, ; the harmonic mean
e^(mean(ln(Ans ; and the geometric mean
End
Ans(1 ;keep the first element in "Ans" and implicitly print it
Notes:
TI-BASIC is a tokenized language. Character count does not equal byte count.
e^(
is this one-byte token.
^-1
is used for this one-byte token.
I opted for writing ^-1
instead because the token looks like ֿ¹
when in a code block.
√(
is this one-byte token.
ΔList(
is this two-byte token.
edited Apr 1 at 13:35
answered Mar 31 at 18:42
TauTau
806313
806313
$begingroup$
I think you can save a parenthesis by putting the geometric mean last.
$endgroup$
– Solomon Ucko
Apr 1 at 11:06
$begingroup$
@SolomonUcko ah, thanks for noticing! Didn't consider that before.
$endgroup$
– Tau
Apr 1 at 13:01
add a comment |
$begingroup$
I think you can save a parenthesis by putting the geometric mean last.
$endgroup$
– Solomon Ucko
Apr 1 at 11:06
$begingroup$
@SolomonUcko ah, thanks for noticing! Didn't consider that before.
$endgroup$
– Tau
Apr 1 at 13:01
$begingroup$
I think you can save a parenthesis by putting the geometric mean last.
$endgroup$
– Solomon Ucko
Apr 1 at 11:06
$begingroup$
I think you can save a parenthesis by putting the geometric mean last.
$endgroup$
– Solomon Ucko
Apr 1 at 11:06
$begingroup$
@SolomonUcko ah, thanks for noticing! Didn't consider that before.
$endgroup$
– Tau
Apr 1 at 13:01
$begingroup$
@SolomonUcko ah, thanks for noticing! Didn't consider that before.
$endgroup$
– Tau
Apr 1 at 13:01
add a comment |
$begingroup$
Java 10, 234 229 214 211 215 206 203 196 180 177 bytes
a->{for(;a[1]-a[0]>4e-9;){double l=a.length,A[]={0,0,0,1};for(var d:a){A[2]+=d/l;A[3]*=Math.pow(d,1/l);A[0]+=1/d;A[1]+=d*d;}A[0]=l/A[0];A[1]=Math.sqrt(A[1]/l);a=A;}return a[0];}
-5 bytes thanks to @PeterCordes.
-15 more bytes thanks to @PeterCordes, inspired by @RobinRyder's R answer.
+4 bytes because I assumed the inputs are pre-ordered.
-27 bytes thanks to @OlivierGrégoire.
Try it online.
Explanation:
a->{ // Method with double-array parameter and double return-type
for(;a[1]-a[0] // Loop as long as the difference between the 2nd and 1st items
>4e-9;){ // is larger than 0.000000004:
double l=a.length, // Set `l` to the amount of values in the array `a`
A[]={0,0,0,1}; // Create an array `A`, filled with the values [0,0,0,1]
for(var d:a){ // Inner loop over the values of `a`:
A[2]+=d/l; // Calculate the sum divided by the length in the third spot
A[3]*=Math.pow(d,1/l);// The product of the power of 1/length in the fourth spot
A[0]+=1/d; // The sum of 1/value in the first spot
A[1]+=d*d; // And the sum of squares in the second spot
} // After the inner loop:
// (the third spot of the array now holds the Arithmetic Mean)
// (the fourth spot of the array now holds the Geometric Mean)
A[0]=l/A[0]; // Divide the length by the first spot
// (this first spot of the array now holds the Harmonic Mean)
A[1]=Math.sqrt(A[1]/l); // Take the square of the second spot divided by the length
// (this second spot of the array now holds the Quadratic Mean)
a=A; // And then replace input `a` with array `A`
} // After the outer loop when all values are approximately equal:
return a[0];} // Return the value in the first spot as result
$endgroup$
$begingroup$
In C you couldf+=Math.abs(d-D)<1e-9;
and get implicit conversion from a boolean compare result to a 0 / 1 integer and thendouble
. Does Java have any compact syntax for that? Or is it possible to dof+=Math.abs(d-D)
and then check that the sum of absolute differences is small enough?
$endgroup$
– Peter Cordes
Apr 2 at 2:28
1
$begingroup$
Yup, for your test cases,f>1e-8
works as a loop condition: 229 bytes.a->{for(double f=1,D,A[],l;f>1e-8;a=A){D=a[0];A=new double[]{f=0,1,0,0};for(var d:a){f+=Math.abs(d-D);A[0]+=d;A[1]*=d;A[2]+=1/d;A[3]+=d*d;}A[0]/=l=a.length;A[1]=Math.pow(A[1],1/l);A[2]=l/A[2];A[3]=Math.sqrt(A[3]/l);}return a[0];}
. With1e-9
, it runs slower (about twice the CPU time), having to do more iterations to get essentially 4 *d-D
down that small. With1e-7
, it's about the same speed as 1e-8. With1e-6
, some of the trailing digits differ for one case.
$endgroup$
– Peter Cordes
Apr 2 at 2:36
1
$begingroup$
@RobinRyder's answer points out that the quadratic mean is always the largest, and harmonic always the smallest, so perhaps you can ditchf
entirely and only checka[3]-a[2]<4e-9
.
$endgroup$
– Peter Cordes
Apr 2 at 7:32
1
$begingroup$
@PeterCordesl==2||
you mean (golfed tol<3|
). But yes, good point; I've added it. :)
$endgroup$
– Kevin Cruijssen
Apr 2 at 8:03
2
$begingroup$
180 bytes by aggregating aggregable reducers.
$endgroup$
– Olivier Grégoire
Apr 2 at 8:45
|
show 12 more comments
$begingroup$
Java 10, 234 229 214 211 215 206 203 196 180 177 bytes
a->{for(;a[1]-a[0]>4e-9;){double l=a.length,A[]={0,0,0,1};for(var d:a){A[2]+=d/l;A[3]*=Math.pow(d,1/l);A[0]+=1/d;A[1]+=d*d;}A[0]=l/A[0];A[1]=Math.sqrt(A[1]/l);a=A;}return a[0];}
-5 bytes thanks to @PeterCordes.
-15 more bytes thanks to @PeterCordes, inspired by @RobinRyder's R answer.
+4 bytes because I assumed the inputs are pre-ordered.
-27 bytes thanks to @OlivierGrégoire.
Try it online.
Explanation:
a->{ // Method with double-array parameter and double return-type
for(;a[1]-a[0] // Loop as long as the difference between the 2nd and 1st items
>4e-9;){ // is larger than 0.000000004:
double l=a.length, // Set `l` to the amount of values in the array `a`
A[]={0,0,0,1}; // Create an array `A`, filled with the values [0,0,0,1]
for(var d:a){ // Inner loop over the values of `a`:
A[2]+=d/l; // Calculate the sum divided by the length in the third spot
A[3]*=Math.pow(d,1/l);// The product of the power of 1/length in the fourth spot
A[0]+=1/d; // The sum of 1/value in the first spot
A[1]+=d*d; // And the sum of squares in the second spot
} // After the inner loop:
// (the third spot of the array now holds the Arithmetic Mean)
// (the fourth spot of the array now holds the Geometric Mean)
A[0]=l/A[0]; // Divide the length by the first spot
// (this first spot of the array now holds the Harmonic Mean)
A[1]=Math.sqrt(A[1]/l); // Take the square of the second spot divided by the length
// (this second spot of the array now holds the Quadratic Mean)
a=A; // And then replace input `a` with array `A`
} // After the outer loop when all values are approximately equal:
return a[0];} // Return the value in the first spot as result
$endgroup$
$begingroup$
In C you couldf+=Math.abs(d-D)<1e-9;
and get implicit conversion from a boolean compare result to a 0 / 1 integer and thendouble
. Does Java have any compact syntax for that? Or is it possible to dof+=Math.abs(d-D)
and then check that the sum of absolute differences is small enough?
$endgroup$
– Peter Cordes
Apr 2 at 2:28
1
$begingroup$
Yup, for your test cases,f>1e-8
works as a loop condition: 229 bytes.a->{for(double f=1,D,A[],l;f>1e-8;a=A){D=a[0];A=new double[]{f=0,1,0,0};for(var d:a){f+=Math.abs(d-D);A[0]+=d;A[1]*=d;A[2]+=1/d;A[3]+=d*d;}A[0]/=l=a.length;A[1]=Math.pow(A[1],1/l);A[2]=l/A[2];A[3]=Math.sqrt(A[3]/l);}return a[0];}
. With1e-9
, it runs slower (about twice the CPU time), having to do more iterations to get essentially 4 *d-D
down that small. With1e-7
, it's about the same speed as 1e-8. With1e-6
, some of the trailing digits differ for one case.
$endgroup$
– Peter Cordes
Apr 2 at 2:36
1
$begingroup$
@RobinRyder's answer points out that the quadratic mean is always the largest, and harmonic always the smallest, so perhaps you can ditchf
entirely and only checka[3]-a[2]<4e-9
.
$endgroup$
– Peter Cordes
Apr 2 at 7:32
1
$begingroup$
@PeterCordesl==2||
you mean (golfed tol<3|
). But yes, good point; I've added it. :)
$endgroup$
– Kevin Cruijssen
Apr 2 at 8:03
2
$begingroup$
180 bytes by aggregating aggregable reducers.
$endgroup$
– Olivier Grégoire
Apr 2 at 8:45
|
show 12 more comments
$begingroup$
Java 10, 234 229 214 211 215 206 203 196 180 177 bytes
a->{for(;a[1]-a[0]>4e-9;){double l=a.length,A[]={0,0,0,1};for(var d:a){A[2]+=d/l;A[3]*=Math.pow(d,1/l);A[0]+=1/d;A[1]+=d*d;}A[0]=l/A[0];A[1]=Math.sqrt(A[1]/l);a=A;}return a[0];}
-5 bytes thanks to @PeterCordes.
-15 more bytes thanks to @PeterCordes, inspired by @RobinRyder's R answer.
+4 bytes because I assumed the inputs are pre-ordered.
-27 bytes thanks to @OlivierGrégoire.
Try it online.
Explanation:
a->{ // Method with double-array parameter and double return-type
for(;a[1]-a[0] // Loop as long as the difference between the 2nd and 1st items
>4e-9;){ // is larger than 0.000000004:
double l=a.length, // Set `l` to the amount of values in the array `a`
A[]={0,0,0,1}; // Create an array `A`, filled with the values [0,0,0,1]
for(var d:a){ // Inner loop over the values of `a`:
A[2]+=d/l; // Calculate the sum divided by the length in the third spot
A[3]*=Math.pow(d,1/l);// The product of the power of 1/length in the fourth spot
A[0]+=1/d; // The sum of 1/value in the first spot
A[1]+=d*d; // And the sum of squares in the second spot
} // After the inner loop:
// (the third spot of the array now holds the Arithmetic Mean)
// (the fourth spot of the array now holds the Geometric Mean)
A[0]=l/A[0]; // Divide the length by the first spot
// (this first spot of the array now holds the Harmonic Mean)
A[1]=Math.sqrt(A[1]/l); // Take the square of the second spot divided by the length
// (this second spot of the array now holds the Quadratic Mean)
a=A; // And then replace input `a` with array `A`
} // After the outer loop when all values are approximately equal:
return a[0];} // Return the value in the first spot as result
$endgroup$
Java 10, 234 229 214 211 215 206 203 196 180 177 bytes
a->{for(;a[1]-a[0]>4e-9;){double l=a.length,A[]={0,0,0,1};for(var d:a){A[2]+=d/l;A[3]*=Math.pow(d,1/l);A[0]+=1/d;A[1]+=d*d;}A[0]=l/A[0];A[1]=Math.sqrt(A[1]/l);a=A;}return a[0];}
-5 bytes thanks to @PeterCordes.
-15 more bytes thanks to @PeterCordes, inspired by @RobinRyder's R answer.
+4 bytes because I assumed the inputs are pre-ordered.
-27 bytes thanks to @OlivierGrégoire.
Try it online.
Explanation:
a->{ // Method with double-array parameter and double return-type
for(;a[1]-a[0] // Loop as long as the difference between the 2nd and 1st items
>4e-9;){ // is larger than 0.000000004:
double l=a.length, // Set `l` to the amount of values in the array `a`
A[]={0,0,0,1}; // Create an array `A`, filled with the values [0,0,0,1]
for(var d:a){ // Inner loop over the values of `a`:
A[2]+=d/l; // Calculate the sum divided by the length in the third spot
A[3]*=Math.pow(d,1/l);// The product of the power of 1/length in the fourth spot
A[0]+=1/d; // The sum of 1/value in the first spot
A[1]+=d*d; // And the sum of squares in the second spot
} // After the inner loop:
// (the third spot of the array now holds the Arithmetic Mean)
// (the fourth spot of the array now holds the Geometric Mean)
A[0]=l/A[0]; // Divide the length by the first spot
// (this first spot of the array now holds the Harmonic Mean)
A[1]=Math.sqrt(A[1]/l); // Take the square of the second spot divided by the length
// (this second spot of the array now holds the Quadratic Mean)
a=A; // And then replace input `a` with array `A`
} // After the outer loop when all values are approximately equal:
return a[0];} // Return the value in the first spot as result
edited Apr 2 at 9:30
answered Apr 1 at 9:12
Kevin CruijssenKevin Cruijssen
42.3k570217
42.3k570217
$begingroup$
In C you couldf+=Math.abs(d-D)<1e-9;
and get implicit conversion from a boolean compare result to a 0 / 1 integer and thendouble
. Does Java have any compact syntax for that? Or is it possible to dof+=Math.abs(d-D)
and then check that the sum of absolute differences is small enough?
$endgroup$
– Peter Cordes
Apr 2 at 2:28
1
$begingroup$
Yup, for your test cases,f>1e-8
works as a loop condition: 229 bytes.a->{for(double f=1,D,A[],l;f>1e-8;a=A){D=a[0];A=new double[]{f=0,1,0,0};for(var d:a){f+=Math.abs(d-D);A[0]+=d;A[1]*=d;A[2]+=1/d;A[3]+=d*d;}A[0]/=l=a.length;A[1]=Math.pow(A[1],1/l);A[2]=l/A[2];A[3]=Math.sqrt(A[3]/l);}return a[0];}
. With1e-9
, it runs slower (about twice the CPU time), having to do more iterations to get essentially 4 *d-D
down that small. With1e-7
, it's about the same speed as 1e-8. With1e-6
, some of the trailing digits differ for one case.
$endgroup$
– Peter Cordes
Apr 2 at 2:36
1
$begingroup$
@RobinRyder's answer points out that the quadratic mean is always the largest, and harmonic always the smallest, so perhaps you can ditchf
entirely and only checka[3]-a[2]<4e-9
.
$endgroup$
– Peter Cordes
Apr 2 at 7:32
1
$begingroup$
@PeterCordesl==2||
you mean (golfed tol<3|
). But yes, good point; I've added it. :)
$endgroup$
– Kevin Cruijssen
Apr 2 at 8:03
2
$begingroup$
180 bytes by aggregating aggregable reducers.
$endgroup$
– Olivier Grégoire
Apr 2 at 8:45
|
show 12 more comments
$begingroup$
In C you couldf+=Math.abs(d-D)<1e-9;
and get implicit conversion from a boolean compare result to a 0 / 1 integer and thendouble
. Does Java have any compact syntax for that? Or is it possible to dof+=Math.abs(d-D)
and then check that the sum of absolute differences is small enough?
$endgroup$
– Peter Cordes
Apr 2 at 2:28
1
$begingroup$
Yup, for your test cases,f>1e-8
works as a loop condition: 229 bytes.a->{for(double f=1,D,A[],l;f>1e-8;a=A){D=a[0];A=new double[]{f=0,1,0,0};for(var d:a){f+=Math.abs(d-D);A[0]+=d;A[1]*=d;A[2]+=1/d;A[3]+=d*d;}A[0]/=l=a.length;A[1]=Math.pow(A[1],1/l);A[2]=l/A[2];A[3]=Math.sqrt(A[3]/l);}return a[0];}
. With1e-9
, it runs slower (about twice the CPU time), having to do more iterations to get essentially 4 *d-D
down that small. With1e-7
, it's about the same speed as 1e-8. With1e-6
, some of the trailing digits differ for one case.
$endgroup$
– Peter Cordes
Apr 2 at 2:36
1
$begingroup$
@RobinRyder's answer points out that the quadratic mean is always the largest, and harmonic always the smallest, so perhaps you can ditchf
entirely and only checka[3]-a[2]<4e-9
.
$endgroup$
– Peter Cordes
Apr 2 at 7:32
1
$begingroup$
@PeterCordesl==2||
you mean (golfed tol<3|
). But yes, good point; I've added it. :)
$endgroup$
– Kevin Cruijssen
Apr 2 at 8:03
2
$begingroup$
180 bytes by aggregating aggregable reducers.
$endgroup$
– Olivier Grégoire
Apr 2 at 8:45
$begingroup$
In C you could
f+=Math.abs(d-D)<1e-9;
and get implicit conversion from a boolean compare result to a 0 / 1 integer and then double
. Does Java have any compact syntax for that? Or is it possible to do f+=Math.abs(d-D)
and then check that the sum of absolute differences is small enough?$endgroup$
– Peter Cordes
Apr 2 at 2:28
$begingroup$
In C you could
f+=Math.abs(d-D)<1e-9;
and get implicit conversion from a boolean compare result to a 0 / 1 integer and then double
. Does Java have any compact syntax for that? Or is it possible to do f+=Math.abs(d-D)
and then check that the sum of absolute differences is small enough?$endgroup$
– Peter Cordes
Apr 2 at 2:28
1
1
$begingroup$
Yup, for your test cases,
f>1e-8
works as a loop condition: 229 bytes. a->{for(double f=1,D,A[],l;f>1e-8;a=A){D=a[0];A=new double[]{f=0,1,0,0};for(var d:a){f+=Math.abs(d-D);A[0]+=d;A[1]*=d;A[2]+=1/d;A[3]+=d*d;}A[0]/=l=a.length;A[1]=Math.pow(A[1],1/l);A[2]=l/A[2];A[3]=Math.sqrt(A[3]/l);}return a[0];}
. With 1e-9
, it runs slower (about twice the CPU time), having to do more iterations to get essentially 4 * d-D
down that small. With 1e-7
, it's about the same speed as 1e-8. With 1e-6
, some of the trailing digits differ for one case.$endgroup$
– Peter Cordes
Apr 2 at 2:36
$begingroup$
Yup, for your test cases,
f>1e-8
works as a loop condition: 229 bytes. a->{for(double f=1,D,A[],l;f>1e-8;a=A){D=a[0];A=new double[]{f=0,1,0,0};for(var d:a){f+=Math.abs(d-D);A[0]+=d;A[1]*=d;A[2]+=1/d;A[3]+=d*d;}A[0]/=l=a.length;A[1]=Math.pow(A[1],1/l);A[2]=l/A[2];A[3]=Math.sqrt(A[3]/l);}return a[0];}
. With 1e-9
, it runs slower (about twice the CPU time), having to do more iterations to get essentially 4 * d-D
down that small. With 1e-7
, it's about the same speed as 1e-8. With 1e-6
, some of the trailing digits differ for one case.$endgroup$
– Peter Cordes
Apr 2 at 2:36
1
1
$begingroup$
@RobinRyder's answer points out that the quadratic mean is always the largest, and harmonic always the smallest, so perhaps you can ditch
f
entirely and only check a[3]-a[2]<4e-9
.$endgroup$
– Peter Cordes
Apr 2 at 7:32
$begingroup$
@RobinRyder's answer points out that the quadratic mean is always the largest, and harmonic always the smallest, so perhaps you can ditch
f
entirely and only check a[3]-a[2]<4e-9
.$endgroup$
– Peter Cordes
Apr 2 at 7:32
1
1
$begingroup$
@PeterCordes
l==2||
you mean (golfed to l<3|
). But yes, good point; I've added it. :)$endgroup$
– Kevin Cruijssen
Apr 2 at 8:03
$begingroup$
@PeterCordes
l==2||
you mean (golfed to l<3|
). But yes, good point; I've added it. :)$endgroup$
– Kevin Cruijssen
Apr 2 at 8:03
2
2
$begingroup$
180 bytes by aggregating aggregable reducers.
$endgroup$
– Olivier Grégoire
Apr 2 at 8:45
$begingroup$
180 bytes by aggregating aggregable reducers.
$endgroup$
– Olivier Grégoire
Apr 2 at 8:45
|
show 12 more comments
$begingroup$
Charcoal, 40 bytes
W‹⌊θ⌈θ≔⟦∕ΣθLθXΠθ∕¹Lθ∕LθΣ∕¹θ₂∕ΣXθ²Lθ⟧θI⊟θ
Try it online! Link is to verbose version of code. Takes input as an array of numbers. Explanation:
W‹⌊θ⌈θ
Repeat while the array contains different values...
≔⟦....⟧θ
... replace the array with a list of values:
∕ΣθLθ
... the mean...
XΠθ∕¹Lθ
... the geometric mean...
∕LθΣ∕¹θ
... the harmonic mean...
₂∕ΣXθ²Lθ
... and the root mean square.
I⊟θ
Cast an element of the array to string and implicitly print it.
$endgroup$
add a comment |
$begingroup$
Charcoal, 40 bytes
W‹⌊θ⌈θ≔⟦∕ΣθLθXΠθ∕¹Lθ∕LθΣ∕¹θ₂∕ΣXθ²Lθ⟧θI⊟θ
Try it online! Link is to verbose version of code. Takes input as an array of numbers. Explanation:
W‹⌊θ⌈θ
Repeat while the array contains different values...
≔⟦....⟧θ
... replace the array with a list of values:
∕ΣθLθ
... the mean...
XΠθ∕¹Lθ
... the geometric mean...
∕LθΣ∕¹θ
... the harmonic mean...
₂∕ΣXθ²Lθ
... and the root mean square.
I⊟θ
Cast an element of the array to string and implicitly print it.
$endgroup$
add a comment |
$begingroup$
Charcoal, 40 bytes
W‹⌊θ⌈θ≔⟦∕ΣθLθXΠθ∕¹Lθ∕LθΣ∕¹θ₂∕ΣXθ²Lθ⟧θI⊟θ
Try it online! Link is to verbose version of code. Takes input as an array of numbers. Explanation:
W‹⌊θ⌈θ
Repeat while the array contains different values...
≔⟦....⟧θ
... replace the array with a list of values:
∕ΣθLθ
... the mean...
XΠθ∕¹Lθ
... the geometric mean...
∕LθΣ∕¹θ
... the harmonic mean...
₂∕ΣXθ²Lθ
... and the root mean square.
I⊟θ
Cast an element of the array to string and implicitly print it.
$endgroup$
Charcoal, 40 bytes
W‹⌊θ⌈θ≔⟦∕ΣθLθXΠθ∕¹Lθ∕LθΣ∕¹θ₂∕ΣXθ²Lθ⟧θI⊟θ
Try it online! Link is to verbose version of code. Takes input as an array of numbers. Explanation:
W‹⌊θ⌈θ
Repeat while the array contains different values...
≔⟦....⟧θ
... replace the array with a list of values:
∕ΣθLθ
... the mean...
XΠθ∕¹Lθ
... the geometric mean...
∕LθΣ∕¹θ
... the harmonic mean...
₂∕ΣXθ²Lθ
... and the root mean square.
I⊟θ
Cast an element of the array to string and implicitly print it.
answered Mar 31 at 9:21
NeilNeil
82.6k745179
82.6k745179
add a comment |
add a comment |
$begingroup$
Jelly, 24 bytes
;µ*Ɱ-;ؽ¤©Æm*İɼ;P½½ƊµÐLḢ
Try it online!
$endgroup$
add a comment |
$begingroup$
Jelly, 24 bytes
;µ*Ɱ-;ؽ¤©Æm*İɼ;P½½ƊµÐLḢ
Try it online!
$endgroup$
add a comment |
$begingroup$
Jelly, 24 bytes
;µ*Ɱ-;ؽ¤©Æm*İɼ;P½½ƊµÐLḢ
Try it online!
$endgroup$
Jelly, 24 bytes
;µ*Ɱ-;ؽ¤©Æm*İɼ;P½½ƊµÐLḢ
Try it online!
answered Mar 31 at 15:42
Erik the OutgolferErik the Outgolfer
33k429106
33k429106
add a comment |
add a comment |
$begingroup$
PowerShell, 182 180 183 bytes
$f={$a=$c=$d=$n=0
$b=1
$m=[math]
$args|%{$n++
$a+=$_
$b*=$_
$c+=1/$_
$d+=+$_*$_}
$p=($a/$n),$m::pow($b,1/$n),($n/$c),$m::sqrt($d/$n)|%{$m::round($_,9)}
if($p-ne$p[0]){$p=&$f @p}$p[0]}
Try it online!
$endgroup$
$begingroup$
171 bytes
$endgroup$
– mazzy
Apr 1 at 12:10
$begingroup$
@mazzy, impressively :)
$endgroup$
– Andrei Odegov
Apr 1 at 14:14
add a comment |
$begingroup$
PowerShell, 182 180 183 bytes
$f={$a=$c=$d=$n=0
$b=1
$m=[math]
$args|%{$n++
$a+=$_
$b*=$_
$c+=1/$_
$d+=+$_*$_}
$p=($a/$n),$m::pow($b,1/$n),($n/$c),$m::sqrt($d/$n)|%{$m::round($_,9)}
if($p-ne$p[0]){$p=&$f @p}$p[0]}
Try it online!
$endgroup$
$begingroup$
171 bytes
$endgroup$
– mazzy
Apr 1 at 12:10
$begingroup$
@mazzy, impressively :)
$endgroup$
– Andrei Odegov
Apr 1 at 14:14
add a comment |
$begingroup$
PowerShell, 182 180 183 bytes
$f={$a=$c=$d=$n=0
$b=1
$m=[math]
$args|%{$n++
$a+=$_
$b*=$_
$c+=1/$_
$d+=+$_*$_}
$p=($a/$n),$m::pow($b,1/$n),($n/$c),$m::sqrt($d/$n)|%{$m::round($_,9)}
if($p-ne$p[0]){$p=&$f @p}$p[0]}
Try it online!
$endgroup$
PowerShell, 182 180 183 bytes
$f={$a=$c=$d=$n=0
$b=1
$m=[math]
$args|%{$n++
$a+=$_
$b*=$_
$c+=1/$_
$d+=+$_*$_}
$p=($a/$n),$m::pow($b,1/$n),($n/$c),$m::sqrt($d/$n)|%{$m::round($_,9)}
if($p-ne$p[0]){$p=&$f @p}$p[0]}
Try it online!
edited Apr 1 at 19:19
answered Mar 31 at 17:38
Andrei OdegovAndrei Odegov
51926
51926
$begingroup$
171 bytes
$endgroup$
– mazzy
Apr 1 at 12:10
$begingroup$
@mazzy, impressively :)
$endgroup$
– Andrei Odegov
Apr 1 at 14:14
add a comment |
$begingroup$
171 bytes
$endgroup$
– mazzy
Apr 1 at 12:10
$begingroup$
@mazzy, impressively :)
$endgroup$
– Andrei Odegov
Apr 1 at 14:14
$begingroup$
171 bytes
$endgroup$
– mazzy
Apr 1 at 12:10
$begingroup$
171 bytes
$endgroup$
– mazzy
Apr 1 at 12:10
$begingroup$
@mazzy, impressively :)
$endgroup$
– Andrei Odegov
Apr 1 at 14:14
$begingroup$
@mazzy, impressively :)
$endgroup$
– Andrei Odegov
Apr 1 at 14:14
add a comment |
$begingroup$
Jelly, 25 24 bytes
Wẋ4¹ÆlÆeƭ²½ƭİ4ƭÆm$€⁺µÐLḢ
Try it online!
Explanation
µÐL | Repeat until unchanged:
W | Wrap as a list
ẋ4 | Copy list 4 times
⁺ | Do twice:
$€ | For each copy of the list:
4ƭ | One of these functions, cycling between them:
¹ | Identity
ÆlÆeƭ | Alternate between log and exp
²½ƭ | Alternate between square and square root
İ | Reciprocal
Æm | Then take the mean
Ḣ| Finally take the first item
$endgroup$
$begingroup$
I am fairly bad at Jelly, but could something similar toP*İL
work for the geometric mean?
$endgroup$
– someone
Mar 31 at 15:08
$begingroup$
@someone it would need to beP*Lİ$
so wouldn’t save bytes. It would mean I could bringÆm
back down a line without costing bytes, but I quite like the fact each one currently has an arithmetic mean at its core.
$endgroup$
– Nick Kennedy
Mar 31 at 15:19
add a comment |
$begingroup$
Jelly, 25 24 bytes
Wẋ4¹ÆlÆeƭ²½ƭİ4ƭÆm$€⁺µÐLḢ
Try it online!
Explanation
µÐL | Repeat until unchanged:
W | Wrap as a list
ẋ4 | Copy list 4 times
⁺ | Do twice:
$€ | For each copy of the list:
4ƭ | One of these functions, cycling between them:
¹ | Identity
ÆlÆeƭ | Alternate between log and exp
²½ƭ | Alternate between square and square root
İ | Reciprocal
Æm | Then take the mean
Ḣ| Finally take the first item
$endgroup$
$begingroup$
I am fairly bad at Jelly, but could something similar toP*İL
work for the geometric mean?
$endgroup$
– someone
Mar 31 at 15:08
$begingroup$
@someone it would need to beP*Lİ$
so wouldn’t save bytes. It would mean I could bringÆm
back down a line without costing bytes, but I quite like the fact each one currently has an arithmetic mean at its core.
$endgroup$
– Nick Kennedy
Mar 31 at 15:19
add a comment |
$begingroup$
Jelly, 25 24 bytes
Wẋ4¹ÆlÆeƭ²½ƭİ4ƭÆm$€⁺µÐLḢ
Try it online!
Explanation
µÐL | Repeat until unchanged:
W | Wrap as a list
ẋ4 | Copy list 4 times
⁺ | Do twice:
$€ | For each copy of the list:
4ƭ | One of these functions, cycling between them:
¹ | Identity
ÆlÆeƭ | Alternate between log and exp
²½ƭ | Alternate between square and square root
İ | Reciprocal
Æm | Then take the mean
Ḣ| Finally take the first item
$endgroup$
Jelly, 25 24 bytes
Wẋ4¹ÆlÆeƭ²½ƭİ4ƭÆm$€⁺µÐLḢ
Try it online!
Explanation
µÐL | Repeat until unchanged:
W | Wrap as a list
ẋ4 | Copy list 4 times
⁺ | Do twice:
$€ | For each copy of the list:
4ƭ | One of these functions, cycling between them:
¹ | Identity
ÆlÆeƭ | Alternate between log and exp
²½ƭ | Alternate between square and square root
İ | Reciprocal
Æm | Then take the mean
Ḣ| Finally take the first item
edited Mar 31 at 23:58
answered Mar 31 at 14:29
Nick KennedyNick Kennedy
1,32649
1,32649
$begingroup$
I am fairly bad at Jelly, but could something similar toP*İL
work for the geometric mean?
$endgroup$
– someone
Mar 31 at 15:08
$begingroup$
@someone it would need to beP*Lİ$
so wouldn’t save bytes. It would mean I could bringÆm
back down a line without costing bytes, but I quite like the fact each one currently has an arithmetic mean at its core.
$endgroup$
– Nick Kennedy
Mar 31 at 15:19
add a comment |
$begingroup$
I am fairly bad at Jelly, but could something similar toP*İL
work for the geometric mean?
$endgroup$
– someone
Mar 31 at 15:08
$begingroup$
@someone it would need to beP*Lİ$
so wouldn’t save bytes. It would mean I could bringÆm
back down a line without costing bytes, but I quite like the fact each one currently has an arithmetic mean at its core.
$endgroup$
– Nick Kennedy
Mar 31 at 15:19
$begingroup$
I am fairly bad at Jelly, but could something similar to
P*İL
work for the geometric mean?$endgroup$
– someone
Mar 31 at 15:08
$begingroup$
I am fairly bad at Jelly, but could something similar to
P*İL
work for the geometric mean?$endgroup$
– someone
Mar 31 at 15:08
$begingroup$
@someone it would need to be
P*Lİ$
so wouldn’t save bytes. It would mean I could bring Æm
back down a line without costing bytes, but I quite like the fact each one currently has an arithmetic mean at its core.$endgroup$
– Nick Kennedy
Mar 31 at 15:19
$begingroup$
@someone it would need to be
P*Lİ$
so wouldn’t save bytes. It would mean I could bring Æm
back down a line without costing bytes, but I quite like the fact each one currently has an arithmetic mean at its core.$endgroup$
– Nick Kennedy
Mar 31 at 15:19
add a comment |
$begingroup$
05AB1E (legacy), 26 24 bytes
Δ©ÅA®PYzmY®zO/®nÅAt)4V}н
Try it online or see the steps of all test cases.
Explanation:
Δ # Loop until the list no longer changes:
© # Store the current list in the register (without popping)
# (which is the implicit input-list in the first iteration)
# Arithmetic mean:
ÅA # Builtin to calculate the arithmetic mean of this list
# Geometric mean:
®P # Take the product of the list from the register
Yz # And 1/`Y` (NOTE: `Y` is 2 by default)
m # Then take the result of the product to the power of 1/`Y`
# (NOTE: a^(1/4) is the same as ⁴√a)
# Harmonic mean:
®z # Get 1/x for each value x in the list from the register
O # Take the sum of that
Y / # And then divide `Y` by that sum
# Quadratic mean:
®n # Take 2 to the power of each number x in the list from the register
ÅA # Calculate the arithmetic mean of this list
t # And then take the square-root of that
) # Wrap all four results into a list
4V # And set variable `Y` now to 4 instead of 2
}н # After the list no longer changes: pop and push its first value
# (which is output implicitly as result)
$endgroup$
add a comment |
$begingroup$
05AB1E (legacy), 26 24 bytes
Δ©ÅA®PYzmY®zO/®nÅAt)4V}н
Try it online or see the steps of all test cases.
Explanation:
Δ # Loop until the list no longer changes:
© # Store the current list in the register (without popping)
# (which is the implicit input-list in the first iteration)
# Arithmetic mean:
ÅA # Builtin to calculate the arithmetic mean of this list
# Geometric mean:
®P # Take the product of the list from the register
Yz # And 1/`Y` (NOTE: `Y` is 2 by default)
m # Then take the result of the product to the power of 1/`Y`
# (NOTE: a^(1/4) is the same as ⁴√a)
# Harmonic mean:
®z # Get 1/x for each value x in the list from the register
O # Take the sum of that
Y / # And then divide `Y` by that sum
# Quadratic mean:
®n # Take 2 to the power of each number x in the list from the register
ÅA # Calculate the arithmetic mean of this list
t # And then take the square-root of that
) # Wrap all four results into a list
4V # And set variable `Y` now to 4 instead of 2
}н # After the list no longer changes: pop and push its first value
# (which is output implicitly as result)
$endgroup$
add a comment |
$begingroup$
05AB1E (legacy), 26 24 bytes
Δ©ÅA®PYzmY®zO/®nÅAt)4V}н
Try it online or see the steps of all test cases.
Explanation:
Δ # Loop until the list no longer changes:
© # Store the current list in the register (without popping)
# (which is the implicit input-list in the first iteration)
# Arithmetic mean:
ÅA # Builtin to calculate the arithmetic mean of this list
# Geometric mean:
®P # Take the product of the list from the register
Yz # And 1/`Y` (NOTE: `Y` is 2 by default)
m # Then take the result of the product to the power of 1/`Y`
# (NOTE: a^(1/4) is the same as ⁴√a)
# Harmonic mean:
®z # Get 1/x for each value x in the list from the register
O # Take the sum of that
Y / # And then divide `Y` by that sum
# Quadratic mean:
®n # Take 2 to the power of each number x in the list from the register
ÅA # Calculate the arithmetic mean of this list
t # And then take the square-root of that
) # Wrap all four results into a list
4V # And set variable `Y` now to 4 instead of 2
}н # After the list no longer changes: pop and push its first value
# (which is output implicitly as result)
$endgroup$
05AB1E (legacy), 26 24 bytes
Δ©ÅA®PYzmY®zO/®nÅAt)4V}н
Try it online or see the steps of all test cases.
Explanation:
Δ # Loop until the list no longer changes:
© # Store the current list in the register (without popping)
# (which is the implicit input-list in the first iteration)
# Arithmetic mean:
ÅA # Builtin to calculate the arithmetic mean of this list
# Geometric mean:
®P # Take the product of the list from the register
Yz # And 1/`Y` (NOTE: `Y` is 2 by default)
m # Then take the result of the product to the power of 1/`Y`
# (NOTE: a^(1/4) is the same as ⁴√a)
# Harmonic mean:
®z # Get 1/x for each value x in the list from the register
O # Take the sum of that
Y / # And then divide `Y` by that sum
# Quadratic mean:
®n # Take 2 to the power of each number x in the list from the register
ÅA # Calculate the arithmetic mean of this list
t # And then take the square-root of that
) # Wrap all four results into a list
4V # And set variable `Y` now to 4 instead of 2
}н # After the list no longer changes: pop and push its first value
# (which is output implicitly as result)
edited Apr 1 at 12:28
answered Apr 1 at 7:46
Kevin CruijssenKevin Cruijssen
42.3k570217
42.3k570217
add a comment |
add a comment |
$begingroup$
C#, 173 bytes
double m(int n,params double[]a)=>(n<1?a[0]:m(n-1,a.Sum()/a.Length,Math.Pow(a.Aggregate((t,x)=>t*x),1.0/a.Length),a.Length/a.Sum(x=>1/x),Math.Sqrt(a.Sum(x=>x*x)/a.Length)));
Try it online!
$endgroup$
2
$begingroup$
This seems to really on a variable that must be passed. Also, you have to includeusing System
andusing System.Linq
in your byte count, since they are required for the program to run. You can change your compiler to the C# Visual Interactive Compiler, which doesn't need those imports. Also,1.0
->1d
$endgroup$
– Embodiment of Ignorance
Mar 31 at 20:59
add a comment |
$begingroup$
C#, 173 bytes
double m(int n,params double[]a)=>(n<1?a[0]:m(n-1,a.Sum()/a.Length,Math.Pow(a.Aggregate((t,x)=>t*x),1.0/a.Length),a.Length/a.Sum(x=>1/x),Math.Sqrt(a.Sum(x=>x*x)/a.Length)));
Try it online!
$endgroup$
2
$begingroup$
This seems to really on a variable that must be passed. Also, you have to includeusing System
andusing System.Linq
in your byte count, since they are required for the program to run. You can change your compiler to the C# Visual Interactive Compiler, which doesn't need those imports. Also,1.0
->1d
$endgroup$
– Embodiment of Ignorance
Mar 31 at 20:59
add a comment |
$begingroup$
C#, 173 bytes
double m(int n,params double[]a)=>(n<1?a[0]:m(n-1,a.Sum()/a.Length,Math.Pow(a.Aggregate((t,x)=>t*x),1.0/a.Length),a.Length/a.Sum(x=>1/x),Math.Sqrt(a.Sum(x=>x*x)/a.Length)));
Try it online!
$endgroup$
C#, 173 bytes
double m(int n,params double[]a)=>(n<1?a[0]:m(n-1,a.Sum()/a.Length,Math.Pow(a.Aggregate((t,x)=>t*x),1.0/a.Length),a.Length/a.Sum(x=>1/x),Math.Sqrt(a.Sum(x=>x*x)/a.Length)));
Try it online!
answered Mar 31 at 19:39
aviatoraviator
1112
1112
2
$begingroup$
This seems to really on a variable that must be passed. Also, you have to includeusing System
andusing System.Linq
in your byte count, since they are required for the program to run. You can change your compiler to the C# Visual Interactive Compiler, which doesn't need those imports. Also,1.0
->1d
$endgroup$
– Embodiment of Ignorance
Mar 31 at 20:59
add a comment |
2
$begingroup$
This seems to really on a variable that must be passed. Also, you have to includeusing System
andusing System.Linq
in your byte count, since they are required for the program to run. You can change your compiler to the C# Visual Interactive Compiler, which doesn't need those imports. Also,1.0
->1d
$endgroup$
– Embodiment of Ignorance
Mar 31 at 20:59
2
2
$begingroup$
This seems to really on a variable that must be passed. Also, you have to include
using System
and using System.Linq
in your byte count, since they are required for the program to run. You can change your compiler to the C# Visual Interactive Compiler, which doesn't need those imports. Also, 1.0
-> 1d
$endgroup$
– Embodiment of Ignorance
Mar 31 at 20:59
$begingroup$
This seems to really on a variable that must be passed. Also, you have to include
using System
and using System.Linq
in your byte count, since they are required for the program to run. You can change your compiler to the C# Visual Interactive Compiler, which doesn't need those imports. Also, 1.0
-> 1d
$endgroup$
– Embodiment of Ignorance
Mar 31 at 20:59
add a comment |
$begingroup$
Clean, 124 bytes
import StdEnv
f=avg o limit o iteratel=let n=toReal(length l)in[avg l,prod l^(1.0/n),n/sum[1.0/x\x<-l],avg[x*x\x<-l]^0.5]
Try it online!
Performs the operation until the result stops changing.
Hurray for limited precision floating point!
$endgroup$
add a comment |
$begingroup$
Clean, 124 bytes
import StdEnv
f=avg o limit o iteratel=let n=toReal(length l)in[avg l,prod l^(1.0/n),n/sum[1.0/x\x<-l],avg[x*x\x<-l]^0.5]
Try it online!
Performs the operation until the result stops changing.
Hurray for limited precision floating point!
$endgroup$
add a comment |
$begingroup$
Clean, 124 bytes
import StdEnv
f=avg o limit o iteratel=let n=toReal(length l)in[avg l,prod l^(1.0/n),n/sum[1.0/x\x<-l],avg[x*x\x<-l]^0.5]
Try it online!
Performs the operation until the result stops changing.
Hurray for limited precision floating point!
$endgroup$
Clean, 124 bytes
import StdEnv
f=avg o limit o iteratel=let n=toReal(length l)in[avg l,prod l^(1.0/n),n/sum[1.0/x\x<-l],avg[x*x\x<-l]^0.5]
Try it online!
Performs the operation until the result stops changing.
Hurray for limited precision floating point!
answered Mar 31 at 21:59
ΟurousΟurous
7,41611136
7,41611136
add a comment |
add a comment |
$begingroup$
Pyth, 32 bytes
h.Wt{H[.OZ@*FZJlZcJscL1Z@.O^R2Z2
Try it online here, or verify all the test cases (bar two, see note below) at once here. Accepts input as a list.
There seem to be some issues with rounding, as certain inputs don't converge correctly when they otherwise should. In particular, test case 0.01 100
gets stuck at values [6.748305820749738, 6.748305820749738, 6.748305820749739, 6.748305820749738]
, and test case 1.61 2.41
gets stuck at [1.9896543776640825, 1.9896543776640825, 1.9896543776640827, 1.9896543776640825]
- note in both cases that the 3rd mean (harmonic mean) differs from the others.
I'm not sure if this problem invalidates my entry, but I'm posting it anyway as it should work. If this isn't acceptable, it can be fixed by instering .RRT
before the [
, to round each of the means to 10 decimal places, as seen in this test suite.
h.Wt{H[.OZ@*FZJlZcJscL1Z@.O^R2Z2)Q Implicit: Q=eval(input())
Trailing )Q inferred
.W Q Funcitonal while: While condition is true, call inner. Starting value Q
t{H Condition function: current input H
{H Deduplicate H
t Discard first value
Empty list is falsey, so while is terminated when means converge
[.OZ@*FZJlZcJscL1Z@.O^R2Z2) Inner function: current input Z
JlZ Take length of Z, store in J
.OZ (1) Arithmetic mean of Z
*FZ Product of Z
@ J (2) Jth root of the above
L Z Map each element of Z...
c 1 ... to its reciprocal
s Sum the above
cJ (3) J / the above
R Z Map each element of Z...
^ 2 ... to its square
.O Arithmetic mean of the above
@ 2 (4) Square root of the above
[ ) Wrap results (1), (2), (3), and (4) in a list
This is used as the input for the next iteration of the loop
h Take the first element of the result, implicit print
$endgroup$
add a comment |
$begingroup$
Pyth, 32 bytes
h.Wt{H[.OZ@*FZJlZcJscL1Z@.O^R2Z2
Try it online here, or verify all the test cases (bar two, see note below) at once here. Accepts input as a list.
There seem to be some issues with rounding, as certain inputs don't converge correctly when they otherwise should. In particular, test case 0.01 100
gets stuck at values [6.748305820749738, 6.748305820749738, 6.748305820749739, 6.748305820749738]
, and test case 1.61 2.41
gets stuck at [1.9896543776640825, 1.9896543776640825, 1.9896543776640827, 1.9896543776640825]
- note in both cases that the 3rd mean (harmonic mean) differs from the others.
I'm not sure if this problem invalidates my entry, but I'm posting it anyway as it should work. If this isn't acceptable, it can be fixed by instering .RRT
before the [
, to round each of the means to 10 decimal places, as seen in this test suite.
h.Wt{H[.OZ@*FZJlZcJscL1Z@.O^R2Z2)Q Implicit: Q=eval(input())
Trailing )Q inferred
.W Q Funcitonal while: While condition is true, call inner. Starting value Q
t{H Condition function: current input H
{H Deduplicate H
t Discard first value
Empty list is falsey, so while is terminated when means converge
[.OZ@*FZJlZcJscL1Z@.O^R2Z2) Inner function: current input Z
JlZ Take length of Z, store in J
.OZ (1) Arithmetic mean of Z
*FZ Product of Z
@ J (2) Jth root of the above
L Z Map each element of Z...
c 1 ... to its reciprocal
s Sum the above
cJ (3) J / the above
R Z Map each element of Z...
^ 2 ... to its square
.O Arithmetic mean of the above
@ 2 (4) Square root of the above
[ ) Wrap results (1), (2), (3), and (4) in a list
This is used as the input for the next iteration of the loop
h Take the first element of the result, implicit print
$endgroup$
add a comment |
$begingroup$
Pyth, 32 bytes
h.Wt{H[.OZ@*FZJlZcJscL1Z@.O^R2Z2
Try it online here, or verify all the test cases (bar two, see note below) at once here. Accepts input as a list.
There seem to be some issues with rounding, as certain inputs don't converge correctly when they otherwise should. In particular, test case 0.01 100
gets stuck at values [6.748305820749738, 6.748305820749738, 6.748305820749739, 6.748305820749738]
, and test case 1.61 2.41
gets stuck at [1.9896543776640825, 1.9896543776640825, 1.9896543776640827, 1.9896543776640825]
- note in both cases that the 3rd mean (harmonic mean) differs from the others.
I'm not sure if this problem invalidates my entry, but I'm posting it anyway as it should work. If this isn't acceptable, it can be fixed by instering .RRT
before the [
, to round each of the means to 10 decimal places, as seen in this test suite.
h.Wt{H[.OZ@*FZJlZcJscL1Z@.O^R2Z2)Q Implicit: Q=eval(input())
Trailing )Q inferred
.W Q Funcitonal while: While condition is true, call inner. Starting value Q
t{H Condition function: current input H
{H Deduplicate H
t Discard first value
Empty list is falsey, so while is terminated when means converge
[.OZ@*FZJlZcJscL1Z@.O^R2Z2) Inner function: current input Z
JlZ Take length of Z, store in J
.OZ (1) Arithmetic mean of Z
*FZ Product of Z
@ J (2) Jth root of the above
L Z Map each element of Z...
c 1 ... to its reciprocal
s Sum the above
cJ (3) J / the above
R Z Map each element of Z...
^ 2 ... to its square
.O Arithmetic mean of the above
@ 2 (4) Square root of the above
[ ) Wrap results (1), (2), (3), and (4) in a list
This is used as the input for the next iteration of the loop
h Take the first element of the result, implicit print
$endgroup$
Pyth, 32 bytes
h.Wt{H[.OZ@*FZJlZcJscL1Z@.O^R2Z2
Try it online here, or verify all the test cases (bar two, see note below) at once here. Accepts input as a list.
There seem to be some issues with rounding, as certain inputs don't converge correctly when they otherwise should. In particular, test case 0.01 100
gets stuck at values [6.748305820749738, 6.748305820749738, 6.748305820749739, 6.748305820749738]
, and test case 1.61 2.41
gets stuck at [1.9896543776640825, 1.9896543776640825, 1.9896543776640827, 1.9896543776640825]
- note in both cases that the 3rd mean (harmonic mean) differs from the others.
I'm not sure if this problem invalidates my entry, but I'm posting it anyway as it should work. If this isn't acceptable, it can be fixed by instering .RRT
before the [
, to round each of the means to 10 decimal places, as seen in this test suite.
h.Wt{H[.OZ@*FZJlZcJscL1Z@.O^R2Z2)Q Implicit: Q=eval(input())
Trailing )Q inferred
.W Q Funcitonal while: While condition is true, call inner. Starting value Q
t{H Condition function: current input H
{H Deduplicate H
t Discard first value
Empty list is falsey, so while is terminated when means converge
[.OZ@*FZJlZcJscL1Z@.O^R2Z2) Inner function: current input Z
JlZ Take length of Z, store in J
.OZ (1) Arithmetic mean of Z
*FZ Product of Z
@ J (2) Jth root of the above
L Z Map each element of Z...
c 1 ... to its reciprocal
s Sum the above
cJ (3) J / the above
R Z Map each element of Z...
^ 2 ... to its square
.O Arithmetic mean of the above
@ 2 (4) Square root of the above
[ ) Wrap results (1), (2), (3), and (4) in a list
This is used as the input for the next iteration of the loop
h Take the first element of the result, implicit print
answered Apr 1 at 9:00
SokSok
4,187925
4,187925
add a comment |
add a comment |
$begingroup$
Japt v2.0a0 -g
, 42 38 bytes
â ÊÉ?ß[Ux²÷(V=UÊ)¬Ux÷V U×qV V÷Ux!÷1]:U
There's got to be a shorter way... This is a monstrosity! Saved 4 bytes thanks to @Shaggy!
Try it
$endgroup$
$begingroup$
38 bytes. But I agree, there must be a shorter way!
$endgroup$
– Shaggy
Apr 1 at 17:12
add a comment |
$begingroup$
Japt v2.0a0 -g
, 42 38 bytes
â ÊÉ?ß[Ux²÷(V=UÊ)¬Ux÷V U×qV V÷Ux!÷1]:U
There's got to be a shorter way... This is a monstrosity! Saved 4 bytes thanks to @Shaggy!
Try it
$endgroup$
$begingroup$
38 bytes. But I agree, there must be a shorter way!
$endgroup$
– Shaggy
Apr 1 at 17:12
add a comment |
$begingroup$
Japt v2.0a0 -g
, 42 38 bytes
â ÊÉ?ß[Ux²÷(V=UÊ)¬Ux÷V U×qV V÷Ux!÷1]:U
There's got to be a shorter way... This is a monstrosity! Saved 4 bytes thanks to @Shaggy!
Try it
$endgroup$
Japt v2.0a0 -g
, 42 38 bytes
â ÊÉ?ß[Ux²÷(V=UÊ)¬Ux÷V U×qV V÷Ux!÷1]:U
There's got to be a shorter way... This is a monstrosity! Saved 4 bytes thanks to @Shaggy!
Try it
edited Apr 1 at 18:36
answered Mar 31 at 21:10
Embodiment of IgnoranceEmbodiment of Ignorance
2,828127
2,828127
$begingroup$
38 bytes. But I agree, there must be a shorter way!
$endgroup$
– Shaggy
Apr 1 at 17:12
add a comment |
$begingroup$
38 bytes. But I agree, there must be a shorter way!
$endgroup$
– Shaggy
Apr 1 at 17:12
$begingroup$
38 bytes. But I agree, there must be a shorter way!
$endgroup$
– Shaggy
Apr 1 at 17:12
$begingroup$
38 bytes. But I agree, there must be a shorter way!
$endgroup$
– Shaggy
Apr 1 at 17:12
add a comment |
$begingroup$
C# (Visual C# Interactive Compiler), 177 bytes
double f(double[]g)=>g.All(c=>Math.Abs(c-g[0])<1e-9)?g[0]:f(new[]{g.Sum()/(z=g.Length),Math.Pow(g.Aggregate((a,b)=>a*b),1d/z),z/g.Sum(x=>1/x),Math.Sqrt(g.Sum(x=>x*x)/z)});int z;
Thanks to @KevinCruijjsen for pointing out that using floating point precision was causing problems! Would be 163 bytes if doubles were perfectly precise.
Try it online!
$endgroup$
$begingroup$
The last two test cases give aStackOverflowException
due to floating point precision. Instead ofc==g[0]
you could do something likeMath.Abs(c-g[0])<1e-9
. Try it online.
$endgroup$
– Kevin Cruijssen
Apr 1 at 9:01
$begingroup$
@KevinCruijssen Thanks, it's such a pain dealing with floating point numbers
$endgroup$
– Embodiment of Ignorance
Apr 2 at 15:32
add a comment |
$begingroup$
C# (Visual C# Interactive Compiler), 177 bytes
double f(double[]g)=>g.All(c=>Math.Abs(c-g[0])<1e-9)?g[0]:f(new[]{g.Sum()/(z=g.Length),Math.Pow(g.Aggregate((a,b)=>a*b),1d/z),z/g.Sum(x=>1/x),Math.Sqrt(g.Sum(x=>x*x)/z)});int z;
Thanks to @KevinCruijjsen for pointing out that using floating point precision was causing problems! Would be 163 bytes if doubles were perfectly precise.
Try it online!
$endgroup$
$begingroup$
The last two test cases give aStackOverflowException
due to floating point precision. Instead ofc==g[0]
you could do something likeMath.Abs(c-g[0])<1e-9
. Try it online.
$endgroup$
– Kevin Cruijssen
Apr 1 at 9:01
$begingroup$
@KevinCruijssen Thanks, it's such a pain dealing with floating point numbers
$endgroup$
– Embodiment of Ignorance
Apr 2 at 15:32
add a comment |
$begingroup$
C# (Visual C# Interactive Compiler), 177 bytes
double f(double[]g)=>g.All(c=>Math.Abs(c-g[0])<1e-9)?g[0]:f(new[]{g.Sum()/(z=g.Length),Math.Pow(g.Aggregate((a,b)=>a*b),1d/z),z/g.Sum(x=>1/x),Math.Sqrt(g.Sum(x=>x*x)/z)});int z;
Thanks to @KevinCruijjsen for pointing out that using floating point precision was causing problems! Would be 163 bytes if doubles were perfectly precise.
Try it online!
$endgroup$
C# (Visual C# Interactive Compiler), 177 bytes
double f(double[]g)=>g.All(c=>Math.Abs(c-g[0])<1e-9)?g[0]:f(new[]{g.Sum()/(z=g.Length),Math.Pow(g.Aggregate((a,b)=>a*b),1d/z),z/g.Sum(x=>1/x),Math.Sqrt(g.Sum(x=>x*x)/z)});int z;
Thanks to @KevinCruijjsen for pointing out that using floating point precision was causing problems! Would be 163 bytes if doubles were perfectly precise.
Try it online!
edited Apr 2 at 15:32
answered Mar 31 at 21:04
Embodiment of IgnoranceEmbodiment of Ignorance
2,828127
2,828127
$begingroup$
The last two test cases give aStackOverflowException
due to floating point precision. Instead ofc==g[0]
you could do something likeMath.Abs(c-g[0])<1e-9
. Try it online.
$endgroup$
– Kevin Cruijssen
Apr 1 at 9:01
$begingroup$
@KevinCruijssen Thanks, it's such a pain dealing with floating point numbers
$endgroup$
– Embodiment of Ignorance
Apr 2 at 15:32
add a comment |
$begingroup$
The last two test cases give aStackOverflowException
due to floating point precision. Instead ofc==g[0]
you could do something likeMath.Abs(c-g[0])<1e-9
. Try it online.
$endgroup$
– Kevin Cruijssen
Apr 1 at 9:01
$begingroup$
@KevinCruijssen Thanks, it's such a pain dealing with floating point numbers
$endgroup$
– Embodiment of Ignorance
Apr 2 at 15:32
$begingroup$
The last two test cases give a
StackOverflowException
due to floating point precision. Instead of c==g[0]
you could do something like Math.Abs(c-g[0])<1e-9
. Try it online.$endgroup$
– Kevin Cruijssen
Apr 1 at 9:01
$begingroup$
The last two test cases give a
StackOverflowException
due to floating point precision. Instead of c==g[0]
you could do something like Math.Abs(c-g[0])<1e-9
. Try it online.$endgroup$
– Kevin Cruijssen
Apr 1 at 9:01
$begingroup$
@KevinCruijssen Thanks, it's such a pain dealing with floating point numbers
$endgroup$
– Embodiment of Ignorance
Apr 2 at 15:32
$begingroup$
@KevinCruijssen Thanks, it's such a pain dealing with floating point numbers
$endgroup$
– Embodiment of Ignorance
Apr 2 at 15:32
add a comment |
$begingroup$
x86 machine code (SIMD 4x float using 128-bit SSE1&AVX) 94 bytes
x86 machine code (SIMD 4x double using 256-bit AVX) 123 bytes
float
passes the test cases in the question, but with a loop-exit threshold small enough to make that happen, it's easy for it to get stuck in an infinite loop with random inputs.
SSE1 packed-single-precision instructions are 3 bytes long, but SSE2 and simple AVX instructions are 4 bytes long. (Scalar-single instructions like sqrtss
are also 4 bytes long, which is why I use sqrtps
even though I only care about the low element. It's not even slower than sqrtss on modern hardware). I used AVX for non-destructive destination to save 2 bytes vs. movaps+op.
In the double version we can still do a couple movlhps
to copy 64-bit chunks (because often we only care about the low element of a horizontal sum). Horizontal sum of a 256-bit SIMD vector also requires an extra vextractf128
to get the high half, vs. the slow but small 2x haddps
strategy for float. The double
version also needs 2x 8-byte constants, instead of 2x 4-byte. Overall it comes out at close to 4/3 the size of the float
version.
mean(a,b) = mean(a,a,b,b)
for all 4 of these means, so we can simply duplicate the input up to 4 elements and never have to implement length=2. Thus we can hardcode geometric mean as 4th-root = sqrt(sqrt), for example. And we only need one FP constant, 4.0
.
We have a single SIMD vector of all 4 [a_i, b_i, c_i, d_i]
. From that, we calculate the 4 means as scalars in separate registers, and shuffle them back together for the next iteration. (Horizontal operations on SIMD vectors are inconvenient, but we need to do the same thing for all 4 elements in enough cases that it balances out. I started on an x87 version of this, but it was getting very long and not fun.)
The loop-exit condition of }while(quadratic - harmonic > 4e-5)
(or a smaller constant for double
) is based on @RobinRyder's R answer, and Kevin Cruijssen's Java answer: quadratic mean is always the largest magnitude, and harmonic mean is always the smallest (ignoring rounding errors). So we can check the delta between those two to detect convergence. We return the arithmetic mean as the scalar result. It's usually between those two and is probably the least susceptible to rounding errors.
Float version: callable as float meanmean_float_avx(__m128);
with the arg and return value in xmm0. (So x86-64 System V, or Windows x64 vectorcall, but not x64 fastcall.) Or declare the return-type as __m128
so you can get at the quadratic and harmonic mean for testing.
Letting this take 2 separate float
args in xmm0 and xmm1 would cost 1 extra byte: we'd need a shufps
with an imm8 (instead of just unpcklps xmm0,xmm0
) to shuffle together and duplicate 2 inputs.
40 address align 32
41 code bytes global meanmean_float_avx
42 meanmean_float_avx:
43 00000000 B9[52000000] mov ecx, .arith_mean ; allows 2-byte call reg, and a base for loading constants
44 00000005 C4E2791861FC vbroadcastss xmm4, [rcx-4] ; float 4.0
45
46 ;; mean(a,b) = mean(a,b,a,b) for all 4 types of mean
47 ;; so we only ever have to do the length=4 case
48 0000000B 0F14C0 unpcklps xmm0,xmm0 ; [b,a] => [b,b,a,a]
49
50 ; do{ ... } while(quadratic - harmonic > threshold);
51 .loop:
52 ;;; XMM3 = geometric mean: not based on addition. (Transform to log would be hard. AVX512ER has exp with 23-bit accuracy, but not log. vgetexp = floor(lofg2(x)), so that's no good.)
53 ;; sqrt once *first*, making magnitudes closer to 1.0 to reduce rounding error. Numbers are all positive so this is safe.
54 ;; both sqrts first was better behaved, I think.
55 0000000E 0F51D8 sqrtps xmm3, xmm0 ; xmm3 = 4th root(x)
56 00000011 F30F16EB movshdup xmm5, xmm3 ; bring odd elements down to even
57 00000015 0F59EB mulps xmm5, xmm3
58 00000018 0F12DD movhlps xmm3, xmm5 ; high half -> low
59 0000001B 0F59DD mulps xmm3, xmm5 ; xmm3[0] = hproduct(sqrt(xmm))
60 ; sqrtps xmm3, xmm3 ; sqrt(hprod(sqrt)) = 4th root(hprod)
61 ; common final step done after interleaving with quadratic mean
62
63 ;;; XMM2 = quadratic mean = max of the means
64 0000001E C5F859E8 vmulps xmm5, xmm0,xmm0
65 00000022 FFD1 call rcx ; arith mean of squares
66 00000024 0F14EB unpcklps xmm5, xmm3 ; [quad^2, geo^2, ?, ?]
67 00000027 0F51D5 sqrtps xmm2, xmm5 ; [quad, geo, ?, ?]
68
69 ;;; XMM1 = harmonic mean = min of the means
70 0000002A C5D85EE8 vdivps xmm5, xmm4, xmm0 ; 4/x
71 0000002E FFD1 call rcx ; arithmetic mean (under inversion)
72 00000030 C5D85ECD vdivps xmm1, xmm4, xmm5 ; 4/. (the factor of 4 cancels out)
73
74 ;;; XMM5 = arithmetic mean
75 00000034 0F28E8 movaps xmm5, xmm0
76 00000037 FFD1 call rcx
77
78 00000039 0F14E9 unpcklps xmm5, xmm1 ; [arith, harm, ?,?]
79 0000003C C5D014C2 vunpcklps xmm0, xmm5,xmm2 ; x = [arith, harm, quad, geo]
80
81 00000040 0F5CD1 subps xmm2, xmm1 ; largest - smallest mean: guaranteed non-negative
82 00000043 0F2E51F8 ucomiss xmm2, [rcx-8] ; quad-harm > convergence_threshold
83 00000047 73C5 jae .loop
84
85 ; return with the arithmetic mean in the low element of xmm0 = scalar return value
86 00000049 C3 ret
87
88 ;;; "constant pool" between the main function and the helper, like ARM literal pools
89 0000004A ACC52738 .fpconst_threshold: dd 4e-5 ; 4.3e-5 is the highest we can go and still pass the main test cases
90 0000004E 00008040 .fpconst_4: dd 4.0
91 .arith_mean: ; returns XMM5 = hsum(xmm5)/4.
92 00000052 C5D37CED vhaddps xmm5, xmm5 ; slow but small
93 00000056 C5D37CED vhaddps xmm5, xmm5
94 0000005A 0F5EEC divps xmm5, xmm4 ; divide before/after summing doesn't matter mathematically or numerically; divisor is a power of 2
95 0000005D C3 ret
96 0000005E 5E000000 .size: dd $ - meanmean_float_avx
0x5e = 94 bytes
(NASM listing created with nasm -felf64 mean-mean.asm -l/dev/stdout | cut -b -34,$((34+6))-
. Strip the listing part and recover the source with cut -b 34- > mean-mean.asm
)
SIMD horizontal sum and divide by 4 (i.e. arithmetic mean) is implemented in a separate function that we call
(with a function pointer to amortize the cost of the address). With 4/x
before/after, or x^2
before and sqrt after, we get the harmonic mean and quadratic mean. (It was painful to write these div
instructions instead of multiplying by an exactly-representable 0.25
.)
Geometric mean is implemented separately with multiply and chained sqrt. Or with one sqrt first to reduce exponent magnitude and maybe help numerical precision. log is not available, only floor(log2(x))
via AVX512 vgetexpps/pd
. Exp is sort of available via AVX512ER (Xeon Phi only), but with only 2^-23 precision.
Mixing 128-bit AVX instructions and legacy SSE is not a performance problem. Mixing 256-bit AVX with legacy SSE can be on Haswell, but on Skylake it just potentially creates a potential false dependency for SSE instructions. I think my double
version avoids any unnecessary loop-carried dep chains, and bottlenecks
on div/sqrt latency/throughput.
Double version:
108 global meanmean_double_avx
109 meanmean_double_avx:
110 00000080 B9[E8000000] mov ecx, .arith_mean
111 00000085 C4E27D1961F8 vbroadcastsd ymm4, [rcx-8] ; float 4.0
112
113 ;; mean(a,b) = mean(a,b,a,b) for all 4 types of mean
114 ;; so we only ever have to do the length=4 case
115 0000008B C4E37D18C001 vinsertf128 ymm0, ymm0, xmm0, 1 ; [b,a] => [b,a,b,a]
116
117 .loop:
118 ;;; XMM3 = geometric mean: not based on addition.
119 00000091 C5FD51D8 vsqrtpd ymm3, ymm0 ; sqrt first to get magnitude closer to 1.0 for better(?) numerical precision
120 00000095 C4E37D19DD01 vextractf128 xmm5, ymm3, 1 ; extract high lane
121 0000009B C5D159EB vmulpd xmm5, xmm3
122 0000009F 0F12DD movhlps xmm3, xmm5 ; extract high half
123 000000A2 F20F59DD mulsd xmm3, xmm5 ; xmm3 = hproduct(sqrt(xmm0))
124 ; sqrtsd xmm3, xmm3 ; xmm3 = 4th root = geomean(xmm0) ;deferred until quadratic
125
126 ;;; XMM2 = quadratic mean = max of the means
127 000000A6 C5FD59E8 vmulpd ymm5, ymm0,ymm0
128 000000AA FFD1 call rcx ; arith mean of squares
129 000000AC 0F16EB movlhps xmm5, xmm3 ; [quad^2, geo^2]
130 000000AF 660F51D5 sqrtpd xmm2, xmm5 ; [quad , geo]
131
132 ;;; XMM1 = harmonic mean = min of the means
133 000000B3 C5DD5EE8 vdivpd ymm5, ymm4, ymm0 ; 4/x
134 000000B7 FFD1 call rcx ; arithmetic mean under inversion
135 000000B9 C5DB5ECD vdivsd xmm1, xmm4, xmm5 ; 4/. (the factor of 4 cancels out)
136
137 ;;; XMM5 = arithmetic mean
138 000000BD C5FC28E8 vmovaps ymm5, ymm0
139 000000C1 FFD1 call rcx
140
141 000000C3 0F16E9 movlhps xmm5, xmm1 ; [arith, harm]
142 000000C6 C4E35518C201 vinsertf128 ymm0, ymm5, xmm2, 1 ; x = [arith, harm, quad, geo]
143
144 000000CC C5EB5CD1 vsubsd xmm2, xmm1 ; largest - smallest mean: guaranteed non-negative
145 000000D0 660F2E51F0 ucomisd xmm2, [rcx-16] ; quad - harm > threshold
146 000000D5 77BA ja .loop
147
148 ; vzeroupper ; not needed for correctness, only performance
149 ; return with the arithmetic mean in the low element of xmm0 = scalar return value
150 000000D7 C3 ret
151
152 ; "literal pool" between the function
153 000000D8 95D626E80B2E113E .fpconst_threshold: dq 1e-9
154 000000E0 0000000000001040 .fpconst_4: dq 4.0 ; TODO: golf these zeros? vpbroadcastb and convert?
155 .arith_mean: ; returns YMM5 = hsum(ymm5)/4.
156 000000E8 C4E37D19EF01 vextractf128 xmm7, ymm5, 1
157 000000EE C5D158EF vaddpd xmm5, xmm7
158 000000F2 C5D17CED vhaddpd xmm5, xmm5 ; slow but small
159 000000F6 C5D35EEC vdivsd xmm5, xmm4 ; only low element matters
160 000000FA C3 ret
161 000000FB 7B000000 .size: dd $ - meanmean_double_avx
0x7b = 123 bytes
C test harness
#include <immintrin.h>
#include <stdio.h>
#include <math.h>
static const struct ab_avg {
double a,b;
double mean;
} testcases[] = {
{1, 1, 1},
{1, 2, 1.45568889},
{100, 200, 145.568889},
{2.71, 3.14, 2.92103713},
{0.57, 1.78, 1.0848205},
{1.61, 2.41, 1.98965438},
{0.01, 100, 6.7483058},
};
// see asm comments for order of arith, harm, quad, geo
__m128 meanmean_float_avx(__m128); // or float ...
__m256d meanmean_double_avx(__m128d); // or double ...
int main(void) {
int len = sizeof(testcases) / sizeof(testcases[0]);
for(int i=0 ; i<len ; i++) {
const struct ab_avg *p = &testcases[i];
#if 1
__m128 arg = _mm_set_ps(0,0, p->b, p->a);
double res = meanmean_float_avx(arg)[0];
#else
__m128d arg = _mm_loadu_pd(&p->a);
double res = meanmean_double_avx(arg)[0];
#endif
double allowed_diff = (p->b - p->a) / 100000.0;
double delta = fabs(p->mean - res);
if (delta > 1e-3 || delta > allowed_diff) {
printf("%f %f => %.9f but we got %.9f. delta = %g allowed=%gn",
p->a, p->b, p->mean, res, p->mean - res, allowed_diff);
}
}
while(1) {
double a = drand48(), b = drand48(); // range= [0..1)
if (a>b) {
double tmp=a;
a=b;
b=tmp; // sorted
}
// a *= 0.00000001;
// b *= 123156;
// a += 1<<11; b += (1<<12)+1; // float version gets stuck inflooping on 2048.04, 4097.18 at fpthreshold = 4e-5
// a *= 1<<11 ; b *= 1<<11; // scaling to large magnitude makes sum of squares loses more precision
//a += 1<<11; b+= 1<<11; // adding to large magnitude is hard for everything, catastrophic cancellation
#if 1
printf("testing float %g, %gn", a, b);
__m128 arg = _mm_set_ps(0,0, b, a);
__m128 res = meanmean_float_avx(arg);
double quad = res[2], harm = res[1]; // same order as double... for now
#else
printf("testing double %g, %gn", a, b);
__m128d arg = _mm_set_pd(b, a);
__m256d res = meanmean_double_avx(arg);
double quad = res[2], harm = res[1];
#endif
double delta = fabs(quad - harm);
double allowed_diff = (b - a) / 100000.0; // calculated in double even for the float case.
// TODO: use the double res as a reference for float res
// instead of just checking quadratic vs. harmonic mean
if (delta > 1e-3 || delta > allowed_diff) {
printf("%g %g we got q=%g, h=%g, a=%g. delta = %g, allowed=%gn",
a, b, quad, harm, res[0], quad-harm, allowed_diff);
}
}
}
Build with:
nasm -felf64 mean-mean.asm &&
gcc -no-pie -fno-pie -g -O2 -march=native mean-mean.c mean-mean.o
Obviously you need a CPU with AVX support, or an emulator like Intel SDE. To compile on a host without native AVX support, use -march=sandybridge
or -mavx
Run: passes the hard-coded test cases, but for the float version, random test cases often fail the (b-a)/10000
threshold set in the question.
$ ./a.out
(note: empty output before the first "testing float" means clean pass on the constant test cases)
testing float 3.90799e-14, 0.000985395
3.90799e-14 0.000985395 we got q=3.20062e-10, h=3.58723e-05, a=2.50934e-05. delta = -3.5872e-05, allowed=9.85395e-09
testing float 0.041631, 0.176643
testing float 0.0913306, 0.364602
testing float 0.0922976, 0.487217
testing float 0.454433, 0.52675
0.454433 0.52675 we got q=0.48992, h=0.489927, a=0.489925. delta = -6.79493e-06, allowed=7.23169e-07
testing float 0.233178, 0.831292
testing float 0.56806, 0.931731
testing float 0.0508319, 0.556094
testing float 0.0189148, 0.767051
0.0189148 0.767051 we got q=0.210471, h=0.210484, a=0.21048. delta = -1.37389e-05, allowed=7.48136e-06
testing float 0.25236, 0.298197
0.25236 0.298197 we got q=0.274796, h=0.274803, a=0.274801. delta = -6.19888e-06, allowed=4.58374e-07
testing float 0.531557, 0.875981
testing float 0.515431, 0.920261
testing float 0.18842, 0.810429
testing float 0.570614, 0.886314
testing float 0.0767746, 0.815274
testing float 0.118352, 0.984891
0.118352 0.984891 we got q=0.427845, h=0.427872, a=0.427863. delta = -2.66135e-05, allowed=8.66539e-06
testing float 0.784484, 0.893906
0.784484 0.893906 we got q=0.838297, h=0.838304, a=0.838302. delta = -7.09295e-06, allowed=1.09422e-06
FP errors are enough that quad-harm comes out less than zero for some inputs.
Or with a += 1<<11; b += (1<<12)+1;
uncommented:
testing float 2048, 4097
testing float 2048.04, 4097.18
^C (stuck in an infinite loop).
None of these problems happen with double
. Comment out the printf
before each test to see that the output is empty (nothing from the if(delta too high)
block).
TODO: use the double
version as a reference for the float
version, instead of just looking at how they converging with quad-harm.
$endgroup$
add a comment |
$begingroup$
x86 machine code (SIMD 4x float using 128-bit SSE1&AVX) 94 bytes
x86 machine code (SIMD 4x double using 256-bit AVX) 123 bytes
float
passes the test cases in the question, but with a loop-exit threshold small enough to make that happen, it's easy for it to get stuck in an infinite loop with random inputs.
SSE1 packed-single-precision instructions are 3 bytes long, but SSE2 and simple AVX instructions are 4 bytes long. (Scalar-single instructions like sqrtss
are also 4 bytes long, which is why I use sqrtps
even though I only care about the low element. It's not even slower than sqrtss on modern hardware). I used AVX for non-destructive destination to save 2 bytes vs. movaps+op.
In the double version we can still do a couple movlhps
to copy 64-bit chunks (because often we only care about the low element of a horizontal sum). Horizontal sum of a 256-bit SIMD vector also requires an extra vextractf128
to get the high half, vs. the slow but small 2x haddps
strategy for float. The double
version also needs 2x 8-byte constants, instead of 2x 4-byte. Overall it comes out at close to 4/3 the size of the float
version.
mean(a,b) = mean(a,a,b,b)
for all 4 of these means, so we can simply duplicate the input up to 4 elements and never have to implement length=2. Thus we can hardcode geometric mean as 4th-root = sqrt(sqrt), for example. And we only need one FP constant, 4.0
.
We have a single SIMD vector of all 4 [a_i, b_i, c_i, d_i]
. From that, we calculate the 4 means as scalars in separate registers, and shuffle them back together for the next iteration. (Horizontal operations on SIMD vectors are inconvenient, but we need to do the same thing for all 4 elements in enough cases that it balances out. I started on an x87 version of this, but it was getting very long and not fun.)
The loop-exit condition of }while(quadratic - harmonic > 4e-5)
(or a smaller constant for double
) is based on @RobinRyder's R answer, and Kevin Cruijssen's Java answer: quadratic mean is always the largest magnitude, and harmonic mean is always the smallest (ignoring rounding errors). So we can check the delta between those two to detect convergence. We return the arithmetic mean as the scalar result. It's usually between those two and is probably the least susceptible to rounding errors.
Float version: callable as float meanmean_float_avx(__m128);
with the arg and return value in xmm0. (So x86-64 System V, or Windows x64 vectorcall, but not x64 fastcall.) Or declare the return-type as __m128
so you can get at the quadratic and harmonic mean for testing.
Letting this take 2 separate float
args in xmm0 and xmm1 would cost 1 extra byte: we'd need a shufps
with an imm8 (instead of just unpcklps xmm0,xmm0
) to shuffle together and duplicate 2 inputs.
40 address align 32
41 code bytes global meanmean_float_avx
42 meanmean_float_avx:
43 00000000 B9[52000000] mov ecx, .arith_mean ; allows 2-byte call reg, and a base for loading constants
44 00000005 C4E2791861FC vbroadcastss xmm4, [rcx-4] ; float 4.0
45
46 ;; mean(a,b) = mean(a,b,a,b) for all 4 types of mean
47 ;; so we only ever have to do the length=4 case
48 0000000B 0F14C0 unpcklps xmm0,xmm0 ; [b,a] => [b,b,a,a]
49
50 ; do{ ... } while(quadratic - harmonic > threshold);
51 .loop:
52 ;;; XMM3 = geometric mean: not based on addition. (Transform to log would be hard. AVX512ER has exp with 23-bit accuracy, but not log. vgetexp = floor(lofg2(x)), so that's no good.)
53 ;; sqrt once *first*, making magnitudes closer to 1.0 to reduce rounding error. Numbers are all positive so this is safe.
54 ;; both sqrts first was better behaved, I think.
55 0000000E 0F51D8 sqrtps xmm3, xmm0 ; xmm3 = 4th root(x)
56 00000011 F30F16EB movshdup xmm5, xmm3 ; bring odd elements down to even
57 00000015 0F59EB mulps xmm5, xmm3
58 00000018 0F12DD movhlps xmm3, xmm5 ; high half -> low
59 0000001B 0F59DD mulps xmm3, xmm5 ; xmm3[0] = hproduct(sqrt(xmm))
60 ; sqrtps xmm3, xmm3 ; sqrt(hprod(sqrt)) = 4th root(hprod)
61 ; common final step done after interleaving with quadratic mean
62
63 ;;; XMM2 = quadratic mean = max of the means
64 0000001E C5F859E8 vmulps xmm5, xmm0,xmm0
65 00000022 FFD1 call rcx ; arith mean of squares
66 00000024 0F14EB unpcklps xmm5, xmm3 ; [quad^2, geo^2, ?, ?]
67 00000027 0F51D5 sqrtps xmm2, xmm5 ; [quad, geo, ?, ?]
68
69 ;;; XMM1 = harmonic mean = min of the means
70 0000002A C5D85EE8 vdivps xmm5, xmm4, xmm0 ; 4/x
71 0000002E FFD1 call rcx ; arithmetic mean (under inversion)
72 00000030 C5D85ECD vdivps xmm1, xmm4, xmm5 ; 4/. (the factor of 4 cancels out)
73
74 ;;; XMM5 = arithmetic mean
75 00000034 0F28E8 movaps xmm5, xmm0
76 00000037 FFD1 call rcx
77
78 00000039 0F14E9 unpcklps xmm5, xmm1 ; [arith, harm, ?,?]
79 0000003C C5D014C2 vunpcklps xmm0, xmm5,xmm2 ; x = [arith, harm, quad, geo]
80
81 00000040 0F5CD1 subps xmm2, xmm1 ; largest - smallest mean: guaranteed non-negative
82 00000043 0F2E51F8 ucomiss xmm2, [rcx-8] ; quad-harm > convergence_threshold
83 00000047 73C5 jae .loop
84
85 ; return with the arithmetic mean in the low element of xmm0 = scalar return value
86 00000049 C3 ret
87
88 ;;; "constant pool" between the main function and the helper, like ARM literal pools
89 0000004A ACC52738 .fpconst_threshold: dd 4e-5 ; 4.3e-5 is the highest we can go and still pass the main test cases
90 0000004E 00008040 .fpconst_4: dd 4.0
91 .arith_mean: ; returns XMM5 = hsum(xmm5)/4.
92 00000052 C5D37CED vhaddps xmm5, xmm5 ; slow but small
93 00000056 C5D37CED vhaddps xmm5, xmm5
94 0000005A 0F5EEC divps xmm5, xmm4 ; divide before/after summing doesn't matter mathematically or numerically; divisor is a power of 2
95 0000005D C3 ret
96 0000005E 5E000000 .size: dd $ - meanmean_float_avx
0x5e = 94 bytes
(NASM listing created with nasm -felf64 mean-mean.asm -l/dev/stdout | cut -b -34,$((34+6))-
. Strip the listing part and recover the source with cut -b 34- > mean-mean.asm
)
SIMD horizontal sum and divide by 4 (i.e. arithmetic mean) is implemented in a separate function that we call
(with a function pointer to amortize the cost of the address). With 4/x
before/after, or x^2
before and sqrt after, we get the harmonic mean and quadratic mean. (It was painful to write these div
instructions instead of multiplying by an exactly-representable 0.25
.)
Geometric mean is implemented separately with multiply and chained sqrt. Or with one sqrt first to reduce exponent magnitude and maybe help numerical precision. log is not available, only floor(log2(x))
via AVX512 vgetexpps/pd
. Exp is sort of available via AVX512ER (Xeon Phi only), but with only 2^-23 precision.
Mixing 128-bit AVX instructions and legacy SSE is not a performance problem. Mixing 256-bit AVX with legacy SSE can be on Haswell, but on Skylake it just potentially creates a potential false dependency for SSE instructions. I think my double
version avoids any unnecessary loop-carried dep chains, and bottlenecks
on div/sqrt latency/throughput.
Double version:
108 global meanmean_double_avx
109 meanmean_double_avx:
110 00000080 B9[E8000000] mov ecx, .arith_mean
111 00000085 C4E27D1961F8 vbroadcastsd ymm4, [rcx-8] ; float 4.0
112
113 ;; mean(a,b) = mean(a,b,a,b) for all 4 types of mean
114 ;; so we only ever have to do the length=4 case
115 0000008B C4E37D18C001 vinsertf128 ymm0, ymm0, xmm0, 1 ; [b,a] => [b,a,b,a]
116
117 .loop:
118 ;;; XMM3 = geometric mean: not based on addition.
119 00000091 C5FD51D8 vsqrtpd ymm3, ymm0 ; sqrt first to get magnitude closer to 1.0 for better(?) numerical precision
120 00000095 C4E37D19DD01 vextractf128 xmm5, ymm3, 1 ; extract high lane
121 0000009B C5D159EB vmulpd xmm5, xmm3
122 0000009F 0F12DD movhlps xmm3, xmm5 ; extract high half
123 000000A2 F20F59DD mulsd xmm3, xmm5 ; xmm3 = hproduct(sqrt(xmm0))
124 ; sqrtsd xmm3, xmm3 ; xmm3 = 4th root = geomean(xmm0) ;deferred until quadratic
125
126 ;;; XMM2 = quadratic mean = max of the means
127 000000A6 C5FD59E8 vmulpd ymm5, ymm0,ymm0
128 000000AA FFD1 call rcx ; arith mean of squares
129 000000AC 0F16EB movlhps xmm5, xmm3 ; [quad^2, geo^2]
130 000000AF 660F51D5 sqrtpd xmm2, xmm5 ; [quad , geo]
131
132 ;;; XMM1 = harmonic mean = min of the means
133 000000B3 C5DD5EE8 vdivpd ymm5, ymm4, ymm0 ; 4/x
134 000000B7 FFD1 call rcx ; arithmetic mean under inversion
135 000000B9 C5DB5ECD vdivsd xmm1, xmm4, xmm5 ; 4/. (the factor of 4 cancels out)
136
137 ;;; XMM5 = arithmetic mean
138 000000BD C5FC28E8 vmovaps ymm5, ymm0
139 000000C1 FFD1 call rcx
140
141 000000C3 0F16E9 movlhps xmm5, xmm1 ; [arith, harm]
142 000000C6 C4E35518C201 vinsertf128 ymm0, ymm5, xmm2, 1 ; x = [arith, harm, quad, geo]
143
144 000000CC C5EB5CD1 vsubsd xmm2, xmm1 ; largest - smallest mean: guaranteed non-negative
145 000000D0 660F2E51F0 ucomisd xmm2, [rcx-16] ; quad - harm > threshold
146 000000D5 77BA ja .loop
147
148 ; vzeroupper ; not needed for correctness, only performance
149 ; return with the arithmetic mean in the low element of xmm0 = scalar return value
150 000000D7 C3 ret
151
152 ; "literal pool" between the function
153 000000D8 95D626E80B2E113E .fpconst_threshold: dq 1e-9
154 000000E0 0000000000001040 .fpconst_4: dq 4.0 ; TODO: golf these zeros? vpbroadcastb and convert?
155 .arith_mean: ; returns YMM5 = hsum(ymm5)/4.
156 000000E8 C4E37D19EF01 vextractf128 xmm7, ymm5, 1
157 000000EE C5D158EF vaddpd xmm5, xmm7
158 000000F2 C5D17CED vhaddpd xmm5, xmm5 ; slow but small
159 000000F6 C5D35EEC vdivsd xmm5, xmm4 ; only low element matters
160 000000FA C3 ret
161 000000FB 7B000000 .size: dd $ - meanmean_double_avx
0x7b = 123 bytes
C test harness
#include <immintrin.h>
#include <stdio.h>
#include <math.h>
static const struct ab_avg {
double a,b;
double mean;
} testcases[] = {
{1, 1, 1},
{1, 2, 1.45568889},
{100, 200, 145.568889},
{2.71, 3.14, 2.92103713},
{0.57, 1.78, 1.0848205},
{1.61, 2.41, 1.98965438},
{0.01, 100, 6.7483058},
};
// see asm comments for order of arith, harm, quad, geo
__m128 meanmean_float_avx(__m128); // or float ...
__m256d meanmean_double_avx(__m128d); // or double ...
int main(void) {
int len = sizeof(testcases) / sizeof(testcases[0]);
for(int i=0 ; i<len ; i++) {
const struct ab_avg *p = &testcases[i];
#if 1
__m128 arg = _mm_set_ps(0,0, p->b, p->a);
double res = meanmean_float_avx(arg)[0];
#else
__m128d arg = _mm_loadu_pd(&p->a);
double res = meanmean_double_avx(arg)[0];
#endif
double allowed_diff = (p->b - p->a) / 100000.0;
double delta = fabs(p->mean - res);
if (delta > 1e-3 || delta > allowed_diff) {
printf("%f %f => %.9f but we got %.9f. delta = %g allowed=%gn",
p->a, p->b, p->mean, res, p->mean - res, allowed_diff);
}
}
while(1) {
double a = drand48(), b = drand48(); // range= [0..1)
if (a>b) {
double tmp=a;
a=b;
b=tmp; // sorted
}
// a *= 0.00000001;
// b *= 123156;
// a += 1<<11; b += (1<<12)+1; // float version gets stuck inflooping on 2048.04, 4097.18 at fpthreshold = 4e-5
// a *= 1<<11 ; b *= 1<<11; // scaling to large magnitude makes sum of squares loses more precision
//a += 1<<11; b+= 1<<11; // adding to large magnitude is hard for everything, catastrophic cancellation
#if 1
printf("testing float %g, %gn", a, b);
__m128 arg = _mm_set_ps(0,0, b, a);
__m128 res = meanmean_float_avx(arg);
double quad = res[2], harm = res[1]; // same order as double... for now
#else
printf("testing double %g, %gn", a, b);
__m128d arg = _mm_set_pd(b, a);
__m256d res = meanmean_double_avx(arg);
double quad = res[2], harm = res[1];
#endif
double delta = fabs(quad - harm);
double allowed_diff = (b - a) / 100000.0; // calculated in double even for the float case.
// TODO: use the double res as a reference for float res
// instead of just checking quadratic vs. harmonic mean
if (delta > 1e-3 || delta > allowed_diff) {
printf("%g %g we got q=%g, h=%g, a=%g. delta = %g, allowed=%gn",
a, b, quad, harm, res[0], quad-harm, allowed_diff);
}
}
}
Build with:
nasm -felf64 mean-mean.asm &&
gcc -no-pie -fno-pie -g -O2 -march=native mean-mean.c mean-mean.o
Obviously you need a CPU with AVX support, or an emulator like Intel SDE. To compile on a host without native AVX support, use -march=sandybridge
or -mavx
Run: passes the hard-coded test cases, but for the float version, random test cases often fail the (b-a)/10000
threshold set in the question.
$ ./a.out
(note: empty output before the first "testing float" means clean pass on the constant test cases)
testing float 3.90799e-14, 0.000985395
3.90799e-14 0.000985395 we got q=3.20062e-10, h=3.58723e-05, a=2.50934e-05. delta = -3.5872e-05, allowed=9.85395e-09
testing float 0.041631, 0.176643
testing float 0.0913306, 0.364602
testing float 0.0922976, 0.487217
testing float 0.454433, 0.52675
0.454433 0.52675 we got q=0.48992, h=0.489927, a=0.489925. delta = -6.79493e-06, allowed=7.23169e-07
testing float 0.233178, 0.831292
testing float 0.56806, 0.931731
testing float 0.0508319, 0.556094
testing float 0.0189148, 0.767051
0.0189148 0.767051 we got q=0.210471, h=0.210484, a=0.21048. delta = -1.37389e-05, allowed=7.48136e-06
testing float 0.25236, 0.298197
0.25236 0.298197 we got q=0.274796, h=0.274803, a=0.274801. delta = -6.19888e-06, allowed=4.58374e-07
testing float 0.531557, 0.875981
testing float 0.515431, 0.920261
testing float 0.18842, 0.810429
testing float 0.570614, 0.886314
testing float 0.0767746, 0.815274
testing float 0.118352, 0.984891
0.118352 0.984891 we got q=0.427845, h=0.427872, a=0.427863. delta = -2.66135e-05, allowed=8.66539e-06
testing float 0.784484, 0.893906
0.784484 0.893906 we got q=0.838297, h=0.838304, a=0.838302. delta = -7.09295e-06, allowed=1.09422e-06
FP errors are enough that quad-harm comes out less than zero for some inputs.
Or with a += 1<<11; b += (1<<12)+1;
uncommented:
testing float 2048, 4097
testing float 2048.04, 4097.18
^C (stuck in an infinite loop).
None of these problems happen with double
. Comment out the printf
before each test to see that the output is empty (nothing from the if(delta too high)
block).
TODO: use the double
version as a reference for the float
version, instead of just looking at how they converging with quad-harm.
$endgroup$
add a comment |
$begingroup$
x86 machine code (SIMD 4x float using 128-bit SSE1&AVX) 94 bytes
x86 machine code (SIMD 4x double using 256-bit AVX) 123 bytes
float
passes the test cases in the question, but with a loop-exit threshold small enough to make that happen, it's easy for it to get stuck in an infinite loop with random inputs.
SSE1 packed-single-precision instructions are 3 bytes long, but SSE2 and simple AVX instructions are 4 bytes long. (Scalar-single instructions like sqrtss
are also 4 bytes long, which is why I use sqrtps
even though I only care about the low element. It's not even slower than sqrtss on modern hardware). I used AVX for non-destructive destination to save 2 bytes vs. movaps+op.
In the double version we can still do a couple movlhps
to copy 64-bit chunks (because often we only care about the low element of a horizontal sum). Horizontal sum of a 256-bit SIMD vector also requires an extra vextractf128
to get the high half, vs. the slow but small 2x haddps
strategy for float. The double
version also needs 2x 8-byte constants, instead of 2x 4-byte. Overall it comes out at close to 4/3 the size of the float
version.
mean(a,b) = mean(a,a,b,b)
for all 4 of these means, so we can simply duplicate the input up to 4 elements and never have to implement length=2. Thus we can hardcode geometric mean as 4th-root = sqrt(sqrt), for example. And we only need one FP constant, 4.0
.
We have a single SIMD vector of all 4 [a_i, b_i, c_i, d_i]
. From that, we calculate the 4 means as scalars in separate registers, and shuffle them back together for the next iteration. (Horizontal operations on SIMD vectors are inconvenient, but we need to do the same thing for all 4 elements in enough cases that it balances out. I started on an x87 version of this, but it was getting very long and not fun.)
The loop-exit condition of }while(quadratic - harmonic > 4e-5)
(or a smaller constant for double
) is based on @RobinRyder's R answer, and Kevin Cruijssen's Java answer: quadratic mean is always the largest magnitude, and harmonic mean is always the smallest (ignoring rounding errors). So we can check the delta between those two to detect convergence. We return the arithmetic mean as the scalar result. It's usually between those two and is probably the least susceptible to rounding errors.
Float version: callable as float meanmean_float_avx(__m128);
with the arg and return value in xmm0. (So x86-64 System V, or Windows x64 vectorcall, but not x64 fastcall.) Or declare the return-type as __m128
so you can get at the quadratic and harmonic mean for testing.
Letting this take 2 separate float
args in xmm0 and xmm1 would cost 1 extra byte: we'd need a shufps
with an imm8 (instead of just unpcklps xmm0,xmm0
) to shuffle together and duplicate 2 inputs.
40 address align 32
41 code bytes global meanmean_float_avx
42 meanmean_float_avx:
43 00000000 B9[52000000] mov ecx, .arith_mean ; allows 2-byte call reg, and a base for loading constants
44 00000005 C4E2791861FC vbroadcastss xmm4, [rcx-4] ; float 4.0
45
46 ;; mean(a,b) = mean(a,b,a,b) for all 4 types of mean
47 ;; so we only ever have to do the length=4 case
48 0000000B 0F14C0 unpcklps xmm0,xmm0 ; [b,a] => [b,b,a,a]
49
50 ; do{ ... } while(quadratic - harmonic > threshold);
51 .loop:
52 ;;; XMM3 = geometric mean: not based on addition. (Transform to log would be hard. AVX512ER has exp with 23-bit accuracy, but not log. vgetexp = floor(lofg2(x)), so that's no good.)
53 ;; sqrt once *first*, making magnitudes closer to 1.0 to reduce rounding error. Numbers are all positive so this is safe.
54 ;; both sqrts first was better behaved, I think.
55 0000000E 0F51D8 sqrtps xmm3, xmm0 ; xmm3 = 4th root(x)
56 00000011 F30F16EB movshdup xmm5, xmm3 ; bring odd elements down to even
57 00000015 0F59EB mulps xmm5, xmm3
58 00000018 0F12DD movhlps xmm3, xmm5 ; high half -> low
59 0000001B 0F59DD mulps xmm3, xmm5 ; xmm3[0] = hproduct(sqrt(xmm))
60 ; sqrtps xmm3, xmm3 ; sqrt(hprod(sqrt)) = 4th root(hprod)
61 ; common final step done after interleaving with quadratic mean
62
63 ;;; XMM2 = quadratic mean = max of the means
64 0000001E C5F859E8 vmulps xmm5, xmm0,xmm0
65 00000022 FFD1 call rcx ; arith mean of squares
66 00000024 0F14EB unpcklps xmm5, xmm3 ; [quad^2, geo^2, ?, ?]
67 00000027 0F51D5 sqrtps xmm2, xmm5 ; [quad, geo, ?, ?]
68
69 ;;; XMM1 = harmonic mean = min of the means
70 0000002A C5D85EE8 vdivps xmm5, xmm4, xmm0 ; 4/x
71 0000002E FFD1 call rcx ; arithmetic mean (under inversion)
72 00000030 C5D85ECD vdivps xmm1, xmm4, xmm5 ; 4/. (the factor of 4 cancels out)
73
74 ;;; XMM5 = arithmetic mean
75 00000034 0F28E8 movaps xmm5, xmm0
76 00000037 FFD1 call rcx
77
78 00000039 0F14E9 unpcklps xmm5, xmm1 ; [arith, harm, ?,?]
79 0000003C C5D014C2 vunpcklps xmm0, xmm5,xmm2 ; x = [arith, harm, quad, geo]
80
81 00000040 0F5CD1 subps xmm2, xmm1 ; largest - smallest mean: guaranteed non-negative
82 00000043 0F2E51F8 ucomiss xmm2, [rcx-8] ; quad-harm > convergence_threshold
83 00000047 73C5 jae .loop
84
85 ; return with the arithmetic mean in the low element of xmm0 = scalar return value
86 00000049 C3 ret
87
88 ;;; "constant pool" between the main function and the helper, like ARM literal pools
89 0000004A ACC52738 .fpconst_threshold: dd 4e-5 ; 4.3e-5 is the highest we can go and still pass the main test cases
90 0000004E 00008040 .fpconst_4: dd 4.0
91 .arith_mean: ; returns XMM5 = hsum(xmm5)/4.
92 00000052 C5D37CED vhaddps xmm5, xmm5 ; slow but small
93 00000056 C5D37CED vhaddps xmm5, xmm5
94 0000005A 0F5EEC divps xmm5, xmm4 ; divide before/after summing doesn't matter mathematically or numerically; divisor is a power of 2
95 0000005D C3 ret
96 0000005E 5E000000 .size: dd $ - meanmean_float_avx
0x5e = 94 bytes
(NASM listing created with nasm -felf64 mean-mean.asm -l/dev/stdout | cut -b -34,$((34+6))-
. Strip the listing part and recover the source with cut -b 34- > mean-mean.asm
)
SIMD horizontal sum and divide by 4 (i.e. arithmetic mean) is implemented in a separate function that we call
(with a function pointer to amortize the cost of the address). With 4/x
before/after, or x^2
before and sqrt after, we get the harmonic mean and quadratic mean. (It was painful to write these div
instructions instead of multiplying by an exactly-representable 0.25
.)
Geometric mean is implemented separately with multiply and chained sqrt. Or with one sqrt first to reduce exponent magnitude and maybe help numerical precision. log is not available, only floor(log2(x))
via AVX512 vgetexpps/pd
. Exp is sort of available via AVX512ER (Xeon Phi only), but with only 2^-23 precision.
Mixing 128-bit AVX instructions and legacy SSE is not a performance problem. Mixing 256-bit AVX with legacy SSE can be on Haswell, but on Skylake it just potentially creates a potential false dependency for SSE instructions. I think my double
version avoids any unnecessary loop-carried dep chains, and bottlenecks
on div/sqrt latency/throughput.
Double version:
108 global meanmean_double_avx
109 meanmean_double_avx:
110 00000080 B9[E8000000] mov ecx, .arith_mean
111 00000085 C4E27D1961F8 vbroadcastsd ymm4, [rcx-8] ; float 4.0
112
113 ;; mean(a,b) = mean(a,b,a,b) for all 4 types of mean
114 ;; so we only ever have to do the length=4 case
115 0000008B C4E37D18C001 vinsertf128 ymm0, ymm0, xmm0, 1 ; [b,a] => [b,a,b,a]
116
117 .loop:
118 ;;; XMM3 = geometric mean: not based on addition.
119 00000091 C5FD51D8 vsqrtpd ymm3, ymm0 ; sqrt first to get magnitude closer to 1.0 for better(?) numerical precision
120 00000095 C4E37D19DD01 vextractf128 xmm5, ymm3, 1 ; extract high lane
121 0000009B C5D159EB vmulpd xmm5, xmm3
122 0000009F 0F12DD movhlps xmm3, xmm5 ; extract high half
123 000000A2 F20F59DD mulsd xmm3, xmm5 ; xmm3 = hproduct(sqrt(xmm0))
124 ; sqrtsd xmm3, xmm3 ; xmm3 = 4th root = geomean(xmm0) ;deferred until quadratic
125
126 ;;; XMM2 = quadratic mean = max of the means
127 000000A6 C5FD59E8 vmulpd ymm5, ymm0,ymm0
128 000000AA FFD1 call rcx ; arith mean of squares
129 000000AC 0F16EB movlhps xmm5, xmm3 ; [quad^2, geo^2]
130 000000AF 660F51D5 sqrtpd xmm2, xmm5 ; [quad , geo]
131
132 ;;; XMM1 = harmonic mean = min of the means
133 000000B3 C5DD5EE8 vdivpd ymm5, ymm4, ymm0 ; 4/x
134 000000B7 FFD1 call rcx ; arithmetic mean under inversion
135 000000B9 C5DB5ECD vdivsd xmm1, xmm4, xmm5 ; 4/. (the factor of 4 cancels out)
136
137 ;;; XMM5 = arithmetic mean
138 000000BD C5FC28E8 vmovaps ymm5, ymm0
139 000000C1 FFD1 call rcx
140
141 000000C3 0F16E9 movlhps xmm5, xmm1 ; [arith, harm]
142 000000C6 C4E35518C201 vinsertf128 ymm0, ymm5, xmm2, 1 ; x = [arith, harm, quad, geo]
143
144 000000CC C5EB5CD1 vsubsd xmm2, xmm1 ; largest - smallest mean: guaranteed non-negative
145 000000D0 660F2E51F0 ucomisd xmm2, [rcx-16] ; quad - harm > threshold
146 000000D5 77BA ja .loop
147
148 ; vzeroupper ; not needed for correctness, only performance
149 ; return with the arithmetic mean in the low element of xmm0 = scalar return value
150 000000D7 C3 ret
151
152 ; "literal pool" between the function
153 000000D8 95D626E80B2E113E .fpconst_threshold: dq 1e-9
154 000000E0 0000000000001040 .fpconst_4: dq 4.0 ; TODO: golf these zeros? vpbroadcastb and convert?
155 .arith_mean: ; returns YMM5 = hsum(ymm5)/4.
156 000000E8 C4E37D19EF01 vextractf128 xmm7, ymm5, 1
157 000000EE C5D158EF vaddpd xmm5, xmm7
158 000000F2 C5D17CED vhaddpd xmm5, xmm5 ; slow but small
159 000000F6 C5D35EEC vdivsd xmm5, xmm4 ; only low element matters
160 000000FA C3 ret
161 000000FB 7B000000 .size: dd $ - meanmean_double_avx
0x7b = 123 bytes
C test harness
#include <immintrin.h>
#include <stdio.h>
#include <math.h>
static const struct ab_avg {
double a,b;
double mean;
} testcases[] = {
{1, 1, 1},
{1, 2, 1.45568889},
{100, 200, 145.568889},
{2.71, 3.14, 2.92103713},
{0.57, 1.78, 1.0848205},
{1.61, 2.41, 1.98965438},
{0.01, 100, 6.7483058},
};
// see asm comments for order of arith, harm, quad, geo
__m128 meanmean_float_avx(__m128); // or float ...
__m256d meanmean_double_avx(__m128d); // or double ...
int main(void) {
int len = sizeof(testcases) / sizeof(testcases[0]);
for(int i=0 ; i<len ; i++) {
const struct ab_avg *p = &testcases[i];
#if 1
__m128 arg = _mm_set_ps(0,0, p->b, p->a);
double res = meanmean_float_avx(arg)[0];
#else
__m128d arg = _mm_loadu_pd(&p->a);
double res = meanmean_double_avx(arg)[0];
#endif
double allowed_diff = (p->b - p->a) / 100000.0;
double delta = fabs(p->mean - res);
if (delta > 1e-3 || delta > allowed_diff) {
printf("%f %f => %.9f but we got %.9f. delta = %g allowed=%gn",
p->a, p->b, p->mean, res, p->mean - res, allowed_diff);
}
}
while(1) {
double a = drand48(), b = drand48(); // range= [0..1)
if (a>b) {
double tmp=a;
a=b;
b=tmp; // sorted
}
// a *= 0.00000001;
// b *= 123156;
// a += 1<<11; b += (1<<12)+1; // float version gets stuck inflooping on 2048.04, 4097.18 at fpthreshold = 4e-5
// a *= 1<<11 ; b *= 1<<11; // scaling to large magnitude makes sum of squares loses more precision
//a += 1<<11; b+= 1<<11; // adding to large magnitude is hard for everything, catastrophic cancellation
#if 1
printf("testing float %g, %gn", a, b);
__m128 arg = _mm_set_ps(0,0, b, a);
__m128 res = meanmean_float_avx(arg);
double quad = res[2], harm = res[1]; // same order as double... for now
#else
printf("testing double %g, %gn", a, b);
__m128d arg = _mm_set_pd(b, a);
__m256d res = meanmean_double_avx(arg);
double quad = res[2], harm = res[1];
#endif
double delta = fabs(quad - harm);
double allowed_diff = (b - a) / 100000.0; // calculated in double even for the float case.
// TODO: use the double res as a reference for float res
// instead of just checking quadratic vs. harmonic mean
if (delta > 1e-3 || delta > allowed_diff) {
printf("%g %g we got q=%g, h=%g, a=%g. delta = %g, allowed=%gn",
a, b, quad, harm, res[0], quad-harm, allowed_diff);
}
}
}
Build with:
nasm -felf64 mean-mean.asm &&
gcc -no-pie -fno-pie -g -O2 -march=native mean-mean.c mean-mean.o
Obviously you need a CPU with AVX support, or an emulator like Intel SDE. To compile on a host without native AVX support, use -march=sandybridge
or -mavx
Run: passes the hard-coded test cases, but for the float version, random test cases often fail the (b-a)/10000
threshold set in the question.
$ ./a.out
(note: empty output before the first "testing float" means clean pass on the constant test cases)
testing float 3.90799e-14, 0.000985395
3.90799e-14 0.000985395 we got q=3.20062e-10, h=3.58723e-05, a=2.50934e-05. delta = -3.5872e-05, allowed=9.85395e-09
testing float 0.041631, 0.176643
testing float 0.0913306, 0.364602
testing float 0.0922976, 0.487217
testing float 0.454433, 0.52675
0.454433 0.52675 we got q=0.48992, h=0.489927, a=0.489925. delta = -6.79493e-06, allowed=7.23169e-07
testing float 0.233178, 0.831292
testing float 0.56806, 0.931731
testing float 0.0508319, 0.556094
testing float 0.0189148, 0.767051
0.0189148 0.767051 we got q=0.210471, h=0.210484, a=0.21048. delta = -1.37389e-05, allowed=7.48136e-06
testing float 0.25236, 0.298197
0.25236 0.298197 we got q=0.274796, h=0.274803, a=0.274801. delta = -6.19888e-06, allowed=4.58374e-07
testing float 0.531557, 0.875981
testing float 0.515431, 0.920261
testing float 0.18842, 0.810429
testing float 0.570614, 0.886314
testing float 0.0767746, 0.815274
testing float 0.118352, 0.984891
0.118352 0.984891 we got q=0.427845, h=0.427872, a=0.427863. delta = -2.66135e-05, allowed=8.66539e-06
testing float 0.784484, 0.893906
0.784484 0.893906 we got q=0.838297, h=0.838304, a=0.838302. delta = -7.09295e-06, allowed=1.09422e-06
FP errors are enough that quad-harm comes out less than zero for some inputs.
Or with a += 1<<11; b += (1<<12)+1;
uncommented:
testing float 2048, 4097
testing float 2048.04, 4097.18
^C (stuck in an infinite loop).
None of these problems happen with double
. Comment out the printf
before each test to see that the output is empty (nothing from the if(delta too high)
block).
TODO: use the double
version as a reference for the float
version, instead of just looking at how they converging with quad-harm.
$endgroup$
x86 machine code (SIMD 4x float using 128-bit SSE1&AVX) 94 bytes
x86 machine code (SIMD 4x double using 256-bit AVX) 123 bytes
float
passes the test cases in the question, but with a loop-exit threshold small enough to make that happen, it's easy for it to get stuck in an infinite loop with random inputs.
SSE1 packed-single-precision instructions are 3 bytes long, but SSE2 and simple AVX instructions are 4 bytes long. (Scalar-single instructions like sqrtss
are also 4 bytes long, which is why I use sqrtps
even though I only care about the low element. It's not even slower than sqrtss on modern hardware). I used AVX for non-destructive destination to save 2 bytes vs. movaps+op.
In the double version we can still do a couple movlhps
to copy 64-bit chunks (because often we only care about the low element of a horizontal sum). Horizontal sum of a 256-bit SIMD vector also requires an extra vextractf128
to get the high half, vs. the slow but small 2x haddps
strategy for float. The double
version also needs 2x 8-byte constants, instead of 2x 4-byte. Overall it comes out at close to 4/3 the size of the float
version.
mean(a,b) = mean(a,a,b,b)
for all 4 of these means, so we can simply duplicate the input up to 4 elements and never have to implement length=2. Thus we can hardcode geometric mean as 4th-root = sqrt(sqrt), for example. And we only need one FP constant, 4.0
.
We have a single SIMD vector of all 4 [a_i, b_i, c_i, d_i]
. From that, we calculate the 4 means as scalars in separate registers, and shuffle them back together for the next iteration. (Horizontal operations on SIMD vectors are inconvenient, but we need to do the same thing for all 4 elements in enough cases that it balances out. I started on an x87 version of this, but it was getting very long and not fun.)
The loop-exit condition of }while(quadratic - harmonic > 4e-5)
(or a smaller constant for double
) is based on @RobinRyder's R answer, and Kevin Cruijssen's Java answer: quadratic mean is always the largest magnitude, and harmonic mean is always the smallest (ignoring rounding errors). So we can check the delta between those two to detect convergence. We return the arithmetic mean as the scalar result. It's usually between those two and is probably the least susceptible to rounding errors.
Float version: callable as float meanmean_float_avx(__m128);
with the arg and return value in xmm0. (So x86-64 System V, or Windows x64 vectorcall, but not x64 fastcall.) Or declare the return-type as __m128
so you can get at the quadratic and harmonic mean for testing.
Letting this take 2 separate float
args in xmm0 and xmm1 would cost 1 extra byte: we'd need a shufps
with an imm8 (instead of just unpcklps xmm0,xmm0
) to shuffle together and duplicate 2 inputs.
40 address align 32
41 code bytes global meanmean_float_avx
42 meanmean_float_avx:
43 00000000 B9[52000000] mov ecx, .arith_mean ; allows 2-byte call reg, and a base for loading constants
44 00000005 C4E2791861FC vbroadcastss xmm4, [rcx-4] ; float 4.0
45
46 ;; mean(a,b) = mean(a,b,a,b) for all 4 types of mean
47 ;; so we only ever have to do the length=4 case
48 0000000B 0F14C0 unpcklps xmm0,xmm0 ; [b,a] => [b,b,a,a]
49
50 ; do{ ... } while(quadratic - harmonic > threshold);
51 .loop:
52 ;;; XMM3 = geometric mean: not based on addition. (Transform to log would be hard. AVX512ER has exp with 23-bit accuracy, but not log. vgetexp = floor(lofg2(x)), so that's no good.)
53 ;; sqrt once *first*, making magnitudes closer to 1.0 to reduce rounding error. Numbers are all positive so this is safe.
54 ;; both sqrts first was better behaved, I think.
55 0000000E 0F51D8 sqrtps xmm3, xmm0 ; xmm3 = 4th root(x)
56 00000011 F30F16EB movshdup xmm5, xmm3 ; bring odd elements down to even
57 00000015 0F59EB mulps xmm5, xmm3
58 00000018 0F12DD movhlps xmm3, xmm5 ; high half -> low
59 0000001B 0F59DD mulps xmm3, xmm5 ; xmm3[0] = hproduct(sqrt(xmm))
60 ; sqrtps xmm3, xmm3 ; sqrt(hprod(sqrt)) = 4th root(hprod)
61 ; common final step done after interleaving with quadratic mean
62
63 ;;; XMM2 = quadratic mean = max of the means
64 0000001E C5F859E8 vmulps xmm5, xmm0,xmm0
65 00000022 FFD1 call rcx ; arith mean of squares
66 00000024 0F14EB unpcklps xmm5, xmm3 ; [quad^2, geo^2, ?, ?]
67 00000027 0F51D5 sqrtps xmm2, xmm5 ; [quad, geo, ?, ?]
68
69 ;;; XMM1 = harmonic mean = min of the means
70 0000002A C5D85EE8 vdivps xmm5, xmm4, xmm0 ; 4/x
71 0000002E FFD1 call rcx ; arithmetic mean (under inversion)
72 00000030 C5D85ECD vdivps xmm1, xmm4, xmm5 ; 4/. (the factor of 4 cancels out)
73
74 ;;; XMM5 = arithmetic mean
75 00000034 0F28E8 movaps xmm5, xmm0
76 00000037 FFD1 call rcx
77
78 00000039 0F14E9 unpcklps xmm5, xmm1 ; [arith, harm, ?,?]
79 0000003C C5D014C2 vunpcklps xmm0, xmm5,xmm2 ; x = [arith, harm, quad, geo]
80
81 00000040 0F5CD1 subps xmm2, xmm1 ; largest - smallest mean: guaranteed non-negative
82 00000043 0F2E51F8 ucomiss xmm2, [rcx-8] ; quad-harm > convergence_threshold
83 00000047 73C5 jae .loop
84
85 ; return with the arithmetic mean in the low element of xmm0 = scalar return value
86 00000049 C3 ret
87
88 ;;; "constant pool" between the main function and the helper, like ARM literal pools
89 0000004A ACC52738 .fpconst_threshold: dd 4e-5 ; 4.3e-5 is the highest we can go and still pass the main test cases
90 0000004E 00008040 .fpconst_4: dd 4.0
91 .arith_mean: ; returns XMM5 = hsum(xmm5)/4.
92 00000052 C5D37CED vhaddps xmm5, xmm5 ; slow but small
93 00000056 C5D37CED vhaddps xmm5, xmm5
94 0000005A 0F5EEC divps xmm5, xmm4 ; divide before/after summing doesn't matter mathematically or numerically; divisor is a power of 2
95 0000005D C3 ret
96 0000005E 5E000000 .size: dd $ - meanmean_float_avx
0x5e = 94 bytes
(NASM listing created with nasm -felf64 mean-mean.asm -l/dev/stdout | cut -b -34,$((34+6))-
. Strip the listing part and recover the source with cut -b 34- > mean-mean.asm
)
SIMD horizontal sum and divide by 4 (i.e. arithmetic mean) is implemented in a separate function that we call
(with a function pointer to amortize the cost of the address). With 4/x
before/after, or x^2
before and sqrt after, we get the harmonic mean and quadratic mean. (It was painful to write these div
instructions instead of multiplying by an exactly-representable 0.25
.)
Geometric mean is implemented separately with multiply and chained sqrt. Or with one sqrt first to reduce exponent magnitude and maybe help numerical precision. log is not available, only floor(log2(x))
via AVX512 vgetexpps/pd
. Exp is sort of available via AVX512ER (Xeon Phi only), but with only 2^-23 precision.
Mixing 128-bit AVX instructions and legacy SSE is not a performance problem. Mixing 256-bit AVX with legacy SSE can be on Haswell, but on Skylake it just potentially creates a potential false dependency for SSE instructions. I think my double
version avoids any unnecessary loop-carried dep chains, and bottlenecks
on div/sqrt latency/throughput.
Double version:
108 global meanmean_double_avx
109 meanmean_double_avx:
110 00000080 B9[E8000000] mov ecx, .arith_mean
111 00000085 C4E27D1961F8 vbroadcastsd ymm4, [rcx-8] ; float 4.0
112
113 ;; mean(a,b) = mean(a,b,a,b) for all 4 types of mean
114 ;; so we only ever have to do the length=4 case
115 0000008B C4E37D18C001 vinsertf128 ymm0, ymm0, xmm0, 1 ; [b,a] => [b,a,b,a]
116
117 .loop:
118 ;;; XMM3 = geometric mean: not based on addition.
119 00000091 C5FD51D8 vsqrtpd ymm3, ymm0 ; sqrt first to get magnitude closer to 1.0 for better(?) numerical precision
120 00000095 C4E37D19DD01 vextractf128 xmm5, ymm3, 1 ; extract high lane
121 0000009B C5D159EB vmulpd xmm5, xmm3
122 0000009F 0F12DD movhlps xmm3, xmm5 ; extract high half
123 000000A2 F20F59DD mulsd xmm3, xmm5 ; xmm3 = hproduct(sqrt(xmm0))
124 ; sqrtsd xmm3, xmm3 ; xmm3 = 4th root = geomean(xmm0) ;deferred until quadratic
125
126 ;;; XMM2 = quadratic mean = max of the means
127 000000A6 C5FD59E8 vmulpd ymm5, ymm0,ymm0
128 000000AA FFD1 call rcx ; arith mean of squares
129 000000AC 0F16EB movlhps xmm5, xmm3 ; [quad^2, geo^2]
130 000000AF 660F51D5 sqrtpd xmm2, xmm5 ; [quad , geo]
131
132 ;;; XMM1 = harmonic mean = min of the means
133 000000B3 C5DD5EE8 vdivpd ymm5, ymm4, ymm0 ; 4/x
134 000000B7 FFD1 call rcx ; arithmetic mean under inversion
135 000000B9 C5DB5ECD vdivsd xmm1, xmm4, xmm5 ; 4/. (the factor of 4 cancels out)
136
137 ;;; XMM5 = arithmetic mean
138 000000BD C5FC28E8 vmovaps ymm5, ymm0
139 000000C1 FFD1 call rcx
140
141 000000C3 0F16E9 movlhps xmm5, xmm1 ; [arith, harm]
142 000000C6 C4E35518C201 vinsertf128 ymm0, ymm5, xmm2, 1 ; x = [arith, harm, quad, geo]
143
144 000000CC C5EB5CD1 vsubsd xmm2, xmm1 ; largest - smallest mean: guaranteed non-negative
145 000000D0 660F2E51F0 ucomisd xmm2, [rcx-16] ; quad - harm > threshold
146 000000D5 77BA ja .loop
147
148 ; vzeroupper ; not needed for correctness, only performance
149 ; return with the arithmetic mean in the low element of xmm0 = scalar return value
150 000000D7 C3 ret
151
152 ; "literal pool" between the function
153 000000D8 95D626E80B2E113E .fpconst_threshold: dq 1e-9
154 000000E0 0000000000001040 .fpconst_4: dq 4.0 ; TODO: golf these zeros? vpbroadcastb and convert?
155 .arith_mean: ; returns YMM5 = hsum(ymm5)/4.
156 000000E8 C4E37D19EF01 vextractf128 xmm7, ymm5, 1
157 000000EE C5D158EF vaddpd xmm5, xmm7
158 000000F2 C5D17CED vhaddpd xmm5, xmm5 ; slow but small
159 000000F6 C5D35EEC vdivsd xmm5, xmm4 ; only low element matters
160 000000FA C3 ret
161 000000FB 7B000000 .size: dd $ - meanmean_double_avx
0x7b = 123 bytes
C test harness
#include <immintrin.h>
#include <stdio.h>
#include <math.h>
static const struct ab_avg {
double a,b;
double mean;
} testcases[] = {
{1, 1, 1},
{1, 2, 1.45568889},
{100, 200, 145.568889},
{2.71, 3.14, 2.92103713},
{0.57, 1.78, 1.0848205},
{1.61, 2.41, 1.98965438},
{0.01, 100, 6.7483058},
};
// see asm comments for order of arith, harm, quad, geo
__m128 meanmean_float_avx(__m128); // or float ...
__m256d meanmean_double_avx(__m128d); // or double ...
int main(void) {
int len = sizeof(testcases) / sizeof(testcases[0]);
for(int i=0 ; i<len ; i++) {
const struct ab_avg *p = &testcases[i];
#if 1
__m128 arg = _mm_set_ps(0,0, p->b, p->a);
double res = meanmean_float_avx(arg)[0];
#else
__m128d arg = _mm_loadu_pd(&p->a);
double res = meanmean_double_avx(arg)[0];
#endif
double allowed_diff = (p->b - p->a) / 100000.0;
double delta = fabs(p->mean - res);
if (delta > 1e-3 || delta > allowed_diff) {
printf("%f %f => %.9f but we got %.9f. delta = %g allowed=%gn",
p->a, p->b, p->mean, res, p->mean - res, allowed_diff);
}
}
while(1) {
double a = drand48(), b = drand48(); // range= [0..1)
if (a>b) {
double tmp=a;
a=b;
b=tmp; // sorted
}
// a *= 0.00000001;
// b *= 123156;
// a += 1<<11; b += (1<<12)+1; // float version gets stuck inflooping on 2048.04, 4097.18 at fpthreshold = 4e-5
// a *= 1<<11 ; b *= 1<<11; // scaling to large magnitude makes sum of squares loses more precision
//a += 1<<11; b+= 1<<11; // adding to large magnitude is hard for everything, catastrophic cancellation
#if 1
printf("testing float %g, %gn", a, b);
__m128 arg = _mm_set_ps(0,0, b, a);
__m128 res = meanmean_float_avx(arg);
double quad = res[2], harm = res[1]; // same order as double... for now
#else
printf("testing double %g, %gn", a, b);
__m128d arg = _mm_set_pd(b, a);
__m256d res = meanmean_double_avx(arg);
double quad = res[2], harm = res[1];
#endif
double delta = fabs(quad - harm);
double allowed_diff = (b - a) / 100000.0; // calculated in double even for the float case.
// TODO: use the double res as a reference for float res
// instead of just checking quadratic vs. harmonic mean
if (delta > 1e-3 || delta > allowed_diff) {
printf("%g %g we got q=%g, h=%g, a=%g. delta = %g, allowed=%gn",
a, b, quad, harm, res[0], quad-harm, allowed_diff);
}
}
}
Build with:
nasm -felf64 mean-mean.asm &&
gcc -no-pie -fno-pie -g -O2 -march=native mean-mean.c mean-mean.o
Obviously you need a CPU with AVX support, or an emulator like Intel SDE. To compile on a host without native AVX support, use -march=sandybridge
or -mavx
Run: passes the hard-coded test cases, but for the float version, random test cases often fail the (b-a)/10000
threshold set in the question.
$ ./a.out
(note: empty output before the first "testing float" means clean pass on the constant test cases)
testing float 3.90799e-14, 0.000985395
3.90799e-14 0.000985395 we got q=3.20062e-10, h=3.58723e-05, a=2.50934e-05. delta = -3.5872e-05, allowed=9.85395e-09
testing float 0.041631, 0.176643
testing float 0.0913306, 0.364602
testing float 0.0922976, 0.487217
testing float 0.454433, 0.52675
0.454433 0.52675 we got q=0.48992, h=0.489927, a=0.489925. delta = -6.79493e-06, allowed=7.23169e-07
testing float 0.233178, 0.831292
testing float 0.56806, 0.931731
testing float 0.0508319, 0.556094
testing float 0.0189148, 0.767051
0.0189148 0.767051 we got q=0.210471, h=0.210484, a=0.21048. delta = -1.37389e-05, allowed=7.48136e-06
testing float 0.25236, 0.298197
0.25236 0.298197 we got q=0.274796, h=0.274803, a=0.274801. delta = -6.19888e-06, allowed=4.58374e-07
testing float 0.531557, 0.875981
testing float 0.515431, 0.920261
testing float 0.18842, 0.810429
testing float 0.570614, 0.886314
testing float 0.0767746, 0.815274
testing float 0.118352, 0.984891
0.118352 0.984891 we got q=0.427845, h=0.427872, a=0.427863. delta = -2.66135e-05, allowed=8.66539e-06
testing float 0.784484, 0.893906
0.784484 0.893906 we got q=0.838297, h=0.838304, a=0.838302. delta = -7.09295e-06, allowed=1.09422e-06
FP errors are enough that quad-harm comes out less than zero for some inputs.
Or with a += 1<<11; b += (1<<12)+1;
uncommented:
testing float 2048, 4097
testing float 2048.04, 4097.18
^C (stuck in an infinite loop).
None of these problems happen with double
. Comment out the printf
before each test to see that the output is empty (nothing from the if(delta too high)
block).
TODO: use the double
version as a reference for the float
version, instead of just looking at how they converging with quad-harm.
edited Apr 3 at 4:49
answered Apr 2 at 18:33
Peter CordesPeter Cordes
2,1521018
2,1521018
add a comment |
add a comment |
$begingroup$
Javascript - 186 bytes
Takes input as an array of numbers. Uses the mean transformations in J42161217's answer to shorten the code.
Try It Online
f=(v,l=[m=(w,k=0)=>w.map(x=>k+=x)&&k/w.length,w=>1/m(w.map(x=>1/x)),w=>Math.E**m(w.map(x=>Math.log(x))),w=>m(w.map(x=>x**2))**.5].map(x=>x(v)).sort((a,b)=>a-b))=>l[3]-l[0]>1e-5?f(l):l[0]
Explanation
f = (
v,
l=[
m=(w,k=0)=>w.map(x=>k+=x)&&k/w.length, // m = w => arithmetic mean of values in w
w=>1/m(w.map(x=>1/x)), // w => harmonic mean of values in w
w=>Math.E**m(w.map(x=>Math.log(x))), // w => geometric mean of values in w
w=>m(w.map(x=>x**2))**.5 // w => quadratic mean of values in w
].map(x=>x(v)) // get array of each mean using input v, stored in l
.sort((a,b)=>a-b) // sort the outputs
) =>
l[3] - l[0] > 1e-5 ? // is the difference between the largest
// and smallest means > 1/100000?
f(l) : // if yes, get the mean mean of the means
l[0] // if no, arbitrarily return the smallest value
// as close enough
$endgroup$
$begingroup$
I thought I was going to be clever and implement the relationship with logarithms but it looks like you and J42161217 got there first!
$endgroup$
– Pureferret
Apr 3 at 11:08
$begingroup$
@Pureferret I take no credit for that,I blatantly stole it :D
$endgroup$
– asgallant
Apr 3 at 15:58
$begingroup$
you wrote it in JavaScript though!
$endgroup$
– Pureferret
Apr 3 at 15:59
1
$begingroup$
That was the easy part. Golfing it was hard.
$endgroup$
– asgallant
Apr 3 at 16:13
1
$begingroup$
The TIL was not configured correctly. I added a TIL link to the answer.
$endgroup$
– asgallant
Apr 4 at 16:27
|
show 2 more comments
$begingroup$
Javascript - 186 bytes
Takes input as an array of numbers. Uses the mean transformations in J42161217's answer to shorten the code.
Try It Online
f=(v,l=[m=(w,k=0)=>w.map(x=>k+=x)&&k/w.length,w=>1/m(w.map(x=>1/x)),w=>Math.E**m(w.map(x=>Math.log(x))),w=>m(w.map(x=>x**2))**.5].map(x=>x(v)).sort((a,b)=>a-b))=>l[3]-l[0]>1e-5?f(l):l[0]
Explanation
f = (
v,
l=[
m=(w,k=0)=>w.map(x=>k+=x)&&k/w.length, // m = w => arithmetic mean of values in w
w=>1/m(w.map(x=>1/x)), // w => harmonic mean of values in w
w=>Math.E**m(w.map(x=>Math.log(x))), // w => geometric mean of values in w
w=>m(w.map(x=>x**2))**.5 // w => quadratic mean of values in w
].map(x=>x(v)) // get array of each mean using input v, stored in l
.sort((a,b)=>a-b) // sort the outputs
) =>
l[3] - l[0] > 1e-5 ? // is the difference between the largest
// and smallest means > 1/100000?
f(l) : // if yes, get the mean mean of the means
l[0] // if no, arbitrarily return the smallest value
// as close enough
$endgroup$
$begingroup$
I thought I was going to be clever and implement the relationship with logarithms but it looks like you and J42161217 got there first!
$endgroup$
– Pureferret
Apr 3 at 11:08
$begingroup$
@Pureferret I take no credit for that,I blatantly stole it :D
$endgroup$
– asgallant
Apr 3 at 15:58
$begingroup$
you wrote it in JavaScript though!
$endgroup$
– Pureferret
Apr 3 at 15:59
1
$begingroup$
That was the easy part. Golfing it was hard.
$endgroup$
– asgallant
Apr 3 at 16:13
1
$begingroup$
The TIL was not configured correctly. I added a TIL link to the answer.
$endgroup$
– asgallant
Apr 4 at 16:27
|
show 2 more comments
$begingroup$
Javascript - 186 bytes
Takes input as an array of numbers. Uses the mean transformations in J42161217's answer to shorten the code.
Try It Online
f=(v,l=[m=(w,k=0)=>w.map(x=>k+=x)&&k/w.length,w=>1/m(w.map(x=>1/x)),w=>Math.E**m(w.map(x=>Math.log(x))),w=>m(w.map(x=>x**2))**.5].map(x=>x(v)).sort((a,b)=>a-b))=>l[3]-l[0]>1e-5?f(l):l[0]
Explanation
f = (
v,
l=[
m=(w,k=0)=>w.map(x=>k+=x)&&k/w.length, // m = w => arithmetic mean of values in w
w=>1/m(w.map(x=>1/x)), // w => harmonic mean of values in w
w=>Math.E**m(w.map(x=>Math.log(x))), // w => geometric mean of values in w
w=>m(w.map(x=>x**2))**.5 // w => quadratic mean of values in w
].map(x=>x(v)) // get array of each mean using input v, stored in l
.sort((a,b)=>a-b) // sort the outputs
) =>
l[3] - l[0] > 1e-5 ? // is the difference between the largest
// and smallest means > 1/100000?
f(l) : // if yes, get the mean mean of the means
l[0] // if no, arbitrarily return the smallest value
// as close enough
$endgroup$
Javascript - 186 bytes
Takes input as an array of numbers. Uses the mean transformations in J42161217's answer to shorten the code.
Try It Online
f=(v,l=[m=(w,k=0)=>w.map(x=>k+=x)&&k/w.length,w=>1/m(w.map(x=>1/x)),w=>Math.E**m(w.map(x=>Math.log(x))),w=>m(w.map(x=>x**2))**.5].map(x=>x(v)).sort((a,b)=>a-b))=>l[3]-l[0]>1e-5?f(l):l[0]
Explanation
f = (
v,
l=[
m=(w,k=0)=>w.map(x=>k+=x)&&k/w.length, // m = w => arithmetic mean of values in w
w=>1/m(w.map(x=>1/x)), // w => harmonic mean of values in w
w=>Math.E**m(w.map(x=>Math.log(x))), // w => geometric mean of values in w
w=>m(w.map(x=>x**2))**.5 // w => quadratic mean of values in w
].map(x=>x(v)) // get array of each mean using input v, stored in l
.sort((a,b)=>a-b) // sort the outputs
) =>
l[3] - l[0] > 1e-5 ? // is the difference between the largest
// and smallest means > 1/100000?
f(l) : // if yes, get the mean mean of the means
l[0] // if no, arbitrarily return the smallest value
// as close enough
edited Apr 4 at 16:48
answered Apr 1 at 23:19
asgallantasgallant
26915
26915
$begingroup$
I thought I was going to be clever and implement the relationship with logarithms but it looks like you and J42161217 got there first!
$endgroup$
– Pureferret
Apr 3 at 11:08
$begingroup$
@Pureferret I take no credit for that,I blatantly stole it :D
$endgroup$
– asgallant
Apr 3 at 15:58
$begingroup$
you wrote it in JavaScript though!
$endgroup$
– Pureferret
Apr 3 at 15:59
1
$begingroup$
That was the easy part. Golfing it was hard.
$endgroup$
– asgallant
Apr 3 at 16:13
1
$begingroup$
The TIL was not configured correctly. I added a TIL link to the answer.
$endgroup$
– asgallant
Apr 4 at 16:27
|
show 2 more comments
$begingroup$
I thought I was going to be clever and implement the relationship with logarithms but it looks like you and J42161217 got there first!
$endgroup$
– Pureferret
Apr 3 at 11:08
$begingroup$
@Pureferret I take no credit for that,I blatantly stole it :D
$endgroup$
– asgallant
Apr 3 at 15:58
$begingroup$
you wrote it in JavaScript though!
$endgroup$
– Pureferret
Apr 3 at 15:59
1
$begingroup$
That was the easy part. Golfing it was hard.
$endgroup$
– asgallant
Apr 3 at 16:13
1
$begingroup$
The TIL was not configured correctly. I added a TIL link to the answer.
$endgroup$
– asgallant
Apr 4 at 16:27
$begingroup$
I thought I was going to be clever and implement the relationship with logarithms but it looks like you and J42161217 got there first!
$endgroup$
– Pureferret
Apr 3 at 11:08
$begingroup$
I thought I was going to be clever and implement the relationship with logarithms but it looks like you and J42161217 got there first!
$endgroup$
– Pureferret
Apr 3 at 11:08
$begingroup$
@Pureferret I take no credit for that,I blatantly stole it :D
$endgroup$
– asgallant
Apr 3 at 15:58
$begingroup$
@Pureferret I take no credit for that,I blatantly stole it :D
$endgroup$
– asgallant
Apr 3 at 15:58
$begingroup$
you wrote it in JavaScript though!
$endgroup$
– Pureferret
Apr 3 at 15:59
$begingroup$
you wrote it in JavaScript though!
$endgroup$
– Pureferret
Apr 3 at 15:59
1
1
$begingroup$
That was the easy part. Golfing it was hard.
$endgroup$
– asgallant
Apr 3 at 16:13
$begingroup$
That was the easy part. Golfing it was hard.
$endgroup$
– asgallant
Apr 3 at 16:13
1
1
$begingroup$
The TIL was not configured correctly. I added a TIL link to the answer.
$endgroup$
– asgallant
Apr 4 at 16:27
$begingroup$
The TIL was not configured correctly. I added a TIL link to the answer.
$endgroup$
– asgallant
Apr 4 at 16:27
|
show 2 more comments
$begingroup$
Perl 5, 92 72 bytes
sub f{@_=map{$m=$_;sum(map$_**$m/@_,@_)**(1/$m)}1,1e-9,-1,2for 1..9;pop}
Try it online!
...using some dirty tricks.
$endgroup$
add a comment |
$begingroup$
Perl 5, 92 72 bytes
sub f{@_=map{$m=$_;sum(map$_**$m/@_,@_)**(1/$m)}1,1e-9,-1,2for 1..9;pop}
Try it online!
...using some dirty tricks.
$endgroup$
add a comment |
$begingroup$
Perl 5, 92 72 bytes
sub f{@_=map{$m=$_;sum(map$_**$m/@_,@_)**(1/$m)}1,1e-9,-1,2for 1..9;pop}
Try it online!
...using some dirty tricks.
$endgroup$
Perl 5, 92 72 bytes
sub f{@_=map{$m=$_;sum(map$_**$m/@_,@_)**(1/$m)}1,1e-9,-1,2for 1..9;pop}
Try it online!
...using some dirty tricks.
edited Apr 2 at 15:12
answered Apr 1 at 20:39
Kjetil S.Kjetil S.
58925
58925
add a comment |
add a comment |
$begingroup$
SNOBOL4 (CSNOBOL4), 296 bytes
X =INPUT
Y =INPUT
A =(X + Y) / 2.
P =X * Y
G =P ^ 0.5
H =P / A
Q =(2 * (A ^ 2) - P) ^ 0.5
O OUTPUT =EQ(Q,A) Q :S(END)
M =A
N =G
O =H
P =Q
A =(M + N + O + P) / 4
G =(M * N * O * P) ^ 0.25
H =4 / (1 / M + 1 / N + 1 / O + 1 / P)
Q =((M ^ 2 + N ^ 2 + O ^ 2 + P ^ 2) / 4) ^ 0.5 :(O)
END
Try it online!
Straightforward implementation. Uses a trick from my answer to a related question to golf a bit more.
$endgroup$
add a comment |
$begingroup$
SNOBOL4 (CSNOBOL4), 296 bytes
X =INPUT
Y =INPUT
A =(X + Y) / 2.
P =X * Y
G =P ^ 0.5
H =P / A
Q =(2 * (A ^ 2) - P) ^ 0.5
O OUTPUT =EQ(Q,A) Q :S(END)
M =A
N =G
O =H
P =Q
A =(M + N + O + P) / 4
G =(M * N * O * P) ^ 0.25
H =4 / (1 / M + 1 / N + 1 / O + 1 / P)
Q =((M ^ 2 + N ^ 2 + O ^ 2 + P ^ 2) / 4) ^ 0.5 :(O)
END
Try it online!
Straightforward implementation. Uses a trick from my answer to a related question to golf a bit more.
$endgroup$
add a comment |
$begingroup$
SNOBOL4 (CSNOBOL4), 296 bytes
X =INPUT
Y =INPUT
A =(X + Y) / 2.
P =X * Y
G =P ^ 0.5
H =P / A
Q =(2 * (A ^ 2) - P) ^ 0.5
O OUTPUT =EQ(Q,A) Q :S(END)
M =A
N =G
O =H
P =Q
A =(M + N + O + P) / 4
G =(M * N * O * P) ^ 0.25
H =4 / (1 / M + 1 / N + 1 / O + 1 / P)
Q =((M ^ 2 + N ^ 2 + O ^ 2 + P ^ 2) / 4) ^ 0.5 :(O)
END
Try it online!
Straightforward implementation. Uses a trick from my answer to a related question to golf a bit more.
$endgroup$
SNOBOL4 (CSNOBOL4), 296 bytes
X =INPUT
Y =INPUT
A =(X + Y) / 2.
P =X * Y
G =P ^ 0.5
H =P / A
Q =(2 * (A ^ 2) - P) ^ 0.5
O OUTPUT =EQ(Q,A) Q :S(END)
M =A
N =G
O =H
P =Q
A =(M + N + O + P) / 4
G =(M * N * O * P) ^ 0.25
H =4 / (1 / M + 1 / N + 1 / O + 1 / P)
Q =((M ^ 2 + N ^ 2 + O ^ 2 + P ^ 2) / 4) ^ 0.5 :(O)
END
Try it online!
Straightforward implementation. Uses a trick from my answer to a related question to golf a bit more.
answered Apr 2 at 19:31
GiuseppeGiuseppe
17.6k31153
17.6k31153
add a comment |
add a comment |
If this is an answer to a challenge…
…Be sure to follow the challenge specification. However, please refrain from exploiting obvious loopholes. Answers abusing any of the standard loopholes are considered invalid. If you think a specification is unclear or underspecified, comment on the question instead.
…Try to optimize your score. For instance, answers to code-golf challenges should attempt to be as short as possible. You can always include a readable version of the code in addition to the competitive one.
Explanations of your answer make it more interesting to read and are very much encouraged.…Include a short header which indicates the language(s) of your code and its score, as defined by the challenge.
More generally…
…Please make sure to answer the question and provide sufficient detail.
…Avoid asking for help, clarification or responding to other answers (use comments instead).
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fcodegolf.stackexchange.com%2fquestions%2f182460%2fcalculate-the-mean-mean-of-two-numbers%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
$begingroup$
Sandbox link
$endgroup$
– someone
Mar 31 at 7:41
$begingroup$
Somewhat related: Calculate the Arithmetic–Geometric Mean
$endgroup$
– user202729
Mar 31 at 13:38
8
$begingroup$
How precise are we supposed to be?
$endgroup$
– Embodiment of Ignorance
Mar 31 at 17:52
2
$begingroup$
Closely related
$endgroup$
– Giuseppe
Apr 1 at 13:18
1
$begingroup$
Can we assume the first input is always smaller than the second, as in all your test cases? (If not I'll rollback my Java answer.)
$endgroup$
– Kevin Cruijssen
Apr 2 at 7:59