long long int型において小数点を含む比の数値比較ができません。

Question

前提・実現したいこと

C言語においてある比 φ = (1 + √5) / 2　の値を用いて、標準入力から a + b φ (a, b は整数) の形の 2 つの数に対し，１つめの方が大きければ １ を，２ つめの方が大きければ -1 を，どちらも等しければ 0 を標準出力に出力するプログラムを作成するにおいて、
各整数をlong long int型で絶対値は10^10以下で考えて良いという問題なのですが、
比の値が小数点を含んでおり、long long int型で計算を行うと正しい結果が出力されません。
この場合どのように数値を比較したら良いのでしょうか？

入力例　　　　　　出力例
5 -3　　　　　　　　-1
-3 2

---

3 -1 1
-2 2

---

3 -2 0
3 -2

---

1392 -1209 1
-2789 1375

---

34264 -18079 1
-40761 28289

---

発生している問題・エラーメッセージ

long long int型において小数点を含む大きな値の数値比較ができない。

該当のソースコード

#include<stdio.h>
#include<math.h>

struct golden {
  long long int a;
  long long int b;
};


int compare(struct golden x, struct golden y) {
  long long int s, t;

  s = (x.a + x.b * ( 1 + sqrt(5) )/2) * 10000000;
  t = y.a + y.b * ( 1 + sqrt(5) )/2 * 10000000;
  //printf("s=%lld t=%lld\n",s ,t);                                             

 if(s == t){ return 0; }
 if(s > t){ return 1;}
 else { return -1; }
}

int main(void) {
  struct golden x, y;
  int n;
  scanf("%lld %lld",&x.a, &x.b);
  scanf("%lld %lld",&y.a, &y.b);
  n = compare(x, y);
  printf("%d\n",n);
  return 0;
}

試したこと

小数点を繰り上げようと10^nをかけて繰り上げて比較しようとしても大きい値の入力を得た時にオーバーフローしてしまいました。
またdouble型にキャストも考えましたがlong long int型より有効桁が少ないためできませんでした。

補足情報（FW/ツールのバージョンなど）

ここにより詳細な情報を記載してください。

Comments

[cateye]

環境、コンパイラは何でしょう？　デバッガっで追いかけてみましたか？　( 1 + sqrt(5) )/2;の結果はどうでした？

[apeirogon0813]

unixでコンパイラはgccです。( 1 + sqrt(5) )/2の値はint型だと小数点は切り捨てされてしまうので1でした。これを防ぐために10^n　乗を最後に全体にかけて比較したのですがそれだと大きい値を入力した時オーバーフローしてしまうのか誤った結果が出力されてしまいます。

[wwbQzhMkhhgEmhU]

式変形して、整数計算だけにしてはいけないの？最終的には5と比較するようにすればいいような・・・

[ikadzuchi]

式変形をするとたぶんルートをとるために2乗が必要なのでオーバーフローすると思います。

[episteme]

で、long double や sqrtl() は使えないの? 高精度で演算するのが目的なんしょ?

[apeirogon0813]

sqrtl()使うことによって高精度で演算できました！
莫大な感謝でいっぱいです！！

Answer 1

精度を考慮して、極力浮動小数点の演算を排除して行う方法です。
又、浮動小数点演算を行う場合は、精度の範囲内の演算であるため、正しい結果が得られることが担保されます。こちらで、確認した限りでは、正しい結果が得られてます。必要であれば、実装したソースの提供も可能です。

X=a1+b1Φ
Y=a2+b2Φ
とします。(a1,b1は1件目の入力値、a2,b2は2件目の入力値、Φ=(1+sqrt(5))/2とします）
X-Yを計算し、その結果をZとすると
Z=0 ならX=Y
Z>0 ならX>Y
Z<0 ならX<Yとなります。

Z=a1+b1Φ-(a2+b2Φ)
=(a1-a2)+(b1-b2)Φ
となります。

a1-a2=A
b1-b2=Bとすると
Z=A+BΦ
となります。

ここで、BΦは、整数×無理数＝無理数なので、
Z=0となる為には、A=0,B=0であることが唯一の解であり、これ以外は、Zは0になりません。
よって、上記以外の場合は、Z>0又はZ<0のケースです。
上記以外のケースは以下のようになります。
１）A>0かつB>0なら、Z>0になります。
２）A<0かつB<0なら、Z<0になります。
３）上記以外で、B=0の場合は、Z=Aとなり、A>0なら、Z>0になります。A<0ならZ<0になります。
４）上記以外で、A=0の場合は、Z=BΦとなり、B>0なら、Z>0になります。B<0ならZ<0になります。

ここまでは、整数の演算のみで判定が可能です。
以降は、AもBも0でない場合のケースです。

５）A>0かつB<0の場合、Aの絶対値をAabs、Bの絶対値をBabsとすると
Aabs<Babs*Φなら、Z<0となります。左記以外は、Z>0となります。

６）A<0かつB>0の場合、Aの絶対値をAabs、Bの絶対値をBabsとすると
Aabs<Babs*Φなら、Z>0となります。左記以外は、Z<0となります。

７）上記の５)６）の判定ですが、以下のようにします。
両辺を2倍した
2Aabsと2BabsΦを比較しても同じなので、これを比較します。
右辺に注目すると、
2BabsΦ=2Babs*(1+sqrt(5))/2=Babs(1+sqrt(5))=Babs+Babssqrt(5)となります。
両辺からBabsを引くと
2Aabs-BabsとBabssqrt(5)を比較しても同じことになります。
W=2Aabs-Babs
V=Babssqrt(5)
とするとWは整数型の演算のみで結果が取得できます。
正しい演算が誤差なく行われたとして
VはBabssqrt(5)の真の値の小数点以下を切り捨てた整数の値とします。
即ち、V<真の値<V+1となります。(sqrt(5)は無理数なので真の値は整数にはなりません)
W=<Vなら左辺＜右辺となります。
上記以外は（W>Vのケース）左辺＞右辺となります。

　　正しいVの算出は、以下のようになります。
long long int Babs;
long long int V;
V = Babs * sqrt(5);
sqrt(5)は、double型なので、有効桁数は10進数に換算して約15桁、少なく見積もって14桁となります。
http://www.cc.kyoto-su.ac.jp/~yamada/programming/float.html
sqrt(5)=2.23606797749979のうち、
2.2360679774997　までは信頼できる数値です。Babsの最大値ですが10^10以下なので,10^10倍としても
22360679774.9979　の数値となり整数部は信頼できる範囲の数値です。

Comments

[apeirogon0813]

上記のように実装してみたのですが
以下の入力に対する出力が誤っていました。
コードが合っているか確認していただきたいです。
====================入力====================
51630685 -25993565
-113949456 76340590
=================================================
====================出力====================
1←これが−１になってしまった
=================================================

long long int A, B, Aabs, Babs, V, W;
double q;

q = ( 1 + sqrt(5))/2;
//printf("%f\n",q);
A = x.a - y.a;
B = x.b - y.b;
Aabs = abs(A);
Babs = abs(B);
V = Babs * sqrt(5);
W = 2 * Aabs - Babs;

if(A == 0 && B == 0) { return 0; }　// 0 )

if(A > 0 && B > 0) { return 1; } // 1 )
if(A < 0 && B < 0) { return -1; } // 2 )

if(B == 0) { // 3 )
if(A > 0) { return 1; }
else { return -1; } //A=0なら前の条件でreturnされているためelse
}
if(A == 0) { // 4 )
if(B > 0) {return 1; }
else { return -1; }
}

if(A > 0 && B < 0) { // 5 )
if(W <= V) { return -1;}
else { return 1; }
}
if(A < 0 && B > 0) { // 6 )
if(W <= V) { return 1;}
else { return -1; }
}

[tatsu99]

コードの内容はあっています。
absはllabsを使った方が良いですが、今回の件とは、無関係です。
こちらでも、確認してみました。
結果が-1になっています。
途中経過を確認すると、
V = Babs * sqrt(5);で
V=228826127　となっています。
正確には、
2288261269.9999999125973932163773447315
の値が、途中でまるめられ、
double型の状態で、228826127.00000　となっています。
そのため、long long int型のVに設定した値は
228826127　となっています。
本来なら
V=228826126となるべきですが、そのようになっていませんでした。
double型を使用し、整数倍した場合、整数部分は正しいはずと書きましたが、間違っていました。
お詫びして訂正します。
従って、この方法では解決できません。
他に解決する案として、
GNPで提供されている多倍長整数ライブラリを使用する方法がありますが、
https://sehermitage.web.fc2.com/etc/gmp.html
そのようなものの使用は可能なのでしょうか。
（課題の制約事項として、double型とlong long int型のみ使用可能なのですか）
（linuxでなら使用可能です。gccは６４ビット対応が望ましい）
もし、使用可能なら、それを使用した解決案が提供可能です。
（こちらでは既に実装し、確認済みです）

[tatsu99]

連投です。
GNPで提供されている多倍長整数ライブラリを利用すると、無制限の桁数の整数が扱えます。
但し、1点、利用する場合の注意点として、
long int型が設定値及び取得値の型となっています。
long int型が32ビットのGCCでは、4バイト整数（int型）と同じ
long int型が64ビットのGCCでは、8バイト整数（long long int型）と同じ
と同じです。
ちなみに、long long int型はどちらの環境でも8バイト整数となります。
あなたの環境はどちらでしょうか。
printf("long int size=%d\n",sizeof(long int));
を実行して、８が表示されれば64ビット環境、4が表示されれば32ビット環境です。
32ビット環境で、long long int型の値を多倍長整数ライブラリに設定する場合、一旦、文字列にしてから、その文字列を、設定しなければいけません。
64ビット環境なら、直接、long long int型の値を多倍長整数ライブラリに設定可能です。

[apeirogon0813]

なるほど大きい整数とdouble型の掛け算だと丸められてしまうのですね。
おっしゃるとおり、課題ではdouble型とlong long int型のしていです。
学校のOSはCent OSとlinuxのubuntuがありますがターミナルで提出する感じなのでよくわかりません。
gccは64bit環境でした。
ここまできたら課題関係なしに正しい出力を得たいのでよければご教示願います。

Answer 2

こういう分野は苦手な問題ですが、回答依頼いただいたので頑張ってみました。
（効率とか知りません）

ご提示のコードでできない原因は、sとtがlong long intで
小数部が無視されるからですね。
（既に気づいていらっしゃる通りです。）

で、計算部分ですが、素直にそのまま計算しました。
（一応オーバーフロー対策で式変形して乗算を使わないようにしましたが。）
ただし、計算結果を整数部と小数部に分けてそれぞれ変数に保存し、

long long int s_int = ...
double s_dec = ...
long long int t_int = ...
double t_dec = ...

最後にそれらを比較することで提示された例通りの出力を出せました。

コード載せちゃうと一応動く答えを載せることになっちゃうので、
とりあえずこういうことしましたよっていう回答にしておきます。

Comments

[apeirogon0813]

q_int = (int) (( 1 + sqrt(5))/2);
q_dec =( 1 + sqrt(5))/2 - (int) (( 1 + sqrt(5))/2);
s_int = x.a + x.b * q_int;
t_int = y.a + y.b * q_int;
s_dec = x.b * q_dec;
t_dec = y.b * q_dec;
s_int = (s_int + s_dec) * 1000000;
t_int = (t_int + t_dec) * 1000000;
このように一度整数と少数を分けて（少数部での積の繰り上がりで整数部も含んでいるが）再び結合して実行したところやはり、オーバーフローとなってしまいます

[dice142]

なぜ整数と小数に分けたのにまたくっつけてるんでしょう？

[apeirogon0813]

以下のようなプログラムで入力を
5 -3
-3 2
としたとき
s_int = 1, t_int = 0, s_dec = -0.854102, t_dec = 0.236068
となり大きい方の整数部から比較し、判定していくと　
s > t
となり１が出力されますが実際は小数部も足し合わせるとs<tでー１が正しい出力となり
ただ、大きい整数部から比較していくだけでは正しく動作しないので、結合してしまいました。

long long int s_int, t_int;
double s_dec, t_dec, q;

q = ( 1 + sqrt(5))/2;

s_int = x.a + (int)(x.b * q);
t_int = y.a + (int)(y.b * q);

s_dec = x.b * q - (int)(x.b * q);
t_dec = y.b * q - (int)(y.b * q);

printf("s=%lld t=%lld %f %f\n",s_int ,t_int,s_dec,t_dec);

if(s_int == t_int) {　//整数部が等しかったら小数部で比較
if(s_dec == t_dec){
return 0;
} else if(s_dec > t_dec) {
return 1;
} else {
return -1;
}
}

if(s_int > t_int){ return 1;}　//整数部が等しくなかったら
else { return -1; }

[dice142]

あー、おそらく負の値の切り捨てで整数部の計算が正しくできていない気がします。
本来は-123.456のような数値は-124になるはずですが、-123と計算されてしまう感じがあります。
intによる切り捨て動作はコンパイラの仕様なので、やるなら切り捨て処理を
独自で実装したほうが良いかもしれません。

[dice142]

私は両辺に1/φを掛ける形に変えて計算したので起こらなかったのかもしれません。

[apeirogon0813]

以下のように小数部を負の場合は正にして前回のようなミスを改善し
入力
34264 -18079
-40761 28289
は正しく出力されましたが
より大きい入力
51630685 -25993565
-113949456 76340590
に対しては正しく出力されず
デバックしてみたところ
s_int=9572214 t_int=9572213 s_dec=-0.658780 t_dec=0.341220　//少数部をいじる前
s_int=9572214 t_int=9572213 s_dec=0.000000 t_dec=1.000000
となり整数部が誤っており、なおかつ少数部の繰り上げもなぜか繰り上げされませんでした。

[dice142]

「以下のように」とありますが、どこのことでしょうか？

[apeirogon0813]

申し訳ないです添付し忘れていました。
//少数部がマイナスだと結果に支障をきたすためその分両方に足して正にする
if(s_dec < t_dec) {//２つの少数部のうち小さい方を両方に足せば良い

if(s_dec < 0.) { //小さい方で、かつ足すのは影響が出る負の場合
t_dec += fabs(s_dec);
s_dec += fabs(s_dec);

if(t_dec >= 1.) {//負の小数点の絶対値を足したことによって（少数部>1.）となった時、繰り上げ
t_int++;
t_dec -= 1;
}
}
} else if(t_dec < 0.) { //同様
s_dec += fabs(t_dec);
t_dec += fabs(t_dec);

if(s_dec >= 1.) {//繰り上げ
s_int++;
s_dec -= 1;
}
}

[dice142]

貼っていただいたコードを試しましたが、整数部への切り上げは行われているようです。
無理数の計算なので、おそらく誤差で0.999999...となっているのかもしれません。

ちなみに、マイナスのintによる切り捨てが+方向に付与するのであれば、
「マイナスの時は-1して(int)で切り捨て」すれば良いという案だったのですが、
よくよく考えればfloorで事足りる気がします。

Answer 3

小数点を繰り上げようと10^nをかけて繰り上げて比較しようとしても大きい値の入力を得た時にオーバーフローしてしまいました

とありますが、問題が

絶対値は10^10以下で考えて良い

で、long long は最大が 9*10^18 なので 数桁の下駄を履かせるくらいは問題ないように思えます。

Comments

[apeirogon0813]

ではなぜ正しく出力されないんでしょうか、、、

[gaya-K]

t の式のカッコが足りないですが、もとのソースも同じですか？

[ikadzuchi]

10桁の数値と比較するには10桁の下駄が必要で、オーバーフローするのでそれができないという問題なのではありませんか?

[gaya-K]

こういう場合、φにも有効桁数10桁必要なんでしたっけ？
久しく考えていなかったので、てっきり小数点4,5桁もあれば十分かと思ってました（汗）

[apeirogon0813]

すみませんtの式のカッコ書き忘れていました。ソースはしっかり書いてますが正しく出力されませんでした。

Answer 4

ここまできたら課題関係なしに正しい出力を得たいのでよければご教示願います。

GNPのMPライブラリを使う方法です。
問題は、V = Babs * sqrt(5);で正しい結果が得られないので、この箇所をGNP のMPで対応します。
５の平方根は、
https://keisan.casio.jp/exec/system/1260402326
で取得します。
2.23606797749978969640917366873
なので、10^29倍すると
223606797749978969640917366873
になります。
Babs*223606797749978969640917366873 を行い
その結果を10^29でわります。
これで正しいVが得られます。
ソースでは、
V = mult_by_sqrt5(Babs);
で、mult_by_sqrt5の関数を追加し、Babsに223606797749978969640917366873を掛けて、
その積を10^29で割り、その商を返します。
ソース使用時は、"gmp.h"をincludeします。
使い方は
https://sehermitage.web.fc2.com/etc/gmp.html
を参考にして作りました。

リンクするときは -lgmp　を指定します。
ソース名がsample.cの場合
gcc -o sample sample.c -lgmp
のようにします。

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include "gmp.h"
struct golden {
    long long int a;
    long long int b;
};

long long int mult_by_sqrt5(long long int x){
    // sqrt(5)の値は下記のサイトで取得
    // https://keisan.casio.jp/exec/system/1260402326
    // sqrt(5) = 2.23606797749978969640917366873
    char str_sq5[] = "223606797749978969640917366873";
    char str_10_29[]="100000000000000000000000000000";
    mpz_t mp_a,mp_b,mp_c,mp_x,mp_sq5,mp_10_29;
    long long int ans;
    mpz_init(mp_a);
    mpz_init(mp_b);
    mpz_init(mp_c);
    mpz_init(mp_x);
    mpz_init(mp_sq5);
    mpz_init(mp_10_29);
    mpz_init_set_str(mp_sq5,str_sq5,10);
    mpz_init_set_str(mp_10_29,str_10_29,10);
    mpz_set_si(mp_x,x);
    mpz_mul(mp_a,mp_x,mp_sq5);
    mpz_tdiv_q(mp_b,mp_a,mp_10_29);
    ans = mpz_get_si(mp_b);
    printf("\nmp_a=");
    mpz_out_str(stdout, 10, mp_a);
    printf("\nmp_b=");
    mpz_out_str(stdout, 10, mp_b);
    printf("\nans=%lld\n",ans);
    mpz_clear(mp_a);
    mpz_clear(mp_b);
    mpz_clear(mp_c);
    mpz_clear(mp_x);
    mpz_clear(mp_sq5);
    mpz_clear(mp_10_29);
    return ans;
}
int compare(struct golden x, struct golden y) {
    //戻り値の定義
    // x > y の時:1
    // x < y の時:-1
    // x = y の時:0
    long long int A, B, Aabs, Babs, V, W;
    double q;

    q = ( 1 + sqrt(5))/2;
    //printf("%f\n",q);
    A = x.a - y.a;
    B = x.b - y.b;
    Aabs = llabs(A);
    Babs = llabs(B);
    V = Babs * sqrt(5);
    printf("native V=%lld\n",V);
    V = mult_by_sqrt5(Babs);
    printf("GMP V=%lld\n",V);

    W = 2 * Aabs - Babs;

    if(A == 0 && B == 0) { return 0; }// 0 )

    if(A > 0 && B > 0) { return 1; } // 1 )
    if(A < 0 && B < 0) { return -1; } // 2 )

    if(B == 0) { // 3 )
    if(A > 0) { return 1; }
    else { return -1; } //A=0なら前の条件でreturnされているためelse
    }
    if(A == 0) { // 4 )
    if(B > 0) {return 1; }
    else { return -1; }
    }

    if(A > 0 && B < 0) { // 5 )
    if(W <= V) { return -1;}
    else { return 1; }
    }
    if(A < 0 && B > 0) { // 6 )
    if(W <= V) { return 1;}
    else { return -1; }
    }
}

int main(void) {
    struct golden x, y;
    int n;
    scanf("%lld %lld",&x.a, &x.b);
    scanf("%lld %lld",&y.a, &y.b);
    printf("x.a=%lld x.b=%lld y.a=%lld y.b=%lld\n",x.a,x.b,y.a,y.b);
    n = compare(x, y);
    printf("%d\n",n);
    return 0;
}

Answer 5

まず、PCが扱う浮動小数点数というのは２進数の数値となります
そして、あなたが扱う数値は１０進数の数値となります
その１０進数と２進数の数値の違いが、いわゆる誤差、となり、あなたのいう正しく出力されない結果となってしまいます

ということで、そこらへん調べてみましょう

＃あなたのいう0.1は、浮動小数点数では表現できません