{"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

\r\n#r "FSharp.PowerPack.dll";\r\n#r @"D:toolsMathProviderNet 4.5MathProvider.dll";\r\n\r\nlet A = matrix [[3.0; 1.0;];\r\n                [2.0; 5.0;]]\r\n\r\nlet det = MathProvider.LinearAlgebra.det A\r\nlet inv = MathProvider.LinearAlgebra.inv A\r\n<\/pre>\n
\u3053\u3093\u306a\u98a8\u306b\u3059\u308b\u3068\u3001<\/p>\n
\r\nval A : matrix = matrix [[3.0; 1.0]\r\n                         [2.0; 5.0]]\r\nval det : float = 13.0\r\nval inv : matrix = matrix [[0.3846153846; -0.07692307692]\r\n                           [-0.1538461538; 0.2307692308]]\r\n<\/pre>\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
\r\nlet inv ( A' : matrix ) =\r\n    let A = A'.Copy()\r\n    let n = A.NumCols-1\r\n    let B = Matrix.identity A.NumCols\r\n    for k in 0..n do\r\n        \/\/ akk\u30921\u306b\u3059\u308b\r\n        let p = A.[k,k]\r\n        for j in 0..n do\r\n            A.[k,j] <- A.[k,j] \/ p\r\n            B.[k,j] <- B.[k,j] \/ p\r\n        \/\/ k\u884c\u4ee5\u5916\u304b\u3089\u5f15\u304f\r\n        for i in 0..n do\r\n            if i <> k then\r\n                let d = A.[i,k]\r\n                for j in 0..n do\r\n                    A.[i,j] <- A.[i,j] - A.[k,j] * d\r\n                    B.[i,j] <- B.[i,j] - B.[k,j] * d\r\n    B\r\n<\/pre>\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
\r\n\/\/\/ Given A[n,n] find it's inverse.\r\n\/\/\/ This call may fail.\r\nlet Inverse a = \r\n  if HaveService() then LinearAlgebraService.inverse a\r\n                   else LinearAlgebraManaged.Inverse a\r\n\r\nlet Inverse A =\r\n    let (n,m) = matrixDims A\r\n    if n <> m then invalid_arg &quot;Matrix must be square when computing its inverse.&quot; \r\n    let P,L,U = LU A\r\n    (SolveTriangularLinearSystems U (SolveTriangularLinearSystems L ((Matrix.identity n).PermuteRows P) true) false)\r\n\r\n    let LU A =\r\n        let (nA,mA) = matrixDims A\r\n        if nA<>mA then invalid_arg &quot;lu: not square&quot;;\r\n        let L = Matrix.zero nA nA\r\n        let U = Matrix.copy A\r\n        let P = [| 0 .. nA-1 |]\r\n        let abs (x:float) = System.Math.Abs x\r\n        let swapR X i j =                           \/\/  REVIEW should we make this a general method?\r\n            let (nX,mX) = matrixDims X\r\n            let t = X.[i .. i,0 .. ]\r\n            for k=0 to mX-1 do\r\n                X.[i,k] <- X.[j,k]\r\n                X.[j,k] <- t.[0,k]\r\n            done\r\n            \r\n        for i=0 to nA-2 do\r\n            let mutable maxi = i        \/\/  Find the biggest pivot element.\r\n            for k=i+1 to nA-1 do\r\n                if abs(U.[maxi,i]) < abs(U.[k,i]) then maxi <- k\r\n            done\r\n            \/\/let maxi = maxIndex (i+1) (nA-1) (fun k -> abs(U.[k,i]))\r\n            \r\n            if maxi <> i then\r\n                swapR U i maxi     \/\/ Swap rows if needed.\r\n                swapR L i maxi     \/\/ REVIEW can be made more performant.\r\n                let t = P.[i]\r\n                P.[i] <- P.[maxi]\r\n                P.[maxi] <- t\r\n            \r\n            for j=i+1 to nA-1 do\r\n                L.[j,i] <- U.[j,i] \/ U.[i,i]\r\n                for k=i+1 to nA-1 do\r\n                    U.[j,k] <- U.[j,k] - L.[j,i] * U.[i,k]\r\n                done\r\n                U.[j,i] <- 0.0\r\n            done\r\n        done\r\n        (((*P.Length,*)Permutation.ofArray P),L + Matrix.identity nA,U)\r\n\r\n    let SolveTriangularLinearSystems K B isLower =\r\n        if isLower then\r\n            let (nK,mK) = matrixDims K\r\n            let (nB,mB) = matrixDims B\r\n            if nK<>mK || nB<> nK then invalid_arg &quot;Cannot do backward substitution on non-square matrices.&quot;\r\n            let X = Matrix.zero nK mK\r\n            for i=0 to nK-1 do\r\n                for k=0 to mB-1 do\r\n                    let s = ref B.[i,k]\r\n                    for j=0 to i-1 do\r\n                        s := !s - K.[i,j] * X.[j,k]\r\n                    done\r\n                    X.[i,k] <- !s \/ K.[i,i]\r\n                done\r\n            done\r\n            X\r\n        else\r\n            let (nK,mK) = matrixDims K\r\n            let (nB,mB) = matrixDims B\r\n            if nK<>mK || nB<> nK then invalid_arg &quot;Cannot do backward substitution on non-square matrices.&quot;\r\n            let X = Matrix.zero nK mK\r\n            for i=0 to nK-1 do\r\n                for k=0 to mB-1 do\r\n                    let s = ref B.[nK-i-1,k]\r\n                    for j=0 to i-1 do\r\n                        s := !s - K.[nK-i-1,nK-j-1] * X.[nK-j-1,k]\r\n                    done\r\n                    X.[nK-i-1,k] <- !s \/ K.[nK-i-1,nK-i-1]\r\n                done\r\n            done\r\n            X\r\n<\/pre>\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_memberships_contains_paid_content":false,"footnotes":""},"categories":[59],"tags":[],"class_list":["post-5291","post","type-post","status-publish","format-standard","hentry","category-f"],"jetpack_featured_media_url":"","jetpack_sharing_enabled":true,"_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}]}}