概要とか

「放物線とマス目の関係」という問題が先日締め切られました。多数の挑戦をいただいたんですが、初回投稿の正解率5%以下という状況でした。 間違ってしまった方のために、どうすればよかったのかをコードを交えて解説させていただきたいと思います。 まずは問題のおさらいです。今回は以下のような問題でした:

【概要】

放物線の方程式とマス目を指定するので、位置関係がどうなっているのかを計算して下さい。

【詳細】

放物線とマス目の位置関係は下表の5種類があります(以下、y 座標が大きい方向が上、x座標が大きい方向が右、になります):

記号 状況
U 完全に放物線より上にある
u マス目の内部は放物線より上にあるが、マス目と放物線に共有点がある。
S マス目を放物線が分断し、放物線より上にも下にも 0 より大きな面積がある
L 完全に放物線より下にある
l マス目の内部は放物線より下にあるが、マス目と放物線に共有点がある。
【入出力】

入力は
1 0 0 0 -3
のようになっています。
空白区切りになっていて、順に a, b, c, X, Y です。
a, b, c は、放物線y=ax2+bx+cの係数です。
X, Y はマス目の左下隅の座標です。マス目の大きさは、 1×1 です。
つまり、マス目の領域は (x,y) = { |x,y| X≦x≦X+1, Y≦y≦Y+1 } です。

出力は、
L
のような感じです。
U , u , S , L , l のいずれかを出力して下さい。

【補足】

問題の方は、1×1のマス目と放物線の位置関係を、5つに分類するというものでした。
放物線は整数係数で、マス目の位置も整数です。
実際のテストケースは以下の通りでした:

# 入力 期待
0 8 -8 2 0 -1 l
1 110 -198 90 0 0 S
2 40 -40 10 0 1 S
3 10 -5 0 0 2 S
4 10 -10 0 0 2 U
5 -8 8 -2 0 1 U
6 -3 7 4 -2 -6 u
7 1 -10 25 3 4 u
8 -200 5160 -33275 12 6 S
9 4 -4 1 0 0 S
10 10 -18 9 0 0 S
11 -1 3 0 1 1 l

改めて見直してみると L が入っていませんね。ちょっと自分でびっくりしました。

実装

キーになる x 座標が3つ、あるいは2つあります。
まずはマス目の左端と右端です。必ずキーになります。
もうひとつは放物線の軸です。ただし、軸がマス目を通らない場合はキーにはならなくなります。
キーになる x 座標を放物線の式に代入して、出て来る値がマス目の上にあるのか、あるいは下にあるのかを観察すればよいわけです。

です。ruby で書くとこんな感じです:

#二次式クラス
class Poly
  def initialize( a, b, c )
    @a=a; @b=b; @c=c
  end
  def f(x); @a*x*x+@b*x+@c; end
  def axis; -@b/(2*@a); end
end

#条件の評価
def evaluate( q )
  case q.uniq.sort.join(",")
  when "-1";  "U"
  when "-1,0";  "u"
  when "1";  "L"
  when "0,1";  "l"
  else;  "S"
  end
end

#問題を解く
def solve(s)
  a, b, c, x, y = s.split(" ").map(&:to_r)
  poly = Poly.new( a, b, c )
  y0=poly.f(x)
  y1=poly.f(x+1)
  q=[y0<=>y, y0<=>(y+1), y1<=>y, y1<=>y+1]
  if (x..x+1)===poly.axis
    ya = poly.f( poly.axis )
    q+=[ya<=>y, ya<=>(y+1) ]
  end
  evaluate(q)
end

puts( solve( gets.strip ) )

実装のポイントですが、まずは、有理数の利用です。「接する」というデリケートな情報を扱うので、誤差のない計算が必要です。有理数を使うことで、頭を使うポイントを減らすことができます。
次に、条件の整理です。
4個または6個の評価ポイント(x座標が2個または3個。y座標が2個なので、都合4個または6個になります)が、放物線よりも上か下かを <=> 演算子を使って評価しています。

という具合です。
「全部-1」とか「0と1」という条件を簡単に評価するために、 uniq.sort という処理を行っています。

ちょっと長くなりますが、上記のコードをそのまま C言語(C99)に移植するとこんな感じです:

// clang -std=c99 -Wall
#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>

//有理数
struct rational
{
  int64_t n, d;
};

// 有理数の生成
struct rational rational( int64_t n, int64_t d )
{
  return ( struct rational ){ n, d };
}

// 有理数の乗算
struct rational mul( struct rational a, struct rational b )
{
  return rational( a.n*b.n, a.d*b.d );
}

// 有理数の加算
struct rational add( struct rational a, struct rational b )
{
  return rational( a.n*b.d+b.n*a.d, a.d*b.d );
}

// 分母を正の数にする
struct rational normalize( struct rational a )
{
  return ( a.d<0 ) ? rational(-a.n, -a.d ) : a;
}

// 有理数と整数の比較。大小関係によって 1, 2, 4 を返す
int compare( struct rational a0, int b )
{
  struct rational a = normalize(a0);
  int64_t diff = a.n - b*a.d;
  return 1<<(( 0<diff ) - ( diff<0 ) + 1);
}

// a0<b0 なら真
bool lt( struct rational a0, struct rational b0 )
{
  struct rational a = normalize(a0);
  struct rational b = normalize(b0);
  return  a.n*b.d < b.n*a.d;
}

// 有理数の逆数
struct rational inv( struct rational a )
{
  return rational( a.d, a.n );
}

//二次式
struct poly
{
  struct rational a, b, c;
};

// 二次式p に有理数 x を代入した結果の値
struct rational f( struct poly p, struct rational x )
{
  return add(mul(add(mul(p.a,x),p.b),x),p.c);
}

// 二次式の生成関数
struct poly make_poly( int a, int b, int c )
{
  return ( struct poly ){ { a, 1 }, { b, 1 }, { c,1 } };
}

// 放物線の軸
struct rational axis( struct poly p )
{
  return mul(p.b,inv(mul(p.a,rational(-2, 1))));
}

// 条件の評価
char evaluate( int q )
{
  switch(q){
    case 1:  return 'U';
    case 1+2:  return 'u';
    case 4:  return 'L';
    case 4+2:  return 'l';
    default:  return 'S';
  }
}

// 問題を解く
char solve(char const * s)
{
  int a, b, c, x, y;
  sscanf(s, "%d %d %d %d %d", &a, &b, &c, &x, &y );
  struct poly poly = make_poly( a, b, c );
  struct rational y0=f( poly, rational(x, 1));
  struct rational y1=f( poly, rational(x+1, 1));
  int q = compare(y0,y)|compare(y0,y+1)|compare(y1,y)|compare(y1,y+1);
  struct rational ax = axis(poly);
  if (lt(rational(x,1),ax) && lt(ax,rational(x+1,1))){
    struct rational ya = f( poly, ax );
    q |= compare( ya, y ) | compare( ya, y+1 );
  }
  return evaluate(q);
}

int main()
{
  char buffer[ 1024 ];
  fgets( buffer, sizeof(buffer), stdin );
  putchar( solve( buffer ) );
}

行数にして約半分が 有理数の実装です。これでも約分をサボって楽をしているんですが。
計算を工夫することで有理数は回避できると思いますが、今回は愚直に実装することにしました。
入力は最大 100000。その3乗が計算上必要なので、10の15乗。32bit では桁あふれするので有理数の内部表現では int64_t を使います。
ruby では uniq.sort をしましたが、C言語で同等の処理を実装するのは大変です。代わりに、集合としての値を表現するためにビットを使いました。
あとは ruby と同じなんですが、長くなりますね。
compare 関数の実装と使い方がややトリッキーですが、C言語なので仕方ないでしょう。

最後に

多くの方の挑戦、ありがとうございました。
これに懲りず、また拙問に挑戦いただければと思います。
楽しんでいただける問題を出していきたいと思います。