Biarkan x
, y
menjadi dua angka floating-point. Apa cara yang benar untuk menghitung artinya?
Cara naif (x+y)/2
dapat menghasilkan luapan saat x
dan y
terlalu besar. Saya pikir 0.5 * x + 0.5 * y
mungkin lebih baik, tetapi melibatkan dua perkalian (yang mungkin tidak efisien), dan saya tidak yakin apakah itu cukup baik. Apakah ada cara yang lebih baik?
Gagasan lain yang saya mainkan adalah (y/2)(1 + x/y)
jika x<=y
. Tetapi sekali lagi, saya tidak yakin bagaimana menganalisis ini dan membuktikan bahwa itu memenuhi persyaratan saya.
Selain itu, saya perlu jaminan bahwa rata-rata yang dihitung akan >= min(x,y)
dan <= max(x,y)
. Seperti yang ditunjukkan dalam jawaban Don Hatch , mungkin cara yang lebih baik untuk mengajukan pertanyaan ini adalah: Apa implementasi dari rata-rata dua angka yang selalu memberikan hasil akurat yang paling mungkin? Yaitu, jika x
dan y
apakah angka floating-point, bagaimana cara menghitung angka floating-point terdekat (x+y)/2
? Dalam hal ini, rata-rata yang dihitung secara otomatis >= min(x,y)
dan <= max(x,y)
. Lihat jawaban Don Hatch untuk detailnya.
Catatan: Prioritas saya adalah akurasi yang kuat. Efisiensi dapat dihabiskan. Namun, jika ada banyak algoritma yang kuat dan akurat, saya akan memilih yang paling efisien.
sumber
Jawaban:
Saya pikir Akurasi dan Stabilitas Algoritma Numerik Higham membahas bagaimana orang dapat menganalisis jenis masalah ini. Lihat Bab 2, terutama latihan 2.8.
Dalam jawaban ini saya ingin menunjukkan sesuatu yang tidak benar-benar dibahas dalam buku Higham (sepertinya tidak terlalu dikenal luas, dalam hal ini). Jika Anda tertarik untuk membuktikan properti dari algoritma numerik sederhana seperti ini, Anda dapat menggunakan kekuatan pemecah SMT modern ( Satisfiability Modulo Theories ), seperti z3 , menggunakan paket seperti sbv di Haskell. Ini agak lebih mudah daripada menggunakan pensil dan kertas.
Misalkan saya diberi nilai , dan saya ingin tahu apakah z = ( x + y ) / 2 memenuhi x ≤ z ≤ y . Berikut kode Haskell0 ≤ x ≤ y z= ( x +y) / 2 x ≤ z≤ y
akan membiarkan saya melakukan ini secara otomatis . Berikutx ≤ f u n ( x , y) ≤ y x , y 0 ≤ x ≤ y
test1 fun
adalah proposisi bahwa untuk semua pelampung terbatas x , y dengan 0 ≤ x ≤ y .Itu meluap. Misalkan saya sekarang mengambil formula Anda yang lain:z= x / 2 + y/ 2
Tidak berfungsi (karena underflow bertahap: , yang mungkin tidak intuitif karena semua aritmatika menjadi basis-2).( x / 2 ) × 2 ≠ x
Sekarang coba :z= x + ( y- x ) / 2
Bekerja! Ini
Q.E.D.
adalah bukti bahwatest1
properti berlaku untuk semua pelampung seperti yang didefinisikan di atas.Bagaimana dengan yang sama, tetapi terbatas pada (bukan 0 ≤ x ≤ y )?x ≤ y 0 ≤ x ≤ y
Oke, jadi jika meluap, bagaimana dengan z = x + ( y / 2 - x / 2 ) ?y- x z= x + ( y/ 2-x / 2)
Jadi sepertinya di antara rumus yang saya coba di sini, tampaknya berfungsi (dengan bukti juga). Pendekatan SMT solver bagi saya merupakan cara yang jauh lebih cepat untuk menjawab kecurigaan tentang rumus titik-mengambang sederhana daripada melalui analisis kesalahan titik-mengambang dengan pensil dan kertas.x + ( y/ 2-x / 2)
Akhirnya, tujuan akurasi dan stabilitas seringkali bertentangan dengan tujuan kinerja. Untuk kinerja, saya tidak benar-benar melihat bagaimana Anda dapat melakukan lebih baik daripada , terutama karena kompiler masih akan melakukan tugas berat menerjemahkan ini menjadi instruksi mesin untuk Anda.( x + y) / 2
SFloat
SDouble
-ffast-math
PPPS Saya terbawa sedikit melihat hanya pada ekspresi aljabar sederhana tanpa persyaratan. Don Hatch 's rumus secara ketat baik.
sumber
>>> x = -1.; y = 1.+2.**-52; print `2**-53`, `(x+y)/2.`, `x+(y/2.-x/2.)`
Pertama, perhatikan bahwa jika Anda memiliki metode yang memberikan jawaban paling akurat dalam semua kasus, maka itu akan memenuhi kondisi yang Anda butuhkan. (Perhatikan bahwa saya mengatakan sebuah jawaban yang paling akurat daripada yang jawabannya paling akurat, karena mungkin ada dua pemenang.) Bukti: Jika, sebaliknya, Anda punya jawaban yang akurat-sebagai-mungkin bahwa tidak tidak memenuhi kondisi yang diperlukan, yang berarti
answer<min(x,y)<=max(x,y)
(dalam halmin(x,y)
ini jawaban yang lebih baik, kontradiksi), ataumin(x,y)<=max(x,y)<answer
(dalam hal inimax(x,y)
ini jawaban yang lebih baik, kontradiksi).Jadi saya pikir itu berarti pertanyaan Anda bermuara pada menemukan jawaban yang paling akurat. Dengan asumsi aritmetika IEEE754 secara keseluruhan, saya mengusulkan yang berikut:
Argumen saya bahwa ini memberikan jawaban yang paling akurat adalah analisis kasus yang agak membosankan. Ini dia:
Kasus
max(abs(x),abs(y)) >= 1.
:x/2.+y/2.
memanipulasi mantra yang sama dan oleh karena itu memberikan jawaban yang sama persis seperti perhitungan(x+y)/2
akan menghasilkan jika kita mengasumsikan eksponen diperluas untuk mencegah luapan. Jawaban ini mungkin tergantung pada mode pembulatan tetapi dalam hal apa pun itu dijamin oleh IEEE754 untuk menjadi jawaban terbaik (dari kenyataan bahwa yang dihitungx+y
dijamin menjadi perkiraan terbaik untuk matematika x + y, dan pembagian dengan 2 tepat dalam hal ini kasus).Huruf x didenormalkan (dan sebagainya
abs(y)>=1
):answer = x/2. + y/2. = y/2. since abs(x/2.) is so tiny compared to abs(y/2.) = the exact mathematical value of y/2 = a best possible answer.
Subkotak y didenormalkan (dan sebagainya
abs(x)>=1
): analog.max(abs(x),abs(y)) < 1.
:x+y
adalah non-denormalized atau denormalized-and- "even": Meskipun dikomputasix+y
mungkin tidak tepat, itu dijamin oleh IEEE754 untuk menjadi perkiraan terbaik untuk matematika x + y. Dalam hal ini pembagian selanjutnya dengan 2 dalam ekspresi(x+y)/2.
adalah tepat, sehingga jawaban yang dihitung(x+y)/2.
adalah perkiraan terbaik untuk matematika (x + y) / 2.x+y
dan "ganjil": Dalam hal ini tepat salah satu dari x, y juga harus didenormalkan-dan- "ganjil", yang berarti yang lain dari x, y didenormalkan dengan tanda yang berlawanan, sehingga yang dihitungx+y
adalah persis matematis x + y, dan yang dihitung(x+y)/2.
dijamin oleh IEEE754 menjadi perkiraan terbaik untuk matematika (x + y) / 2.sumber
Untuk format titik-mengambang biner IEEE-754, dicontohkan oleh
binary64
perhitungan (presisi ganda), S. Boldo secara resmi membuktikan bahwa algoritma sederhana yang ditunjukkan di bawah ini memberikan rata-rata bulat yang benar.Sylvie Boldo, "Verifikasi formal program menghitung rata-rata titik mengambang." Dalam Konferensi Internasional tentang Metode Teknik Formal , hlm. 17-32. Springer, Cham, 2015. ( konsep online )
Karena pembagian oleh dua akurat dalam aritmatika titik-mengambang biner, asalkan aliran bawah tidak terjadi , tampaknya secara intuitif jelas bahwa dengan memilih salah satu dari dua rumus( x + y) / 2 dan x / 2 + y/ 2 yang sesuai (berdasarkan besarnya input) sekali harus mencapai rata-rata yang dibulatkan secara akurat. Kertas Boldo menunjukkan bahwa untuk IEEE-754 C∈ [ 2- 967, 2970] akan cukup. Orang mungkin memilihC sehingga memberikan kinerja terbaik untuk use case tertentu.
binary64
ada titik perpindahanIni menghasilkan
ISO-C99
kode teladan berikut :Dalam pekerjaan tindak lanjut baru-baru ini, S. Boldo dan rekan penulis menunjukkan bagaimana untuk mencapai hasil terbaik untuk format floating-point desimal IEEE-754 dengan memanfaatkan operasi FMA multiply-add (FMA) dan presisi yang terkenal. blok bangunan penggandaan (TwoSum):
Sylvie Boldo, Florian Faissole, dan Vincent Tourneur, "Algoritma yang Terbukti Secara Resmi untuk Menghitung Rata-Rata yang Benar dari Angka Titik Apung Desimal." Dalam Simposium IEEE 25 tentang Aritmatika Komputer (ARITH 25) , Juni 2018, hlm. 69-75. ( konsep online )
sumber
Meskipun ini mungkin bukan kinerja-efisien super-bijaksana, ada cara yang sangat sederhana untuk (1) memastikan tidak ada angka yang lebih besar dari salah satu
x
atauy
(tidak ada luapan) dan (2) menjaga titik mengambang sebagai "akurat" seperti mungkin (dan (3) , sebagai bonus tambahan, meskipun pengurangan digunakan, tidak ada nilai yang akan disimpan sebagai angka negatif.Bahkan, jika Anda benar - benar ingin mendapatkan akurasi, Anda bahkan tidak perlu melakukan pembagian di tempat; cukup kembalikan nilai-nilai
min(x, y)
dandifference
yang dapat Anda gunakan untuk menyederhanakan secara logis atau memanipulasi nanti.sumber
2,4,9
, tidak sama dengan rata-rata3,9
.x
dany
apakah floating-point, perhitungan Anda menghasilkan floating-point terdekat(x+y)/2
?Konversikan ke precission yang lebih tinggi, tambahkan nilainya di sana dan konversi kembali.
Seharusnya tidak ada overflow dalam precisi yang lebih tinggi dan jika keduanya berada dalam kisaran floating point yang valid, jumlah yang dihitung harus di dalam juga.
Dan itu harus di antara mereka, kasus terburuk hanya setengah dari jumlah yang lebih besar jika precisi tidak mencukupi.
sumber
Secara teoritis,
x/2
dapat dihitung dengan mengurangi 1 dari mantissa.Namun, sebenarnya menerapkan operasi bitwise seperti ini tidak selalu langsung, terutama jika Anda tidak tahu format angka floating point Anda.
Jika Anda dapat melakukan ini, seluruh operasi dikurangi menjadi 3 tambahan / pengurangan, yang seharusnya merupakan peningkatan yang signifikan.
sumber
Saya berpikir sepanjang jalan yang sama dengan @Rand Heath tetapi belum bisa berkomentar, inilah pendapat saya:
x/2
dapat dihitung dengan mengurangi 1 dari eksponen (bukan mantissa, mengurangi 1 dari mantissa mengurangi2^(value_of_exponent-length_of_mantissa)
dari nilai keseluruhan).Tanpa batasan kasus umum, mari kita asumsikan
x < y
. (Jikax > y
, beri label ulang variabel. Jikax = y
,(x+y) / 2
sepele.)(x+y) / 2
menjadix/2 + y/2
, yang dapat dilakukan oleh dua pengurangan integer (oleh satu dari eksponen)x
akan membuatx/2
lebih kecil dari yang dapat diwakili (dengan asumsi mantissa diwakili dengan terkemuka 1).x
,x
mantissa bergeser ke kanan dengan satu (dan tambahkan memimpin 1 implisit, jika ada).x
ke kanan sesuai dengan eksponeny
.x
telah bergeser sepenuhnya. Jika kedua eksponen minimal, yang memimpin akan meluap, yang ok, karena yang meluap seharusnya menjadi yang terdepan lagi.sumber