^{Отказ от ответственности: среднее значение составлено мной}

Определите среднее арифметическое чисел как Определите среднее геометрическое чисел как Определить среднее гармоническое для n чисел как M _ {- 1} (x_1, ..., x_n) = \ frac {n} {\ frac {1 } {x_2} + \ frac {1} {x_2} + ... + \ frac {1} {x_n}} Определить среднеквадратичное значение n чисел как M_2 (x_1, ..., x_n) = \ root \ of {\ frac {x_1 ^ 2 + x_2 ^ 2 + ... + x_n ^ 2} {n}} Среднее значение ( M_M ) определяется следующим образом: Определите четыре последовательности ( a_k, b_k, c_k, d_k ) как$n$

$$M_1(x_1,...,x_n)=\frac{x_1+x_2+...+x_n}{n}$$ $n$ $$M_0(x_1,...,x_n)=\root{n}\of{x_1x_2...x_n}$$ $n$ $$M_{-1}(x_1,...,x_n)=\frac{n}{\frac{1}{x_2} + \frac{1}{x_2} + ... + \frac{1}{x_n}}$$ $n$ $$M_2(x_1,...,x_n)=\root\of{\frac{x_1^2+x_2^2+...+x_n^2}{n}}$$ $M_M$к,Ьк,ск,dKс0=М1(х1,...,хп),$a_k, b_k, c_k, d_k$ $$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)$$ Все четыре последовательности сходятся к тот же номер,$M_M(x_1,x_2,...,x_n)$ .

пример

Среднее значение 1 и 2 рассчитывается следующим образом: начните с

$$a_0 = (1+2)/2 = 1.5, b_0 = \root\of{1 * 2} = \root\of2 \approx 1.4142,\\ c_0 = \frac2{\frac1{1}+\frac1{2}} = \frac4{3} \approx 1.3333, d_0 = \root\of{\frac{1^2+2^2}2} = \root\of{\frac5{2}} \approx 1.5811.$$ Тогда $$a_1 = \frac{1.5+1.4142+1.3333+1.5811}4 \approx 1.4571,\\ b_1 = \root^4\of{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 = \root\of{\frac{1.5^2+1.4142^2+1.3333^2+1.5811^2}4} \approx 1.4601.$$ Дальнейший расчет последовательностей должен быть понятен. Видно, что они сходятся к одному и тому же числу, примерно $1.45568889$ .

Вызов

Для двух положительных действительных чисел $a$ и $b$ ( $a<b$ ) вычислите их среднее значение $M_M(a,b)$ .

Контрольные примеры

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

Заметки

• Ваша программа действительна, если разница между ее выходом и правильным выходом не превышает 1/100000 от абсолютного значения разности между входными числами.
• На выходе должно быть одно число.

Это код-гольф , поэтому выигрывает самый короткий код!

Wolfram Language (Mathematica) , 52 байта

#//.x_:>N@{M@x,E^M@Log@x,1/M[1/x],M[x^2]^.5}&
M=Mean

Попробуйте онлайн!

В моем первом подходе я использовал эти встроенные функции
Mean GeometricMean HarmonicMeanиRootMeanSquare

Вот некоторые замены для сохранения байтов

HarmonicMean-> 1/Mean[1/x] от @Robin Ryder (сохранено 3 байта)
GeometricMean-> E^Mean@Log@xот @A. Рекс (2 байта сохранены)
RootMeanSquare-> Mean[x^2]^.5@A. Рекс (4 байта сохранены)

наконец , мы можем приписать Meanк M(как это было предложено @ovs) и сохранить более 5 байт

R 70 69 67 байт

x=scan();`?`=mean;while(x-?x)x=c((?x^2)^.5,?x,2^?log2(x),1/?1/x);?x

Попробуйте онлайн!

_{-1 байт с лучшей подготовкой.
-2 байта при переключении на базу 2.}

Как и некоторые другие ответы, здесь используется выражение среднего геометрического в качестве среднего арифметического на логарифмической шкале (здесь, в базе 2):

$$M_0(x_1,\ldots, x_n) = 2^{M_1(\log_2 x_1,\ldots,\log_2 x_n)}.$$

Также используется тот факт, что $\forall k, d_k\geq a_k \geq b_k\geq c_k$ , то есть $d_k=\max(a_k, b_k, c_k, d_k)$ . Следовательно, условие $a_k=b_k=c_k=d_k$ эквивалентно $d_k = M_1(a_k, b_k, c_k, d_k)$ , что я и использую в цикле while; это достигается путем злоупотребления синтаксисом,whileкоторый учитывает первый элемент, только когда условие является вектором, отсюда и порядок, в котором хранятся средства. (Обратите внимание, что мы могли бы также использовать$c_k$ вместо этого, так как это минимум из четырех, но мы не могли использовать$a_k$ или$b_k$ в условии.)

Когда мы выходим из цикла while, xэто постоянный вектор. Финал ?xвычисляет его значение, чтобы уменьшить его до скаляра.

J , 34 байта

(31 как выражение без присваивания переменной f)

f=:1{(^.z,%z,*:z,[z=:(+/%#)&.:)^:_

Попробуйте онлайн!

Для функций aи b, a &.: b( «а при Ь» ( связанная с этим проблема )) равносильна (b inv) a b- применить Ь, то а, то обратный Ь. В этом случае среднее геометрическое / гармоническое / квадратичное - это среднее арифметическое «под» логарифмом, инверсией и квадратом соответственно.

TI-BASIC, 42 35 34 байта

-1 байт благодаря @SolomonUcko

While max(ΔList(Ans:{mean(Ans),√(mean(Ans²)),mean(Ans^-1)^-1,e^(mean(ln(Ans:End:Ans(1

Ввод представляет собой список из двух целых чисел в Ans.
Вывод сохраняется Ansи автоматически распечатывается по завершении программы.

Формулы, используемые для геометрических, гармонических и квадратичных средних, основаны на объяснениях пользователя 202729 .

Пример:

{1,2
           {1 2}
prgmCDGFB
     1.455688891
{100,200
       {100 200}
prgmCDGFB
     145.5688891

Пояснение:
(Новые строки были добавлены для пояснения. Они НЕ появляются в коде.)

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

---

Заметки:

TI-BASIC - это токенизированный язык. Количество символов не равно количеству байтов.

e^(это этот один байт маркера.

^-1используется для этого однобайтового токена. Вместо этого
я выбрал запись, ^-1потому что токен выглядит как ֿ¹в блоке кода.

√(это этот один байт маркера.

ΔList(это это два байта маркера.

Ява 10, 234 229 214 211 215 206 203 196 180 177 байт

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 байт благодаря @PeterCordes .
-15 еще байт благодаря @PeterCordes , вдохновленных @RobinRyder R ответа «s .
+4 байта, потому что я предполагал, что входы предварительно упорядочены.
-27 байт благодаря @ OlivierGrégoire .

Попробуйте онлайн.

Объяснение:

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

Древесный уголь , 40 байт

Ｗ‹⌊θ⌈θ≔⟦∕ΣθＬθＸΠθ∕¹Ｌθ∕ＬθΣ∕¹θ₂∕ΣＸθ²Ｌθ⟧θＩ⊟θ

Попробуйте онлайн! Ссылка на подробную версию кода. Принимает ввод как массив чисел. Объяснение:

Ｗ‹⌊θ⌈θ

Повторите, пока массив содержит разные значения ...

≔⟦....⟧θ

... заменить массив списком значений:

∕ΣθＬθ

... Значение...

ＸΠθ∕¹Ｌθ

... среднее геометрическое ...

∕ＬθΣ∕¹θ

... среднее гармоническое ...

₂∕ΣＸθ²Ｌθ

... и корень означает квадрат.

Ｉ⊟θ

Приведите элемент массива к строке и неявно напечатайте его.

Желе , 24 байта

;µ*Ɱ-;ؽ¤©Æm*İɼ;P½½ƊµÐLḢ

Попробуйте онлайн!

PowerShell , 182 180 183 байта

$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]}

Попробуйте онлайн!

05AB1E , 26 24 23 байта

Δ©ÅA®.²ÅAo®zÅAz®nÅAt)}н

Попробуйте онлайн или посмотрите шаги всех тестовых случаев .

-1 байт благодаря @Grimy .

23-байтовая альтернатива для среднего геометрического:

Δ©P®gzm®ÅA®zÅAz®nÅAt)}н

Попробуйте онлайн или посмотрите шаги всех тестовых случаев .

Объяснение:

Δ         # Loop until the list no longer changes:
 ©        #  Store the current list in variable `®` (without popping)
          #  (which is the implicit input-list in the first iteration)
          #  Arithmetic mean:
  ÅA      #   Builtin to calculate the arithmetic mean of the list
          #  Geometric mean:
  ®.²     #   Take the base-2 logarithm of each value in the list `®`
     ÅA   #   Get the arithmetic mean of that list
       o  #   And take 2 to the power of this mean
          #  Harmonic mean:
  ®z      #   Get 1/x for each value x in the list `®`
    ÅA    #   Get the arithmetic mean of that list
      z   #   And calculate 1/y for this mean y
          #  Quadratic mean:
  ®n      #   Take the square of each number x in the list from the register
    ÅA    #   Calculate the arithmetic mean of this list
      t   #   And take the square-root of that mean
  )       #  Wrap all four results into a list
}н        # After the list no longer changes: pop and push its first value
          # (which is output implicitly as result)

Желе , 25 24 байта

Wẋ4¹ÆlÆeƭ²½ƭİ4ƭÆm$€⁺µÐLḢ

Попробуйте онлайн!

объяснение

                    µÐ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

Python 3 , 152 байта

from math import*
s=sum
def f(*a):l=len(a);return 2>len({*a})and{*a}or f(s(a)/l,l/s(map(pow,a,l*[-1])),exp(s(map(log,a))/l),(s(map(pow,a,l*[2]))/l)**.5)

Попробуйте онлайн!

Рекурсивная функция f, будет сходиться к точности с плавающей точкой. В принципе работает для всех списков положительных чисел любого размера, но ограничен рекурсивным пределом Python для некоторых тестовых случаев ошибка округления.

---

В качестве альтернативы можно рассчитать с точностью до 9 десятичных знаков:

Python 3 , 169 байт

from math import*
s=sum
def f(*a):l=len(a);return(2>len({round(i,9)for i in a}))*a[0]or f(s(a)/l,l/s(map(pow,a,l*[-1])),exp(s(map(log,a))/l),(s(map(pow,a,l*[2]))/l)**.5)

Попробуйте онлайн!

C # , 173 байта

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)));

Попробуйте онлайн!

Чисто , 124 байта

import StdEnv
f=avg o limit o iterate\l=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]

Попробуйте онлайн!

Выполняет операцию, пока результат не перестанет меняться.

Ура для ограниченной точности с плавающей точкой!

Pyth, 32 байта

h.Wt{H[.OZ@*FZJlZcJscL1Z@.O^R2Z2

Попробуйте это онлайн здесь , или проверьте все тесты (план два, см. Примечание ниже) сразу здесь . Принимает ввод в виде списка.

Кажется, есть некоторые проблемы с округлением, поскольку некоторые входные данные не сходятся правильно, когда они должны были бы. В частности, тестовый пример 0.01 100застревает в значениях [6.748305820749738, 6.748305820749738, 6.748305820749739, 6.748305820749738], а тестовый пример 1.61 2.41застревает в [1.9896543776640825, 1.9896543776640825, 1.9896543776640827, 1.9896543776640825]- обратите внимание, что в обоих случаях 3-е среднее (среднее гармоническое) отличается от других.

Я не уверен, что эта проблема делает мою запись недействительной, но я все равно публикую ее, так как она должна работать. Если это неприемлемо, это можно исправить, введя .RRTперед, округляя [каждое из средств до 10 десятичных знаков, как показано в этом наборе тестов .

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

Japt v2.0a0 -g, 42 38 байт

â ÊÉ?ß[Ux²÷(V=UÊ)¬Ux÷V U×qV V÷Ux!÷1]:U

Должен быть более короткий путь ... Это чудовище! Сохранено 4 байта благодаря @Shaggy!

Попытайся

C # (интерактивный компилятор Visual C #) , 177 байт

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;

Спасибо @KevinCruijjsen за указание на то, что использование точности с плавающей точкой вызывает проблемы! Было бы 163 байта, если бы двойники были совершенно точными.

Попробуйте онлайн!

машинный код x86 (SIMD 4x float с использованием 128-битного SSE1 и AVX) 94 байта

машинный код x86 (SIMD 4x double с использованием 256-битного AVX) 123 байта

float проходит тестовые случаи в вопросе, но с порогом выхода из цикла, достаточно маленьким, чтобы это произошло, легко застрять в бесконечном цикле со случайными входами.

Инструкции SSE1 упакованной одинарной точности имеют длину 3 байта, а инструкции SSE2 и простые AVX имеют длину 4 байта. (Отдельные скалярные инструкции, такие как sqrtss, также имеют длину 4 байта, поэтому я использую их, sqrtpsхотя мне важен только низкий элемент. Он даже медленнее, чем sqrtss на современном оборудовании). Я использовал AVX для неразрушающего назначения, чтобы сохранить 2 байта против movaps + op.
В двойной версии мы все еще можем сделать пару movlhpsдля копирования 64-битных кусков (потому что часто мы заботимся только о младшем элементе горизонтальной суммы). Горизонтальная сумма 256-битного вектора SIMD также требует дополнительного, vextractf128чтобы получить верхнюю половину, по сравнению с медленной, но небольшой 2-кратной haddpsстратегией для float . doubleВерсия также нуждается в 2x 8-байтовых константах вместо 2x 4-байтовых. В целом он выходит примерно на 4/3 размера floatверсии.

mean(a,b) = mean(a,a,b,b)для всех 4-х из этих средств , поэтому мы можем просто дублировать ввод до 4-х элементов и никогда не реализовывать length = 2. Таким образом, мы можем жестко закодировать геометрическое среднее как, например, 4-й корень = sqrt (sqrt). И нам нужна только одна константа FP 4.0.

У нас есть один вектор SIMD из всех 4 [a_i, b_i, c_i, d_i]. Исходя из этого, мы вычисляем 4 средних как скаляры в отдельных регистрах и перемешиваем их вместе для следующей итерации. (Горизонтальные операции над векторами SIMD неудобны, но мы должны сделать то же самое для всех 4 элементов в достаточном количестве случаев, чтобы он уравновешивался. Я начал на этой версии для x87, но она становилась очень длинной и не увлекательной.)

Условие петлевой выход }while(quadratic - harmonic > 4e-5)(или меньше константа для double) основана на @ R ответа RobinRyder в и Java ответ Кевин Cruijssen в : среднее квадратичное всегда по величине величины и среднее гармоническое всегда наименьшее ( без учета ошибок округления). Таким образом, мы можем проверить дельту между этими двумя, чтобы обнаружить сходимость. Мы возвращаем среднее арифметическое как скалярный результат. Это обычно между этими двумя и, вероятно, наименее подвержены ошибкам округления.

Версия с плавающей запятой: может вызываться как float meanmean_float_avx(__m128);с аргументом arg и возвращаемым значением в xmm0. (Таким образом, x86-64 System V, или Windows x64 vectorcall, но не x64 fastcall.) Или объявите тип возврата как __m128таковой, чтобы вы могли получить квадратичное и гармоническое среднее для тестирования.

Если для этого потребуются 2 отдельных floatаргумента в xmm0 и xmm1, это будет стоить 1 дополнительный байт: нам понадобится a shufpsс imm8 (вместо просто unpcklps xmm0,xmm0), чтобы перемешать и дублировать 2 входа.

    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 создан с помощью nasm -felf64 mean-mean.asm -l/dev/stdout | cut -b -34,$((34+6))-. Удалите часть списка и восстановите источник с помощью cut -b 34- > mean-mean.asm)

SIMD горизонтальная сумма и деление на 4 (т.е. среднее арифметическое) реализована в отдельной функции, которую мы call(с указателем функции, чтобы амортизировать стоимость адреса). С 4/xдо / после или x^2до и после, мы получаем среднее гармоническое и квадратичное среднее. (Было больно писать эти divинструкции вместо умножения на точно представимое 0.25.)

Среднее геометрическое реализовано отдельно с умноженным и цепным квадратами. Или с одним sqrt первым, чтобы уменьшить величину экспоненты и, возможно, помочь численной точности. floor(log2(x))Журнал недоступен, только через AVX512 vgetexpps/pd. Exp является своего рода доступным через AVX512ER (только Xeon Phi), но только с точностью 2 ^ -23.

Смешивание 128-битных инструкций AVX и устаревших SSE не является проблемой производительности. Смешивание 256-битного AVX с устаревшим SSE может быть на Haswell, но на Skylake это просто потенциально создает потенциальную ложную зависимость для инструкций SSE. Я думаю, что моя doubleверсия избегает любых ненужных переносимых по цепочке цепочек депов и узких мест в div / sqrt латентности / пропускной способности.

Двойная версия:

   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

---

Испытательный жгут

#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=%g\n",
                   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, %g\n", 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, %g\n", 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=%g\n",
                   a, b, quad, harm, res[0], quad-harm, allowed_diff);
        }
    }

}

Сборка с:

nasm -felf64 mean-mean.asm &&
gcc -no-pie -fno-pie -g -O2 -march=native mean-mean.c mean-mean.o

Очевидно, вам нужен процессор с поддержкой AVX или эмулятор, такой как Intel SDE. Чтобы скомпилировать на хосте без встроенной поддержки AVX, используйте -march=sandybridgeили-mavx

Выполнить: проходит жестко заданные тестовые примеры, но для плавающей версии случайные тестовые примеры часто не соответствуют (b-a)/10000пороговому значению, заданному в вопросе.

$ ./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 достаточны для того, чтобы квад-вред оказался меньше нуля для некоторых входных данных.

Или с a += 1<<11; b += (1<<12)+1;комментариями

testing float 2048, 4097
testing float 2048.04, 4097.18
^C  (stuck in an infinite loop).

Ни одна из этих проблем не происходит с double. Закомментируйте printfперед каждым тестом, чтобы увидеть, что вывод пуст (ничего из if(delta too high)блока).

TODO: используйте doubleверсию как ссылку на floatверсию, вместо того, чтобы просто посмотреть, как они сходятся с квад-вредом.

Javascript - 186 байт

Принимает ввод как массив чисел. Использует средние преобразования в ответе J42161217, чтобы сократить код.

Попробуйте онлайн

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]

объяснение

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

Perl 5 , 92 72 байта

sub f{@_=map{$m=$_;sum(map$_**$m/@_,@_)**(1/$m)}1,1e-9,-1,2for 1..9;pop}

Попробуйте онлайн!

... используя некоторые грязные трюки.

SNOBOL4 (CSNOBOL4) , 296 байт

	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

Попробуйте онлайн!

Простая реализация. Использую трюк из моего ответа на связанный вопрос, чтобы немного больше в гольфе.