leif-ibsen / BigInt

Arbitrary-precision integer arithmetic in Swift
MIT License
44 stars 18 forks source link

Believe I have found a large division that fails. #14

Closed wjamesjr closed 1 year ago

wjamesjr commented 1 year ago

`import Foundation

let b = BInt("4782202788054612029528392986600059097414971724022365008513345109918378950942662970278927686112707894586824720981524256319306585052676834087480834429433264797425893247623688331021633208954847354805799943341309825989013743806187109581043148680813778321530496715601563282624414040398143207622036272190408590790537203475256105564071579263867875240985573356522656108542128577321057879052328865035355873615679363655889925711574420153832091752422843046918811427400662135559303516853703976812686385750376227787949580582081831261725701003498206512329872677233489510953469375683037038373999696771585788905639115522613405495707184524158219208223766442059014593330657009722153962376853423770486138578089775621301167811299166407361746606697808186757966914671246073712904200588408923186387737887675292886953797066980967406053530122853539036965490224784924649007954898678503314655546475504501686187354866964374552614120640782949622452027788962138602665933147687696322089504278791624651519312327831756553779377194524673395819281486668576384019590720179413349582970319393884388810494546040342087536563628332152073181614300721769371426238517540520845214665313301183551962591849558938499025348780376716477073930634436840084468255937443451690315999349137664638968972614199015304906547819056227171224947070739716300953775743441307920501863532234466545645695774331885044978250148663467372130392099894852145190998232878772486650513010816769902892518719250066947215706536216248696237476619548294280997466823267797934620158926088915926052113157320399149350072309099187966099088189643497513687928960839375289097070057698545052161834612441440300447910931798424330069494272576133664129571446005941039050449109761662280133298161365994610326139093159039351518411873726699434670218627794470163865731797258607504515467706568302430378820423176019440925413294968918946753234527062907457039327403799514654681482620563976230609510390120177018153378953914211667001350128690497063627989286291690626363431460266484993616898632404483747215192933193062784093302610061072941730610833278036184381530158256539585838588799897714285126286720260761560779987177142728853496833739616878241477395730921177557832273916603145310867588335248530358959095181725126608131193726511914438842826945650698812108368961639111997274610216346844001820682880739430341057410931474128548720115358991355261589298277633257899790971197558786178367745055254357923445045308508288059601997158296253020130655356530996489598877891324128600243475505599610596315196496323542179888792209304804935537197649858327858787700693018393896295993380379016269866179084654381355939982578243521489624917379652576570974212281051937534650720271599262478255593167843763594436777566708623653325250511141415367412037101277051776800861048051216123083433710443499239892012310615341329522047064721158043794002611244331208099097510130494454107820850957462322893961198397453714743859503326095545071195512880003809279")! let a = b * b let c = a / b if c != b { print("failed") } else { print("passed") } exit(0)

`

Please try the above code on your machine. I am using your divide routines in a package that I am developing and the above test fails. I tested with your package and get the same failure. The product b * b is correct, but the division a / b generates a wrong quotient. I apologize for the large size of the number but it resulted from some Mersenne prime test code that I was exercising when I found it. Sorry that I can't provide a smaller number that fails. I am using an Intel Mac but I don't think that should make a difference.

wjamesjr commented 1 year ago

OK, I have spent the morning trying to find another big integer that fails like the above. I literally tried 10s of millions of random numbers with limb count from 60 to 152 (this number) and none fail. Then I tried 10s of millions of random numbers with the same bit count (9689) and again no failures. I then tried the above number shifted right a bit at a time and it fails with bit shifts (0-2137, 2265, 2393, 2521, 2649, 2777, 2905, 3033, 3161, 3289, 3417, 3545, 3673, 3801, 3929, 4057, 4185, 4313, 4441, 4569, 4697, 4825, 4953, 5081, 5209, 5337, 5465, 5593, and 5721.

The number shifted 5721 bits is:

3069183075406921794879924958191326267980083732617257090967156575624265074992632621684985264011687141177463521361517911553021245099194483198628542584486348563466378125452177821735146139462477045607889665505906919596910519383246040059165099934889611077055111196533329299014397179377993006736606277358858761491911450868432133382295467939853925371563407906143030316759860088351677624708599821169550833590036545299179057356555493184183038259794167997381510402990078867498040076511934068262113082029079685194092585200557212952087870250306201570762147399263274328246438930899678576215096246269986660643586905677861972798413323900862279393294098586674450767550167631066398651798197094996096391423300063966130361456553414900428770512722081658104476364858583339960519510499984305492480977699199899564863730624237563221068841591071268000190139910882859705395243310564643593059495182047032409956901498079495621606539723949410789406025669661263490770568454509838497464363714458645954467110174994932482618754701038967174720237992450773625322184776742653465690932102840101388417452307232873785183975126764865419953046185303713655389674631027509028941278698713684303566707778368015911611508474550527223037689855

Smaller but still very mysterious.

I have this sinking feeling that I have a processor issue behind this. I have been putting off a hardware upgrade until M3 macs ship, but this might push be to upgrade now.

Have you been able to duplicate?

wjamesjr commented 1 year ago

I think I have a clue now. Printing the first number in hex and all bits are set.

wjamesjr commented 1 year ago

I spent some time looking over the burnikel code. Knowing that the code works 99.999999% of the time I looked for ways it could have a one-off type of problem. I stumbled on a change that makes the problem go away. Not sure if the fix is the proper one but this is what I found. Changing:

if A1.compare(B1) < 0 {

with:

if A1.compare(B1) <= 0 {

in the 6th line of function DI3n2n() makes the problem go away. I am testing with random values now with success as well. Let me know what you find and think. Thanks,

leif-ibsen commented 1 year ago

Thanks for your investigations. Being on vacation, I just saw your issue now. I will look into it.

leif-ibsen commented 1 year ago

I ran your test and it failed as expected. I also ran it on an iMac where it failed as well.

There is a workaround:

In the BurnikelZiegler.swift file change line 12 to (say)

static let BZ_DIV_LIMIT = 60000000000

Then division will be performed by the basecase Knuth algorithm D, which works (although not as fast)

wjamesjr commented 1 year ago

Thanks for the confirmation! Not urgent, I just wanted to make sure that I did not have an issue. Waiting on a fix is not a problem.

leif-ibsen commented 1 year ago

The bug is fixed in the new release 1.13.1

wjamesjr commented 1 year ago

I integrated the new function in my code and tested and all tests passed. Did you happen to performance test the new function against the old? My timing is saying that the new function is definitely slower. There are several factors though. I am running on an old Intel MacBook versus your M1 I believe. The crossover points for BZDiv and Knuth, and base multiply, Karatsuba, and ToomCook are quite different for me. I have slightly different functions for a number of functions including add, subtract, shift, etc. I liked the code organization of the new BZDiv function better than the old, I think it is easier to understand and follow but for now, I am going to stick with a modified version of the old function.

leif-ibsen commented 1 year ago

1) The old version was coded at the Limb level, the new version is coded at the BInt level which makes it esier to follow the code logic which closely follows the BZ algorithms 1, 2 and 3

2) I didn't do a real comparison of the old performance vs the new, but I did notice that the new is a bit slower which is to be expected when you move from Limbs to BInt's

3) The best crossover point between Knuth and BZ can be debated, I think. I let it stay at 60 Limbs because then I didn't had to change any documentation. Maybe it should be higher, say 100 - 200 Limbs

4) Knuth vs BZ: I found at the 1000 Limb level BZ 3 - 4 times faster than Knuth, at the 10000 Limb level about 10 times faster and at the 100000 Limb level about 30-50 times faster

5) Although I have an iMac M1, i did all development and performance comparisons on an Intel MacBook

6) I included your two failure examples in the test suite