コラッツ予想 (Collatz conjecture) その三

前回(↓)の続き。
nazodane.hatenadiary.org

前回のコードを適当にベクトライズするとこうなる。

// gcc rcollatz.c -O3 -march=haswell
// clang rcollatz.c -O3 -mavx
#include <stdint.h>
#include <stdio.h>
//#include <assert.h>
void loop(uint64_t x, int l){
//  assert(x&1);
  uint64_t xa[64];
  for (int i = 0; i < 63; i++) xa[i] = x;
  for (int i = 0; i < 63; i++) xa[i] = xa[i] << (i+1); // shift + add vs. fma
  for (int i = 0; i < 63; i++) xa[i] = xa[i] + 0xaaaaaaaaaaaaaaab * ((1 << (i+1)) - 1);
  for (int i = 0; i < 63; i++) {
    if (0x5555555555555555 >= xa[i] /*暫定→*/ && xa[i] < 0xffff){
      printf("l=%d::%0lx\n", l, xa[i]);
      if (l) loop(0xaaaaaaaaaaaaaaab * (xa[i] ^ 1), l-1);
    }
  }
}

int main(void){
  uint64_t s = 1, x = 0;
  printf("l=ss::1\n");
  for(int i = 1; i < 16; i++){
    s <<= 2; s++;
    printf("l=s::%0lx\n", s);

    x <<= 2; x += (0xaaaaaaaaaaaaaaab << 2);

    loop(x, 99);
  }
  return 0;
}

以下のコードはシフトより上と下に分けることが出来る。

bool ge_shift_add(uint32_t x){
  return (CMPVAL >= (x<<SHIFT) + ADDVAL);
}

そのためxのシフトを無くしてこう変形できるはず。

bool ge_shift_add2(uint32_t x){
  return ((((CMPVAL>>SHIFT) - !!((CMPVAL & ((1<<SHIFT)-1)) < (ADDVAL & ((1<<SHIFT)-1))))&(0xffffffff>>SHIFT)) >= ((x + (ADDVAL>>SHIFT))&(0xffffffff>>SHIFT)));
}

ADDVAL>>SHIFT側をpythonで計算するとこんな感じ。

>>> print [(0xaaaaaaaaaaaaaaab * ((1 << (i+1)) - 1) &0xffffffffffffffff) >> (i+1) for i in range(63)]
[6148914691236517205L, 0L, 1537228672809129301L, 0L, 384307168202282325L, 0L, 96076792050570581L, 0L, 24019198012642645L, 0L, 6004799503160661L, 0L, 1501199875790165L, 0L, 375299968947541L, 0L, 93824992236885L, 0L, 23456248059221L, 0L, 5864062014805L, 0L, 1466015503701L, 0L, 366503875925L, 0L, 91625968981L, 0L, 22906492245L, 0L, 5726623061L, 0L, 1431655765L, 0L, 357913941L, 0L, 89478485L, 0L, 22369621L, 0L, 5592405L, 0L, 1398101L, 0L, 349525L, 0L, 87381L, 0L, 21845L, 0L, 5461L, 0L, 1365L, 0L, 341L, 0L, 85L, 0L, 21L, 0L, 5L, 0L, 1L]

偶数が0、奇数が[(4**i)/3 for i in range(1, 33)]の逆順になってる。

((((CMPVAL>>SHIFT) - !!((CMPVAL & ((1<<SHIFT)-1)) < (ADDVAL & ((1<<SHIFT)-1))))&(0xffffffff>>SHIFT)) 側をpythonで計算するとこんな感じ。

>>> print [(lambda shift=i+1, mask=(1<<(i+1))-1:((0x5555555555555555 >> shift) - ((0x5555555555555555&mask) < ((0xaaaaaaaaaaaaaaab*mask) & mask)))&(0xffffffffffffffff>>shift))() for i in range(63)]
[3074457345618258602L, 1537228672809129301L, 768614336404564650L, 384307168202282325L, 192153584101141162L, 96076792050570581L, 48038396025285290L, 24019198012642645L, 12009599006321322L, 6004799503160661L, 3002399751580330L, 1501199875790165L, 750599937895082L, 375299968947541L, 187649984473770L, 93824992236885L, 46912496118442L, 23456248059221L, 11728124029610L, 5864062014805L, 2932031007402L, 1466015503701L, 733007751850L, 366503875925L, 183251937962L, 91625968981L, 45812984490L, 22906492245L, 11453246122L, 5726623061L, 2863311530L, 1431655765L, 715827882L, 357913941L, 178956970L, 89478485L, 44739242L, 22369621L, 11184810L, 5592405L, 2796202L, 1398101L, 699050L, 349525L, 174762L, 87381L, 43690L, 21845L, 10922L, 5461L, 2730L, 1365L, 682L, 341L, 170L, 85L, 42L, 21L, 10L, 5L, 2L, 1L, 0L]

つまり、[(2**(i+2)-3-(-1)**i)/6 for i in range(63)] の逆順になってる。

うーん、このあともシフト計算無くして行ける気がするけど、式が長くなって頭こんがらがっていくので今日はここまで。


別式の追加

コラッツ予想 (Collatz conjecture) その二

前回(↓)の続き。
nazodane.hatenadiary.org

コラッツ問題の計算で、最後の手前の数字は「3の倍数かつ1を足すと2の累乗数になる数」になるはず(偶数処理を遅延した場合はそれに2の累乗数を掛ける必要あるけど、今回は遅延しないで考えることとする)。

1を足すと2の累乗数になる数というのは二進法で1が連続してる数のこと(e.g. 0b1、0b11、0b111、...)であり、
その中で3の倍数になるのは1の連続が偶数個の場合となっている (以下の検証はPython)。

>>> print [(2**x-1) for x in range(1, 32)]
[1, 3, 7, 15, 31, 63, 127, 255, 511, 1023, 2047, 4095, 8191, 16383, 32767, 65535, 131071, 262143, 524287, 1048575, 2097151, 4194303, 8388607, 16777215, 33554431, 67108863, 134217727, 268435455, 536870911, 1073741823, 2147483647]
>>> print [(2**x-1)%3 for x in range(1, 32)] # 0の時に割り切れる
[1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1]
>>> print [(2**(x*2)-1)%3 for x in range(1, 32)]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

3を掛ける前はどうなるかというと以下のようになる。

>>> print [(2**(x*2)-1)/3 for x in range(1, 32)]
[1, 5, 21, 85, 341, 1365, 5461, 21845, 87381, 349525, 1398101, 5592405, 22369621, 89478485, 357913941, 1431655765, 5726623061, 22906492245, 91625968981, 366503875925, 1466015503701, 5864062014805, 23456248059221, 93824992236885, 375299968947541, 1501199875790165, 6004799503160661, 24019198012642645, 96076792050570581, 384307168202282325, 1537228672809129301]
>>> print [bin((2**(x*2)-1)/3) for x in range(1, 32)]
['0b1', '0b101', '0b10101', '0b1010101', '0b101010101', '0b10101010101', '0b1010101010101', '0b101010101010101', '0b10101010101010101', '0b1010101010101010101', '0b101010101010101010101', '0b10101010101010101010101', '0b1010101010101010101010101', '0b101010101010101010101010101', '0b10101010101010101010101010101', '0b1010101010101010101010101010101', '0b101010101010101010101010101010101', '0b10101010101010101010101010101010101', '0b1010101010101010101010101010101010101', '0b101010101010101010101010101010101010101', '0b10101010101010101010101010101010101010101', '0b1010101010101010101010101010101010101010101', '0b101010101010101010101010101010101010101010101', '0b10101010101010101010101010101010101010101010101', '0b1010101010101010101010101010101010101010101010101', '0b101010101010101010101010101010101010101010101010101', '0b10101010101010101010101010101010101010101010101010101', '0b1010101010101010101010101010101010101010101010101010101', '0b101010101010101010101010101010101010101010101010101010101', '0b10101010101010101010101010101010101010101010101010101010101', '0b1010101010101010101010101010101010101010101010101010101010101']

ここから2**nを掛けて(右側に0を増やす)、1を引いて(LSBをクリアしてその下のビットを全て立てる)、3の倍数を作っては3で割るというのを繰り返していけば良いはず。
前二つの操作は合わせると0b1をクリアして右側に1を増やしていくということになる。

3で割り切れるかどうかは以下みたいに変形できるらしい(gcc -O3の結果より)。

#include <stdint.h>
uint64_t is_divable3(uint64_t x){
  return (x % 3 == 0);
}

uint64_t is_divable3_opt(uint64_t x){
  return 0x5555555555555555 >= 0xaaaaaaaaaaaaaaab * x;
}

上記を全て合わせると逆から計算するこんなコードが出来上がったけど上手く動いているかは未確認。

#include <stdint.h>
#include <stdio.h>
//#include <assert.h>
void loop(uint64_t x, int l){
//  assert(x&1);
  x = 0xaaaaaaaaaaaaaaab * (x ^ 1);
  for (int i = 1; i < 64; i ++){
    x <<= 1; x += 0xaaaaaaaaaaaaaaab;
    if (0x5555555555555555 >= x /*暫定→*/ && x < 0xffff){
      printf("l=%d::%0lx\n", l, x);
      if (l) loop(x, l-1);
    }
  }
}
int main(void){
  uint64_t s = 1;
  printf("l=ss::1\n");
  for(int i = 1; i < 16; i++){
    s <<= 2; s++;
    printf("l=s::%0lx\n", s);
    loop(s, 99);
  }
  return 0;
}


偶数のケースと1のケースを修正。


同じ値でループするケースを削除。


重複無くしたはず。

コラッツ予想 (Collatz conjecture)

1億円欲しいからコラッツ予想に手を出してみる。

science.srad.jp
ja.wikipedia.org

まず偶数処理をまとめたコラッツ問題の計算プログラムを適当に書く。

// gcc collatz.c -march=haswell
// ./a.out 13179928405231
// Copyright: Nazodane
// License: GPLv2 or later
#include <stdio.h>
#include <x86intrin.h>
int main(int argc, char *argv[]){
  u_int64_t n = 0;

  if(argc < 2){
    printf("./collatz num");
    return 0;
  }
  sscanf(argv[1], "%lu", &n);

  u_int64_t step = _tzcnt_u64(n);
  u_int64_t odd = n >> step; // 偶数 / 2 を繰り返す
//  printf("%lu\n", step);
  while(odd > 1){
    u_int64_t even = ((odd+1) + (odd<<1)); // 奇数 * 3 + 1
    u_int64_t cnt = _tzcnt_u64(even);
    odd = even >> cnt; // 偶数 / 2 を繰り返す
    step += 1 + cnt;
//  printf("%lu\n", step);
  }

  printf("step = %lu\n", step);
}

オーバーフローしやすくなるのを無視して偶数処理を遅延させる。

// gcc collatz.c -march=haswell
// ./a.out 13179928405231
// Copyright: Nazodane
// License: GPLv2 or later
#include <stdio.h>
#include <x86intrin.h>
int main(int argc, char *argv[]){
  u_int64_t n = 0;

  if(argc < 2){
    printf("./collatz num");
    return 0;
  }
  sscanf(argv[1], "%lu", &n);

  u_int64_t step = 0;
  u_int64_t x = n&-n;
  while(n != x){ // 奇数処理
    n = ((n+x) + (n<<1));
    x = n&-n;
    step++;
  }
  step += _tzcnt_u64(x); // 偶数処理
  printf("step = %lu\n", step);
}

適当に多倍長化してみる。

// gcc collatz.c -march=haswell
// ./a.out 13179928405231
// Copyright: Nazodane
// License: GPLv2 or later
#include <stdio.h>
#include <x86intrin.h>

#define N 10000
int main(int argc, char *argv[]){
  u_int64_t n[N] = {0};

  if(argc < 2){
    printf("./collatz num");
    return 0;
  }

  for(int a = 1; a < argc; a++)
    sscanf(argv[a], "%lx", &n[argc - a - 1]);

  u_int64_t step = 0;
  u_int64_t idx = 0;
  while(idx < N && !n[idx]) idx++; // TODO: overflow check
  u_int64_t max = N-1;
  while(max > 0 && !n[max]) max--;
//  printf("idx:max = %lu:%lu\n",idx,max);

  u_int64_t x = n[idx]&-n[idx];
  while(max != idx || n[idx] != x){
    unsigned __int128 t = n[idx];
    t = (t+x) + (t<<1);
    n[idx] = t;
    u_int64_t ovfl = t >> 64;
    if(idx == max && ovfl) max++; // TODO: overflow check

    for (int i = idx+1; i < max+1; i++){
      t = n[i];
      t = (t+ovfl) + (t<<1);
      n[i] = t;
      ovfl = t >> 64;
      if(i == max && ovfl) max++; // TODO: overflow check
    }

    while(idx < max && !n[idx]) idx++;
    x = n[idx]&-n[idx];
//    printf("idx:max = %lu:%lu\n",idx,max);
    step++;
  }
  step += _tzcnt_u64(x) + idx*64; // 偶数処理
  printf("step = %lu\n", step);
}

下にある例を試してみる。
math.stackexchange.com


The best 500-digit number I found used 25858 steps, and is:
16910⋅2^1646−1=

52907,40618,32017,54152,63867,78240,01344,94071,10672,87955,14970,95443,32426,37432,94723,15343,69589,15667,94197,77289,84251,30595,66986,13736,01999,30240,57853,88392,15872,66353,93648,16021,76528,02318,47325,39772,11544,04475,63943,10264,65156,58431,35122,30262,85026,21429,61128,92928,64147,40168,11368,14718,14844,43105,54823,58463,32284,16963,01649,12670,32815,41060,21677,81365,13557,54204,91127,46553,29004,21933,21888,08910,76179,94524,51718,32033,27205,55270,66366,57606,00582,53081,68267,14042,35673,44288,21454,35891,44482,68697,66804,20865,40074,26362,70517,10800,03191,26585,23809,38239

16進数に変換して16文字で区切ってコマンドラインに入れる。

$ ./a.out 10837FFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF
step = 25858

上手く動いてるみたいだけど、さてここからどうしよう?

真面目に闘病中

今は別の界隈に居たりしますが、長年の病気の原因がやっとこさ分かって治る見込みが出来てきたのでここにも書いておこうかなと思いました。

 

ついでに昔の日記を読み返してみたけど当時のこと思い出して何とも言えない気分になってきてアレ(^^;

文体が安定しないのは今もな気がする…

 

治っても治らなくてもまた報告します。

[その他] Anitubeの広告配信元

怪しい記事を見かけたので、こちらの環境で解析したAnitubeの広告配信経路をここに書いておきます。なおこれは一例であって、環境や機会により広告関連スクリプトの流れも変わるはずです。


とりあえず、以下のページのソースコードを開きます。
http://www.anitube.se/video/#####/XYZ

すると、ページのヘッダ、サイド、フッタ、動画の下の広告表示部分に以下のようなHTMLコードが埋めこまれています。これが、広告表示スクリプトを読み込むための起点となります。ここでは蘭Adglare社のスクリプトを読み込んでいます。

		<!-- Advert Start -->
			<div id="pageAdvert" align="center">
			<script async src='//anigrupo.engine.adglare.net/?#########'></script><span id=#########></span>
			</div>
		<!-- Advert End -->

その後、以下のURLのJavascriptが読み込まれ実行されます。(なお、URLの後ろの#########はzIDとのこと。)

//anigrupo.engine.adglare.net/?#########

スクリプトの内容は大幅に省略しますが、以下のような感じで、別のJavascriptを読み込むためのものです。

(前略)
AdGlare.loadJS('//anigrupo.engine.adglare.net/?#########&t=1&tt=##########-########');

loadJS関数でURLが色々加工された後、最終的に以下のURLにあるJavascriptが読み込まれ実行されます。

http://anigrupo.engine.adglare.net/?#########&t=1&tt=##########-########&winID=#####&screen=####x####&framed=0&vb=1&crIDsLoaded=&referer=http%3A%2F%2Fwww.anitube.se%2Fvideo%2F#####%2FXYZ

上記URLのJavascriptの内容はこんな感じ。

(前略)
var AdGlareSettings_###### = { zID: #########, width: ###, height: ###, impURL: "//anigrupo.engine.adglare.net/imp?data=(略)", isThirdParty: true, autoClickTracking: true, isSync: false, alwaysOnScreen: false, alwaysOnScreenData: '', animationType: 'none', allowClose: true, closeButtonURL: '//anigrupo.cdn.adglare.net/inventory/close_button.png', isBase64Encoded: true, zoneHTML: "(BASE64っぽい文字列)", (略)}

(略)document.write('<span id=zone'+this.settings.zID+'>'+this.settings.zoneHTML+'</span>'(略)

BASE64っぽい文字列をデコードすると、HTMLの断片が現れます。

<!--  ad tags Size: 728x90 ZoneId:#######-->
<script type="text/javascript" src="http://js.sprout-ad.com/t/###/###/x#######.js"></script>

そして、そこに書かれたスクリプトが読み込まれます。このURL (*.sprout-ad.com) は日本のトラストリッジ社のものですね。

内容は下のような感じ。

(略)
var di={(略),vd:"sprout-ad.genieesspv.jp",(略)};
(略)
u+="//"+di.vd+"/yie/ld/jsk";
u+=(略)

そして、ここで以下のURLが生成されます。このURL (*.genieesspv.jp) は日本のジーニー社のものですね。

http://sprout-ad.genieesspv.jp/yie/ld/jsk?zoneid=#######&cb=###########&charset=UTF-8&loc=http%3A%2F%2Fwww.anitube.se%2Fvideo%2F#####%2FXYZ&referer=http%3A%2F%2Fwww.anitube.se%2Fvideo%2F#####%2FXYZ&sw=####&sh=####&topframe=0

生成されたURLからスクリプトが更に読み込まれます。内容は下のような感じ。

gen_tag = "%3cscript%20type%3d%27text%2fjavascript%27(略)%2fscript%3e";document.write(decodeURIComponent(gen_tag));

gen_tagのエスケープを解いてやると、ついに広告本体のHTMLが出てきます (謎なiframeも付属してましたが解説では省略)。今回は最終的な広告のURL (*.gsspat.jp) がジーニー社のものでした。

<script type='text/javascript'>()<a href=\"http:\/\/rt.gsspat.jp\/c?c=(略)&amp;y=1&amp;p=XXXXXXXXXXX&amp;do=(略)&amp;vs=(略)\" style=\"border:none;margin:0;padding:0;\" target=\"_top\" >\\n');(略)document.write(\"<img src=\\\"\" + scheme + \"\/\/rt.gsspat.jp\/b?p=XXXXXXXXXXX&amp;y=1&amp;v=(略)\\\" height='1' width='1' style='display: none;'>\");(略)</script>

ということで、今回調べた時は蘭Adglare社→日トラストリッジ社→日ジーニー社という流れでしたが、別の会社の広告が出てくることも確認していますし、経路はまだまだあるはずです。
それと、これだけでは契約の流れがどうなってるのかまでは分かりませんね。更に中間業者が挟まっている気がします。

# …にしても凄く長いタライ回しだなぁ…。

[カラーサイエンス][Python] 正しく変換したアニメカラーデータ

ねこまたやさんのアニメカラー測定データがどうにも怪しい (実測データと目視データが異なる) ので調べてみた。
同サイトのD65 XYZ値は100倍になっているようで、まずここでハマる。測定機器の凸版 CS CM1000は調べても詳細が出てこないが、下記を見る限りD50光源らしい?


このツールでの測色はD50にしないと正常なプロファイルが取れません。
光源の切替えはCS-CM1000付属のソフト上で行います。

ピンと来たので、測定データのD65 XYZ値をXYZ Scaling変換でD50 XYZ値に一旦戻してから、より高精度なBradford変換でD65 XYZ値に変換してみたら、どうにも其れっぽい色が出てきた。値のズレはXYZ Scaling変換にも一因があるようだ。

import numpy as np
def fix_color(xyz):
	xyz_mat = np.matrix(xyz).T

	# http://technorgb.blogspot.jp/2015/08/blog-post_22.html
	# http://www.brucelindbloom.com/index.html?Eqn_ChromAdapt.html
	xyz_to_lms_mat = np.matrix([
		[ 0.8951,  0.2664, -0.1614],
		[-0.7502,  1.7135,  0.0367],
		[ 0.0389, -0.0685,  1.0296]
	])

	d65_mat = np.matrix([0.95047, 1.0, 1.08883]).T
	d50_mat = np.matrix([0.96422, 1.0, 0.82521]).T

	d65_xyz_to_d50_xyz_mat = np.diagflat(d50_mat.A / d65_mat.A) # XYZ Scaling method
	d50_xyz_to_d65_xyz_mat = xyz_to_lms_mat.I * np.diagflat((xyz_to_lms_mat * d65_mat) / (xyz_to_lms_mat * d50_mat)) * xyz_to_lms_mat # Bradford method
	xyz_mat = d50_xyz_to_d65_xyz_mat * (d65_xyz_to_d50_xyz_mat * xyz_mat)	

	return xyz_mat.A1

sRGBへの変換 (Numpy使いきれてないのがバレる)

import numpy as np

# https://www.w3.org/Graphics/Color/srgb
def linearsrgb_to_srgb(lin):
#	return [0 if c < 0 else 12.92 * c if c <= 0.0031308 else 1.055 * np.power(c, (1/2.4)) - 0.055 for c in lin]
	return [12.92 * c if abs(c) <= 0.0031308 else math.copysign(1.055 * np.power(abs(c), (1/2.4)) - 0.055, c) for c in lin]

def d65xyz_to_srgb(xyz):
	d65xyz_to_linearsrgb_mat = np.matrix([
		[ 3.2406255, -1.537208 , -0.4986286],
		[-0.9689307,  1.8757561,  0.0415175],
		[ 0.0557101, -0.2040211,  1.0569959]
	])

	xyz_mat = np.matrix(xyz).T

	lin_srgb = d65xyz_to_linearsrgb_mat * xyz_mat

	srgb = linearsrgb_to_srgb(lin_srgb)
	return [i*255 for i in srgb]

書き出し

import numpy as np
import csv

with open("test.html", "w") as htmlfile:
	htmlfile.write("<html><head></head><body>")
	# http://www.nekomataya.info/teck_info/taiyo_color/TaiyoChart.txt
	with open("TaiyoChart.txt", "rb") as cvsfile:
		cl = 0
		chartreader = csv.reader(cvsfile, delimiter='\t')
		for row in chartreader:
			if cl < 2:
				cl = cl + 1
				continue
			if len(row) < 9:
				continue
			cname = row[2]
			xyz   = map(float,[row[7], row[8], row[9]])
			xyz_str = ",".join(map(str,xyz))
			srgb  = d65xyz_to_srgb(fix_color(np.array(xyz)/100))
			srgb_str = ",".join([str(int(round(c))) for c in srgb])
			o_srgb = [row[12], row[13], row[14]]
			o_srgb_str = ",".join(o_srgb)
			m_srgb = [row[4], row[5], row[6]]
			m_srgb_str = ",".join(m_srgb)
			print("%s = XYZ(%s), sRGB(%s), o_sRGB(%s)"%(cname, xyz_str, srgb_str, o_srgb_str))
			htmlfile.write("<div style='display: inline-block; height 5em; width: 10em;'>\
							<p style='margin: 0'>%s</p>\
							<p style='background-color: rgb(%s); margin: 0'>XYZ=%s<br>sRGB=%s</p>\
							<p style='background-color: rgb(%s); margin: 0'>o_sRGB=%s</p>\
							<p style='background-color: rgb(%s); margin: 0'>m_sRPG=%s</p></div>"%\
							(cname, srgb_str, xyz_str, srgb_str, o_srgb_str, o_srgb_str, m_srgb_str, m_srgb_str))
	htmlfile.write("</body></html>")


若干表現変更。

暗い色はsRGB変換よりもgamma 2.2変換の方が、アニメカラー測定データにある目視カラーデータに近くなるように見える (X7より後はsRGB変換の方が近い)。うーむ、目視カラーデータの方も怪しい気がするような。

AbemaTVの仕様とHLSの暗号化の弱さ

AbemaTVの仕様について気になったので調べてみた (研究目的です念の為)。
AbemaTVはPCへの動画配信において、配信プロトコルにHLSを使用しているようだ。HLSはMPEG-DASHと異なりDRMが使えず (厳密にはMac環境のFairplayなどの例外もあるが) 、AbemaTVでは鍵の生成に若干工夫を行ってるのみのようだ。


まず、APIを使ってチャンネル一覧をダウンロード。

$ curl https://api.abema.io/v1/channels
{"channels":[{"id":"abema-news","name":"AbemaNewsチャンネル","playback":{"hls":"https://linear-abematv.akamaized.net/channel/abema-news/playlist.m3u8"}},{"id":"abema-special","name":"AbemaSPECIALチャンネル","playback":{"hls":"https://linear-abematv.akamaized.net/channel/abema-special/playlist.m3u8"}},(後略)

次に画質一覧をダウンロード。

$ curl https://linear-abematv.akamaized.net/channel/abema-news/playlist.m3u8
#EXTM3U
#EXT-X-STREAM-INF:PROGRAM-ID=1,BANDWIDTH=184000
180/playlist.m3u8
#EXT-X-STREAM-INF:PROGRAM-ID=1,BANDWIDTH=300000
240/playlist.m3u8
#EXT-X-STREAM-INF:PROGRAM-ID=1,BANDWIDTH=900000
360/playlist.m3u8
#EXT-X-STREAM-INF:PROGRAM-ID=1,BANDWIDTH=1400000
480/playlist.m3u8
#EXT-X-STREAM-INF:PROGRAM-ID=1,BANDWIDTH=2200000
720/playlist.m3u8
#EXT-X-STREAM-INF:PROGRAM-ID=1,BANDWIDTH=4200000
1080/playlist.m3u8

映像のフラグメント一覧をダウンロード。

$ curl https://linear-abematv.akamaized.net/channel/abema-news/1080/playlist.m3u8
#EXTM3U
#EXT-X-VERSION:3
#EXT-X-TARGETDURATION:6
#EXT-X-MEDIA-SEQUENCE:951004
#EXT-X-DISCONTINUITY-SEQUENCE:16994
#EXT-X-KEY:METHOD=AES-128,URI="abematv-license://XXXXXXXXXXXXXXXXXXXXXX",IV=0x000000000000000000000000000000
#EXTINF:5.005000,
https://linear-abematv.akamaized.net/tsnews/abema-news/h264/1080/xxxxxxxxxxxxxxxxxxxxxx.ts
#EXTINF:5.005000,
https://linear-abematv.akamaized.net/tsnews/abema-news/h264/1080/xxxxxxxxxxxxxxxxxxxxxx.ts
#EXTINF:5.005000,
https://linear-abematv.akamaized.net/tsnews/abema-news/h264/1080/xxxxxxxxxxxxxxxxxxxxxx.ts
#EXTINF:5.005000,
https://linear-abematv.akamaized.net/tsnews/abema-news/h264/1080/xxxxxxxxxxxxxxxxxxxxxx.ts

さて、映像はAES-128方式で暗号化されているようだ。暗号の鍵には初期化ベクトル(IV)とURIが指定されているが、URIに使われているabematv-licenseスキーマとは何だろう。仕組みは良く分からないが、Chromeの通信ログを見ると、スキーマの後ろの部分 (XXXXXXXXXXXXXXXXXXXXXX) と何処かにあるトークンを使って、とあるURLにアクセスしているようだ。

トークンはローカルストレージにあるものと同じようなので、Chromeのコンソールからwindow.localStorage["abm_mediaToken"]と打つと手に入る。このトークンとスキーマの後ろの部分を使って、ライセンスキーの種を手に入れる。

$ curl https://license.abema.io/abematv-hls?t=トークン --data '{"lt":"スキーマの後ろの部分","kv":"wd","kg":166}'
{"cid":"abema-news","k":"XXXXXXXXXXXXXXXXXXXXXXX"}

さて、どうやってライセンスキーの種 (k) をキーに変換するのだろう? 調べた所、遅延ロードされた以下のJavascriptがこの変換を処理しているようだ。
https://abema.tv/xhrp.js

若干難読化されているけれども、肝心の部分はそのままだし、コードインジェクションもし放題なので割と何とかなる。
キー計算の表面部分のロジックはこんな感じ。

function _0x569113(cid, uid, k){
  var _k = k.substring(0,k.length-1);
  var c = k.charAt(k.length-1);
  return c=='5'?_0x1e2ccc(cid, uid, _k):
         c=='4'?_0xa25b8f(cid, uid, _k):
                _0x2782e2(cid, uid, _k);
}
var _0x5ee3af=_0x569113(cid, window.localStorage["abm_userId"], k)

キー計算の中心部分は解読していないけれど、alert(_0x5ee3af);をインジェクションしてコードを実行するだけでキーが手に入る。
手に入ったキーは、バイナリ化してkey.binとして保存しておく。あとはそのキーを使って再生するだけ。

$ wget -N https://linear-abematv.akamaized.net/channel/abema-news/1080/playlist.m3u8 \
&& sed -i 's/URI=.*\,/URI=\"key.bin\",/g' playlist.m3u8 \
&& ffplay playlist.m3u8 -protocol_whitelist file,http,https,tcp,tls,crypto

フラグメント毎にしか再生できないので実用性には欠けるけれども、何にせよHLSが弱いことは証明できたので良いかな。
今後、AbemaTVでも強固なDRM付きのMPEG-DASHが導入されていくらしいので期待。



追記。何となく「簡単さ」が伝わってないようで残念なので、Python + Selenium WebDriverで自動化した方法を書いておきます。といっても大したものでもないですが。

from selenium import webdriver
from time import sleep
import requests
import re
import os

if __name__ == '__main__':
    browser = webdriver.Chrome(executable_path = "/usr/lib/chromium-browser/chromedriver")
    browser.get("https://abema.tv/now-on-air/abema-news")
    sleep(1)
    js = requests.get("https://abema.tv/xhrp.js").text
    mod_js = re.sub("(_0x31a687=.*?);", "\\1;window.key=_0x31a687;", js)
    browser.execute_script(mod_js)
    sleep(1)
    key = browser.execute_script("return window.key;")
    print(key)
    browser.close()

    f = open("key.bin", "wb")
    f.write(bytearray(key))
    f.close()

    pl = requests.get("https://linear-abematv.akamaized.net/channel/abema-news/1080/playlist.m3u8").text
    mod_pl = re.sub('URI=.*?\,', 'URI=\"key.bin\",', pl)

    f = open("playlist.m3u8", "w")
    f.write(mod_pl)
    f.close()

    os.system("ffplay playlist.m3u8 -protocol_whitelist file,http,https,tcp,tls,crypto")

なお、これはHLSの弱さを伝えるための単なる技術デモであり、フラグメント毎にしか再生できないため実用的ではなく、研究目的以外での利用は想定していません。

再追記。普通に独自スキーマXMLHttpRequestするだけで良いと聞いたのでテストコード。

from selenium import webdriver
from time import sleep
import requests
import re
import os

if __name__ == '__main__':
    browser = webdriver.Chrome(executable_path = "/usr/lib/chromium-browser/chromedriver")
    browser.get("https://abema.tv/now-on-air/abema-news")
    sleep(2)

    pl = requests.get("https://linear-abematv.akamaized.net/channel/abema-news/1080/playlist.m3u8").text

    key_url = re.search(u'URI=\"(.*?)\"\,',pl).group(1)

    browser.execute_script('''
var xhr = new XMLHttpRequest();
xhr.onreadystatechange = function() {
    if (xhr.readyState == 4 && xhr.status == 200)
        window.key = new Uint8Array(xhr.response)
}
xhr.open("GET", "%s");
xhr.send();
'''%key_url)

    sleep(1)
    key = browser.execute_script("return window.key;")
    browser.close()

    f = open("key.bin", "wb")
    f.write(bytearray(key))
    f.close()

    mod_pl = re.sub('URI=.*?\,', 'URI=\"key.bin\",', pl)

    f = open("playlist.m3u8", "w")
    f.write(mod_pl)
    f.close()

    os.system("ffplay playlist.m3u8 -protocol_whitelist file,http,https,tcp,tls,crypto")

ありゃ、本当だ。色々難しく考えすぎてたようです。Javascriptへのcode injectionは不要だし、他のサイトにも使えそう。

ついでに複数フラグメントについても調べてみたら、単にプレーヤー側がリロードを繰り返すだけとのこと。ちょっと信じられないので、とりあえずPython3でプロキシを書いてみた。

from http.server import HTTPServer, SimpleHTTPRequestHandler
import requests
import re

class MyHandler(SimpleHTTPRequestHandler):

    def do_GET(self):
        if self.path == "/key.bin":
            f = open("key.bin", "rb")
            body = f.read()
            f.close()
        else:
            pl = requests.get("https://linear-abematv.akamaized.net/channel/abema-news/1080/playlist.m3u8").content
            body = re.sub(b'URI=.*?\,', b'URI=\"key.bin\",', pl)
        self.send_response(200)
        self.send_header('Content-type', 'application/x-mpegURL')
        self.send_header('Content-length', len(body))
        self.end_headers()
        self.wfile.write(body)

httpd = HTTPServer(('localhost', 8000), MyHandler)
httpd.serve_forever()

…本当ですね。思った以上にザルだった。MPEG-DASHへの一本化が待たれます。

なお、これらコードはエラーチェックが適当ですし、何故かたまに途切れりします。研究目的以外での利用は想定していません。また、暗号化においてHLSを使用することも推奨しません。