bn_gcd_ext_stein returns different Bezout coefficients

[guidovranken]

#include <relic_conf.h>
#include <relic.h>

static void print(const bn_t bn) {
    const size_t size = bn_size_str(bn, 10);
    char* s = malloc(size);
    bn_write_str(s, size, bn, 10);
    printf("%s\n", s);
    free(s);
}

int main(void)
{
    if ( core_init() != RLC_OK ) abort();

    bn_t A, B, R1, R2, R3;

    bn_null(A); bn_new(A);
    bn_null(B); bn_new(B);
    bn_null(R1); bn_new(R1);
    bn_null(R2); bn_new(R2);
    bn_null(R3); bn_new(R3);

    const char* s1 = "5164793035168374968207600271600905026892845810428681291073785203688775362247239268";
    const char* s2 = "5283805872177215917744890733324222647467627372282454631805636520118426898077005012726054741991100315546033159766903297571835109551593675061360850733933177617502325701667760312233485684363633142493422170270643448379568945292725497736661586140491201415317126501617867591215350104665119015015817521615759243407495199353028315105172061139602974617671400272636810437246755816823338394035369836196755181835281965585364403724003277355729041318640395005014628477314706822320341755899800556000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000014757395258967641291";

    bn_read_str(A, s1, strlen(s1), 10);
    bn_read_str(B, s2, strlen(s2), 10);

    bn_gcd_ext(R1, R2, R3, A, B);
    //bn_gcd_ext_stein(R1, R2, R3, A, B);
    printf("X: ");
    print(R2);
    printf("Y: ");
    print(R3);

    bn_free(A);
    bn_free(B);
    bn_free(R1);
    bn_free(R2);
    bn_free(R3);
    return 0;
}

This prints:

X: -2334493656996182823275386098025722837267847670342583447812770965002487040035891667645931211575012067474869815141001371309742222798923582281991903361433816857284793171172202132868762546512738353438924652504161395593771842537881698987155001154022766039701782361067593977767534348745405658716857775521475982325570619648183136384350990434617045215586772640313287885173302080311963162792321219673260978643314314323397819418813873750161820384441780750523021354675280480157084050036814376229406937663406239565105028476082113614911424952113945189610229378030339835936836889085628174260786552146600075710792422081703915494050345069709444888014533732600032661190639204349489820858872040677181673791693570976920590298172331705518969454255362551314741029642246293337165725861217003547776706675706235517044007999616836174040141420777918911034678327710831516964376848368134412431426270231023223035329042517545495058026078980422
Y: 2281911347990235805690167191576998187378956241207477141540806199256473656659570667

But if you comment out bn_gcd_ext and uncomment bn_gcd_ext_stein, it prints:

X: 2949312215181033094469504635298499810199779701939871183992865555115939858041113345080123530416088248071163344625901926262092886752670092779368947372499360760217532530495558179364723137850894789054497517766482052785797102754843798749506584986468435375615344140550273613447815755919713356298959746094283261081924579704845178720821070704985929402084627632323522552073453736511375231243048616523494203191967651261966584305189403605567220934198614254491607122639426342163257705862986179770593062336593760434894971523917886385088575047886054810389770621969660164063163110914371825739213447853399924289207577918296084505949654930290555111985466267399967338809360795650510179141127959322818326208306429023079409701827668294481030545744637448685258970357753706662834274138782996452223293324293764482955992000383163825959858579222081088965321672289168483035623151631865587568573729768976776964670957482469262337232888660869
Y: -2882881687178139162517433080023906839513889569221204149532979004432301705587668601

Both of these are correct in that they satisfy the equation X*A + Y*B == GCD(A,B).

While Bezout coefficients for a given GCD(A,B) are not unique, it is my understanding that the particular values returned by the extended GCD algorithm are unique.

Specifically, if my reading of Wikipedia is correct, they must satisfy:

X < (B / (2 * GCD(A,B)))
Y < (A / (2 * GCD(A,B)))

The result returned by bn_gcd_ext_stein does not satisfy these constraints. Python:

A = 5164793035168374968207600271600905026892845810428681291073785203688775362247239268
B = 5283805872177215917744890733324222647467627372282454631805636520118426898077005012726054741991100315546033159766903297571835109551593675061360850733933177617502325701667760312233485684363633142493422170270643448379568945292725497736661586140491201415317126501617867591215350104665119015015817521615759243407495199353028315105172061139602974617671400272636810437246755816823338394035369836196755181835281965585364403724003277355729041318640395005014628477314706822320341755899800556000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000014757395258967641291

GCD = 1
GCDx2 = GCD * 2

# bn_gcd_ext_stein
X1 = 2949312215181033094469504635298499810199779701939871183992865555115939858041113345080123530416088248071163344625901926262092886752670092779368947372499360760217532530495558179364723137850894789054497517766482052785797102754843798749506584986468435375615344140550273613447815755919713356298959746094283261081924579704845178720821070704985929402084627632323522552073453736511375231243048616523494203191967651261966584305189403605567220934198614254491607122639426342163257705862986179770593062336593760434894971523917886385088575047886054810389770621969660164063163110914371825739213447853399924289207577918296084505949654930290555111985466267399967338809360795650510179141127959322818326208306429023079409701827668294481030545744637448685258970357753706662834274138782996452223293324293764482955992000383163825959858579222081088965321672289168483035623151631865587568573729768976776964670957482469262337232888660869
Y1 = -2882881687178139162517433080023906839513889569221204149532979004432301705587668601

# bn_gcd_ext
X2 = -2334493656996182823275386098025722837267847670342583447812770965002487040035891667645931211575012067474869815141001371309742222798923582281991903361433816857284793171172202132868762546512738353438924652504161395593771842537881698987155001154022766039701782361067593977767534348745405658716857775521475982325570619648183136384350990434617045215586772640313287885173302080311963162792321219673260978643314314323397819418813873750161820384441780750523021354675280480157084050036814376229406937663406239565105028476082113614911424952113945189610229378030339835936836889085628174260786552146600075710792422081703915494050345069709444888014533732600032661190639204349489820858872040677181673791693570976920590298172331705518969454255362551314741029642246293337165725861217003547776706675706235517044007999616836174040141420777918911034678327710831516964376848368134412431426270231023223035329042517545495058026078980422
Y2 = 2281911347990235805690167191576998187378956241207477141540806199256473656659570667

assert(A*X1 + B*Y1 == GCD)
assert(A*X2 + B*Y2 == GCD)

assert(X2 < B / GCDx2)
assert(Y2 < A / GCDx2)

# Fails
assert(X1 < B / GCDx2)
assert(Y1 < A / GCDx2)

Can you confirm this is a bug?

[dfaranha]

Agreed it's a bug!

I never really cared about the bounds and it should be fairly easy to fix. Thanks! :)

[guidovranken]

It now hangs:

#include <relic_conf.h>
#include <relic.h>

int main(void)
{
    if ( core_init() != RLC_OK ) abort();

    bn_t A, B, R1, R2, R3;

    bn_null(A); bn_new(A);
    bn_null(B); bn_new(B);
    bn_null(R1); bn_new(R1);
    bn_null(R2); bn_new(R2);
    bn_null(R3); bn_new(R3);

    const char* s1 = "555555555555555550000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000080000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000081624404648562574776480988416386318723821493777533412583444266365419065804535170339197542765321752482359143484437390596115585268050595800018570555555555555555555555555555500000000000000000000000000000000000000000000000000000";
    const char* s2 = "66666666666666666666666666660000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001123911187445684340258072669751633781322533717977849594343489334935505860661000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000036363310000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000031696932928800000000000000000000000000000000000000";

    bn_read_str(A, s1, strlen(s1), 10);
    bn_read_str(B, s2, strlen(s2), 10);

    bn_gcd_ext_stein(R1, R2, R3, A, B);

    bn_free(A);
    bn_free(B);
    bn_free(R1);
    bn_free(R2);
    bn_free(R3);
    return 0;
}

[dfaranha]

I pushed a not very robust heuristic to fix this later case in 94e5653 to see how it goes. Notice that I also did a minor refactoring of the configuration and now the function is called bn_gcd_ext_binar to match the rest of the library.

[guidovranken]

This one still hangs:

#include <relic_conf.h>
#include <relic.h>

int main(void)
{
    if ( core_init() != RLC_OK ) abort();

    bn_t A, B, R1, R2, R3;

    bn_null(A); bn_new(A);
    bn_null(B); bn_new(B);
    bn_null(R1); bn_new(R1);
    bn_null(R2); bn_new(R2);
    bn_null(R3); bn_new(R3);

    const char* s1 = "7000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000008000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000";
    const char* s2 = "800000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000";

    bn_read_str(A, s1, strlen(s1), 10);
    bn_read_str(B, s2, strlen(s2), 10);

    bn_gcd_ext_binar(R1, R2, R3, A, B);

    bn_free(A);
    bn_free(B);
    bn_free(R1);
    bn_free(R2);
    bn_free(R3);
    return 0;
}

[dfaranha]

I just pushed another, more principled, attempt.

[guidovranken]

I can't find any more issues.
I've enabled the function in my fuzzer: guidovranken/cryptofuzz@947c568
If there is anything else, OSS-Fuzz will notify us.