プログラミングのお題スレ Part22
お題:数値が入力されるのでその数値に最も近い回分数を出力せよ 回分数とは回分になっている数(負数含まず)のことである 最も近い回分数が2つある場合は2つとも出力せよ 入力 0 出力 0 入力 17 出力 22 入力 100 出力 99 出力 101 >>234 Rust fn foo(n: usize) -> (usize, Option<usize>) { let n2b = |n: usize| { let mut o = Some(n); iter::from_fn(|| { let n = o.take()?; o = (n >= 10).then(|| n / 10); Some((n % 10) as i8) }).collect::<Vec<i8>>() }; let b2n = |b: &[i8]| b.iter().rev().fold(0_usize, |n, b| n * 10 + *b as usize); let pal = |b: &mut [i8]| { let len = b.len() / 2; let (l, u) = b.split_at_mut(len); iter::zip(l, u.iter().rev()).for_each(|(l, u)| *l = *u); }; let inc = |b: &mut [i8]| { let len = b.len() / 2; let mut c = 1; b[len..].iter_mut().for_each(|b| { *b += c; if *b > 9 { *b = 0; c = 1; } else { c = 0; }}); }; let dec = |b: &mut [i8]| { let len = b.len() / 2; let mut c = 1; b[len..].iter_mut().for_each(|b| { *b -= c; if *b < 0 { *b = 9; c = 1; } else { c = 0; }}); }; let fix = |b: &mut [i8]| { if b.last() == Some(&0) { if b.len() & 1 == 0 { b[(b.len() - 1) / 2] = 9; } true } else { false } }; let mut b = n2b(n); pal(&mut b); let n1 = b2n(&b); match n.cmp(&n1) { Ordering::Equal => return (n, None), Ordering::Greater => inc(&mut b), Ordering::Less => dec(&mut b), } if fix(&mut b) { b.pop(); } pal(&mut b); let n2 = b2n(&b); match n.abs_diff(n1).cmp(&n.abs_diff(n2)) { Ordering::Less => (n1, None), Ordering::Greater => (n2, None), Ordering::Equal => (n1, Some(n2)), } } >>234 >>237 は入力が1〜9のとき出力が正しくなかった。function内の1行目に if ($n -le 9) {return $n} を 挿入すると修正される。 Rでは添字の開始値は1で添字0では空のデータが返るので、入力が1〜9のときの場合分けは不要。 []演算子と+演算子を文字列でも使えるように再定義した。 https://ideone.com/PP5swB Dでは添字範囲指定は半開区間なので、入力が1〜9のときの場合分けは不要。 https://ideone.com/hvNBia >>234 Perl5 for $n (0,17,100,123459321) { my %a; for (0..$n) { $i = $n - $_; $a{$i} = $i if 0 <= $i and $i =~ /^((\d)(?1)\2|\d?)$/; $j = $n + $_; $a{$j} = $j if $j =~ /^((\d)(?1)\2|\d?)$/; last if keys %a; } @a = keys %a; print "$n -> @a\n"; } ※見やすくするためインデントを全角スペースに置換してあります。 実行結果 $ perl 22_234_palindromic_number.pl 0 -> 0 17 -> 22 100 -> 99 101 123459321 -> 123464321 123454321 >>241 last if keys %a; } @a = keys %a; は last if @a = keys %a; } とコンパクトに書けるんだった、まぁいいや >>234 Perl5、小さい方の検索は0で止まるので負の値を避ける必要はなかった、書き直し。 $r = qr/^((\d)(?1)\2|\d?)$/; for $n (0,17,100,123459321) { my %a; for (0..$n) { $a{$n - $_} = 1 if ($n - $_) =~ $r; $a{$n + $_} = 1 if ($n + $_) =~ $r; last if @a = keys %a; } print "$n -> @a\n"; } >>234 Kotlin 何か画期的なアルゴリズムを使ったわけではなく、むしろほとんど何も考えずただ作られただけのプログラム。 https://paiza.io/projects/S5qsLnHz_pZD3um9jYRg_Q >>234 Python3 def f(k): s = str(k) return s == s[::-1] for n in [0, 17, 100, 123459321]: l = set() for i in range(n + 1): if f(n - i): l.add(n - i) if f(n + i): l.add(n + i) if l: print(n, l) break ※見易くするためインデントは全角空白に置換してあります 実行結果 $ python3 22_234_palindromic_number..py 0 {0} 17 {22} 100 {99, 101} 123459321 {123454321, 123464321} >>251 >[[0, [0]], [17, [11]], [100, [99, 101]]] 17 は、22 だよ >>235 しらみ潰しで失格 >>236 しらみ潰しで失格 >>240 しらみ潰しで失格 >>243 しらみ潰しで失格 >>244 しらみ潰しで失格 >>245 しらみ潰しで失格 >>246 しらみ潰しで失格 >>247 しらみ潰しで失格 >>248 しらみ潰しで失格 >>250 しらみ潰しで失格 >>251 しらみ潰しで失格 >>234 ruby https://ideone.com/N0w91j f = -> n { (0..n).lazy.map {|i| [n - i, n + i].select {|x| x.to_s.reverse.to_i == x}}.find(&:any?).uniq } >>252 (`・ω・´)ゞ 誤:a - 1, a + 1 正:a - 1, b + 1 >>234 dart https://ideone.com/e23wRv void main() { var rev = (n) => int.parse(n.toString().split('').reversed.join()); var f = (n) => Iterable.generate(n + 1).map((i) => [n - i, n + i].where((x) => x == rev(x))).firstWhere((a) => a.isNotEmpty).toSet().toList(); print([0, 17, 100].map((n) => [n, f(n)])); } しらみ潰しとは例えば1から順番に見つかるまで全てを試していく最悪な方法を指す 今回の場合だと与えた数から順番に見つかるまで±1を続けて全てを試していって探すプログラムが該当する >>262 私は >>248 だけれども、解法としてはそれしかないと思いますね >>249 と >>238 は しらみ潰しではなく きちんとプログラミングして算出しているようにみえますね >>264 私は解答は提出していないが、ざっくりと自分が思いついた方法 まず、以下のような操作を考える A. 1234という入力に対して1234321を返す B. 1234という入力に対して12344321を返す ここで、xという入力に対してA,Bが返す数をA(x),B(x)と表すことにする 次に、与えられた数の桁数で場合分け (1)与えられた数字の桁数が奇数の場合 例として5桁の数字を考える N=a*10000+b*1000+c*100+d*10+e*1 (a~eは1桁の自然数, aは0でない) が与えられたとき、 M=a*100+b*10+c*1 とすると、N=10000の場合を除いて、Nに最も近い回文数は A(M), A(M+1), A(M-1) の3つの候補に絞られる(厳密にはA(M)とNとの大小比較からA(M±1)の何れかは明らかに候補にならないので2つを考えれば良い) N=10000の場合は9999と10001が答え (2)与えられた数の桁数が偶数の場合 例として6桁の数を考える (1)と同様に N=a*100000+b*10000+c*1000+d*100+e*10+f*1 に対して M=a*100+b*10+c*1 とすると、N=100000の場合を除いて B(M), B(M+1), B(M-1) のどれかがNに最も近い回文数(厳密には以下略) N=100000の場合は99999と100001が答え 十分大きな数に対しては虱潰しに回文判定していくより速く求まる >>259 Iterable.generate(n).map(f)は単に Iterable.generate(n, f)で良かったと判明 >>234 >>249 をC++で書き換え(入力値は64ビット整数の範囲内限定) https://ideone.com/e1AM8A 元々はCで書き、4行目はなし、15行目と24行目はstrrev(s + i);だったが、Windowsのgccでは コンパイルできたのにideoneではできなかったので、仕方なくC++にしてstd::reverseで代用した。 >>234 Rust fn nearest_palindrome_numbers(n: usize) -> Vec<usize> { let mut dd = DecimalDigits::new(n); dd.palindrome_using_upper_half(); let n1 = dd.to_number(); match compare(n, n1) { Equal => return vec![n], Greater => dd.increment_upper_half(), Less => dd.decrement_upper_half(), } if dd.is_most_upper_zero() { return vec![n - 1, n + 1]; } dd.palindrome_using_upper_half(); let n2 = dd.to_number(); match compare_absolute_diff((n, n1), (n, n2)) { Less => return vec![n1], Greater => return vec![n2], Equal => return if n1 < n2 { vec![n1, n2] } else { vec![n2, n1] }, } } >>267 N=17=1*10+7*1のとき、Nは2桁(偶数桁)でM=1*1=1 B(M)=B(1)=11はNより小さいのでB(M-1)は考えなくてよい B(M+1)=B(2)=22なので11,22が答えの候補 11より22のほうが17に近いので22が答え ちょっとNに対するMの説明が足りてなかったけど言葉で上手く言い表せないすみません(上位半分以上かつ最小の桁数を抜き出す、的な) >>234 >>269 の一部でC++の機能をどうせ使ってしまったので、この際、全部をC++流に変えたら C流よりすっきり書けた。 https://ideone.com/38bo2E >>273 > [1000, [1001]] 誤:ps = [p.(s), p.(t.to_i.pred.abs.to_s + u), p.(t.succ + u)] 正:ps = [p.(s), p.(t.to_i.pred.abs.to_s + u), p.(t.succ + u), p.(s.to_i.pred.abs.to_s)] とりあえず雑に修正してみたが? (ノ∀`)アチャー >>262 >>245 のKotlinのプログラムは何も考えてなくて本当に馬鹿正直に±1して一つ一つ検査する方式で作ったんだけど、それでもあなたのテストではダメということになったの? まあ Int (符号付32bit整数) 使ってるからその限界超えたらダメではあるんだけど、そういう問題ではなく? >>276 それはしらみつぶしと言われる駄目プログラミングだよ 例えば求める解法が存在する方程式を解くのに値を±1しながら順に代入して試していくのと同じ >>277 あー。プログラムにバグがあってまともに答えが出ないっていうことではなく何の捻りもないプログラムだからダメっていう感想ね。それならわかる。 こちらもアルゴリズム思い浮かばないけどとりあえず作ってみただけだし。ダメというほどではないが良いとも思えないプログラムなので。 >>234 ruby https://ideone.com/E9VSE3 ・273の[1000, [1001]]バグ修正版 ・275とは違う方法で修正してみたがやっつけ感大 >>234 ruby 2.5.5 https://ideone.com/1zqSr1 ・いわゆる(?)ジェネレータ版 ・「終端を持たない範囲オブジェクト」はRuby 2.6.0から >>234 ruby 2.5.5 https://ideone.com/04fxGM ・ジェネレータ版ちょっとアレンジ ・to_sしてto_iするのをやめた 異なる自然数 a, b (a > b) における a^3 - b^3 を「a, b の三乗差」と呼ぶことにする。 異なる5通りの組(a, b) (c, d) ... (j, k) について三乗差がすべて相等しいとき その組(a, b)...(j, k) および三乗差自体を求めよ 異なる6通りの組で三乗差が相等しい場合があるかも検討せよ >>282 [(1134, 357), (1155, 504), (1246, 805), (2115, 2004), (4746, 4725)] a^3 - b^3 == 1412774811 >>282 aが5000以下限定で>>283 の解はRで5秒以内に求められた。8〜9行目は下三角行列の要素番号から 行番号と列番号を求めているだけで、>>282 を解く特別なアルゴリズムというわけではない。 https://ideone.com/w4X3DC >>282 ruby 2.5.5 https://ideone.com/4GdMaO ・285の無駄なループ回数を削減 ・でも5秒じゃ_ >>282 >>287 は下三角行列の要素番号から列番号を求めるのに2次方程式を解くのが分かりにくかったから、 各列の最下行の要素番号を計算しておいてそれを二分探索するように変更した。 https://ideone.com/Tefgkh n = 25000で実行してみても、n = 5000ときの解のa, bを両方とも2, 3, 4, 5倍した解しか別に 見つからなかった。 >>282 >>289 と似た手順でC++ https://ideone.com/nFRsrK Rではソートする前に重複値だけを抽出しているが、C++のunique関数はソート済みデータにしか 使えないので使っていない。 >>282 下三角行列の各列内の要素は昇順で既に並んでいるのに、>>290 は下三角行列の全要素を ソートして無駄なので、列のマージに変更(要するにマージソートを途中段階から開始) したら少し速くなった。 https://ideone.com/EZSvB3 n = 10000でも5秒以内に終わった。 https://ideone.com/huQiBe >>282 >>292 とは別の方法で>>290 を高速化 https://ideone.com/QxjmyT 1³, 2³, 3³, …, 5000³をD = 5001で割った余りはすべて異なる値になるから、d = a³ − b³を Dで割った余りはどれか1つの値に偏ることなく均等に分布する。dをDで割った余りによりdを 区分すれば、各区分に入る個数はどれも多すぎないのでソートに時間が余りかからない。 Dの値はconstexpr関数によりコンパイラに計算させている。n = 10000, 15000のときは それぞれD = 10002, 15009になる。 >>282 c https://ideone.com/JnJpJW ・291から省メモリ化 旧:unsigned int values[10]; 新:unsinged short values[4]; >>282 C++ https://ideone.com/1c4s5I >>294 はa, bの二重ループ内でa³ − b³をD = 5001で割った余りrにより区分していたが、 rのループ内でa, bを変化させるように変更したら、2次元配列がなくなってすっきりした。 その結果、メモリ使用量が激減し、nが大きい場合でも実行できるようになった。 >>296 の続き n = 1000000, m = 6で実行すると、12通りの解が見つかった。 [6つ組解] 424910390480793: (75978, 23919), (77385, 33768), (83482, 53935), (141705, 134268), (317982, 316575), (596001, 595602) 620174235433536: (86184, 27132), (87780, 38304), (90237, 48573), (94696, 61180), (160740, 152304), (360696, 359100) 1238805803151000: (107487, 14487), (108540, 34170), (110550, 48240), (119260, 77050), (454260, 452250), (851430, 850860) 1384074844012224: (112152, 29844), (125324, 83600), (130050, 93426), (159372, 138624), (224928, 215412), (357447, 353799) 1936290882196125: (127629, 52254), (133320, 75675), (149285, 111620), (228525, 215430), (246510, 235395), (290214, 282339) 4589726535576000: (170172, 69672), (177760, 100900), (185265, 120945), (304700, 287240), (328680, 313860), (386952, 376452) 4961393883468288: (172368, 54264), (175560, 76608), (180474, 97146), (189392, 122360), (321480, 304608), (721392, 718200) 11072598752097792: (224304, 59688), (250648, 167200), (260100, 186852), (318744, 277248), (449856, 430824), (714894, 707598) 36717812284608000: (340344, 139344), (355520, 201800), (370530, 241890), (609400, 574480), (657360, 627720), (773904, 752904) 52279853819295375: (382887, 156762), (399960, 227025), (447855, 334860), (685575, 646290), (739530, 706185), (870642, 847017) [7つ組解] 15490327057569000: (249281, 6281), (255258, 104508), (266640, 151350), (298570, 223240), (457050, 430860), (493020, 470790), (580428, 564678) 123922616460552000: (498562, 12562), (510516, 209016), (533280, 302700), (555795, 362835), (597140, 446480), (914100, 861720), (986040, 941580) 6つ組解の(2, 7), (4, 8), (5, 10), (6, 9)番目は各括弧内で自然数比になっている。 6つ組解の5番目の2倍は7つ組解の1番目のうちの6組を構成している。 >>282 C++ https://ideone.com/tM1cuo >>296 のrのループ内でa³ − b³をD2 = 5003で割った余りr2により区分し、それぞれの区分ごとに 解を探すようにしたら速くなった。ただし、nが大きい場合にはかえって遅くなる。 >>297 出題者です。 すごいです。ありがとうございます。私の手元ではまだ6通り解、7通り解のひとつも入手できていないので、参考になりました 私のアルゴリズムは効率が悪いようですね >>282 C++ https://ideone.com/LEU7EV >>298 でnを大きくするにつれ>>296 に対する高速化効果が薄れていくのは、ABをvectorでなく 配列にしたらある程度改善された。n = 5000のときの実行時間は>>296 の半分以下になった。 ただし、n = 1000000まで大きくすると、296よりやっぱり遅くなる。 >>299 どんなプログラムを書いたのか見せて。 >>282 C++ https://ideone.com/PG6UiY >>300 の実行時間を分析すると、最も時間が掛かっているのは46〜と47行目だと判明した。 そこで配列ABの第1次元と第2次元を入れ替えてみると、n = 5000では変わらないが、 1万, 2万, 5万, 10万, 20万では35%前後高速になった。これは、改良前には第2次元の添字が 小さい要素に書き込みが集中しているため、改良後のように第1次元に入れ替えた方が 纏まったメモリ領域に書き込みが集中しキャッシュの効きが良くなるからだと考えられる。 一方、n = 100万で高速化しないのは、書き込み集中領域が大きすぎるからだろう。 https://ideone.com/6RzW0n n = 100万の場合にはr2の値によってデータを多数の列へ振り分けるのをやめ、列を1つにして、 その内部でr2の値により2種類に区分し、それぞれの内部で2種類にさらに区分し、…と再帰的に 区分していけば(要するにクイックソートの変形版)、1つの配列内での要素のスワップだけで済み、 キャッシュの効きが改善されるとの予想通り、n = 100万で実行速度は>>296 より25%速くなった。 (原理的には>>300 より非効率なのでn = 5000では>>300 より当然遅い) >>301 ハズレが多いから2passは効果ある? >>282 C++ https://ideone.com/xQD1W8 関数mainのループで配列A, B, Pに書き込まずdにだけ書き込むようにし、関数FindDuplicatesで dの添字Pではなくdそのものをソートするように変えて、n = 1000000の場合に>>301 より10%高速化。 関数PrintSolutionでa, bをmainでと同じ方法で再計算するのは非効率だが、PrintSolutionは僅か12回しか 呼ばれないため、全体の実行時間への影響は無視できる。 お題 460円 580円 600円 の3種類の商品があります これらを組み合わせて合計10個買ったら5360円になりました 組み合わせを求めるプログラムを書いてください ちなみに答えの一つは ・600円×2 ・580円×4 ・460円×4 だそうです https://rio2016.5ch.net/test/read.cgi/cigaret/1706726196/56-57 >>304 面倒なので全て460円を引くと A=0円 B=120円 C=140円 10個で760円という問題 面倒なのでさらに20で割ると A=0円 B=6 C=7円 10個で38円という問題 つまり唯一奇数のCは偶数個が確定 Cが6個以上だと42円以上でオーバーしてNG Cが4個だと28円で残り10円をA,Bで作れないからNG Cが2個だと14円で残り24円はBが4個で残り4個がA Cが0個だと0円で残り38円をA,Bで作れないからNG つまり解は(A,B,C)=(4,4,2)しかない >>304 Rで全探索でなくちゃんと解くと https://ideone.com/F44pCL 解が複数ある場合と全くない場合の例として、600円を540円と520円に変更したときの出力も載せた。 2pass案は多少工夫したらかなり速い n ␣␣m ␣296␣ ␣301-1 ␣301-2 ␣303␣ ␣2pass 5k␣␣5 ␣ 0.5s ␣ 0.1s ␣ 0.5s ␣ 0.4s ␣ 0.1s 25k ␣5 ␣12.7s ␣ 2.5s ␣13.9s ␣11.1s ␣ 1.7s 100k␣5 ␣3m52s ␣49.3s ␣4m13s ␣3m26s ␣38.9s 1M* ␣6 ␣8h23m ␣2h50m ␣8h51m ␣6h43m ␣1h11m *n=100万は1万サンプルの部分ループ500k≦r<510kから100倍 >>301 の296と301-2の比較記述と違う傾向があるのはキャッシュ階層の違いだと思う 2passは301-1に近いけど1pass目でのランダムアクセスサイズを落としながらも 誤判定率を低く抑える(0.2%~2%)工夫をするのがお楽しみだと思う >>304 a = (600, 580, 460) m = min(a) h = set() def buy(b, yen): if yen < m: return for i in range(0, len(a)): v = a[i] if yen >= v: b[i] += 1 if yen == v: h.add(str(b)) else: buy(b, yen - v) b[i] -= 1 buy([0, 0, 0], 5360) for s in h: print(s) >>304 Haskell [(a, b, c) | a <- [0..20], b <- [0..20], c <- [0..20], a * 460 + b * 580 + c * 600 == 5360] output: [(0,2,7),(4,4,2)] 合計10個という条件忘れてた。 [(a, b, c) | a <- [0..20], b <- [0..20], c <- [0..20], a + b + c == 10, a * 460 + b * 580 + c * 600 == 5360] output: [(4,4,2)] 大幅に手直しした 特に前回数値が一部出てこない状態になっていたので色々と手動で最適化した 新しいアイディアを思いつかない限りはシングルスレッドでの限界に近いと思う n m 301-1 303 2pass 2pass' 5k 5 0.1s 0.4s 0.1s 0.1s 25k 5 2.5s 11.1s 2.3s* 1.7s 100k 5 49.3s 3m26s 38.9s 27.7s 1M* 6 2h50m 6h43m 1h11m 48m10s 2M* 6 17h06m 28h27m 5h47m 3h13m Max* 6 35h51m 51h23m 11h09m 5h47m *前回>>309 2pass n=25kの再計測値 *n=1Mは部分ループ500k<=r<510kから100倍 *n=2Mは部分ループ500k<=r<505kから400倍 *Max:=2642245は3乗がUINT64に収まる最大 *n=Maxは部分ループ500k<=r<500k+3785から2642245/3785倍 ヒント含みの数値がこちら n D1 D2 D3 = 5000 5001 5003 5009 false_positive = 23 / 5001 = 0.46% total_t_pass1 = 64.220 ms 2.568 ns/iter total_t_pass2 = 0.044 ms 0.381 ns/iter real 0m0.097s n D1 D2 D3 = 25000 25003 25005 25006 false_positive = 171 / 25003 = 0.68% total_t_pass1 = 1654.681 ms 2.647 ns/iter total_t_pass2 = 1.407 ms 0.329 ns/iter real 0m1.709s false_positive = 2211 / 100005 = 2.21% total_t_pass1 = 27338.298 ms 2.734 ns/iter total_t_pass2 = 78.402 ms 0.355 ns/iter real 0m27.692s n D1 D2 D3 = 1000000 1000002 1000009 1000015 false_positive = 18 / 10000 = 0.18% total_t_pass1 = 28674.338 ms 2.867 ns/iter total_t_pass2 = 5.642 ms 0.313 ns/iter real 0m28.897s n D1 D2 D3 = 2000000 2000003 2000013 2000015 false_positive = 13 / 5000 = 0.26% total_t_pass1 = 28777.424 ms 2.878 ns/iter total_t_pass2 = 8.620 ms 0.332 ns/iter real 0m29.015s n D1 D2 D3 = 2642245 2642246 2642253 2642258 false_positive = 315 / 3785 = 8.32% total_t_pass1 = 29210.857 ms 2.921 ns/iter total_t_pass2 = 336.864 ms 0.405 ns/iter real 0m29.800s >>304 ぶっちゃけ、他の言語の人と同じっぽくないので心配なんだが…。 自分なりにHaskellで全探索じゃないバージョン書いてみた。 Haskell [(a, b, c) | a <- [0..10], b <- [0..10 - a], c <- [0..10 - (a + b)], a * 460 + b * 580 + c * 600 == 5360, a + b + c == 10] 答えは同じ[(4,4,2)]。 少し一般化して、N個の商品があり、i番目の商品はA_i円です 合計M個購入し、価格の合計がS円であるような購入の仕方を998244353で割った余りを求めてください だとO(N M S)より小さい計算量で解けるのかな 2個、4個、8個…みたいにメモ化すればMはlogMにできるかもしれんね 空間がlogM倍されそうだが >>319 https://ideone.com/YCi6qe ヘロンが作ったもう1つの式である平方根を加算と除算の繰り返しで求める式も使用。 sqrt関数を呼び出すより実行形式ファイルサイズがほんの少しだけ小さくなる。 >>321 無限ループにならない? 機械イプシロン(DBL_EPSILON)とか気になる 解答例 #include <stdio.h> #include <math.h> void heron(double, double, double); int main(void) { double a, b, c; printf("3辺a, b, cを入力せよ "); scanf("%lf,%lf,%lf", &a, &b, &c); heron(a, b, c); } void heron(double x, double y, double z) // heronの定義 { double s, t; s = (x+y+z)*0.5; t = s*(s-x)*(s-y)*(s-z); printf("3角形の面積は S=%g\n", sqrt(t)); return; } >>322 p < 0 のとき(= 三角形を作れない場合)は浮動小数点数の特性に関係なく無限ループになる。 sqrt(p) と同様にNANを返すには、if (p < 0) return 0 / (p - p); を追加すれば良い。 p > 0 のときは無限ループにならないはず。以下が検証プログラム。 https://ideone.com/mzemEM x = sqrt(p), y = p / x とすると、浮動小数点数の特性により x == y とならない場合は存在する。 このとき、xとyの仮数部を整数と見なした値(以降では「仮数整数」と呼ぶ)の差は1なので、 z = (x + y) / 2 はxとyのうち仮数整数が偶数の方に一致する。zを新たなxとして代入しyとzを 再計算すれば、今度はxの仮数整数が偶数なのでzはxに必ず一致し、>>321 の収束判定条件が成立する。 具体例で見ると、p = 2 のときはxの仮数整数が奇数なので x != z となるが、zを新たなxとして代入し 再計算すれば x == z が成立する。桁上がりが起こる p = 3.9999999999999996 のときも、同様に 再計算で x == z が成立する。p = 3 のときはxの仮数整数が偶数なので x == z が成立し再計算は不要。 >>325 解説ありがとう 俺には理解できないレベルだと分かりましたw 俺なら収束の自信が無くてDBL_EPSILONを使った判定と ループ回数上限を組み合わせて実装しそうだ ヘロンの公式はそのままだと、数値計算での安定性が良くないらしいぞ 解決策は、Wikipediaの英語版の方に… tps://en.wikipedia.org/wiki/Heron%27s_formula#Numerical_stability の長さ この公式で計算される面積は、理論的には正しい値です。しかし、実際には、以下の理由で誤差が生じる可能性があります。 数値計算の誤差: 計算機で数値を扱う場合、有限桁しか扱えないため、丸め誤差が生じます。特に、辺の長さの値が大きく異なる三角形の場合、この誤差が顕著になります。 四捨五入誤差: 計算結果を小数点以下n桁まで表示する場合、n桁目以降の数字を切り捨てます。この四捨五入誤差も、面積の誤差に影響を与えます。 by Gemini >>327 に載っている参考文献 William M. Kahan, ‘Miscalculating Area and Angles of a Needle-like Triangle’ http://www.cs.berkeley.edu/ ~wkahan/Triangle.pdf のTable 1の問題がパソコン等でのC++プログラムでも再現されるか試してみた。 https://ideone.com/r4toUS Table 1とは違い、Accurate Δが概ね正確な場合にHeron's Δ'が大きく懸け離れた不正確な値に なってしまうことはなく、ほぼ同じ値になり差はごく僅かしかない。Table 1のような不安定性は Table 1の計算に使われたプログラマブル関数電卓に特有の問題で、パソコン等のプログラムでは 再現されない。(パソコン等のdoubleの方が精度が高いので当然と言えば当然だが) 一方、(a, b, c) = (5278.64055, 94721.35941, 99999.99996)の場合は、逆にHeron's Δ' = 0が 正確なのにAccurate Δ = 9.53674324543714が大きく懸け離れた不正確な値になってしまう 重大な欠点がある。これは、Accurate Δの式の根号内の第2因数c - (a - b)が正確には0なのに 3.63797880709171e-12と計算されてしまい、この誤差が他の因数との乗算により増幅されるから。 Heron's Δ'の式の根号内の第4因数s - cは0と計算されるので問題ない。 double向けの入力値(a, b, c) = (31622.77777777662, 0.000000000023, 31622.77777777661)を 作れば、Heron's Δ' = 2.30085990753844e-07, Accurate Δ = 3.20111707955507e-07となり、 相対差は確かに大きくなるが、200ビットで計算したほぼ正確な値3.27490470056059e-07から 見れば両方とも不正確だから、Accurate Δの利点はない。 だから、パソコン等のプログラムでは改良版の式を使う必要がないどころか使うべきではなく、 ヘロンの公式をそのまま使う方が良い。 >(a, b, c) = (5278.64055, 94721.35941, 99999.99996) >c - (a - b)が正確には0なのに3.63797880709171e-12と計算されてしまい この例に限らず、たいていの場合a,b,cはdoubleでexactに格納されて無くて この例では「c - (a - b)が正確には0」なのをチョイスしただけでは? パソロジカルな三角形をパラメトライズして面積を積分する検証はどう? 数式計算での正確な値 Heronで面積計算した時の数値積分 Accurateで面積計算した時の数値積分 を比べるのがフェアかなぁと > 200ビットで計算したほぼ正確な値3.27490470056059e-07 この例だけ見るとAccurate Δの方が優れているように見えるので >>331 の様なチェリーピックはどちらの計算式でも出来るので平均的に近似が近い方が精度的に優れているかと >>331 その通り。そして、(a, b, c) = (10000.1, 10000.2, 20000.3)とすれば、正しい面積は0なのに Heron's Δ' = 2.69745899635295とAccurate Δ = 1.34872949817647は両方とも大間違いになる。 この場合のようにHeron's Δ'での問題がAccurate Δで改善されないだけでなく、>>331 の引用の 場合のようにHeron's Δ'では結果的に問題ないのにAccurate Δでは新たな問題が生じてしまうのは、 参考文献の11ページで述べられた An algorithm stood convicted of numerical instability if it could be replaced by a new algorithm at least about as fast and accurate as the old for all data, and good for all data for which the old algorithm was bad. すべてのデータに対して旧アルゴリズムと少なくとも同じくらい高速かつ正確であり、 かつ旧アルゴリズムが悪くなるすべてのデータに対して良くなる新アルゴリズムによって 置き換えることができるとしたら、旧アルゴリズムは数値的に不安定と判定される。 という判定条件を満たさないから、Accurate Δは改良版としての適性を欠く。 >>333 その例では有効桁数がHeron's Δ'は0桁、Accurate Δは1桁しかなく、どちらの品質も絶対的に 劣悪で、それらの間の相対的な優劣に大した意味はない。 そもそも針のように異様に細長い三角形が重箱の隅をつつくような話で、普通はそんな場合は 想定しなくても良く、ヘロンの公式で充分。そこを敢えてつつくなら、ヘロンの公式だけでなく 改良式もぼろが出てしまうだけ。 read.cgi ver 07.5.4 2024/05/19 Walang Kapalit ★ | Donguri System Team 5ちゃんねる