複天一流:どんな手を使ってでも問題解決を図るブログ

宮本武蔵の五輪書の教えに従い、どんな手を使ってでも問題解決を図るブログです(特に、科学、数学、工学の問題についてですが)

東大数学2024問題6 (part 11):細かい場合分の作業に突入....

前回のあらすじ

$g(x) = xk(x), k(x)=x^2+ax+b$が与えられており、$x\rightarrow n\in\mathbb{Z}$の場合に($\mathbb{Z}$は整数の集合)、$g(n)$が素数となるかどうか($g(n)\in\mathbb{P}?$)を調べる問題をずっとやっているが、前回は、$a,b\in\mathbb{Z}$を固定した場合に、 $k(p)=1$かつ$k(-p')=-1$を満たす$p,p'\in\mathbb{P}$を4つ見つけることはできないことを証明した。ということで、考察すべき残りの状況は以下の5つとなった。

  • (a) ケース3が2つの整数解を与え、ケース1とケース2が残りの2つの整数解を与えて4つとなる場合。
  • (b) ケース3が2つの整数解を与え、ケース4が一つの整数解、ケース1あるいはケース2が残りの整数解を与えて4つとなる場合。
  • (c) ケース3をケース4と言い換えた場合の、(a)の場合。
  • (d) ケース3をケース4と言い換えた場合の、(b)の場合。
  • (e) ケース1とケース2が1つずつの整数解を与え(合計2つ)、さらにケース3が一つ、ケース4が一つの解を与える場合。

本日は(a)と(b)の場合について調べてみたい。

(a)の場合とは?

ケース3というのは、$x=p\in\mathbb{P}$, かつ$k(p)=1$となる場合のことである。ただし、$\mathbb{P}$は素数の集合を表す。したがって、$g(p) = p$となるから与えられた条件を満たすというわけである。

前に数値的に探した時は、$(a,b)=(-5,7)$という具体例が一つ見つかった。この場合、$k(x)=x^2-5x+7$であるから、$k(2)=4-10+7=1$、および$k(3)=9-15+7=1$となり、$g(2)=2, g(3)=3$となって与えられた条件を満たす。

一方でケース4というのは、$x=-p'$(ただし$p'\in\mathbb{P}$)の場合に、$k(-p')=-1$となる場合のことである。この場合は、$g(-p')=(-p')(-1)=p'$となって与えられた条件を満たすというわけであるが、ケース3とケース4を同時に実現するような$(a,b)$は存在しないことを前回証明したのである。

問題を解く観点からすると蛇足に近い形になってしまうが、ケース4についても条件を満たす2つの異なる素数$p,p'$がみつかるかやってみよう。ケース4を満たすための「必要条件」、$k(-p')=-1$が$p'$についての2次方程式となり、異なる2つの正の実数解を与える条件、は \begin{equation} b < \frac{a^2}{4}-1, \quad a > 0, \quad b>-1 \end{equation} であった。下の図で紫色の影で塗られているのが該当する領域である(ちなみに緑色の影で塗られているのがケース3において2つの正の実数解を持てる条件に対応する領域)。

この図で確認できる格子点を探し、上の図に重ねてみる。 境界上の格子点は条件からは外さなければならない(重解となるため、解の個数が不足してしまう)が、一応計算だけはやっておくことにする。

前回作ったpythonのプログラム(discriminant.py)を再利用することができる。

まずは境界上にある格子点から計算してみる。データの見方は

通し番号, a, b, 判別式, 平方数かどうか、整数解1、整数解2

である。今考えているのはケース4なので、整数解はすべて負の整数(素数にマイナスをかけたもの)になってないといけない。

   1  0 -1 0 True 0.0 0.0
     2  1 -1 1 True 0.0 -1.0
     3  2 -1 4 True 0.0 -2.0
     4  2 0 0 True -1.0 -1.0
     5  3 -1 9 True 0.0 -3.0
     6  4 -1 16 True 0.0 -4.0
     7  4 3 0 True -2.0 -2.0
     8  5 -1 25 True 0.0 -5.0
     9  Programme completed...

整数解が出ていたとしても、ケース4の場合は(絶対値が)素数でないといけないから、0.0および1.0になってしまったものを除外しないといけない。また、重解のものも外れる。ということで、予想通り境界上の格子点に条件を満たすものは存在していないことが確認できた。

次は、ちゃんと境界内部に存在する格子点(のいくつか)を確認してみよう。

     1  3 1 1 True -1.0 -2.0
     2  4 2 4 True -1.0 -3.0
     3  5 3 9 True -1.0 -4.0
     4  5 5 1 True -2.0 -3.0
     5  Programme completed..

すでに説明したおおり、整数解として-1.0が含まれている場合は除外する。したがって、4行目の$(a,b)=(5,5)$が条件に該当するものであることがわかる。検算してみよう。

\begin{equation} k(-p)=p^2-5p+5 \end{equation} に対して、$k(-2)=4-10+5=-1$, また$k(-3) =9-15+5=-1$となって条件が満たされていることが確認できた。

もちろん、この解はケース3の領域に含まれていないからだめである。といっても、実数解が2つ出るところまでは大丈夫であるが、二つの正の解が出てこないという点で失敗する。確かめてみよう。

$k(p)=p^2+5p+5=1$となる$p$が素数になるかどうかである。$p$についての方程式を解くと、 \begin{equation} p = \frac{-5\pm\sqrt{25-16}}{2} = \frac{-5\pm 3}{2} = -4, -1 < 0 \end{equation} となってダメである(といっても整数解が出た部分までは「すごい」けれど、素数でもないし、正数でもないからダメはダメである)。

ケース1の場合

ケース1とは、$x\rightarrow n=1$となり、$k(1)=1+a+b$が素数になる場合である。素数になるかどうかの条件式を与えるのは面倒なので、まずは素数の必要条件である「正の数」という条件を課して、領域を狭めてみたい。$k(1)>0$であるから、 \begin{equation} b > -a -1 \end{equation} という不等式が条件に対応する。

ケース2の場合

ケース2は$x\rightarrow n=-1$となり、$k(-1) = 1-a+b$が素数に負符号をかけた値となる場合、つまり$k(-1)=-p$、ただし$p\in\mathbb{P}$の場合である。この場合も「必要条件」として$k(-1)<0$を課すことができるから、 \begin{equation} b < a-1 \end{equation} という不等式の形の条件を得ることができる。

ケース1とケース2が同時に成立すること

ケース4、あるいはケース3で見つかった$(a,b)$に対し、$k(1)$と$k(-1)$が同時に(異なる)素数となれば4つの整数解が見つかってしまい、証明は破綻する。まずは$k(1)$と$k(-1)$が同時に正数となる領域を調べてみよう。これは共通の切片$b=-1$をもつ直線グラフ(傾きは$\pm 1$)の「右側」の領域となる。正確な場所は下図の黄色で塗られた領域である。

黄色の領域は緑の領域とは絶対に共通領域を持たないことがわかる。つまり、ケース3で2つの解が見つかりつつ、ケース1とケース2で残りの2つの解が提供されることはあり得ないことが証明されたことになる。

つまり(a)の場合、4つの整数解が見つかることはないのである。

(c)の場合

(c)の場合、というのは「ケース3」を「ケース4」と言い換えた場合の(a)の場合である。この可能性はまだ否定されていない。つまり、ケース4で2つの解が見つかり、かつケース1とケース2で残りの二つの解が補充される場合である。この状況をわかりやすく領域を塗り直すと次のようになる。 紫と黄色に塗られた領域の和がそれである。

紫色の領域

この領域は閉じた領域であり、目で見て数えてしまうことが可能で、$(a,b)=(3,0), (3,1), (4,0),(4,1),(4,2)$であるが、これらの点はすでに条件を満たさないことを計算して確認済みである(境界上の点も不適であることをすでに確かめた)。つまり、可能性があるとするならば、黄色の領域に限られる。

黄色の領域

黄色の領域も、$a=5$の格子点、つまり$(a,b)=(a,b)$、ただし$b=0,1,2,3,4$、については条件を満たさないことを確認済みである($b=5$の場合のみが条件をかろうじて満たしたが、黄色の領域の外側になってしまった)。

したがって、$a\ge 6$の格子点で、$0\le b\le a-1$を満たすものだけを調べればよい。

東大数学2024問題6 (part 10):ケース3とケース4が同時に成立して4つの整数解が見つかる場合

前回のあらすじ

前回は、ケース3.の場合について詳しくみてみた。(a,b)=(-5,7)が候補となったが、$f(2)=2, f(3)=3$だけが題意を満たす場合であることがわかり、4つの素数を見つけるどころか、3つの素数が与えられた式を満たす状況にはなっていないことが確認できた。

ただ、具体例が一つ見つかったことで、何を証明してよいかぼんやりと見えてきた。今回は、ケース3.が2つの素数を生み出す場合、ケース4とは相入れないことを示したい。つまり、ケース3が2つの整数、ケース4が2つの整数をそれぞれ解として与えることで、全体として4つの整数が見つかるということは絶対にない、ことを示したいと思う。

ケース4の状況

この場合は、素数$p$に対し,$x=-p$によって$g(-p)$が素数になる場合に相当する。このとき、$g(-p)=(-p)k(-p) =(-p)(-1)=p$が成り立つ。したがって、条件式は \begin{equation} k(-p) = -1 \rightarrow p^2 -ap + b = -1 \rightarrow p^2-ap + b+1 =0 \end{equation} 解の公式により、 \begin{equation} p=\frac{a\pm\sqrt{a^2-4(b+1)}}{2} \end{equation}

したがって、まずは判別式が正となる条件が得られる。それは \begin{equation} a^2-4(b+1) > 0 \rightarrow b < \frac{a^2}{4}-1 \end{equation} である。

まずはこの条件をケース3が2つの整数を与える条件と重ねてみることにする。 色が濃くなっている場所が該当領域である。前回の考察で(a,b)=(-5,7)が候補となったが、ケース4が複素数になってしまったことを思い出そう。上のグラフでは今回求めた放物線の上の領域に(-5,7)が位置しており、複素数になってしまった理由がはっきりわかる。

上の領域をわかりやすく表示すると、次のようになる。

次に2つの解が両方ともに正となる条件を足す。 \begin{equation} -a < \sqrt{a^2-4(b+1)} < a \end{equation} この条件は$a>0$の時だけ成立する。そのとき、左側の不等式は自明な関係式となるので、右側のみを考慮すると \begin{equation} b > -1 \end{equation} を得る。

すでに$a>0$の条件が見つかった段階で、上のグラフとの共通領域は消滅してしまっている。つまり、ケース3とケース4を組み合わせて、4つの整数解を与えるような(a,b)が存在しないことが証明できた。

考察すべき、残りの状況

となると、考察すべき状況として残ったのは次のような場合である。

  • (a) ケース3が2つの整数解を与え、ケース1とケース2が残りの2つの整数解を与えて4つとなる場合。
  • (b) ケース3が2つの整数解を与え、ケース4が一つの整数解、ケース1あるいはケース2が残りの整数解を与えて4つとなる場合。
  • (c) ケース3をケース4と言い換えた場合の、(a)の場合。
  • (d) ケース3をケース4と言い換えた場合の、(b)の場合。
  • (e) ケース1とケース2が1つずつの整数解を与え(合計2つ)、さらにケース3が一つ、ケース4が一つの解を与える場合。

以上である。

東大数学2024問題6 (part 9):gnuplotで見つけた格子点の確認

前回のあらすじ

東大入試問題の数学2024問題6の後半にとりかかった。まずは二次方程式の解が2つとも正となる条件について調べることにした。条件を満たす領域がわかり、数値的にいくつかの格子点を(手で)抜き出してみたら、次のような感じとなた。

赤点が条件を満たす格子点

これが果たして、素数を生み出すかどうか検算してみる。

判別式を計算してみる

ところで判別式は英語でDiscriminantというそうだ。だから$D=b^2-4ac$とDで表すのであろう。

さて、計算結果をじっと睨んで、手で抜き出した格子点をデータファイルにしてみた。これをpythonのプログラムに読み込ませて、判別式を計算し、ちゃんと素数が出てくるか調べてみたい。ということで、プログラム名はdiscriminant.pyとする。

import math

def discr(a,b):
    return a**2 - 4*(b-1)

while True:
    try:
        ab_data = input().split()
        a = int(ab_data[0])
        b = int(ab_data[1])
        d = discr(a,b)
        sqd = math.sqrt(d)
        if sqd.is_integer() == True:
            sol0 = (-a + sqd)/2.
            sol1 = (-a - sqd)/2.
            print(a,b,d,sqd.is_integer(),sol0,sol1)
    except EOFError:
        break
print("Programme completed...")

真ん中あたりにあるif文は、判別式が「平方数」になっているかを確認している。つまり、うまい具合に平方根が開けて、整数となれば、解が素数になる可能性が出てくるのである。もし平方根が開けなければ、解は「実数」になってしまって題意を満たすことはできない。その判定にis_integer()という関数を利用したというわけである。

さて、格子点の中でうまい具合に平方数になるものは意外にたくさんあって、全部で12個見つかった。

% python3 discriminant.py < ab-data.dat|cat -n
     1  -5 1 25 True 5.0 0.0
     2  -4 1 16 True 4.0 0.0
     3  -3 1 9 True 3.0 0.0
     4  -2 1 4 True 2.0 0.0
     5  -1 1 1 True 1.0 0.0
     6  0 1 0 True 0.0 0.0
     7  -2 2 0 True 1.0 1.0
     8  -3 3 1 True 2.0 1.0
     9  -4 4 4 True 3.0 1.0
    10  -4 5 0 True 2.0 2.0
    11  -5 5 9 True 4.0 1.0
    12  -5 7 1 True 3.0 2.0
    13  Programme completed...

右端の2つの整数が2次方程式の解であるが、これが素数となっているのは12番目のケースのみである。0や1は素数の定義により除外されるからである。つまり、我々の探索範囲で条件を満たしたのは、$(a,b)=(-5,7)$のときに発生する、$f(2)=2, f(3)=3$という2つの素数の解であった。

素数の解は4つあるだろうか?

さて、(a,b)=(-5,7)の場合、ケース3.により2つの解が見つかったことになる。果たして、この(a,b)に対し、ケース1.およびケース2.そしてケース4.は解を2つ提供できるだろうか?まずは似たような状況としてケース4.が2つの素数解を提供できるか確認してみよう。

この場合、$x=-p<0$に対して, $k(p)=p^2+5p+7=-1$とならねばならない。つまり、 \begin{equation} p^2+5p+8 =0 \end{equation} の解が素数になっていなくてはならないのだが、解の公式で解いてみると \begin{equation} p = \frac{-5\pm\sqrt{25-32}}{2} \end{equation} となってしまい、これは「複素数」である。実数ですらないので素数には到底なり得ないのは明らかである。ということで、あえなく唯一の望みは消えてしまったのである。

ケース4.がだめなら、ケース1.とケース2.を満たす素数が2つ見つかればよい。やってみると \begin{equation} k(1)=1-5+7 =3, \quad k(-1) = 1+5+7 =13 \end{equation} となった。$k(1)=3$はよいのだが、2次方程式の解の一つとしてすでに見つかっているため、「だぶり」である。もうひとつのケースについては$k(-1)=13$となってダメである。これはケース2.なので$k(-1)<0$とならなければならないからである。13は素数であるから、かなりいい線まで肉薄できたのだが、正負の関係で失敗である。ということで、今回候補となった(a,b)については、問題6の前半で見つけたような素数の解が3つあるケースではなく、素数の解が2つある場合に過ぎないことが確認された。

東大数学2024問題6 (part 8): gnuplotで条件に合致する領域を塗ってみる

前回のあらすじ

しばらく前になるが、問題6の前半部分を解析的、および数値的に解いた。後半部分は、前半部分の応用であるからあまりおもしろくない、と思っていたのだが、gnuplotの練習の観点から少し面白そうなネタありそうだったので、取り上げてみることにする。

後半部分

後半部分は、考える関数が「若干」一般化され、整数$a,b$を用いて \begin{equation}g(x)=x^3+ax^2+bx=x(x^2+ax+b)\end{equation}となる。問題の前半は$a=10, b=20$という特別な場合を考えていたのであった。後半でも、$x\rightarrow n\in \mathbb{Z}$としたときに、$g(n)$が素数になる場合について考えさせる問題である。$g(n)$が素数になるケースは、問題前半で検討したように、4つの場合が存在する。ただし、以下では因数分解したときの二次関数の部分を$k(x)=x^2+ax+b$と表すことにする。

  1. $x=1$, かつ$k(1)>0$が素数になる場合
  2. $x=-1$, かつ$k(-1)<$が負の符号がついた素数($-p$)になる場合
  3. $p$を素数とする。$x=p>0$に対して$k(p)=p^2+ap+b =1$となる場合
  4. $x=-p <0$に対して$k(-p)=p^2-ap+b=-1$となる場合。

この4つの場合の組み合わせを丁寧に追いかければ証明はできるはずである。とはいえ、いろいろな組み合わせが発生するので、取りこぼしがないように気をつけねばならない。

最初のケース

まずは3.のケースを考えたい。$k(p)=1$となるような素数$p$が見つかるかどうかである。これは$p$についての2次方程式を解くことになるので、一般解は解の公式を用いて \begin{equation} p = \frac{-a\pm\sqrt{a^2-4(b-1)}}{2} \end{equation} と「形式的」にかける。もちろん、このままでは複素数だったり、無理数を一般的に許すことになるので、制限を課さねばならない。

絶対にまずいのが、複素数になってしまう場合である。これを避けるために判別式(つまり平方根の中身)が正値あるいは零値をもつことを課さねばならない。すなわち最初の条件式は \begin{equation} a^2-4(b-1) \ge 0 \rightarrow b \le \frac{a^2}{4}+1 \end{equation} である。この条件式を満たす$(a,b)$ならば, 素数$p$が複素数になってしまうことはない(笑)。

次は無理数にならない条件などを考えたいところだが、ちょっと面倒なので、まずは素数が正数であるという定義から攻めることにしよう。ここで考えるのは、2次方程式の2つの解が両方とも正値をとるための条件である。片方だけ正値ということはあるかもしれないが、まずは両方ともに正値になる場合について考えてみる。 \begin{equation} -a\pm\sqrt{a^2-4(b-1)} > 0 \rightarrow -a > \sqrt{a^2-4(b-1)} > a \end{equation} $a>0$の場合は成立しないので、$a=-a'<0$の場合だけを考える。すなわち \begin{equation} a' > \sqrt{a'^2-4(b-1)} > -a' \end{equation} 右側の不等式は自明な関係式なので、条件式としては利用できない。左側の不等式を自乗して整理すると \begin{equation} b > 1 \end{equation} という条件式が手に入る。

gnuplotの出番

ここで、gnuplotの出番である。3.のケースが2つの素数を許す「必要条件」として次の3つが浮かび上がった。 \begin{align} a &<0 \\ b & > 1\\ b &\le \frac{a^2}{4}+1 \end{align} 我々が興味あるのは、この領域を図示するにはどうすればよいかというgnuplotの問題である。まずはそれぞれの領域を塗りつぶし、重ね塗りすることで対応する領域を見つけ出したい。

いろいろ検索すると、弘前大学の相対論グループが発表しているgnuplotマニュアルが役立ちそうである。

まずは上の条件式のなかに含まれる放物線$b=a^2/4+1$を描いてみよう。今回は、座標軸も描き入れてみよう。

unset key
unset arrow

set size square
set grid

f(x)=x**2/4.+1

xmax=5; xmin=-5
ymax=8; ymin=-2
set arrow from xmin,0 to xmax,0 lw 2
set arrow from 0,ymin to 0,ymax lw 2

plot [xmin:xmax][ymin:ymax] f(x) lw 2

この放物線の下側の領域が3.の場合の2つの解が実数となる領域に対応する。塗りつぶすには次のようにする。

plot [xmin:xmax][ymin:ymax] f(x) w filledcurves y=ymin fc "#114422"

この「塗り絵」はちょっと重苦しいので、透明度を上げて印象を軽くしよう。

set style fill transparent solid 0.2
...
plot [xmin:xmax][ymin:ymax] f(x) w filledcurves y=ymin fc "#114422", f(x) w l lw 3

を付け足すとよい。 境界線の放物線グラフを際立たせるために「2度描き」していることに注意。

次は2次方程式の解が2つとも正になるための2つ目の条件$b>1$を取り込もう。$1<b<\text{ymax}$の部分を塗りつぶせば良いので、

plot [xmin:xmax][ymin:ymax] \
     f(x) w filledc y=ymin fc "#114422",\
     1 w filledc y=ymax fc "#114422",\
     f(x) lw 2 lt 2, \
     1 lw 2 lt 2

となる。重ね塗りされた部分が正しい領域になる。

最後に$a<0$を取り込む。描画領域を負の領域と正の領域に分ける必要があり、このためにはサンプリングのテクニックを使う。サンプリングはplot sampleで実現できる。具体的には次のようにする。

plot [xmin:xmax][ymin:ymax] sample \
     [xmin:xmax] f(x) w filledc y=ymin fc "#114422",\
     [xmin:xmax] 1 w filledc y=ymax fc "#114422",\
     [xmin:00] ymax w filledc y=ymin fc "#114422",\
     [xmin:xmax] f(x) lw 2 lt 2, \
     [xmin:xmax] 1 lw 2 lt 2

4行目の定義域が「部分的」、つまり負の領域だけが指定されている。つまりsampleによって定義域を分割(sampling)することができるのである。表示結果は次のような感じになる。流石に3重塗りだとよくわからなくなってくる。

2つの関数グラフに挟み込間れた部分だけを塗る

gnuplotは「挟み込んだ領域」を塗りつぶす命令をもっているのでそれを利用する。上の場合は、挟み込む領域が放物線と$y=ymin$だったので、放物線と$y=1$に変えればよいだけである。

plot [xmin:xmax][ymin:ymax]  f(x) w filledc y=1 fc "#114422",\
     f(x) lw 2, 1 lw 2 lt 2

最後にsampleのテクニックを盛り込んで完了である。随分簡単に描ける。

plot [xmin:xmax][ymin:ymax] sample \
     [xmin:0] f(x) w filledc y=1 fc "#114422", f(x) lw 2, 1 lw 2 lt 2

ちょっとだけ数学の問題に戻ってみる

さて、問題を解くという観点に少しだけ戻ってみよう。

塗りつぶされた領域のどこでも良いわけではなく、整数の格子点だけが許されることに注意しなくてはならない。このグラフに含まれている部分だけを抜き出せば、(a,b)=(-5,0),(-4, 0), (-3,0), (-2, 0), (-1,0), (0,0); (-5,1),(-4,1),(-3,1),(-2,1); (-5, 2), (-4,2),(-3,2),(-2,2); (-5,3),(-4,3), (-3,3); (-5,4),(-4,4), (-5,5),(-5,4),(-5,6),(-5,7)が候補となる。これらの値を$k(n)$に代入したときに、素数となるようケースが見つかるだろうか?(つづく)

画像処理にチャレンジ(part 4): jpeg(JFIF)のヘッダー情報

前回はjpegのヘッダー情報を読み込むための関数を書いた。そして、JPEGを表す「魔法数」、つまりSOIを表す整数である「0xff 0xd8」を読み取ることに成功した。今回は、この関数を用いて更なるヘッダー情報を読み込んでみる。

次のヘッダーを読み込む(APP0)

最初の2バイトが魔法数になっているのは、JPEGBMPも同じだった。たぶん、pngやmp3なんかも同様の魔法数があるのであろう(ちょっと検索してみたらどうも2バイトとは限らないようである)。JPEGの場合、魔法数の次から始まるのが種々の「セグメント(情報の塊)」で、それぞれのセグメントは「マーカー」という2バイトの整数値によって区切られている。したがって、最初の読み込みの後は再度2バイトを読み込んで、最初のマーカーを識別することになる。

多くのJPEG形式の場合、この最初のマーカーはFF Enという値になってようだが、これはAPPnと呼ばれるマーカーである。JPEG画像を生成したアプリケーションが付加したヘッダー情報が続くことを示すマーカーである。nは0-15つまり16進数の一桁分に相当する整数である。したがってAPPnセグメントは16種類ある。

もっとも頻発するのがn=0の場合、つまりAPP0である。マーカー値はFF E0になっていて、その値に続いてJPEGの拡張版であるJFIF形式の画像データが始まる。

APP1セグメントはある場合とない場合がある

ときどきAPP1が付随しているJPEG画像がネットなどで見つかるかもしれない。APP1とはExifという形式で記録された画像のヘッダー情報である。たとえば、どんなデジカメを使って撮影したか、どこで撮影したか、などといった情報が記録される。よく犯罪で利用されるのはこの部分である。個人情報などがここに記録される場合があるので、インターネットに画像をアップロードして不特定多数に公開するときは、この部分を削ってからアップしたほうがセキュリティが高まる。「削る」といっても、バイナリファイルをエディタで編集するのは難しいので、自作のプログラムを使って削除するのが手っ取り早いだろう(ここでの研究が役立つはずである)。

ファイルを開けておけば、次々と読み続けることができる。

JPEGのセグメントはマーカー値の次に、2バイトのセグメント長が来る。この情報はどれだけファイルを読み進めるか判断するときに役立つ。C言語のファイルポインタは、現在の読み取り位置がファイルのどの位置まできているか記録してくれるので、ファイルを開けっぱなしにしておけば(つまりfclose()を実行しなければ)、順番にデータを読み続けていくことが可能である。

ohtani.jpgの場合はちょっと古いver.1.01のJFIFであった

ohtani.jpgの場合は、JPEGの魔法数ff d8が来た後、ff e0というマーカーが置かれていた。つまり、APP0であり、JFIF形式であることが宣言されていた。JFIFのバージョンを調べてみると1.01であった。

コード

それではAPP0までを丸ごと読み取って、そのバージョンを表示させるコードを書いてみる。

#include <stdio.h>

int *readJpegHeader(FILE *fp, int byteSize, int *header_p){
    int i;
    for(i = 0; i < byteSize; ++i){
        header_p[i] = fgetc(fp);
    }
    return header_p;
}

int main(){
    int i, byteSize;
    int header_info[256];
    int *header_p;
    FILE *fp;
    if( (fp = fopen("ohtani.jpg","rb")) == NULL){
        printf("File not found.\n");
        return 1;
    }
    printf("File found.\n");
    header_p = readJpegHeader(fp, byteSize = 2, header_info);
    if( header_p[0] == 0xff && header_p[1] == 0xd8){
        printf("JPEG file found: success.\n");
        header_p = readJpegHeader(fp, byteSize = 2, header_info);
    }else{
        printf("Not a JPEG: error.\n");
        return 1;
    }

    if( header_p[0] == 0xff && header_p[1] == 0xe0){
        printf("Marker APP0 found.\n");
        header_p = readJpegHeader(fp, byteSize = 2, header_info);
        byteSize = header_p[1] * 256 + header_p[0];
    }else{
        printf("Marker APP0 not found. Maybe elsewhere\n");
        return 1;
    }

    header_p = readJpegHeader(fp, byteSize, header_info);
    for(i = 0; i < 5; ++i)
        printf("%c ",header_p[i]);
    printf("\n");
    printf("Version: %3.2f\n",header_p[5] + header_p[6]/100.);

    printf("-- next marker: %x %x\n",header_p[byteSize - 2], header_p[byteSize -1]);

    return 0;
}

マーカーの次に置かれているセグメント長は、マーカーを除いた部分のバイト長を意味する。とすると、最後から7行目にある、自作関数を最後に呼び出したところで、読み出しサイズを(バイト長を除いた)byteSize - 2ではなく、単にbyteSizeとしているのを不自然に感じるかもしれない。これはあえてそうすること次のマーカーを読み取っているのである。最後のprintf文でそれを確認しており、マーカーがff dbつまり、APP0の次にDQTセグメントが続くことが確認できる。

このプログラムをコンパイルしてから実行すると、次のようになるはずである。

% a.out    

File found.
JPEG file found: success.
Marker APP0 found.
Size of APP0: 16
J F I F  
Version: 1.01
-- next marker: ff db

さらにその先は?

このような感じで、次々とマーカーを識別しつつ、対応するセグメント長を計算しながらヘッダー情報を先へ先へと読み進めていけば良い。最後に、マーカーがff daとなったところが最後のセグメントである。大抵は12バイト長であるが、それを読み切った先が、JPEG画像データの本体となるのである。

ohtani.jpgの場合、画像データの本体の最初の値はf9 a1である。これはマーカーではない。画像データそのものである。

画像処理にチャレンジ(part 3): jpegのヘッダー情報を読み込む関数をつくる

前回のあらすじ

前回はJPEGのヘッダー情報を読み取ってみた。最初の「魔法数」に相当するSOIがff d8であること、つづくff e0がAPP0に相当し、JFIF形式の画像(JPEGの拡張版)となっていることなどが読み込めた。

しかし、1バイトずつfgetc()関数で読み取るのは、なんだか素人っぽくて恥ずかしい感じもしないでもない。そこでC言語らしく「自作の関数」でも書いて「どうだ!」とやってみたい。ということで、今回はそれをつかってJPEG(JFIF形式)のヘッダー情報を読み取ってみたい。

自作の関数Ver.0

今回も読み取るのは、LA Dodgersの大谷選手の画像である(MLBのHPより拝借したもの)。この画像は"ohtani.jpg"という名前で保存した。まずは基本部分のコードを書き込もう(main.cと名付ける)。

#include <stdio.h>

int main(){
    FILE *fp;
    if ( (fp = fopen("ohtani.jpg","rb")) == NULL) {
        printf("File not found.\n");
        return 1;
    }
    printf("File found.\n");
    fclose(fp);
    return 0;
}

この状態ですでにコンパイルすることができるが、当然ながら実行しても何も起きない(ohtani.jpgがあるかないかは判断してくれる)。

これまでは、main()のなかで直接fgetc()を利用していたのだが、関数を呼び出してそこで読み出すようにしてみる。関数名はreadJpegHeader()としよう。 読み出したヘッダー情報は「整数」なので整数型の関数となる。また、ファイルの中身にアクセスするためには、ファイルポインタ*fpを引き取る必要もあるだろう。 すなわち、次のような形になることが予想される。

int readJpegHeader(FILE *fp){
    int headerInfo;
    headerInfo = fgetc(fp);
    return headerInfo;
}

この部分をmain.cに書き込んでもいいし、readJpegHeader.cという別ファイルに書き込んで後でまとめてコンパイルしてもよいし、いろいろあるだろう。まだ関数が短いので、まずはmain.cの先頭に書き加えてみる。

#include <stdio.h>

int readJpegHeader(FILE *fp){
    int headerInfo;
    headerInfo = fgetc(fp);
    return headerInfo;
}

int main(){
    FILE *fp;
    if ( (fp = fopen("ohtani.jpg","rb")) == NULL) {
        printf("File not found.\n");
        return 1;
    }
    printf("File found.\n");
    header_info = readJpegHeader(fp);
    printf("%x\n",header_info);
    fclose(fp);
    return 0;
}

コンパイルは通る。実行すると、ohtani.jpgの最初の1バイトの情報が出力される。それは当然ながら0xffである。

jpegか、それともbmpかどうかを判定するには、最初の2バイトを読み出す必要がある。そこで任意の長さのデータの読み出しを可能にするために、関数の引数にbyteSizeという整数値を加えることにする。そうすると自作関数は次のように書き変わる。が、すぐに問題にぶち当たる。まずはコードを書いてみよう。

int readJpegHeader(FILE *fp, int byteSize){
    int i, headerInfo;
    for( i=0; i < byteSize; ++i){
        headerInfo = fgetc(fp);
    }
    return headerInfo;
}

このままコンパイルしてもエラーは出ないが、出力は2バイト目の0xd8だけとなる。もちろん、先ほどの実行結果から最初の2バイトが0xff 0xd8となったことは推測できる。しかし、1回の実行で2バイト分の情報を一括して出力させたい。1文字ずつ関数を呼ぶという手も考えられるが、それではfgetcを1文字ずつ呼び出すのと大差ないから、全くの無駄な解決法である。コンピュータというのは「繰り返し」をやらせるための機械であるから、これができないと手作業と大差ないことになってしまう。こういう場合に利用すべきものが配列である。そこで、headerInfoを整数型の配列に拡張してみる。

int *readJpegHeader(FILE *fp, int byteSize){
    int i, headerInfo[byteSize];
    for( i=0; i < byteSize; ++i){
        headerInfo[i] = fgetc(fp);
    }
    return headerInfo;
}

さて、このコードはコンパイルできるだろうか? まず、配列自体の受け渡しがC言語ではできないので、ポインタだけを自作関数からmain()関数に渡さねばならない。そのためには関数の型はint *つまり整数のポインタを返す関数として定義しなくてはならない。またmain()でこのポインタを受け取るためには、main()でポインタを用意しておく必要がある。ということで、全体として次のような形になる。

#include <stdio.h>

int *readJpegHeader(FILE *fp, int byteSize){
  int i;
  int headerInfo[byteSize];
  for( i=0; i < byteSize; ++i){
    headerInfo[i] = fgetc(fp);
  }
  return headerInfo;
}

int main(){
  int i;
  int *header_info, byteSize;
  FILE *fp;
  if( (fp = fopen("ohtani.jpg","rb")) == NULL){
    printf("File not found.\n");
    return 1;
  }
  printf("File found.\n");
  header_info = readJpegHeader(fp, byteSize=2);
  for( i=0; i < byteSize; ++i){
    printf("%d %x\n",i, header_info[i]);
  }
  fclose(fp);
  return 0;
}

コンパイルするとwarningが出るが、一応実行形式a.outは生成された。しかし実行してみると予想外の結果となってしまった。

%a.out 

File found.
0 ff
1 1

最後の行が"1 1"となっている。ここはJPEGのSOIとなるように"1 d8"となってほしいところである。自作関数で読み取りに失敗が起きているのであろうか?関数内で配列の中身を表示させてみる。

int *readJpegHeader(FILE *fp, int byteSize){
  int i;
  int headerInfo[byteSize];
  for( i=0; i < byteSize; ++i){
    headerInfo[i] = fgetc(fp);
  }
  for (i = 0; i < byteSize; ++i){
    printf("%d %x\n",i, headerInfo[i]);
  }
  return headerInfo;
}

コンパイルするとwarningは出るが、先ほど同じようにa.outは生成されるので、実行してみた。

File found.
0 ff
1 d8
0 ff
1 1

2、3行目が関数内で配列の中身を印字させた結果であり、正しくJPEGのSOIが表示されている。しかし、その配列(のポインタ)をmain()で参照した途端に問題が派生している。問題と言っても最初のffは正しいのだが、2つめのd8が壊れて1に変わってしまったという「部分的な問題」である。

あきらかにWarningを無視しているのが良くない点である。仕方ないので、よく読んで対処することにする。まずはwarningの内容を見てみよう。

warning: address of stack memory associated with local variable 'headerInfo' returned [-Wreturn-stack-address]
  return headerInfo;
         ^~~~~~~~~~
1 warning generated.

局所変数に関するスタックメモリのアドレスが返り値になっているから警告ね」と言っているように思える。だめなんですか?としか言いようがないが、ダメなのである。

局所変数はスタックと呼ばれる「簡易メモリ領域」に一時保存されるが、その内容は保証されないらしい。最初の要素は運良く保存してくれたが、2つ目はダメでした、ということなんであろう。スタックを利用せず、局所変数の値を保存するにはどうすればよいのだろうか?

というか、局所変数の内容を利用しているから警告なのである。とすれば、mainでheaderInfoを配列として宣言し、そのポインタを渡してやったらどうだろう。これなら、スタックを経由しないのではないだろうか?

ということで、次のように書き換えてみた。

#include <stdio.h>

#define MAX_HEADER 14

int *readJpegHeader(FILE *fp, int byteSize, int *headerInfo){
  int i;
  for(i = 0; i < byteSize; ++i){
    headerInfo[i] = fgetc(fp);
  }
  for(i = 0; i < byteSize; ++i){
    printf("%d %x\n",i, headerInfo[i]);
  }
  return headerInfo;
}

int main(){
  int i;
  int header_info[MAX_HEADER], byteSize;
  int *header_p;
  FILE *fp;
  if( (fp = fopen("ohtani.jpg","rb")) == NULL){
    printf("File not found.\n");
    return 1;
  }
  printf("File found.\n");
  printf("%p \n",header_info);
  header_p = readJpegHeader(fp, byteSize=2, header_info);
  printf("%p \n",header_p);
  for( i=0; i < byteSize; ++i){
    printf("%d %x\n",i, header_p[i]);
  }
  fclose(fp);
  return 0;
}

配列header_info[]をmain()の中で定義し、その先頭要素をポインタによって自作関数とやりとりする形式である。これなら局所変数を呼び出さないのでwarningは出ないはずである。コンパイルしてみると、警告は全て消えa.outが生成された!成功である。実行してみると予想通りの結果となり成功である。

File found.
0x16f6db6b0 
0 ff
1 d8
0x16f6db6b0 
0 ff
1 d8

最初にmain()で定義された配列header_infoの先頭アドレスが表示され、それが0x16f6db6b0であることがわかる。自作関数内で正しいJPEGヘッダーデータが読み込まれていることが次の2行"0 ff, 1 d8"の部分でわかる。次の長い16進数はmain()の中で表示させたポインタheader_pの中身である。配列header_infoの先頭アドレスと同じになっていることが確認できた。そして最後の2行である。0 ff, 1 d8となり、JPEG(JFIF)のSOI識別コードであるff d8が表示できた!

ちなみに、main()の中でheader_pの代わりにheader_infoを表示させてみたらどうなるだろうか?つまり、最後の部分であるfclose(fp);の文の上の箇所でheader_pの代わりにheader_infoに変えたらどうなるか?

...(省略)

for( i=0; i < byteSize; ++i){
    printf("%d %x\n",i, header_info[i]);
  }
fclose(fp);

...(省略)

コンパイルは通ってしまって、実行形式を走らせても同じ結果を得ることが確認できた。つまり、目的だった「自作関数からの複数データの受け渡し」に成功したことになる(やった!)。

画像処理にチャレンジ(part 2): jpegのヘッダー情報解読

前回のあらすじ

前回はBMP画像ファイルの読み込みに挑戦し、まずはヘッダーファイルという画像情報の記録に挑戦し、少しだけその「正体」を理解することができた。バイナリファイルの読み書きについての練習もできたし、実りの多い試みであった。

引き続いてBMPファイルの画像データそのものにアクセスし、加工してみようかと思ったのだが、ひとつアイデアが出てきてしまったので、最初にそちらを片付けてから本筋に戻ることにした。そのアイデアとはJPEGのヘッダー情報の読み取りの挑戦である!

もともと、デジカメのデータはjpegで記録することが多いので、画像処理したい対象はjpegpngであった。しかし「圧縮」という複雑な処理が施されているので、画像処理の素人には敷居の高い代物に見えた。しかし、よくよく考えるとヘッダー情報なら「圧縮」されているわけはないので、単なるバイナリファイルの読み取りで事済むはずじゃないだろうか?と思いたったというわけである。また、これまでに書き下したBMPヘッダーファイル読み取りプログラムは、少し書き換えればバイナリファイル一般に応用可能である。つまり、jpegに踏み込める上に、コードは再利用できるということで「一石二鳥」というわけである。

最初のマジックナンバー

jpegファイルはバイナリファイルなので、ファイルサイズ以内の読み取りサイズ(私のコード中でbyteSizeと定義されているもの)を指定してファイルをひらけば、なんらかのデータは読み取れるはずである。BMPの場合は最初14バイトがBITMAPFILEHEADERという「構造体」になっていたが、JPEGはどうなんだろうか?14バイト以上の画像データをどこかで拾ってきて、最初の14バイトを兎にも角にも読み取ってみよう。

拾ってきたのは、LA Dodgersの大谷選手の紹介HPに利用されていた画像である。今回のサンプルプログラムで読み取りたいのは、この画像のデータである。

MLBのLAドジャースのHPより拝借...
まずは何も考えずに最初の14バイトを(先日bmp用に書いたコードを使って)読み取ってみた。結果は次の通りである。

File found.
255 ff 
216 d8 
255 ff 
224 e0 
0 0 
16 10 
74 4a 
70 46 
73 49 
70 46 
0 0 
1 1 
1 1 
0 0 

読めた!読めたぞ!

最初の2バイトに"ff d8"というのがあって、これがBMPファイルのマジックナンバー"B M (0x42 0x4d)"に似た感じがある。続く"ff e0"もなんらかの記号のような印象を受ける。この数字を使って検索をかけてみると案の定であった。ただし、JPEGの場合には「魔法数」とは呼ばず、SOI(Start of Image)と呼ぶそうである。

JPEGのファイルであることを示す最初のマジックナンバー(魔法数)は、先頭2バイトの"ff d8"であった。つまり、最初の2バイトを読み取ると、BMPか、JPEGか識別できるというわけだ。これは後で画像処理アプリを開発するときに利用できそうな情報である。

JPEGには終わりのマークもある(EOI)

検索結果を読み進めると、BMPと違って、JPEGにはファイルの終わりを示すマークもあることがわかった。それはEOI(End of Image)と呼ばれる"0xff 0xd9"であるようだ。1バイトずつ読み取っていけば、最後にEOIが出てきて、そこでJPEGファイルが終わるということなんだろうか?この実験は後でやってみることにしよう。

ff e0の意味

JPEGの「魔法数」に相当するSOI、つまり0xff 0xd8の数字の次に記録されていたのが、0xff 0xe0であった。検索結果のなかにこの意味の解説も書かれていた。

これはAPPnと呼ばれるデータセグメントの先頭で、JPEGを作成したアプリケーションが付与したヘッダー情報になっているそうである。nの値は0から15まで(つまり16進数)の範囲をもっている。今回の大谷の画像で見つかったのはe0だからAPPnの最初のカテゴリーであるという意味であろう。

ff e0に続く次の2バイトはApp0セグメントの長さを(バイト単位で)示す情報で、大谷の画像では0x00 0x10つまり、16byteがApp0に割かれていることがわかる。次の5バイトは10 4a 46 49 46 00となっているが、これはASCIIコードになっているそうである。文字形式で印字してみると、JFIF\0と印字された。これは厳密な意味での”JPEG”画像データではなく、拡張されたJPEGという意味らしい。とはいえ、その拡張子はjpg, jpegを取ることが許されていて、緩い意味でのJPEG画像ファイルとして認識されているという。(これは知らなかったので勉強になった。)

JFIF\0の宣言文字に続くのは0と1の値であるが、これはなんらかの画像情報のデータであろう。

最初の試みとしてはこんなところであろう

ということで、jpegファイルもヘッダー部分は単なるバイナリデータであることが確認でき、私の未熟な技能でも読み取りが可能であることがわかった。手も足も出なかったJPEGファイルに一太刀浴びせることができたのは嬉しい限りである。JPEGのヘッダーファイルはBMPのそれよりも長く複雑であるそうだが、各種のマーク記号によって区切られており、セグメントごとに読み取ればなんとか処理できそうである。次はそちら方面をやってみようかと思っている。たとえば、JPEGのファイルサイズとか、撮影日時なんかのヘッダー情報が手に入れば、それはかなり実用面に近づくと思うのである。

おまけ(ファイルの終端)

JPEGの終端がff d9で終わっているか、まずは手っ取り早くxddコマンドを使って調べてみた。結果は次のとおりで、予想通りであった。