読者です 読者をやめる 読者になる 読者になる

ゆるふわブログ

東京大学理科 I 類 2 年の学生です.プログラミング,大学の勉強,日常生活で感じたことをゆるふわに書いていけたらなと思います.技術的に拙いところがあっても温かい目で見守っていただければ幸いです.

Desmos で遊んでみた

Desmos っていう手軽にグラフが書けるサイトを見つけたのでちょっと遊んでみました.

Square Wave (Fourier Transform)

(-1)^n を加えるだけで円の動きに味が出る!

Cube

Desmos で 3 次元の表現はできないだろうかと思って少しやってみた.線分でやるとなんか重かったので,直線でやってみました.(見にくい)
 (a,b) を再生ボタン押して適当にアニメーションさせると,消失点がカメラの画面への射影  (x_c,y_c) と一致して不変なのがよくわかります.
これは立方体の奥に向かう直線の方向ベクトルが (0,0,1) だから.



ちなみにこの記事で使ってたのを見て知りました.グラフアートってすごい…

corollary2525.hatenablog.com

Lorentz 変換と双曲線

1 年生の夏学期,「基礎方程式とその意味」という講義で時期尚早にもやった Lorentz 変換.
導出の過程は追うことができましたが,直感的な理解ができずしっくりと来ていませんでした.
天下り的に双曲線関数を使うときれいになるということも,どのような背景があるのかさっぱりでした.

今日から,統計と多変量解析の本を読み終わってひと段落ついた (それについての記事はまた後日) ので,途中で挫折した相対論を勉強しようと本を読み始めました.
そして,このサイトを見つけてしまいました.

d.hatena.ne.jp

非常に美しく・わかりやすく Lorentz 変換を説明しています.
また,この資料でさらに理解が深まりました.

https://www.kaijo.ed.jp/wp-content/uploads/2016/02/2013summer-5_2.pdf

内積空間と絡めて説明しています.Lorentz 変換における不変量 (の絶対値) が Minkowski ノルムの 2 乗であることがわかります.

この感動を是非伝えるべく記事にまとめねば (笑) と思い,僕なりにまとめてみました.

github.com

数式混じりの記事は LaTeX にしたほうがいいという結論に至った.

LaTeX 練習

最近教習所に通ってあとはグダグダしててあんまりプログラミングしてないです()

今日あるの知らずにノー勉で突入した仮免学科試験が合格していたことがわかりました.
47/50 だったけど…

春休みはその他学習塾で個別指導のバイトをしていて,受験生が卒業してしまったので今持ってる生徒が高 2 の子 1 人なんですね.暇だし,1 人相手なら余裕あるし,何かできないかなと…

そこで,バイトの生徒のために LaTeX の練習がてら数学のプリントを作りました.TikZ パッケージで図も入れられるようになった.

github.com

数式のフォントは数学ガールのあれです.

何気に GitHub 使ったの初めてです… 使い方よくわかってない…

Haskell で ABC 001 を解く

qiita.com

この記事を読んだからやろうと思った.深い意味はない.

AtCoder Beginner Contest 001 - AtCoder Beginner Contest 001 | AtCoder

A - 積雪深差

main :: IO ()
main = do
  h1 <- readLn
  h2 <- readLn

  print $ h1-h2

B - 視程の通報

import Control.Applicative
import Text.Printf

main :: IO ()
main = do
  m <- (/1000) <$> readLn :: IO Double

  printf "%02.0f\n" $ cond m

cond :: Double -> Double
cond m
  | m < 0.1            = 0
  | m >= 0.1 && m <= 5 = m*10
  | m >= 6 && m <= 30  = m+50
  | m >= 35 && m <= 70 = (m-30)/5+80
  | m > 70             = 89

MultiWayIf 拡張が Compile Error で切れそうだった.

C - 風力観測

import Control.Applicative ((<$>))

main :: IO ()
main = do
  [deg, dis] <- map read . words <$> getLine
  let ansW = w $ roundAt 1 $ dis/60

  putStrLn $ if ansW == 0
    then "C 0"
    else dir (deg/10) ++ " " ++ show ansW

roundAt :: Int -> Double -> Double
roundAt d x
  | dec >= 0.5 = fromIntegral (int+1)/10^d
  | otherwise  = fromIntegral int/10^d
  where
    (int,dec) = properFraction $ x*10^d

dir :: Double -> String
dir deg
  | null cand = "N"
  | otherwise = snd . head $ cand
  where
    cand = filter (\(i,_) -> 11.25+i*22.5 <= deg && deg < 33.75+i*22.5) $ zip [0..] ["NNE","NE","ENE","E","ESE","SE","SSE","S","SSW","SW","WSW","W","WNW","NW","NNW"]

w :: Double -> Int
w dis = head $ filter (\i -> fst (lims!!i) <= dis && dis <= snd (lims!!i)) [0..12]
  where
    lims = [(0.0,0.2),(0.3,1.5),(1.6,3.3),(3.4,5.4),(5.5,7.9),(8.0,10.7),(10.8,13.8),(13.9,17.1),(17.2,20.7),(20.8,24.4),(24.5,28.4),(28.5,32.6),(32.7,200.0)]

Haskell の round は HALF_EVEN とかいうそうで偶数に丸めるらしく,挙動が望んだものと違ったので自分で実装せざるを得なかった.このサイトを参考にしました.
oropon.hatenablog.com

D - 感雨時刻の整理

import Control.Applicative ((<$>))
import Control.Monad (replicateM)
import Data.List (sort, splitAt)

type Section = (String, String)

main :: IO ()
main = do
  n <- readLn
  secs <- fmap sort $ replicateM n $ do
    (s,_:e) <- splitAt 4 <$> getLine
    return (roundTF s, roundTS e)

  mapM_ (\(s,e) -> putStrLn $ s ++ "-" ++ e) . reverse $ foldl f [] secs

  where
    f :: [Section] -> Section -> [Section]
    f [] s = [s]
    f (x:xs) s
      | x<&>s     = x<|>s:xs
      | otherwise = s:x:xs

show02d :: Int -> String
show02d i
  | i < 10    = "0" ++ show i
  | otherwise = show i

roundTF :: String -> String
roundTF t
  | m >= 60   = show02d (h+1) ++ show02d (m-60)
  | otherwise = show02d h ++ show02d m
  where
    roundI m = m`div`5*5

    h = read $ take 2 t
    m = roundI . read $ drop 2 t

roundTS :: String -> String
roundTS t
  | m >= 60   = show02d (h+1) ++ show02d (m-60)
  | otherwise = show02d h ++ show02d m
  where
    roundI m
      | m `mod` 5 /= 0 = (m`div`5+1)*5
      | otherwise      = m

    h = read $ take 2 t
    m = roundI . read $ drop 2 t

(<&>) :: Section -> Section -> Bool
(<&>) (a,b) (c,d) = c <= b

(<|>) :: Section -> Section -> Section
(<|>) (a,b) (c,d) = (min a c, max b d)

Data.List.Split が Compiler Error で splitOn "-" ってできなくて切れそうだった.

Haskell で競プロの Straightforward な問題をやるといい練習になる (?)

補足 - Text.Printf.printf の仕組み

printf はどうして可変長引数を取れるんでしょうか.
その秘密は PrintfType 型クラスにありました.
PrintfType 型クラスは,

class PrintfType r where
  printf :: String -> r

と定義されていて,PrintfType のインスタンスには,

instance PrintfType String
instance PrintfType (IO ())
instance (PrintfArg a, PrintfType r) => PrintfType (a -> r)

があります.最初の 2 つが printf の型の終端にあたります.
3 つ目のように,再帰的にインスタンス宣言することで,printf :: String -> a -> a -> … -> a -> (IO () or String) となります.よって,printf は IO アクションとしてだけでなく,C 言語でいう sprintf のように出力先が文字列でも使えます.
以下のサイトを参考にしました.
stackoverflow.com

連成振動と 2 次超曲面

テスト無事終わりました!
今日は振動・波動論のテスト中に気づいたことをまとめてみようと思います.

振動・波動論のテストの大問 1 は連成振動の問題でした.3 連成振子.3 つの質量 M のおもりが互いにバネで接続されているやつです.ここでは一般に N 連成振子を考えます.N 個のおもりの位置を  x_n\ (n=1,2,3,\cdots) とすると,運動方程式は,

 \begin{align*}
  M\ddot{x}_n &= -k(x_n-x_{n-1})+k(x_{n+1}-x_n)\\
  \ \\
  M \ddot{\mathbf{x}} &= -K\mathbf{x} \\
  \ \\
  \mathbf{x} &= \begin{pmatrix}
    x_1 \\
    x_2 \\
    \vdots \\
    x_N
  \end{pmatrix} \\
  \ \\
  K &= \begin{pmatrix}
    2k & -k & 0 & \cdots & \cdots & \cdots & 0 \\
    -k & 2k & -k & 0 & \cdots & \cdots & 0 \\
    0 & -k & 2k & -k & 0 & \cdots & 0 \\
    \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
    0 & \cdots & \cdots & 0 & -k & 2k & -k \\
    0 & \cdots & \cdots & \cdots & 0 & -k & 2k
  \end{pmatrix}
\end{align*}

となります. x_0 = x_{N+1} = 0 という形式的な変数を用意すれば,第 1 式のように綺麗に表されます.それを行列とベクトルでまとめたのが第 2 式です.規則性のおかげで, K は対称行列になります.よって,固有ベクトルを正規直交基底にしたものを束ねた直交行列で対角化できます.この直交行列を  P とすると,

 \begin{align*}
  P^{-1} K P &= \begin{pmatrix}
    \lambda_1 & 0 & \cdots & 0 \\
    0 & \lambda_2 & \cdots & 0 \\
    \vdots & \vdots & \ddots & \vdots \\
    0 & 0 & \cdots & \lambda_N
  \end{pmatrix} \\
  \ \\
  M\ddot{\mathbf{x}} &= -P\begin{pmatrix}
    \lambda_1 & 0 & \cdots & 0 \\
    0 & \lambda_2 & \cdots & 0 \\
    \vdots & \vdots & \ddots & \vdots \\
    0 & 0 & \cdots & \lambda_N
  \end{pmatrix} P^{-1} \mathbf{x} \\
  \ \\
  MP^{-1}\ddot{\mathbf{x}} &= -\begin{pmatrix}
    \lambda_1 & 0 & \cdots & 0 \\
    0 & \lambda_2 & \cdots & 0 \\
    \vdots & \vdots & \ddots & \vdots \\
    0 & 0 & \cdots & \lambda_N
  \end{pmatrix} P^{-1} \mathbf{x}
\end{align*}
よって,
 \begin{equation*}
  \mathbf{Q} = \begin{pmatrix}
    Q_1 \\
    Q_2 \\
    \vdots \\
    Q_N
  \end{pmatrix} = P^{-1} \begin{pmatrix}
    x_1 \\
    x_2 \\
    \vdots \\
    x_N
  \end{pmatrix} = P^{-1} \mathbf{x}
\end{equation*}
と基準座標を導入すると,
 \begin{align*}
  M\ddot{\mathbf{Q}} &= -\begin{pmatrix}
    \lambda_1 & 0 & \cdots & 0 \\
    0 & \lambda_2 & \cdots & 0 \\
    \vdots & \vdots & \ddots & \vdots \\
    0 & 0 & \cdots & \lambda_N
  \end{pmatrix} \mathbf{Q} \\
  \ \\
  M\ddot{Q}_n &= -\lambda_n Q_n
\end{align*}
となり,各方程式は単振動型の微分方程式になります.よって, Q_n は解けるので,基準座標の定義を遡って  \mathbf{x} = P \mathbf{Q} と解が求められます.つまり,各  x_n は角振動数  \sqrt{\frac{\lambda_1}{M}},\ \sqrt{\frac{\lambda_2}{M}},\cdots, \sqrt{\frac{\lambda_N}{M}} の単振動の重ね合わせ (線形結合) で表されることがわかります.

つまり,1 つの固有値が 1 つの基準振動数  \omega_n = \sqrt{\frac{\lambda_n}{M}} に対応しているわけですが,固有ベクトルは対応する基準振動の振幅比を表します.
これは, P が直交行列なので  P^{-1} = {}^t\!P であることから,固有ベクトル \mathbf{v_n} とすると,

 \begin{align*}
  P &= \begin{pmatrix}
    \mathbf{v_1} & \mathbf{v_2} & \cdots & \mathbf{v_N}
  \end{pmatrix} \\
  \ \\
  \mathbf{Q} &= P^{-1}\mathbf{x} \\
  &= {}^t\!P\mathbf{x} \\
  &= \begin{pmatrix}
    \mathbf{v_1} \\
    \mathbf{v_2} \\
    \vdots \\
    \mathbf{v_N}
  \end{pmatrix} \begin{pmatrix}
    x_1 \\
    x_2 \\
    \vdots \\
    x_N
  \end{pmatrix}
\end{align*}

となり,基準座標  Q_n \mathbf{v_n} = (v_{n1}, v_{n2}, \cdots, v_{nN}) の各成分で  x_n を線形結合した  \displaystyle\sum_{i=1}^N v_{ni}x_i で表されるからです.

で,テストで,基準座標で全エネルギー  E を表しなさいって問題が出たんですね.まぁ,基準振動のエネルギーの和だろうから,

 \displaystyle \begin{equation*}
  E = \sum_{n=1}^N \frac{1}{2}M \left(\dot{Q}_n^2 + \omega_n^2 Q_n^2\right)
\end{equation*}

で終わりにしようと思ってたんですが,もうちょっとちゃんと書いといたほうがいいのかなと思って.プリントに書いてあったことを思い出してたんですが… プリントには,連成振動の全エネルギーをまず  x_n の方で表してから, Q_n に持って行くみたいな流れだったんです.で, x_n で全エネルギーを表すと,

 \begin{equation*}
  E = \frac{1}{2}M{}^t\!\dot{\mathbf{x}}\dot{\mathbf{x}} + \frac{1}{2} {}^t\!\mathbf{x}K\mathbf{x}
\end{equation*}

って書いてあったんですね! 弾性エネルギーが 2 次超曲面の形で書けちゃうんです!
(タイトル回収)
すごく綺麗だと思いました. \mathbf{x} を一つのデータと捉えた時,2 次曲面はまさにその "2 次" の形なのだなぁと.

なんでこうなるのかちょっと考えてみました.
 x_n運動方程式 \dot{x}_n をかけると,力 × 速度で仕事率になるので,

 \displaystyle \begin{align*}
  M \ddot{x}_n \dot{x}_n &= +kx_{n-1}\dot{x}_n -2kx_n\dot{x}_n +kx_{n+1}\dot{x}_n \\
  \frac{d}{dt}\left(\frac{1}{2}M\dot{x}_n^2\right) &= +kx_{n-1}\dot{x}_n -\frac{d}{dt}\left(\frac{1}{2}2kx_n^2\right) +kx_{n+1}\dot{x}_n
\end{align*}

これを n が 1 から N まで足し合わせると, x_0 = x_{N+1} = 0 は消えて,積の微分から  \displaystyle\dot{x}_nx_{n+1} + x_n\dot{x}_{n+1} = \frac{d}{dt}\left(x_nx_{n+1}\right) とまとめられるので,

 \displaystyle \begin{align*}
  \frac{d}{dt}\left(\sum_{n=1}^N \frac{1}{2}M\dot{x}_n^2\right) &= \frac{d}{dt}\left(\sum_{n=1}^N \left(k x_{n}x_{n+1}-\frac{1}{2}2kx_n^2\right)\right) \\
  &= -\frac{d}{dt}\left(\frac{1}{2}{}^t\!\mathbf{x}K\mathbf{x}\right) \\
  \frac{d}{dt}\left(\frac{1}{2}M{}^t\!\dot{\mathbf{x}}\dot{\mathbf{x}} + \frac{1}{2} {}^t\!\mathbf{x}K\mathbf{x}\right) &= 0
\end{align*}

となります. K は対称行列ですから,i < j を満たす (i,j) 成分 (対角成分の上側) と i > j を満たす (i,j) 成分 (対角成分の下側) が同じです.よって, \dfrac{1}{2} {}^t\!\mathbf{x}K\mathbf{x} \dfrac{1}{2} によって 1 つ分が残ります.また,対角成分は 1 つだけなので, \dfrac{1}{2} によって,半分になります.第 1 式の右辺を見ると, \dfrac{d}{dt} の中で, x_nx_{n+1} には  \dfrac{1}{2} がかかっていませんが, x_n^2 の項にはかかっています.よって,これは 2 次超曲面の形で書けるのです*1

これを  \mathbf{x} = P\mathbf{Q} を代入して基準座標の式に変えます.

 \displaystyle \begin{align*}
  \frac{1}{2}M{}^t\!\dot{\mathbf{x}}\dot{\mathbf{x}} + \frac{1}{2} {}^t\!\mathbf{x}K\mathbf{x} &= \frac{1}{2}M{}^t\!{(P\dot{\mathbf{Q}})}P\dot{\mathbf{Q}} + \frac{1}{2} {}^t\!{(P\mathbf{Q})}KP\mathbf{Q} \\
  &= \frac{1}{2}M{}^t\!\dot{\mathbf{Q}}P^{-1}P\dot{\mathbf{Q}} + \frac{1}{2} {}^t\!\mathbf{Q}P^{-1}KP\mathbf{Q} \\
  \ \\
  &= \frac{1}{2}M{}^t\!\dot{\mathbf{Q}}\dot{\mathbf{Q}} + \frac{1}{2} {}^t\!\mathbf{Q}\begin{pmatrix}
    \lambda_1 & 0 & \cdots & 0 \\
    0 & \lambda_2 & \cdots & 0 \\
    \vdots & \vdots & \ddots & \vdots \\
    0 & 0 & \cdots & \lambda_N
  \end{pmatrix}\mathbf{Q}
\end{align*}

これで,形は少し違いますが,上にあげた基準座標でのエネルギー  E の表現が導けました.
ここで注目したいのが 2 次曲面がどうなったかです.上の方で説明したように,対角化を連立微分方程式を解くのに利用していると捉えることもできますが, 2 次超曲面によって弾性エネルギーを表すと,対角化は 2 次曲面  {}^t\!\mathbf{x}K\mathbf{x} を標準形  \displaystyle \sum_{n=1}^N \lambda_n Q_n^2 に変形していると捉えることができるのです!対角化の応用例というと,この連成振動の話と,微分積分学線形代数学を大学初年次で習っている身としては,Hessian による 2 変数関数の極値判定をあげたくなりますが,あれも, f_x(x,y) = f_y(x,y) = 0 を満たす停留点の周りで 2 次近似をして,2 次曲面の標準形を求め, f(x,y) = x^2+y^2,\ f(x,y) = -x^2-y^2,\ f(x,y) = x^2-y^2 の場合に帰着していると直感的には捉えることができますから*2,この話と繋がっていたのです.一見バラバラだった対角化の応用例が,実は同じ話だったと自分の中で繋がったのでとてもうれしかったです.微分積分学線形代数学,振動・波動論,今期取っていた 3 つの授業で暗に同じことを繰り返し教わっていたのだなぁとしみじみ思いました.

*1:2 次超曲面では, K の (i,j) 成分が  x_ix_j の係数に対応します.

*2:つまり,2 次超曲面の係数行列が  \begin{pmatrix} f_{xx} & f_{xy} \\ f_{yx} & f_{yy} \end{pmatrix} となるということです.

ブロック崩し (by Haskell)

Haskell から JavaScript に変換してくれる Haste というライブラリを使ってとりあえずブロック崩しを作って見ました.
忙しいので詳細は後日書きます.
以下ソースコード

{-# LANGUAGE LambdaCase #-}
import Haste
import Haste.DOM
import Haste.Events
import Haste.Graphics.Canvas
import Lens.Micro
import Control.Applicative
import Control.Arrow ((***))
import Control.Monad
import Data.Array
import Data.IORef
import System.Random

data BlkState = NoDamage | Hit | Destroy
  deriving (Eq, Enum)

data ScrState = Prepare | Game | Clear
  deriving Eq

data State = State {
  scrRef :: IORef ScrState,
  keyFRef :: IORef (Bool, Bool),
  bxRef :: IORef Double,
  bpRef :: IORef Point,
  bvRef :: IORef Point,
  missRef :: IORef Int,
  blocks :: IORef (Array (Int, Int) BlkState)
}

bnw, bnh :: Int
(bnw,bnh) = (12, 12)

main :: IO ()
main = do
  elemById "canvas" >>= \case
    Just cid -> do
      w <- read <$> getProp cid "width"
      h <- read <$> getProp cid "height"
      let (bw, bh) = (floor (w/fi bnw), floor (h/fi (2*bnh)))

      scrRef <- newIORef Prepare
      keyFRef <- newIORef (False, False)
      bxRef <- newIORef (-50)
      bpRef <- newIORef (0, h/2-70)
      theta <- randomRIO (-pi/3, pi/3) :: IO Double
      bvRef <- newIORef (5*cos (theta-pi/2), 5*sin (theta-pi/2))
      missRef <- newIORef 0
      arr <- newIORef $ listArray ((0,0),(bh,bw)) $ repeat NoDamage

      documentBody `onEvent` KeyDown $ \(KeyData code _ _ _ _) -> do
        when (code==37) $ modifyIORef keyFRef (&_1.~True)        -- Left
        when (code==39) $ modifyIORef keyFRef (&_2.~True)        -- Right
        scr <- readIORef scrRef
        when (code==13 && scr /= Clear) $ writeIORef scrRef Game -- Enter
      documentBody `onEvent` KeyUp $ \(KeyData code _ _ _ _) -> do
        when (code==37) $ modifyIORef keyFRef (&_1.~False) -- Left
        when (code==39) $ modifyIORef keyFRef (&_2.~False) -- Right
      
      getCanvas cid >>= \case
        Just cvs -> loop cvs (w,h) (bw,bh) (State scrRef keyFRef bxRef bpRef bvRef missRef arr)
        Nothing  -> error "Canvas could not be found!"
    Nothing  -> error "Canvas ID could not be found!"

loop :: Canvas -> Vector -> (Int, Int) -> State -> IO ()
loop cvs (w,h) (bw,bh) (State scrRef keyFRef bxRef bpRef bvRef missRef arr) = do
  let org = (w/2,h/2)

  scr <- readIORef scrRef
  bx <- readIORef bxRef

  when (scr /= Prepare) $ do
    (leftF, rightF) <- readIORef keyFRef
    when (leftF && bx > -w/2) $ modifyIORef bxRef (subtract 5)
    when (rightF && bx+fi bw-10 < w/2) $ modifyIORef bxRef (+5)

    bv <- readIORef bvRef
    modifyIORef bpRef ((+fst bv)***(+snd bv))

  (bpx, bpy) <- readIORef bpRef
  miss <- readIORef missRef
  blocks <- readIORef arr 

  render cvs $ translate org $ do
    color (RGB 40 40 40) $ fill $ rect (bx,h/2-50) (bx+fi bw-10,h/2-50+fi bh-5)
    color (RGB 40 40 40) $ fill $ circle (bpx,bpy) 10
    color (RGB 40 40 40) $ font "30px ヒラギノ角ゴ" $ text (-w/2+10,-h/2+38) $ "Miss: " ++ show miss
    forM_ [0..bnh-1] $ \i -> do
      forM_ [0..bnw-1] $ \j -> do
        let alpha = case blocks!(i,j) of
                      NoDamage -> 1.0
                      Hit -> 0.2
                      Destroy -> 0

        opacity alpha $ do
          color (hsv ((fi bnh-1-fi i+fi j)*360/fi bnh) 180 240) $ fill $ rect ((w-fi bnw*fi bw)/2+fi bw*fi j+5-w/2, 60+fi i*fi bh-h/2) ((w-fi bnw*fi bw)/2+fi bw*fi j+5-w/2+fi bw-10, 60+fi i*fi bh+fi bh-5-h/2)
          let (rx,ry) = ((w-fi bnw*fi bw)/2+fi bw*fi j+5-w/2+2, 60+fi i*fi bh-h/2+2)
          color (hsv ((fi bnh-1-fi i+fi j)*360/fi bnh) 120 240) $ lineWidth 5 $ stroke $ path [(rx,ry), (rx+fi bw-10-4,ry), (rx+fi bw-10-4,ry+fi bh-5-4), (rx,ry+fi bh-5-4), (rx,ry-2.5)]

    when (scr /= Game) $ do
      opacity 0.8 $ color (RGB 40 40 40) $ fill $ rect (-w/2,-h/2) (w/2,h/2)
      color (RGB 255 255 255) $ font "60px ヒラギノ角ゴ" $ do
        when (scr == Prepare) $ text (-450/2,0) "Enter To Start!"
        when (scr == Clear) $ text (-240/2,0) "Cleared!"

  when (bpx < -w/2+10 || bpx > w/2-10) $ modifyIORef bvRef (&_1%~negate)
  when (bpy < -h/2+10 || bpy > h/2-10) $ modifyIORef bvRef (&_2%~negate) 
  when (bpy > h/2-10) $ modifyIORef missRef (+1)
  reflect (bw,bh) (bx,h/2-50) (bpx,bpy) bpRef bvRef

  forM_ [0..bnh-1] $ \i -> do
    forM_ [0..bnw-1] $ \j -> when (blocks!(i,j) /= Destroy) $ do
      refF <- reflect (bw,bh) ((w-fi bnw*fi bw)/2+fi bw*fi j+5-w/2, 60+fi i*fi bh-h/2) (bpx,bpy) bpRef bvRef

      when refF $ writeIORef arr $ blocks // [((i,j), succ $ blocks!(i,j))]

  when (all (==Destroy) $ elems blocks) $ writeIORef scrRef Clear

  setTimer (Once 10) (loop cvs (w,h) (bw,bh) (State scrRef keyFRef bxRef bpRef bvRef missRef arr))
  return ()

reflect :: (Int, Int) -> Point -> Point -> IORef Point -> IORef Point -> IO Bool
reflect (bw,bh) (bx,by) (x,y) bpRef bvRef = do
  if (y+10 >= by && y-10 <= by+fi bh-5 && x+10 >= bx && x-10 <= bx+fi bw-10)
    then do
      modifyIORef bvRef (&_2%~negate)
      when (y+10 < by+fi bh-5) $ modifyIORef bpRef (&_2%~(subtract 1))
      when (y-10 > by) $ modifyIORef bpRef (&_2%~(+1))
      when (y+10 < by || y-10 > by+fi bh-5) $ modifyIORef bvRef (&_1%~negate)
      return True
    else return False

fi :: (Integral a, Num b) => a -> b
fi = fromIntegral

hsv :: Double -> Double -> Double -> Color
hsv h s v = case hi of
  0 -> fRGB v k m
  1 -> fRGB n v m
  2 -> fRGB m v k
  3 -> fRGB m n v
  4 -> fRGB k m v
  5 -> fRGB v m n
  where
    h' = h-(fromIntegral $ floor $ h/360)*360
    hi = fromIntegral $ floor $ h' / 60
    f = h' / 60 - hi
    m = v * (1-s/255)
    n = v * (1-s/255*f)
    k = v * (1-s/255*(1-f))
    fRGB :: Double -> Double -> Double -> Color
    fRGB r g b = RGB (floor r) (floor g) (floor b)

練習あるのみと思って下手くそなりに書いてみたわけだけど,絶対もっとうまく書けるよね… IORef 使う以外にないのかなぁ…
Haste サンプルが少なくて使い方よくわかってないけど,とりあえず描画とキーの取得さえできればゲームはできる (?)

f:id:Ysmr_Ry:20170119094626g:plain

sinh 型の置換積分

A セメスター期末試験まであと 1 週間,今日で大学 1 年生の全講義が終了しました.感慨深いものがあります.
テスト勉強で忙しいので手軽に書ける記事を.

電磁気の過去問を同じクラスの強い人が解いて答えを上げてくれたのですが,そこで,あまり見る頻度の少ない置換積分をしていたのでメモっとこうと.
直線電流の周りの磁場を求める積分です.アンペールの法則使えば一瞬とか言いっこなし.

 \displaystyle\begin{equation*}
  \frac{\mu_0 I}{4\pi}\int_{-\infty}^\infty \frac{R}{\sqrt{R^2+z^2}^3}\,dz
\end{equation*}
僕は  \displaystyle z=R\tan\theta,\,dz = \frac{R}{\cos^2\theta}\,d\theta と置換していました.
 \begin{align*}
  \frac{\mu_0 I}{4\pi}\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \frac{R}{\sqrt{R^2(1+\tan^2 \theta)}^3}\, \frac{R}{\cos^2 \theta}d\theta &= \frac{\mu_0 I}{4\pi R}\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \cos \theta\,d\theta \\
&= \frac{\mu_0 I}{4\pi R} \left[\sin \theta\right]_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \\
&= \frac{\mu_0 I}{2\pi R}
\end{align*}
しかし,  z=R\sinh t と置換することもできます.
 \begin{align*}
  \frac{\mu_0 I}{4\pi}\int_{-\infty}^\infty \frac{R}{\sqrt{R^2(1+\sinh^2 t)}^3}R\cosh t\,dt &= \frac{\mu_0 I}{4\pi R}\int_{-\infty}^\infty \frac{1}{\cosh^2 t}\,dt \\
&= \frac{\mu_0 I}{4\pi R} \left[\tanh t\right]_{-\infty}^\infty \\
&= \frac{\mu_0 I}{2\pi R}
\end{align*}
これは  \displaystyle \left(\sinh^{-1} \frac{t}{a}\right)' = \frac{1}{\sqrt{a^2+t^2}} に由来するものです.
式変形の途中双曲線関数の公式  \cosh^2 t-\sinh^2 t=1 を利用しました.
また, \displaystyle\lim_{t \to \pm\infty} \tanh t = \lim_{t \to \pm\infty} \frac{e^t-e^{-t}}{e^t+e^{-t}} = \pm1,\, (\tanh t)' = \frac{1}{\cosh^2 t} です.

以下が参考になるかも.

mathtrain.jp
mathtrain.jp