{"id":5291,"date":"2014-01-06T15:10:26","date_gmt":"2014-01-05T21:10:26","guid":{"rendered":"http:\/\/www.moonmile.net\/blog\/?p=5291"},"modified":"2014-01-06T15:17:44","modified_gmt":"2014-01-06T06:17:44","slug":"f-%e9%80%86%e8%a1%8c%e5%88%97%e3%82%92%e8%a8%88%e7%ae%97%e3%81%99%e3%82%8b","status":"publish","type":"post","link":"http:\/\/www.moonmile.net\/blog\/archives\/5291","title":{"rendered":"[F#] \u9006\u884c\u5217\u3092\u8a08\u7b97\u3059\u308b"},"content":{"rendered":"

F#\u3067\u6570\u5024\u30fb\u7dda\u5f62\u4ee3\u6570\u8a08\u7b97\u3092\u3059\u308b\u305f\u3081\u306e\u30e9\u30a4\u30d6\u30e9\u30ea\u7d39\u4ecb\uff08F# PowerPack, F# MathProvider\uff09
\nhttp:\/\/d.hatena.ne.jp\/teramonagi\/20111215\/1323874810#20111215fn6<\/a><\/p>\n

\u3092\u898b\u308b\u3068\u3001F# PowerPack<\/a> \u3068\u00a0F# MathProvider<\/a> \u3092\u4f7f\u3046\u3068\u3001\u9006\u884c\u5217\u304c\u7c21\u5358\u306b\u89e3\u3051\u307e\u3059\u3002<\/p>\n

[code lang=”csharp”]
\n#r "FSharp.PowerPack.dll";
\n#r @"D:toolsMathProviderNet 4.5MathProvider.dll";<\/p>\n

let A = matrix [[3.0; 1.0;];
\n [2.0; 5.0;]]<\/p>\n

let det = MathProvider.LinearAlgebra.det A
\nlet inv = MathProvider.LinearAlgebra.inv A
\n[\/code]<\/p>\n

\u3053\u3093\u306a\u98a8\u306b\u3059\u308b\u3068\u3001<\/p>\n

[code lang=”csharp”]
\nval A : matrix = matrix [[3.0; 1.0]
\n [2.0; 5.0]]
\nval det : float = 13.0
\nval inv : matrix = matrix [[0.3846153846; -0.07692307692]
\n [-0.1538461538; 0.2307692308]]
\n[\/code]<\/p>\n

\u3068\u51fa\u307e\u3059\u3002\u3081\u3067\u305f\u3057\u3081\u3067\u305f\u3057…\u306a\u306e\u306f\u3001\u5f8c\u304b\u3089\u77e5\u3063\u305f\u308f\u3051\u3067\u3001\u5b9f\u306f F# MathProvider \u306f\u5fc5\u305a Blas.dll \u3068 lapack.dll \u304c\u5fc5\u8981\u306a\u306e\u304b\u3068\u601d\u3063\u3066\u8e8a\u8e87\u3057\u3066\u305f\u3093\u3067\u3059\u3088\u306d\u3002\u3067\u304d\u308b\u3053\u3068\u306a\u3089\u3070\u3001F# \u3060\u3051\u3067\u4f5c\u308a\u305f\u3044\u3001\u3068\u3044\u3046\u304b\u3001\u30b9\u30c8\u30a2\u30a2\u30d7\u30ea\u3068\u304b\u81ea\u524d\u3067\u9ad8\u901f\u5316\uff08\u6709\u9650\u8981\u7d20\u6cd5\u30671\u4e07\u7bc0\u70b9\u3068\u304b\uff09\u3092\u8003\u3048\u308b\u3068\u3001\u4e2d\u8eab\u3092\u7406\u89e3\u3057\u3066\u304a\u304d\u305f\u3044\u3001\u3068\u3044\u3046\u9858\u671b\u304c\u3042\u3063\u305f\u306e\u3067…\u8eca\u8f2a\u306e\u518d\u767a\u660e\u3092\u3057\u3066\u3057\u307e\u3044\u307e\u3057\u305f\u3002<\/p>\n

[code lang=”csharp”]
\nlet inv ( A’ : matrix ) =
\n let A = A’.Copy()
\n let n = A.NumCols-1
\n let B = Matrix.identity A.NumCols
\n for k in 0..n do
\n \/\/ akk\u30921\u306b\u3059\u308b
\n let p = A.[k,k]
\n for j in 0..n do
\n A.[k,j] <- A.[k,j] \/ p
\n B.[k,j] <- B.[k,j] \/ p
\n \/\/ k\u884c\u4ee5\u5916\u304b\u3089\u5f15\u304f
\n for i in 0..n do
\n if i <> k then
\n let d = A.[i,k]
\n for j in 0..n do
\n A.[i,j] <- A.[i,j] – A.[k,j] * d
\n B.[i,j] <- B.[i,j] – B.[k,j] * d
\n B
\n[\/code]<\/p>\n

\u9006\u884c\u5217\u306e\u8a08\u7b97
\nhttp:\/\/www.asahi-net.or.jp\/~uc3k-ymd\/Lesson\/Section03\/invmat.html<\/a><\/p>\n

\u3092\u53c2\u8003\u306b\u3057\u3066\u3079\u305f\u3079\u305f\u306a\u624b\u8a08\u7b97\u306e\u65b9\u6cd5\u3092F#\u306b\u76f4\u3057\u3066\u3044\u307e\u3059\u3002\u7d50\u679c\u304c0\u306b\u306a\u308b\u3068\u3053\u308d\u3082\u8a08\u7b97\u3057\u3066\u3057\u307e\u3063\u3066\u7121\u99c4\u304c\u591a\u3044\u306e\u3067\u3001\u9ad8\u901f\u5316\u306e\u3057\u7532\u6590\u304c\u3042\u308b…\u3068\u3044\u3046\u304b\u3001Fortran \u3068\u304b C\u8a00\u8a9e\u3068\u304b\u306e\u30b3\u30fc\u30c9\u3082\u3042\u3061\u3053\u3061\u306b\u3042\u308b\u306e\u3067\u3001\u305d\u308c\u3089\u3092\u30b3\u30f3\u30d0\u30fc\u30c8\u3057\u305f\u307b\u3046\u304c\u826f\u3044\u304b\u3082\u3002
\n\u305f\u3060\u3057\u3001\u6709\u9650\u8981\u7d20\u6cd5\u3067\u4f7f\u3046\u9006\u884c\u5217\u306e\u5834\u5408\u306b\u306f\u3001\u5bfe\u79f0\u6027\u3092\u3046\u307e\u304f\u5229\u7528\u3057\u305f\u308a\u8a08\u7b97\u91cf\u3092\u6e1b\u3089\u3059\u3053\u3068\u304c\u3067\u304d\u308b\u306e\u3067\u3001\u4e00\u822c\u7684\u306a\u9006\u884c\u5217\u3088\u308a\u3082\u624b\u3092\u52a0\u3048\u305f\u307b\u3046\u304c\u5b9f\u7528\u7684\u304b\u3082\u3057\u308c\u307e\u305b\u3093\u3002\u304c\u3001ADVENTURE<\/a> \u306e\u3088\u3046\u306b\u5927\u898f\u6a21\u306a\u69cb\u9020\uff08\uff11\u5104\u7bc0\u70b9\u306a\u3069\uff09\u3068\u306f\u9055\u3063\u3066\u3001\u7bc0\u70b9\u6570\u306e\u30aa\u30fc\u30c0\u30fc\u304c\u9055\u3048\u3070\u3001\u6700\u8fd1\u306ePC\u306a\u3089\u3070\u3055\u3057\u3066\u9ad8\u901f\u5316\u305b\u305a\u3068\u3082\u7d50\u69cb\u306a\u3068\u3053\u308d\u307e\u3067\u3044\u3051\u308b\u306e\u3067\u306f\uff1f\u3068\u8003\u3048\u3066\u3044\u307e\u3059\u3002\u3053\u308c\u306f\u5148\u884c\u304d\u8a66\u3057\u3066\u307f\u305f\u3044\u3068\u3053\u308d\u3067\u3059\u3002<\/p>\n

\u5b9f\u306f\u3053\u306eF#\u306e\u30b3\u30fc\u30c9\u3001\u534a\u65e5\u307b\u3069\u304b\u304b\u3063\u305f\u3093\u3067\u3059\u3088\u306d\u3002\u306a\u304b\u306a\u304b\u6570\u5024\u304c\u5408\u308f\u306a\u304f\u3066\u82e6\u52b4\u3057\u305f\u306e\u306f\u3001F#\u306e\u6587\u6cd5\u306b\u99b4\u308c\u3066\u3044\u306a\u304b\u3063\u305f\u305f\u3081\u3082\u3042\u308b\u306e\u3067\u3059\u304c\u3001\u65e2\u5b58\u306e Fortran \u30b3\u30fc\u30c9\u3068\u306f\u5225\u306e\u5f62\u3067\u5b9f\u88c5\u3057\u3066\u624b\u8a08\u7b97\u305d\u306e\u307e\u307e\u3092\u79fb\u3057\u305f\u305f\u3081\u3067\u3059\u3002\u307e\u3042\u3001\u52c9\u5f37\u306b\u306a\u3063\u305f\u304b\u3089\u3044\u3044\u304b\u3002\u6b63\u6708\u3060\u3057\u3002<\/p>\n

\u3053\u308c\u304c\u300c\u3067\u304d\u305f\u30fc\u3063!!!\u300d\u5f8c\u306b\u6c17\u4ed8\u3044\u305f\u306e\u304c\u3001F# MathProvider \u306e\u5b9f\u88c5\u3067\u3059\u3002Blas.dll \u3092\u4f7f\u308f\u306a\u3044\u5834\u5408\u306b\u306f\u3001\u81ea\u524d\u306e F# \u30b3\u30fc\u30c9\u3092\u4f7f\u3046\u3088\u3046\u306b\u3067\u304d\u3066\u3044\u308b\u3093\u3067\u3059\u306d\u3002\u306a\u308b\u307b\u3069\u3002
\n\u4ee5\u4e0b\u304c\u629c\u7c8b\u3001LU\u5206\u89e3\u3068\u4e0a\u4e09\u89d2\u3092\u4f7f\u3063\u3066\u3044\u307e\u3059\u3002<\/p>\n

[code lang=”csharp”]
\n\/\/\/ Given A[n,n] find it’s inverse.
\n\/\/\/ This call may fail.
\nlet Inverse a =
\n if HaveService() then LinearAlgebraService.inverse a
\n else LinearAlgebraManaged.Inverse a<\/p>\n

let Inverse A =
\n let (n,m) = matrixDims A
\n if n <> m then invalid_arg "Matrix must be square when computing its inverse."
\n let P,L,U = LU A
\n (SolveTriangularLinearSystems U (SolveTriangularLinearSystems L ((Matrix.identity n).PermuteRows P) true) false)<\/p>\n

let LU A =
\n let (nA,mA) = matrixDims A
\n if nA<>mA then invalid_arg "lu: not square";
\n let L = Matrix.zero nA nA
\n let U = Matrix.copy A
\n let P = [| 0 .. nA-1 |]
\n let abs (x:float) = System.Math.Abs x
\n let swapR X i j = \/\/ REVIEW should we make this a general method?
\n let (nX,mX) = matrixDims X
\n let t = X.[i .. i,0 .. ]
\n for k=0 to mX-1 do
\n X.[i,k] <- X.[j,k]
\n X.[j,k] <- t.[0,k]
\n done<\/p>\n

for i=0 to nA-2 do
\n let mutable maxi = i \/\/ Find the biggest pivot element.
\n for k=i+1 to nA-1 do
\n if abs(U.[maxi,i]) < abs(U.[k,i]) then maxi <- k
\n done
\n \/\/let maxi = maxIndex (i+1) (nA-1) (fun k -> abs(U.[k,i]))<\/p>\n

if maxi <> i then
\n swapR U i maxi \/\/ Swap rows if needed.
\n swapR L i maxi \/\/ REVIEW can be made more performant.
\n let t = P.[i]
\n P.[i] <- P.[maxi]
\n P.[maxi] <- t<\/p>\n

for j=i+1 to nA-1 do
\n L.[j,i] <- U.[j,i] \/ U.[i,i]
\n for k=i+1 to nA-1 do
\n U.[j,k] <- U.[j,k] – L.[j,i] * U.[i,k]
\n done
\n U.[j,i] <- 0.0
\n done
\n done
\n (((*P.Length,*)Permutation.ofArray P),L + Matrix.identity nA,U)<\/p>\n

let SolveTriangularLinearSystems K B isLower =
\n if isLower then
\n let (nK,mK) = matrixDims K
\n let (nB,mB) = matrixDims B
\n if nK<>mK || nB<> nK then invalid_arg "Cannot do backward substitution on non-square matrices."
\n let X = Matrix.zero nK mK
\n for i=0 to nK-1 do
\n for k=0 to mB-1 do
\n let s = ref B.[i,k]
\n for j=0 to i-1 do
\n s := !s – K.[i,j] * X.[j,k]
\n done
\n X.[i,k] <- !s \/ K.[i,i]
\n done
\n done
\n X
\n else
\n let (nK,mK) = matrixDims K
\n let (nB,mB) = matrixDims B
\n if nK<>mK || nB<> nK then invalid_arg "Cannot do backward substitution on non-square matrices."
\n let X = Matrix.zero nK mK
\n for i=0 to nK-1 do
\n for k=0 to mB-1 do
\n let s = ref B.[nK-i-1,k]
\n for j=0 to i-1 do
\n s := !s – K.[nK-i-1,nK-j-1] * X.[nK-j-1,k]
\n done
\n X.[nK-i-1,k] <- !s \/ K.[nK-i-1,nK-i-1]
\n done
\n done
\n X
\n[\/code]<\/p>\n

\u307e\u3042\u3001\u6570\u5f0f\u306f\u7406\u89e3\u3057\u3066\u4f7f\u3063\u305f\u307b\u3046\u304c\u3044\u3044\u306e\u3067\u3001\u3057\u3070\u3089\u304f\u306f\u81ea\u524d\u3067\u5b9f\u88c5\u3068\u3044\u3046\u3053\u3068\u3067\u3002<\/p>\n","protected":false},"excerpt":{"rendered":"

F#\u3067\u6570\u5024\u30fb\u7dda\u5f62\u4ee3\u6570\u8a08\u7b97\u3092\u3059\u308b\u305f\u3081\u306e\u30e9\u30a4\u30d6\u30e9\u30ea\u7d39\u4ecb\uff08F# PowerPack, F# MathProvider\uff09 http:\/\/d.hatena.ne.jp\/teramonagi\/20111215\/1323874810# … \u7d9a\u304d\u3092\u8aad\u3080 →<\/span><\/a><\/p>\n","protected":false},"author":2,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"_jetpack_newsletter_access":"","_jetpack_dont_email_post_to_subs":false,"_jetpack_newsletter_tier_id":0,"_jetpack_memberships_contains_paywalled_content":false,"_jetpack_hide_featured_image":false,"_jetpack_featured_image_duplicate_dismissed":0,"_jetpack_feature_clip_id":0,"_jetpack_memberships_contains_paid_content":false,"footnotes":"","jetpack_post_was_ever_published":false},"categories":[59],"tags":[],"class_list":["post-5291","post","type-post","status-publish","format-standard","hentry","category-f"],"jetpack_sharing_enabled":true,"jetpack_featured_media_url":"","_links":{"self":[{"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/posts\/5291","targetHints":{"allow":["GET"]}}],"collection":[{"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/comments?post=5291"}],"version-history":[{"count":4,"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/posts\/5291\/revisions"}],"predecessor-version":[{"id":5293,"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/posts\/5291\/revisions\/5293"}],"wp:attachment":[{"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/media?parent=5291"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/categories?post=5291"},{"taxonomy":"post_tag","embeddable":true,"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/tags?post=5291"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}