BigInteger.cs 102 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380238123822383238423852386238723882389239023912392239323942395239623972398239924002401240224032404240524062407240824092410241124122413241424152416241724182419242024212422242324242425242624272428242924302431243224332434243524362437243824392440244124422443244424452446244724482449245024512452245324542455245624572458245924602461246224632464246524662467246824692470247124722473247424752476247724782479248024812482248324842485248624872488248924902491249224932494249524962497249824992500250125022503250425052506250725082509251025112512251325142515251625172518251925202521252225232524252525262527252825292530253125322533253425352536253725382539254025412542254325442545254625472548254925502551255225532554255525562557255825592560256125622563256425652566256725682569257025712572257325742575257625772578257925802581258225832584258525862587258825892590259125922593259425952596259725982599260026012602260326042605260626072608260926102611261226132614261526162617261826192620262126222623262426252626262726282629263026312632263326342635263626372638263926402641264226432644264526462647264826492650265126522653265426552656265726582659266026612662266326642665266626672668266926702671267226732674267526762677267826792680268126822683268426852686268726882689269026912692269326942695269626972698269927002701270227032704270527062707270827092710271127122713271427152716271727182719272027212722272327242725272627272728272927302731273227332734273527362737273827392740274127422743274427452746274727482749275027512752275327542755275627572758275927602761276227632764276527662767276827692770277127722773277427752776277727782779278027812782278327842785278627872788278927902791279227932794279527962797279827992800280128022803280428052806280728082809281028112812281328142815281628172818281928202821282228232824282528262827282828292830283128322833283428352836283728382839284028412842284328442845284628472848284928502851285228532854285528562857285828592860286128622863286428652866286728682869287028712872287328742875287628772878287928802881288228832884288528862887288828892890289128922893289428952896289728982899290029012902290329042905290629072908290929102911291229132914291529162917291829192920292129222923292429252926292729282929293029312932293329342935293629372938293929402941294229432944294529462947294829492950295129522953295429552956295729582959296029612962296329642965296629672968296929702971297229732974297529762977297829792980298129822983298429852986298729882989299029912992299329942995299629972998299930003001300230033004300530063007300830093010301130123013301430153016301730183019302030213022302330243025302630273028302930303031303230333034303530363037303830393040304130423043304430453046304730483049305030513052305330543055305630573058305930603061306230633064306530663067306830693070307130723073307430753076307730783079308030813082308330843085308630873088308930903091309230933094309530963097309830993100310131023103310431053106310731083109311031113112311331143115311631173118311931203121312231233124312531263127312831293130313131323133313431353136313731383139314031413142314331443145314631473148314931503151315231533154315531563157315831593160316131623163316431653166316731683169317031713172317331743175317631773178317931803181318231833184318531863187318831893190319131923193319431953196319731983199320032013202320332043205320632073208320932103211321232133214321532163217321832193220322132223223322432253226322732283229323032313232323332343235323632373238323932403241324232433244324532463247324832493250325132523253325432553256325732583259326032613262326332643265326632673268326932703271327232733274327532763277327832793280328132823283328432853286328732883289329032913292329332943295329632973298329933003301330233033304330533063307330833093310331133123313331433153316331733183319332033213322332333243325332633273328332933303331333233333334333533363337333833393340334133423343334433453346334733483349335033513352335333543355335633573358335933603361336233633364336533663367336833693370337133723373337433753376337733783379338033813382338333843385338633873388338933903391339233933394339533963397339833993400340134023403340434053406340734083409341034113412341334143415341634173418341934203421342234233424342534263427342834293430343134323433343434353436343734383439344034413442344334443445344634473448344934503451345234533454345534563457345834593460346134623463346434653466346734683469347034713472347334743475347634773478347934803481348234833484348534863487348834893490349134923493349434953496349734983499350035013502350335043505350635073508350935103511351235133514351535163517351835193520352135223523352435253526352735283529353035313532353335343535353635373538353935403541354235433544354535463547354835493550355135523553355435553556355735583559356035613562356335643565356635673568356935703571357235733574357535763577357835793580358135823583358435853586358735883589359035913592359335943595359635973598359936003601360236033604
  1. #if !BESTHTTP_DISABLE_ALTERNATE_SSL && (!UNITY_WEBGL || UNITY_EDITOR)
  2. #pragma warning disable
  3. using System;
  4. using System.Collections;
  5. using System.Diagnostics;
  6. using System.Globalization;
  7. using System.Text;
  8. using BestHTTP.SecureProtocol.Org.BouncyCastle.Security;
  9. using BestHTTP.SecureProtocol.Org.BouncyCastle.Utilities;
  10. namespace BestHTTP.SecureProtocol.Org.BouncyCastle.Math
  11. {
  12. #if !(NETCF_1_0 || NETCF_2_0 || SILVERLIGHT || PORTABLE || NETFX_CORE)
  13. [Serializable]
  14. #endif
  15. public class BigInteger
  16. {
  17. // The first few odd primes
  18. /*
  19. 3 5 7 11 13 17 19 23 29
  20. 31 37 41 43 47 53 59 61 67 71
  21. 73 79 83 89 97 101 103 107 109 113
  22. 127 131 137 139 149 151 157 163 167 173
  23. 179 181 191 193 197 199 211 223 227 229
  24. 233 239 241 251 257 263 269 271 277 281
  25. 283 293 307 311 313 317 331 337 347 349
  26. 353 359 367 373 379 383 389 397 401 409
  27. 419 421 431 433 439 443 449 457 461 463
  28. 467 479 487 491 499 503 509 521 523 541
  29. 547 557 563 569 571 577 587 593 599 601
  30. 607 613 617 619 631 641 643 647 653 659
  31. 661 673 677 683 691 701 709 719 727 733
  32. 739 743 751 757 761 769 773 787 797 809
  33. 811 821 823 827 829 839 853 857 859 863
  34. 877 881 883 887 907 911 919 929 937 941
  35. 947 953 967 971 977 983 991 997 1009
  36. 1013 1019 1021 1031 1033 1039 1049 1051
  37. 1061 1063 1069 1087 1091 1093 1097 1103
  38. 1109 1117 1123 1129 1151 1153 1163 1171
  39. 1181 1187 1193 1201 1213 1217 1223 1229
  40. 1231 1237 1249 1259 1277 1279 1283 1289
  41. */
  42. // Each list has a product < 2^31
  43. internal static readonly int[][] primeLists = new int[][]
  44. {
  45. new int[]{ 3, 5, 7, 11, 13, 17, 19, 23 },
  46. new int[]{ 29, 31, 37, 41, 43 },
  47. new int[]{ 47, 53, 59, 61, 67 },
  48. new int[]{ 71, 73, 79, 83 },
  49. new int[]{ 89, 97, 101, 103 },
  50. new int[]{ 107, 109, 113, 127 },
  51. new int[]{ 131, 137, 139, 149 },
  52. new int[]{ 151, 157, 163, 167 },
  53. new int[]{ 173, 179, 181, 191 },
  54. new int[]{ 193, 197, 199, 211 },
  55. new int[]{ 223, 227, 229 },
  56. new int[]{ 233, 239, 241 },
  57. new int[]{ 251, 257, 263 },
  58. new int[]{ 269, 271, 277 },
  59. new int[]{ 281, 283, 293 },
  60. new int[]{ 307, 311, 313 },
  61. new int[]{ 317, 331, 337 },
  62. new int[]{ 347, 349, 353 },
  63. new int[]{ 359, 367, 373 },
  64. new int[]{ 379, 383, 389 },
  65. new int[]{ 397, 401, 409 },
  66. new int[]{ 419, 421, 431 },
  67. new int[]{ 433, 439, 443 },
  68. new int[]{ 449, 457, 461 },
  69. new int[]{ 463, 467, 479 },
  70. new int[]{ 487, 491, 499 },
  71. new int[]{ 503, 509, 521 },
  72. new int[]{ 523, 541, 547 },
  73. new int[]{ 557, 563, 569 },
  74. new int[]{ 571, 577, 587 },
  75. new int[]{ 593, 599, 601 },
  76. new int[]{ 607, 613, 617 },
  77. new int[]{ 619, 631, 641 },
  78. new int[]{ 643, 647, 653 },
  79. new int[]{ 659, 661, 673 },
  80. new int[]{ 677, 683, 691 },
  81. new int[]{ 701, 709, 719 },
  82. new int[]{ 727, 733, 739 },
  83. new int[]{ 743, 751, 757 },
  84. new int[]{ 761, 769, 773 },
  85. new int[]{ 787, 797, 809 },
  86. new int[]{ 811, 821, 823 },
  87. new int[]{ 827, 829, 839 },
  88. new int[]{ 853, 857, 859 },
  89. new int[]{ 863, 877, 881 },
  90. new int[]{ 883, 887, 907 },
  91. new int[]{ 911, 919, 929 },
  92. new int[]{ 937, 941, 947 },
  93. new int[]{ 953, 967, 971 },
  94. new int[]{ 977, 983, 991 },
  95. new int[]{ 997, 1009, 1013 },
  96. new int[]{ 1019, 1021, 1031 },
  97. new int[]{ 1033, 1039, 1049 },
  98. new int[]{ 1051, 1061, 1063 },
  99. new int[]{ 1069, 1087, 1091 },
  100. new int[]{ 1093, 1097, 1103 },
  101. new int[]{ 1109, 1117, 1123 },
  102. new int[]{ 1129, 1151, 1153 },
  103. new int[]{ 1163, 1171, 1181 },
  104. new int[]{ 1187, 1193, 1201 },
  105. new int[]{ 1213, 1217, 1223 },
  106. new int[]{ 1229, 1231, 1237 },
  107. new int[]{ 1249, 1259, 1277 },
  108. new int[]{ 1279, 1283, 1289 },
  109. };
  110. internal static readonly int[] primeProducts;
  111. private const long IMASK = 0xFFFFFFFFL;
  112. private const ulong UIMASK = 0xFFFFFFFFUL;
  113. private static readonly int[] ZeroMagnitude = new int[0];
  114. private static readonly byte[] ZeroEncoding = new byte[0];
  115. private static readonly BigInteger[] SMALL_CONSTANTS = new BigInteger[17];
  116. public static readonly BigInteger Zero;
  117. public static readonly BigInteger One;
  118. public static readonly BigInteger Two;
  119. public static readonly BigInteger Three;
  120. public static readonly BigInteger Ten;
  121. //private readonly static byte[] BitCountTable =
  122. //{
  123. // 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4,
  124. // 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
  125. // 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
  126. // 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  127. // 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
  128. // 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  129. // 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  130. // 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  131. // 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
  132. // 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  133. // 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  134. // 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  135. // 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  136. // 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  137. // 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  138. // 4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8
  139. //};
  140. private readonly static byte[] BitLengthTable =
  141. {
  142. 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4,
  143. 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
  144. 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
  145. 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
  146. 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
  147. 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
  148. 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
  149. 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
  150. 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
  151. 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
  152. 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
  153. 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
  154. 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
  155. 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
  156. 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
  157. 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8
  158. };
  159. // TODO Parse radix-2 64 bits at a time and radix-8 63 bits at a time
  160. private const int chunk2 = 1, chunk8 = 1, chunk10 = 19, chunk16 = 16;
  161. private static readonly BigInteger radix2, radix2E, radix8, radix8E, radix10, radix10E, radix16, radix16E;
  162. private static readonly SecureRandom RandomSource = new SecureRandom();
  163. /*
  164. * These are the threshold bit-lengths (of an exponent) where we increase the window size.
  165. * They are calculated according to the expected savings in multiplications.
  166. * Some squares will also be saved on average, but we offset these against the extra storage costs.
  167. */
  168. private static readonly int[] ExpWindowThresholds = { 7, 25, 81, 241, 673, 1793, 4609, Int32.MaxValue };
  169. private const int BitsPerByte = 8;
  170. private const int BitsPerInt = 32;
  171. private const int BytesPerInt = 4;
  172. static BigInteger()
  173. {
  174. Zero = new BigInteger(0, ZeroMagnitude, false);
  175. Zero.nBits = 0; Zero.nBitLength = 0;
  176. SMALL_CONSTANTS[0] = Zero;
  177. for (uint i = 1; i < SMALL_CONSTANTS.Length; ++i)
  178. {
  179. SMALL_CONSTANTS[i] = CreateUValueOf(i);
  180. }
  181. One = SMALL_CONSTANTS[1];
  182. Two = SMALL_CONSTANTS[2];
  183. Three = SMALL_CONSTANTS[3];
  184. Ten = SMALL_CONSTANTS[10];
  185. radix2 = ValueOf(2);
  186. radix2E = radix2.Pow(chunk2);
  187. radix8 = ValueOf(8);
  188. radix8E = radix8.Pow(chunk8);
  189. radix10 = ValueOf(10);
  190. radix10E = radix10.Pow(chunk10);
  191. radix16 = ValueOf(16);
  192. radix16E = radix16.Pow(chunk16);
  193. primeProducts = new int[primeLists.Length];
  194. for (int i = 0; i < primeLists.Length; ++i)
  195. {
  196. int[] primeList = primeLists[i];
  197. int product = primeList[0];
  198. for (int j = 1; j < primeList.Length; ++j)
  199. {
  200. product *= primeList[j];
  201. }
  202. primeProducts[i] = product;
  203. }
  204. }
  205. private int[] magnitude; // array of ints with [0] being the most significant
  206. private int sign; // -1 means -ve; +1 means +ve; 0 means 0;
  207. private int nBits = -1; // cache BitCount() value
  208. private int nBitLength = -1; // cache BitLength() value
  209. private int mQuote = 0; // -m^(-1) mod b, b = 2^32 (see Montgomery mult.), 0 when uninitialised
  210. private static int GetByteLength(
  211. int nBits)
  212. {
  213. return (nBits + BitsPerByte - 1) / BitsPerByte;
  214. }
  215. internal static BigInteger Arbitrary(int sizeInBits)
  216. {
  217. return new BigInteger(sizeInBits, RandomSource);
  218. }
  219. private BigInteger(
  220. int signum,
  221. int[] mag,
  222. bool checkMag)
  223. {
  224. if (checkMag)
  225. {
  226. int i = 0;
  227. while (i < mag.Length && mag[i] == 0)
  228. {
  229. ++i;
  230. }
  231. if (i == mag.Length)
  232. {
  233. this.sign = 0;
  234. this.magnitude = ZeroMagnitude;
  235. }
  236. else
  237. {
  238. this.sign = signum;
  239. if (i == 0)
  240. {
  241. this.magnitude = mag;
  242. }
  243. else
  244. {
  245. // strip leading 0 words
  246. this.magnitude = new int[mag.Length - i];
  247. Array.Copy(mag, i, this.magnitude, 0, this.magnitude.Length);
  248. }
  249. }
  250. }
  251. else
  252. {
  253. this.sign = signum;
  254. this.magnitude = mag;
  255. }
  256. }
  257. public BigInteger(
  258. string value)
  259. : this(value, 10)
  260. {
  261. }
  262. public BigInteger(
  263. string str,
  264. int radix)
  265. {
  266. if (str.Length == 0)
  267. throw new FormatException("Zero length BigInteger");
  268. NumberStyles style;
  269. int chunk;
  270. BigInteger r;
  271. BigInteger rE;
  272. switch (radix)
  273. {
  274. case 2:
  275. // Is there anyway to restrict to binary digits?
  276. style = NumberStyles.Integer;
  277. chunk = chunk2;
  278. r = radix2;
  279. rE = radix2E;
  280. break;
  281. case 8:
  282. // Is there anyway to restrict to octal digits?
  283. style = NumberStyles.Integer;
  284. chunk = chunk8;
  285. r = radix8;
  286. rE = radix8E;
  287. break;
  288. case 10:
  289. // This style seems to handle spaces and minus sign already (our processing redundant?)
  290. style = NumberStyles.Integer;
  291. chunk = chunk10;
  292. r = radix10;
  293. rE = radix10E;
  294. break;
  295. case 16:
  296. // TODO Should this be HexNumber?
  297. style = NumberStyles.AllowHexSpecifier;
  298. chunk = chunk16;
  299. r = radix16;
  300. rE = radix16E;
  301. break;
  302. default:
  303. throw new FormatException("Only bases 2, 8, 10, or 16 allowed");
  304. }
  305. int index = 0;
  306. sign = 1;
  307. if (str[0] == '-')
  308. {
  309. if (str.Length == 1)
  310. throw new FormatException("Zero length BigInteger");
  311. sign = -1;
  312. index = 1;
  313. }
  314. // strip leading zeros from the string str
  315. while (index < str.Length && Int32.Parse(str[index].ToString(), style) == 0)
  316. {
  317. index++;
  318. }
  319. if (index >= str.Length)
  320. {
  321. // zero value - we're done
  322. sign = 0;
  323. magnitude = ZeroMagnitude;
  324. return;
  325. }
  326. //////
  327. // could we work out the max number of ints required to store
  328. // str.Length digits in the given base, then allocate that
  329. // storage in one hit?, then Generate the magnitude in one hit too?
  330. //////
  331. BigInteger b = Zero;
  332. int next = index + chunk;
  333. if (next <= str.Length)
  334. {
  335. do
  336. {
  337. string s = str.Substring(index, chunk);
  338. ulong i = ulong.Parse(s, style);
  339. BigInteger bi = CreateUValueOf(i);
  340. switch (radix)
  341. {
  342. case 2:
  343. // TODO Need this because we are parsing in radix 10 above
  344. if (i >= 2)
  345. throw new FormatException("Bad character in radix 2 string: " + s);
  346. // TODO Parse 64 bits at a time
  347. b = b.ShiftLeft(1);
  348. break;
  349. case 8:
  350. // TODO Need this because we are parsing in radix 10 above
  351. if (i >= 8)
  352. throw new FormatException("Bad character in radix 8 string: " + s);
  353. // TODO Parse 63 bits at a time
  354. b = b.ShiftLeft(3);
  355. break;
  356. case 16:
  357. b = b.ShiftLeft(64);
  358. break;
  359. default:
  360. b = b.Multiply(rE);
  361. break;
  362. }
  363. b = b.Add(bi);
  364. index = next;
  365. next += chunk;
  366. }
  367. while (next <= str.Length);
  368. }
  369. if (index < str.Length)
  370. {
  371. string s = str.Substring(index);
  372. ulong i = ulong.Parse(s, style);
  373. BigInteger bi = CreateUValueOf(i);
  374. if (b.sign > 0)
  375. {
  376. if (radix == 2)
  377. {
  378. // NB: Can't reach here since we are parsing one char at a time
  379. Debug.Assert(false);
  380. // TODO Parse all bits at once
  381. // b = b.ShiftLeft(s.Length);
  382. }
  383. else if (radix == 8)
  384. {
  385. // NB: Can't reach here since we are parsing one char at a time
  386. Debug.Assert(false);
  387. // TODO Parse all bits at once
  388. // b = b.ShiftLeft(s.Length * 3);
  389. }
  390. else if (radix == 16)
  391. {
  392. b = b.ShiftLeft(s.Length << 2);
  393. }
  394. else
  395. {
  396. b = b.Multiply(r.Pow(s.Length));
  397. }
  398. b = b.Add(bi);
  399. }
  400. else
  401. {
  402. b = bi;
  403. }
  404. }
  405. // Note: This is the previous (slower) algorithm
  406. // while (index < value.Length)
  407. // {
  408. // char c = value[index];
  409. // string s = c.ToString();
  410. // int i = Int32.Parse(s, style);
  411. //
  412. // b = b.Multiply(r).Add(ValueOf(i));
  413. // index++;
  414. // }
  415. magnitude = b.magnitude;
  416. }
  417. public BigInteger(
  418. byte[] bytes)
  419. : this(bytes, 0, bytes.Length)
  420. {
  421. }
  422. public BigInteger(
  423. byte[] bytes,
  424. int offset,
  425. int length)
  426. {
  427. if (length == 0)
  428. throw new FormatException("Zero length BigInteger");
  429. // TODO Move this processing into MakeMagnitude (provide sign argument)
  430. if ((sbyte)bytes[offset] < 0)
  431. {
  432. this.sign = -1;
  433. int end = offset + length;
  434. int iBval;
  435. // strip leading sign bytes
  436. for (iBval = offset; iBval < end && ((sbyte)bytes[iBval] == -1); iBval++)
  437. {
  438. }
  439. if (iBval >= end)
  440. {
  441. this.magnitude = One.magnitude;
  442. }
  443. else
  444. {
  445. int numBytes = end - iBval;
  446. byte[] inverse = new byte[numBytes];
  447. int index = 0;
  448. while (index < numBytes)
  449. {
  450. inverse[index++] = (byte)~bytes[iBval++];
  451. }
  452. Debug.Assert(iBval == end);
  453. while (inverse[--index] == byte.MaxValue)
  454. {
  455. inverse[index] = byte.MinValue;
  456. }
  457. inverse[index]++;
  458. this.magnitude = MakeMagnitude(inverse, 0, inverse.Length);
  459. }
  460. }
  461. else
  462. {
  463. // strip leading zero bytes and return magnitude bytes
  464. this.magnitude = MakeMagnitude(bytes, offset, length);
  465. this.sign = this.magnitude.Length > 0 ? 1 : 0;
  466. }
  467. }
  468. private static int[] MakeMagnitude(
  469. byte[] bytes,
  470. int offset,
  471. int length)
  472. {
  473. int end = offset + length;
  474. // strip leading zeros
  475. int firstSignificant;
  476. for (firstSignificant = offset; firstSignificant < end
  477. && bytes[firstSignificant] == 0; firstSignificant++)
  478. {
  479. }
  480. if (firstSignificant >= end)
  481. {
  482. return ZeroMagnitude;
  483. }
  484. int nInts = (end - firstSignificant + 3) / BytesPerInt;
  485. int bCount = (end - firstSignificant) % BytesPerInt;
  486. if (bCount == 0)
  487. {
  488. bCount = BytesPerInt;
  489. }
  490. if (nInts < 1)
  491. {
  492. return ZeroMagnitude;
  493. }
  494. int[] mag = new int[nInts];
  495. int v = 0;
  496. int magnitudeIndex = 0;
  497. for (int i = firstSignificant; i < end; ++i)
  498. {
  499. v <<= 8;
  500. v |= bytes[i] & 0xff;
  501. bCount--;
  502. if (bCount <= 0)
  503. {
  504. mag[magnitudeIndex] = v;
  505. magnitudeIndex++;
  506. bCount = BytesPerInt;
  507. v = 0;
  508. }
  509. }
  510. if (magnitudeIndex < mag.Length)
  511. {
  512. mag[magnitudeIndex] = v;
  513. }
  514. return mag;
  515. }
  516. public BigInteger(
  517. int sign,
  518. byte[] bytes)
  519. : this(sign, bytes, 0, bytes.Length)
  520. {
  521. }
  522. public BigInteger(
  523. int sign,
  524. byte[] bytes,
  525. int offset,
  526. int length)
  527. {
  528. if (sign < -1 || sign > 1)
  529. throw new FormatException("Invalid sign value");
  530. if (sign == 0)
  531. {
  532. this.sign = 0;
  533. this.magnitude = ZeroMagnitude;
  534. }
  535. else
  536. {
  537. // copy bytes
  538. this.magnitude = MakeMagnitude(bytes, offset, length);
  539. this.sign = this.magnitude.Length < 1 ? 0 : sign;
  540. }
  541. }
  542. public BigInteger(
  543. int sizeInBits,
  544. Random random)
  545. {
  546. if (sizeInBits < 0)
  547. throw new ArgumentException("sizeInBits must be non-negative");
  548. this.nBits = -1;
  549. this.nBitLength = -1;
  550. if (sizeInBits == 0)
  551. {
  552. this.sign = 0;
  553. this.magnitude = ZeroMagnitude;
  554. return;
  555. }
  556. int nBytes = GetByteLength(sizeInBits);
  557. byte[] b = new byte[nBytes];
  558. random.NextBytes(b);
  559. // strip off any excess bits in the MSB
  560. int xBits = BitsPerByte * nBytes - sizeInBits;
  561. b[0] &= (byte)(255U >> xBits);
  562. this.magnitude = MakeMagnitude(b, 0, b.Length);
  563. this.sign = this.magnitude.Length < 1 ? 0 : 1;
  564. }
  565. public BigInteger(
  566. int bitLength,
  567. int certainty,
  568. Random random)
  569. {
  570. if (bitLength < 2)
  571. throw new ArithmeticException("bitLength < 2");
  572. this.sign = 1;
  573. this.nBitLength = bitLength;
  574. if (bitLength == 2)
  575. {
  576. this.magnitude = random.Next(2) == 0
  577. ? Two.magnitude
  578. : Three.magnitude;
  579. return;
  580. }
  581. int nBytes = GetByteLength(bitLength);
  582. byte[] b = new byte[nBytes];
  583. int xBits = BitsPerByte * nBytes - bitLength;
  584. byte mask = (byte)(255U >> xBits);
  585. byte lead = (byte)(1 << (7 - xBits));
  586. for (;;)
  587. {
  588. random.NextBytes(b);
  589. // strip off any excess bits in the MSB
  590. b[0] &= mask;
  591. // ensure the leading bit is 1 (to meet the strength requirement)
  592. b[0] |= lead;
  593. // ensure the trailing bit is 1 (i.e. must be odd)
  594. b[nBytes - 1] |= 1;
  595. this.magnitude = MakeMagnitude(b, 0, b.Length);
  596. this.nBits = -1;
  597. this.mQuote = 0;
  598. if (certainty < 1)
  599. break;
  600. if (CheckProbablePrime(certainty, random, true))
  601. break;
  602. for (int j = 1; j < (magnitude.Length - 1); ++j)
  603. {
  604. this.magnitude[j] ^= random.Next();
  605. if (CheckProbablePrime(certainty, random, true))
  606. return;
  607. }
  608. }
  609. }
  610. public BigInteger Abs()
  611. {
  612. return sign >= 0 ? this : Negate();
  613. }
  614. /**
  615. * return a = a + b - b preserved.
  616. */
  617. private static int[] AddMagnitudes(
  618. int[] a,
  619. int[] b)
  620. {
  621. int tI = a.Length - 1;
  622. int vI = b.Length - 1;
  623. long m = 0;
  624. while (vI >= 0)
  625. {
  626. m += ((long)(uint)a[tI] + (long)(uint)b[vI--]);
  627. a[tI--] = (int)m;
  628. m = (long)((ulong)m >> 32);
  629. }
  630. if (m != 0)
  631. {
  632. while (tI >= 0 && ++a[tI--] == 0)
  633. {
  634. }
  635. }
  636. return a;
  637. }
  638. public BigInteger Add(
  639. BigInteger value)
  640. {
  641. if (this.sign == 0)
  642. return value;
  643. if (this.sign != value.sign)
  644. {
  645. if (value.sign == 0)
  646. return this;
  647. if (value.sign < 0)
  648. return Subtract(value.Negate());
  649. return value.Subtract(Negate());
  650. }
  651. return AddToMagnitude(value.magnitude);
  652. }
  653. private BigInteger AddToMagnitude(
  654. int[] magToAdd)
  655. {
  656. int[] big, small;
  657. if (this.magnitude.Length < magToAdd.Length)
  658. {
  659. big = magToAdd;
  660. small = this.magnitude;
  661. }
  662. else
  663. {
  664. big = this.magnitude;
  665. small = magToAdd;
  666. }
  667. // Conservatively avoid over-allocation when no overflow possible
  668. uint limit = uint.MaxValue;
  669. if (big.Length == small.Length)
  670. limit -= (uint) small[0];
  671. bool possibleOverflow = (uint) big[0] >= limit;
  672. int[] bigCopy;
  673. if (possibleOverflow)
  674. {
  675. bigCopy = new int[big.Length + 1];
  676. big.CopyTo(bigCopy, 1);
  677. }
  678. else
  679. {
  680. bigCopy = (int[]) big.Clone();
  681. }
  682. bigCopy = AddMagnitudes(bigCopy, small);
  683. return new BigInteger(this.sign, bigCopy, possibleOverflow);
  684. }
  685. public BigInteger And(
  686. BigInteger value)
  687. {
  688. if (this.sign == 0 || value.sign == 0)
  689. {
  690. return Zero;
  691. }
  692. int[] aMag = this.sign > 0
  693. ? this.magnitude
  694. : Add(One).magnitude;
  695. int[] bMag = value.sign > 0
  696. ? value.magnitude
  697. : value.Add(One).magnitude;
  698. bool resultNeg = sign < 0 && value.sign < 0;
  699. int resultLength = System.Math.Max(aMag.Length, bMag.Length);
  700. int[] resultMag = new int[resultLength];
  701. int aStart = resultMag.Length - aMag.Length;
  702. int bStart = resultMag.Length - bMag.Length;
  703. for (int i = 0; i < resultMag.Length; ++i)
  704. {
  705. int aWord = i >= aStart ? aMag[i - aStart] : 0;
  706. int bWord = i >= bStart ? bMag[i - bStart] : 0;
  707. if (this.sign < 0)
  708. {
  709. aWord = ~aWord;
  710. }
  711. if (value.sign < 0)
  712. {
  713. bWord = ~bWord;
  714. }
  715. resultMag[i] = aWord & bWord;
  716. if (resultNeg)
  717. {
  718. resultMag[i] = ~resultMag[i];
  719. }
  720. }
  721. BigInteger result = new BigInteger(1, resultMag, true);
  722. // TODO Optimise this case
  723. if (resultNeg)
  724. {
  725. result = result.Not();
  726. }
  727. return result;
  728. }
  729. public BigInteger AndNot(
  730. BigInteger val)
  731. {
  732. return And(val.Not());
  733. }
  734. public int BitCount
  735. {
  736. get
  737. {
  738. if (nBits == -1)
  739. {
  740. if (sign < 0)
  741. {
  742. // TODO Optimise this case
  743. nBits = Not().BitCount;
  744. }
  745. else
  746. {
  747. int sum = 0;
  748. for (int i = 0; i < magnitude.Length; ++i)
  749. {
  750. sum += BitCnt(magnitude[i]);
  751. }
  752. nBits = sum;
  753. }
  754. }
  755. return nBits;
  756. }
  757. }
  758. public static int BitCnt(int i)
  759. {
  760. uint u = (uint)i;
  761. u = u - ((u >> 1) & 0x55555555);
  762. u = (u & 0x33333333) + ((u >> 2) & 0x33333333);
  763. u = (u + (u >> 4)) & 0x0f0f0f0f;
  764. u += (u >> 8);
  765. u += (u >> 16);
  766. u &= 0x3f;
  767. return (int)u;
  768. }
  769. private static int CalcBitLength(int sign, int indx, int[] mag)
  770. {
  771. for (;;)
  772. {
  773. if (indx >= mag.Length)
  774. return 0;
  775. if (mag[indx] != 0)
  776. break;
  777. ++indx;
  778. }
  779. // bit length for everything after the first int
  780. int bitLength = 32 * ((mag.Length - indx) - 1);
  781. // and determine bitlength of first int
  782. int firstMag = mag[indx];
  783. bitLength += BitLen(firstMag);
  784. // Check for negative powers of two
  785. if (sign < 0 && ((firstMag & -firstMag) == firstMag))
  786. {
  787. do
  788. {
  789. if (++indx >= mag.Length)
  790. {
  791. --bitLength;
  792. break;
  793. }
  794. }
  795. while (mag[indx] == 0);
  796. }
  797. return bitLength;
  798. }
  799. public int BitLength
  800. {
  801. get
  802. {
  803. if (nBitLength == -1)
  804. {
  805. nBitLength = sign == 0
  806. ? 0
  807. : CalcBitLength(sign, 0, magnitude);
  808. }
  809. return nBitLength;
  810. }
  811. }
  812. //
  813. // BitLen(value) is the number of bits in value.
  814. //
  815. internal static int BitLen(int w)
  816. {
  817. uint v = (uint)w;
  818. uint t = v >> 24;
  819. if (t != 0)
  820. return 24 + BitLengthTable[t];
  821. t = v >> 16;
  822. if (t != 0)
  823. return 16 + BitLengthTable[t];
  824. t = v >> 8;
  825. if (t != 0)
  826. return 8 + BitLengthTable[t];
  827. return BitLengthTable[v];
  828. }
  829. private bool QuickPow2Check()
  830. {
  831. return sign > 0 && nBits == 1;
  832. }
  833. public int CompareTo(
  834. object obj)
  835. {
  836. return CompareTo((BigInteger)obj);
  837. }
  838. /**
  839. * unsigned comparison on two arrays - note the arrays may
  840. * start with leading zeros.
  841. */
  842. private static int CompareTo(
  843. int xIndx,
  844. int[] x,
  845. int yIndx,
  846. int[] y)
  847. {
  848. while (xIndx != x.Length && x[xIndx] == 0)
  849. {
  850. xIndx++;
  851. }
  852. while (yIndx != y.Length && y[yIndx] == 0)
  853. {
  854. yIndx++;
  855. }
  856. return CompareNoLeadingZeroes(xIndx, x, yIndx, y);
  857. }
  858. private static int CompareNoLeadingZeroes(
  859. int xIndx,
  860. int[] x,
  861. int yIndx,
  862. int[] y)
  863. {
  864. int diff = (x.Length - y.Length) - (xIndx - yIndx);
  865. if (diff != 0)
  866. {
  867. return diff < 0 ? -1 : 1;
  868. }
  869. // lengths of magnitudes the same, test the magnitude values
  870. while (xIndx < x.Length)
  871. {
  872. uint v1 = (uint)x[xIndx++];
  873. uint v2 = (uint)y[yIndx++];
  874. if (v1 != v2)
  875. return v1 < v2 ? -1 : 1;
  876. }
  877. return 0;
  878. }
  879. public int CompareTo(
  880. BigInteger value)
  881. {
  882. return sign < value.sign ? -1
  883. : sign > value.sign ? 1
  884. : sign == 0 ? 0
  885. : sign * CompareNoLeadingZeroes(0, magnitude, 0, value.magnitude);
  886. }
  887. /**
  888. * return z = x / y - done in place (z value preserved, x contains the
  889. * remainder)
  890. */
  891. private int[] Divide(
  892. int[] x,
  893. int[] y)
  894. {
  895. int xStart = 0;
  896. while (xStart < x.Length && x[xStart] == 0)
  897. {
  898. ++xStart;
  899. }
  900. int yStart = 0;
  901. while (yStart < y.Length && y[yStart] == 0)
  902. {
  903. ++yStart;
  904. }
  905. Debug.Assert(yStart < y.Length);
  906. int xyCmp = CompareNoLeadingZeroes(xStart, x, yStart, y);
  907. int[] count;
  908. if (xyCmp > 0)
  909. {
  910. int yBitLength = CalcBitLength(1, yStart, y);
  911. int xBitLength = CalcBitLength(1, xStart, x);
  912. int shift = xBitLength - yBitLength;
  913. int[] iCount;
  914. int iCountStart = 0;
  915. int[] c;
  916. int cStart = 0;
  917. int cBitLength = yBitLength;
  918. if (shift > 0)
  919. {
  920. // iCount = ShiftLeft(One.magnitude, shift);
  921. iCount = new int[(shift >> 5) + 1];
  922. iCount[0] = 1 << (shift % 32);
  923. c = ShiftLeft(y, shift);
  924. cBitLength += shift;
  925. }
  926. else
  927. {
  928. iCount = new int[] { 1 };
  929. int len = y.Length - yStart;
  930. c = new int[len];
  931. Array.Copy(y, yStart, c, 0, len);
  932. }
  933. count = new int[iCount.Length];
  934. for (;;)
  935. {
  936. if (cBitLength < xBitLength
  937. || CompareNoLeadingZeroes(xStart, x, cStart, c) >= 0)
  938. {
  939. Subtract(xStart, x, cStart, c);
  940. AddMagnitudes(count, iCount);
  941. while (x[xStart] == 0)
  942. {
  943. if (++xStart == x.Length)
  944. return count;
  945. }
  946. //xBitLength = CalcBitLength(xStart, x);
  947. xBitLength = 32 * (x.Length - xStart - 1) + BitLen(x[xStart]);
  948. if (xBitLength <= yBitLength)
  949. {
  950. if (xBitLength < yBitLength)
  951. return count;
  952. xyCmp = CompareNoLeadingZeroes(xStart, x, yStart, y);
  953. if (xyCmp <= 0)
  954. break;
  955. }
  956. }
  957. shift = cBitLength - xBitLength;
  958. // NB: The case where c[cStart] is 1-bit is harmless
  959. if (shift == 1)
  960. {
  961. uint firstC = (uint) c[cStart] >> 1;
  962. uint firstX = (uint) x[xStart];
  963. if (firstC > firstX)
  964. ++shift;
  965. }
  966. if (shift < 2)
  967. {
  968. ShiftRightOneInPlace(cStart, c);
  969. --cBitLength;
  970. ShiftRightOneInPlace(iCountStart, iCount);
  971. }
  972. else
  973. {
  974. ShiftRightInPlace(cStart, c, shift);
  975. cBitLength -= shift;
  976. ShiftRightInPlace(iCountStart, iCount, shift);
  977. }
  978. //cStart = c.Length - ((cBitLength + 31) / 32);
  979. while (c[cStart] == 0)
  980. {
  981. ++cStart;
  982. }
  983. while (iCount[iCountStart] == 0)
  984. {
  985. ++iCountStart;
  986. }
  987. }
  988. }
  989. else
  990. {
  991. count = new int[1];
  992. }
  993. if (xyCmp == 0)
  994. {
  995. AddMagnitudes(count, One.magnitude);
  996. Array.Clear(x, xStart, x.Length - xStart);
  997. }
  998. return count;
  999. }
  1000. public BigInteger Divide(
  1001. BigInteger val)
  1002. {
  1003. if (val.sign == 0)
  1004. throw new ArithmeticException("Division by zero error");
  1005. if (sign == 0)
  1006. return Zero;
  1007. if (val.QuickPow2Check()) // val is power of two
  1008. {
  1009. BigInteger result = this.Abs().ShiftRight(val.Abs().BitLength - 1);
  1010. return val.sign == this.sign ? result : result.Negate();
  1011. }
  1012. int[] mag = (int[]) this.magnitude.Clone();
  1013. return new BigInteger(this.sign * val.sign, Divide(mag, val.magnitude), true);
  1014. }
  1015. public BigInteger[] DivideAndRemainder(
  1016. BigInteger val)
  1017. {
  1018. if (val.sign == 0)
  1019. throw new ArithmeticException("Division by zero error");
  1020. BigInteger[] biggies = new BigInteger[2];
  1021. if (sign == 0)
  1022. {
  1023. biggies[0] = Zero;
  1024. biggies[1] = Zero;
  1025. }
  1026. else if (val.QuickPow2Check()) // val is power of two
  1027. {
  1028. int e = val.Abs().BitLength - 1;
  1029. BigInteger quotient = this.Abs().ShiftRight(e);
  1030. int[] remainder = this.LastNBits(e);
  1031. biggies[0] = val.sign == this.sign ? quotient : quotient.Negate();
  1032. biggies[1] = new BigInteger(this.sign, remainder, true);
  1033. }
  1034. else
  1035. {
  1036. int[] remainder = (int[]) this.magnitude.Clone();
  1037. int[] quotient = Divide(remainder, val.magnitude);
  1038. biggies[0] = new BigInteger(this.sign * val.sign, quotient, true);
  1039. biggies[1] = new BigInteger(this.sign, remainder, true);
  1040. }
  1041. return biggies;
  1042. }
  1043. public override bool Equals(
  1044. object obj)
  1045. {
  1046. if (obj == this)
  1047. return true;
  1048. BigInteger biggie = obj as BigInteger;
  1049. if (biggie == null)
  1050. return false;
  1051. return sign == biggie.sign && IsEqualMagnitude(biggie);
  1052. }
  1053. private bool IsEqualMagnitude(BigInteger x)
  1054. {
  1055. int[] xMag = x.magnitude;
  1056. if (magnitude.Length != x.magnitude.Length)
  1057. return false;
  1058. for (int i = 0; i < magnitude.Length; i++)
  1059. {
  1060. if (magnitude[i] != x.magnitude[i])
  1061. return false;
  1062. }
  1063. return true;
  1064. }
  1065. public BigInteger Gcd(
  1066. BigInteger value)
  1067. {
  1068. if (value.sign == 0)
  1069. return Abs();
  1070. if (sign == 0)
  1071. return value.Abs();
  1072. BigInteger r;
  1073. BigInteger u = this;
  1074. BigInteger v = value;
  1075. while (v.sign != 0)
  1076. {
  1077. r = u.Mod(v);
  1078. u = v;
  1079. v = r;
  1080. }
  1081. return u;
  1082. }
  1083. public override int GetHashCode()
  1084. {
  1085. int hc = magnitude.Length;
  1086. if (magnitude.Length > 0)
  1087. {
  1088. hc ^= magnitude[0];
  1089. if (magnitude.Length > 1)
  1090. {
  1091. hc ^= magnitude[magnitude.Length - 1];
  1092. }
  1093. }
  1094. return sign < 0 ? ~hc : hc;
  1095. }
  1096. // TODO Make public?
  1097. private BigInteger Inc()
  1098. {
  1099. if (this.sign == 0)
  1100. return One;
  1101. if (this.sign < 0)
  1102. return new BigInteger(-1, doSubBigLil(this.magnitude, One.magnitude), true);
  1103. return AddToMagnitude(One.magnitude);
  1104. }
  1105. public int IntValue
  1106. {
  1107. get
  1108. {
  1109. if (sign == 0)
  1110. return 0;
  1111. int n = magnitude.Length;
  1112. int v = magnitude[n - 1];
  1113. return sign < 0 ? -v : v;
  1114. }
  1115. }
  1116. /**
  1117. * return whether or not a BigInteger is probably prime with a
  1118. * probability of 1 - (1/2)**certainty.
  1119. * <p>From Knuth Vol 2, pg 395.</p>
  1120. */
  1121. public bool IsProbablePrime(int certainty)
  1122. {
  1123. return IsProbablePrime(certainty, false);
  1124. }
  1125. internal bool IsProbablePrime(int certainty, bool randomlySelected)
  1126. {
  1127. if (certainty <= 0)
  1128. return true;
  1129. BigInteger n = Abs();
  1130. if (!n.TestBit(0))
  1131. return n.Equals(Two);
  1132. if (n.Equals(One))
  1133. return false;
  1134. return n.CheckProbablePrime(certainty, RandomSource, randomlySelected);
  1135. }
  1136. private bool CheckProbablePrime(int certainty, Random random, bool randomlySelected)
  1137. {
  1138. Debug.Assert(certainty > 0);
  1139. Debug.Assert(CompareTo(Two) > 0);
  1140. Debug.Assert(TestBit(0));
  1141. // Try to reduce the penalty for really small numbers
  1142. int numLists = System.Math.Min(BitLength - 1, primeLists.Length);
  1143. for (int i = 0; i < numLists; ++i)
  1144. {
  1145. int test = Remainder(primeProducts[i]);
  1146. int[] primeList = primeLists[i];
  1147. for (int j = 0; j < primeList.Length; ++j)
  1148. {
  1149. int prime = primeList[j];
  1150. int qRem = test % prime;
  1151. if (qRem == 0)
  1152. {
  1153. // We may find small numbers in the list
  1154. return BitLength < 16 && IntValue == prime;
  1155. }
  1156. }
  1157. }
  1158. // TODO Special case for < 10^16 (RabinMiller fixed list)
  1159. // if (BitLength < 30)
  1160. // {
  1161. // RabinMiller against 2, 3, 5, 7, 11, 13, 23 is sufficient
  1162. // }
  1163. // TODO Is it worth trying to create a hybrid of these two?
  1164. return RabinMillerTest(certainty, random, randomlySelected);
  1165. // return SolovayStrassenTest(certainty, random);
  1166. // bool rbTest = RabinMillerTest(certainty, random);
  1167. // bool ssTest = SolovayStrassenTest(certainty, random);
  1168. //
  1169. // Debug.Assert(rbTest == ssTest);
  1170. //
  1171. // return rbTest;
  1172. }
  1173. public bool RabinMillerTest(int certainty, Random random)
  1174. {
  1175. return RabinMillerTest(certainty, random, false);
  1176. }
  1177. internal bool RabinMillerTest(int certainty, Random random, bool randomlySelected)
  1178. {
  1179. int bits = BitLength;
  1180. Debug.Assert(certainty > 0);
  1181. Debug.Assert(bits > 2);
  1182. Debug.Assert(TestBit(0));
  1183. int iterations = ((certainty - 1) / 2) + 1;
  1184. if (randomlySelected)
  1185. {
  1186. int itersFor100Cert = bits >= 1024 ? 4
  1187. : bits >= 512 ? 8
  1188. : bits >= 256 ? 16
  1189. : 50;
  1190. if (certainty < 100)
  1191. {
  1192. iterations = System.Math.Min(itersFor100Cert, iterations);
  1193. }
  1194. else
  1195. {
  1196. iterations -= 50;
  1197. iterations += itersFor100Cert;
  1198. }
  1199. }
  1200. // let n = 1 + d . 2^s
  1201. BigInteger n = this;
  1202. int s = n.GetLowestSetBitMaskFirst(-1 << 1);
  1203. Debug.Assert(s >= 1);
  1204. BigInteger r = n.ShiftRight(s);
  1205. // NOTE: Avoid conversion to/from Montgomery form and check for R/-R as result instead
  1206. BigInteger montRadix = One.ShiftLeft(32 * n.magnitude.Length).Remainder(n);
  1207. BigInteger minusMontRadix = n.Subtract(montRadix);
  1208. do
  1209. {
  1210. BigInteger a;
  1211. do
  1212. {
  1213. a = new BigInteger(n.BitLength, random);
  1214. }
  1215. while (a.sign == 0 || a.CompareTo(n) >= 0
  1216. || a.IsEqualMagnitude(montRadix) || a.IsEqualMagnitude(minusMontRadix));
  1217. BigInteger y = ModPowMonty(a, r, n, false);
  1218. if (!y.Equals(montRadix))
  1219. {
  1220. int j = 0;
  1221. while (!y.Equals(minusMontRadix))
  1222. {
  1223. if (++j == s)
  1224. return false;
  1225. y = ModPowMonty(y, Two, n, false);
  1226. if (y.Equals(montRadix))
  1227. return false;
  1228. }
  1229. }
  1230. }
  1231. while (--iterations > 0);
  1232. return true;
  1233. }
  1234. // private bool SolovayStrassenTest(
  1235. // int certainty,
  1236. // Random random)
  1237. // {
  1238. // Debug.Assert(certainty > 0);
  1239. // Debug.Assert(CompareTo(Two) > 0);
  1240. // Debug.Assert(TestBit(0));
  1241. //
  1242. // BigInteger n = this;
  1243. // BigInteger nMinusOne = n.Subtract(One);
  1244. // BigInteger e = nMinusOne.ShiftRight(1);
  1245. //
  1246. // do
  1247. // {
  1248. // BigInteger a;
  1249. // do
  1250. // {
  1251. // a = new BigInteger(nBitLength, random);
  1252. // }
  1253. // // NB: Spec says 0 < x < n, but 1 is trivial
  1254. // while (a.CompareTo(One) <= 0 || a.CompareTo(n) >= 0);
  1255. //
  1256. //
  1257. // // TODO Check this is redundant given the way Jacobi() works?
  1258. //// if (!a.Gcd(n).Equals(One))
  1259. //// return false;
  1260. //
  1261. // int x = Jacobi(a, n);
  1262. //
  1263. // if (x == 0)
  1264. // return false;
  1265. //
  1266. // BigInteger check = a.ModPow(e, n);
  1267. //
  1268. // if (x == 1 && !check.Equals(One))
  1269. // return false;
  1270. //
  1271. // if (x == -1 && !check.Equals(nMinusOne))
  1272. // return false;
  1273. //
  1274. // --certainty;
  1275. // }
  1276. // while (certainty > 0);
  1277. //
  1278. // return true;
  1279. // }
  1280. //
  1281. // private static int Jacobi(
  1282. // BigInteger a,
  1283. // BigInteger b)
  1284. // {
  1285. // Debug.Assert(a.sign >= 0);
  1286. // Debug.Assert(b.sign > 0);
  1287. // Debug.Assert(b.TestBit(0));
  1288. // Debug.Assert(a.CompareTo(b) < 0);
  1289. //
  1290. // int totalS = 1;
  1291. // for (;;)
  1292. // {
  1293. // if (a.sign == 0)
  1294. // return 0;
  1295. //
  1296. // if (a.Equals(One))
  1297. // break;
  1298. //
  1299. // int e = a.GetLowestSetBit();
  1300. //
  1301. // int bLsw = b.magnitude[b.magnitude.Length - 1];
  1302. // if ((e & 1) != 0 && ((bLsw & 7) == 3 || (bLsw & 7) == 5))
  1303. // totalS = -totalS;
  1304. //
  1305. // // TODO Confirm this is faster than later a1.Equals(One) test
  1306. // if (a.BitLength == e + 1)
  1307. // break;
  1308. // BigInteger a1 = a.ShiftRight(e);
  1309. //// if (a1.Equals(One))
  1310. //// break;
  1311. //
  1312. // int a1Lsw = a1.magnitude[a1.magnitude.Length - 1];
  1313. // if ((bLsw & 3) == 3 && (a1Lsw & 3) == 3)
  1314. // totalS = -totalS;
  1315. //
  1316. //// a = b.Mod(a1);
  1317. // a = b.Remainder(a1);
  1318. // b = a1;
  1319. // }
  1320. // return totalS;
  1321. // }
  1322. public long LongValue
  1323. {
  1324. get
  1325. {
  1326. if (sign == 0)
  1327. return 0;
  1328. int n = magnitude.Length;
  1329. long v = magnitude[n - 1] & IMASK;
  1330. if (n > 1)
  1331. {
  1332. v |= (magnitude[n - 2] & IMASK) << 32;
  1333. }
  1334. return sign < 0 ? -v : v;
  1335. }
  1336. }
  1337. public BigInteger Max(
  1338. BigInteger value)
  1339. {
  1340. return CompareTo(value) > 0 ? this : value;
  1341. }
  1342. public BigInteger Min(
  1343. BigInteger value)
  1344. {
  1345. return CompareTo(value) < 0 ? this : value;
  1346. }
  1347. public BigInteger Mod(
  1348. BigInteger m)
  1349. {
  1350. if (m.sign < 1)
  1351. throw new ArithmeticException("Modulus must be positive");
  1352. BigInteger biggie = Remainder(m);
  1353. return (biggie.sign >= 0 ? biggie : biggie.Add(m));
  1354. }
  1355. public BigInteger ModInverse(
  1356. BigInteger m)
  1357. {
  1358. if (m.sign < 1)
  1359. throw new ArithmeticException("Modulus must be positive");
  1360. // TODO Too slow at the moment
  1361. // // "Fast Key Exchange with Elliptic Curve Systems" R.Schoeppel
  1362. // if (m.TestBit(0))
  1363. // {
  1364. // //The Almost Inverse Algorithm
  1365. // int k = 0;
  1366. // BigInteger B = One, C = Zero, F = this, G = m, tmp;
  1367. //
  1368. // for (;;)
  1369. // {
  1370. // // While F is even, do F=F/u, C=C*u, k=k+1.
  1371. // int zeroes = F.GetLowestSetBit();
  1372. // if (zeroes > 0)
  1373. // {
  1374. // F = F.ShiftRight(zeroes);
  1375. // C = C.ShiftLeft(zeroes);
  1376. // k += zeroes;
  1377. // }
  1378. //
  1379. // // If F = 1, then return B,k.
  1380. // if (F.Equals(One))
  1381. // {
  1382. // BigInteger half = m.Add(One).ShiftRight(1);
  1383. // BigInteger halfK = half.ModPow(BigInteger.ValueOf(k), m);
  1384. // return B.Multiply(halfK).Mod(m);
  1385. // }
  1386. //
  1387. // if (F.CompareTo(G) < 0)
  1388. // {
  1389. // tmp = G; G = F; F = tmp;
  1390. // tmp = B; B = C; C = tmp;
  1391. // }
  1392. //
  1393. // F = F.Add(G);
  1394. // B = B.Add(C);
  1395. // }
  1396. // }
  1397. if (m.QuickPow2Check())
  1398. {
  1399. return ModInversePow2(m);
  1400. }
  1401. BigInteger d = this.Remainder(m);
  1402. BigInteger x;
  1403. BigInteger gcd = ExtEuclid(d, m, out x);
  1404. if (!gcd.Equals(One))
  1405. throw new ArithmeticException("Numbers not relatively prime.");
  1406. if (x.sign < 0)
  1407. {
  1408. x = x.Add(m);
  1409. }
  1410. return x;
  1411. }
  1412. private BigInteger ModInversePow2(BigInteger m)
  1413. {
  1414. Debug.Assert(m.SignValue > 0);
  1415. Debug.Assert(m.BitCount == 1);
  1416. if (!TestBit(0))
  1417. {
  1418. throw new ArithmeticException("Numbers not relatively prime.");
  1419. }
  1420. int pow = m.BitLength - 1;
  1421. long inv64 = ModInverse64(LongValue);
  1422. if (pow < 64)
  1423. {
  1424. inv64 &= ((1L << pow) - 1);
  1425. }
  1426. BigInteger x = BigInteger.ValueOf(inv64);
  1427. if (pow > 64)
  1428. {
  1429. BigInteger d = this.Remainder(m);
  1430. int bitsCorrect = 64;
  1431. do
  1432. {
  1433. BigInteger t = x.Multiply(d).Remainder(m);
  1434. x = x.Multiply(Two.Subtract(t)).Remainder(m);
  1435. bitsCorrect <<= 1;
  1436. }
  1437. while (bitsCorrect < pow);
  1438. }
  1439. if (x.sign < 0)
  1440. {
  1441. x = x.Add(m);
  1442. }
  1443. return x;
  1444. }
  1445. private static int ModInverse32(int d)
  1446. {
  1447. // Newton's method with initial estimate "correct to 4 bits"
  1448. Debug.Assert((d & 1) != 0);
  1449. int x = d + (((d + 1) & 4) << 1); // d.x == 1 mod 2**4
  1450. Debug.Assert(((d * x) & 15) == 1);
  1451. x *= 2 - d * x; // d.x == 1 mod 2**8
  1452. x *= 2 - d * x; // d.x == 1 mod 2**16
  1453. x *= 2 - d * x; // d.x == 1 mod 2**32
  1454. Debug.Assert(d * x == 1);
  1455. return x;
  1456. }
  1457. private static long ModInverse64(long d)
  1458. {
  1459. // Newton's method with initial estimate "correct to 4 bits"
  1460. Debug.Assert((d & 1L) != 0);
  1461. long x = d + (((d + 1L) & 4L) << 1); // d.x == 1 mod 2**4
  1462. Debug.Assert(((d * x) & 15L) == 1L);
  1463. x *= 2 - d * x; // d.x == 1 mod 2**8
  1464. x *= 2 - d * x; // d.x == 1 mod 2**16
  1465. x *= 2 - d * x; // d.x == 1 mod 2**32
  1466. x *= 2 - d * x; // d.x == 1 mod 2**64
  1467. Debug.Assert(d * x == 1L);
  1468. return x;
  1469. }
  1470. /**
  1471. * Calculate the numbers u1, u2, and u3 such that:
  1472. *
  1473. * u1 * a + u2 * b = u3
  1474. *
  1475. * where u3 is the greatest common divider of a and b.
  1476. * a and b using the extended Euclid algorithm (refer p. 323
  1477. * of The Art of Computer Programming vol 2, 2nd ed).
  1478. * This also seems to have the side effect of calculating
  1479. * some form of multiplicative inverse.
  1480. *
  1481. * @param a First number to calculate gcd for
  1482. * @param b Second number to calculate gcd for
  1483. * @param u1Out the return object for the u1 value
  1484. * @return The greatest common divisor of a and b
  1485. */
  1486. private static BigInteger ExtEuclid(BigInteger a, BigInteger b, out BigInteger u1Out)
  1487. {
  1488. BigInteger u1 = One, v1 = Zero;
  1489. BigInteger u3 = a, v3 = b;
  1490. if (v3.sign > 0)
  1491. {
  1492. for (;;)
  1493. {
  1494. BigInteger[] q = u3.DivideAndRemainder(v3);
  1495. u3 = v3;
  1496. v3 = q[1];
  1497. BigInteger oldU1 = u1;
  1498. u1 = v1;
  1499. if (v3.sign <= 0)
  1500. break;
  1501. v1 = oldU1.Subtract(v1.Multiply(q[0]));
  1502. }
  1503. }
  1504. u1Out = u1;
  1505. return u3;
  1506. }
  1507. private static void ZeroOut(
  1508. int[] x)
  1509. {
  1510. Array.Clear(x, 0, x.Length);
  1511. }
  1512. public BigInteger ModPow(BigInteger e, BigInteger m)
  1513. {
  1514. if (m.sign < 1)
  1515. throw new ArithmeticException("Modulus must be positive");
  1516. if (m.Equals(One))
  1517. return Zero;
  1518. if (e.sign == 0)
  1519. return One;
  1520. if (sign == 0)
  1521. return Zero;
  1522. bool negExp = e.sign < 0;
  1523. if (negExp)
  1524. e = e.Negate();
  1525. BigInteger result = this.Mod(m);
  1526. if (!e.Equals(One))
  1527. {
  1528. if ((m.magnitude[m.magnitude.Length - 1] & 1) == 0)
  1529. {
  1530. result = ModPowBarrett(result, e, m);
  1531. }
  1532. else
  1533. {
  1534. result = ModPowMonty(result, e, m, true);
  1535. }
  1536. }
  1537. if (negExp)
  1538. result = result.ModInverse(m);
  1539. return result;
  1540. }
  1541. private static BigInteger ModPowBarrett(BigInteger b, BigInteger e, BigInteger m)
  1542. {
  1543. int k = m.magnitude.Length;
  1544. BigInteger mr = One.ShiftLeft((k + 1) << 5);
  1545. BigInteger yu = One.ShiftLeft(k << 6).Divide(m);
  1546. // Sliding window from MSW to LSW
  1547. int extraBits = 0, expLength = e.BitLength;
  1548. while (expLength > ExpWindowThresholds[extraBits])
  1549. {
  1550. ++extraBits;
  1551. }
  1552. int numPowers = 1 << extraBits;
  1553. BigInteger[] oddPowers = new BigInteger[numPowers];
  1554. oddPowers[0] = b;
  1555. BigInteger b2 = ReduceBarrett(b.Square(), m, mr, yu);
  1556. for (int i = 1; i < numPowers; ++i)
  1557. {
  1558. oddPowers[i] = ReduceBarrett(oddPowers[i - 1].Multiply(b2), m, mr, yu);
  1559. }
  1560. int[] windowList = GetWindowList(e.magnitude, extraBits);
  1561. Debug.Assert(windowList.Length > 0);
  1562. int window = windowList[0];
  1563. int mult = window & 0xFF, lastZeroes = window >> 8;
  1564. BigInteger y;
  1565. if (mult == 1)
  1566. {
  1567. y = b2;
  1568. --lastZeroes;
  1569. }
  1570. else
  1571. {
  1572. y = oddPowers[mult >> 1];
  1573. }
  1574. int windowPos = 1;
  1575. while ((window = windowList[windowPos++]) != -1)
  1576. {
  1577. mult = window & 0xFF;
  1578. int bits = lastZeroes + BitLengthTable[mult];
  1579. for (int j = 0; j < bits; ++j)
  1580. {
  1581. y = ReduceBarrett(y.Square(), m, mr, yu);
  1582. }
  1583. y = ReduceBarrett(y.Multiply(oddPowers[mult >> 1]), m, mr, yu);
  1584. lastZeroes = window >> 8;
  1585. }
  1586. for (int i = 0; i < lastZeroes; ++i)
  1587. {
  1588. y = ReduceBarrett(y.Square(), m, mr, yu);
  1589. }
  1590. return y;
  1591. }
  1592. private static BigInteger ReduceBarrett(BigInteger x, BigInteger m, BigInteger mr, BigInteger yu)
  1593. {
  1594. int xLen = x.BitLength, mLen = m.BitLength;
  1595. if (xLen < mLen)
  1596. return x;
  1597. if (xLen - mLen > 1)
  1598. {
  1599. int k = m.magnitude.Length;
  1600. BigInteger q1 = x.DivideWords(k - 1);
  1601. BigInteger q2 = q1.Multiply(yu); // TODO Only need partial multiplication here
  1602. BigInteger q3 = q2.DivideWords(k + 1);
  1603. BigInteger r1 = x.RemainderWords(k + 1);
  1604. BigInteger r2 = q3.Multiply(m); // TODO Only need partial multiplication here
  1605. BigInteger r3 = r2.RemainderWords(k + 1);
  1606. x = r1.Subtract(r3);
  1607. if (x.sign < 0)
  1608. {
  1609. x = x.Add(mr);
  1610. }
  1611. }
  1612. while (x.CompareTo(m) >= 0)
  1613. {
  1614. x = x.Subtract(m);
  1615. }
  1616. return x;
  1617. }
  1618. private static BigInteger ModPowMonty(BigInteger b, BigInteger e, BigInteger m, bool convert)
  1619. {
  1620. int n = m.magnitude.Length;
  1621. int powR = 32 * n;
  1622. bool smallMontyModulus = m.BitLength + 2 <= powR;
  1623. uint mDash = (uint)m.GetMQuote();
  1624. // tmp = this * R mod m
  1625. if (convert)
  1626. {
  1627. b = b.ShiftLeft(powR).Remainder(m);
  1628. }
  1629. int[] yAccum = new int[n + 1];
  1630. int[] zVal = b.magnitude;
  1631. Debug.Assert(zVal.Length <= n);
  1632. if (zVal.Length < n)
  1633. {
  1634. int[] tmp = new int[n];
  1635. zVal.CopyTo(tmp, n - zVal.Length);
  1636. zVal = tmp;
  1637. }
  1638. // Sliding window from MSW to LSW
  1639. int extraBits = 0;
  1640. // Filter the common case of small RSA exponents with few bits set
  1641. if (e.magnitude.Length > 1 || e.BitCount > 2)
  1642. {
  1643. int expLength = e.BitLength;
  1644. while (expLength > ExpWindowThresholds[extraBits])
  1645. {
  1646. ++extraBits;
  1647. }
  1648. }
  1649. int numPowers = 1 << extraBits;
  1650. int[][] oddPowers = new int[numPowers][];
  1651. oddPowers[0] = zVal;
  1652. int[] zSquared = Arrays.Clone(zVal);
  1653. SquareMonty(yAccum, zSquared, m.magnitude, mDash, smallMontyModulus);
  1654. for (int i = 1; i < numPowers; ++i)
  1655. {
  1656. oddPowers[i] = Arrays.Clone(oddPowers[i - 1]);
  1657. MultiplyMonty(yAccum, oddPowers[i], zSquared, m.magnitude, mDash, smallMontyModulus);
  1658. }
  1659. int[] windowList = GetWindowList(e.magnitude, extraBits);
  1660. Debug.Assert(windowList.Length > 1);
  1661. int window = windowList[0];
  1662. int mult = window & 0xFF, lastZeroes = window >> 8;
  1663. int[] yVal;
  1664. if (mult == 1)
  1665. {
  1666. yVal = zSquared;
  1667. --lastZeroes;
  1668. }
  1669. else
  1670. {
  1671. yVal = Arrays.Clone(oddPowers[mult >> 1]);
  1672. }
  1673. int windowPos = 1;
  1674. while ((window = windowList[windowPos++]) != -1)
  1675. {
  1676. mult = window & 0xFF;
  1677. int bits = lastZeroes + BitLengthTable[mult];
  1678. for (int j = 0; j < bits; ++j)
  1679. {
  1680. SquareMonty(yAccum, yVal, m.magnitude, mDash, smallMontyModulus);
  1681. }
  1682. MultiplyMonty(yAccum, yVal, oddPowers[mult >> 1], m.magnitude, mDash, smallMontyModulus);
  1683. lastZeroes = window >> 8;
  1684. }
  1685. for (int i = 0; i < lastZeroes; ++i)
  1686. {
  1687. SquareMonty(yAccum, yVal, m.magnitude, mDash, smallMontyModulus);
  1688. }
  1689. if (convert)
  1690. {
  1691. // Return y * R^(-1) mod m
  1692. MontgomeryReduce(yVal, m.magnitude, mDash);
  1693. }
  1694. else if (smallMontyModulus && CompareTo(0, yVal, 0, m.magnitude) >= 0)
  1695. {
  1696. Subtract(0, yVal, 0, m.magnitude);
  1697. }
  1698. return new BigInteger(1, yVal, true);
  1699. }
  1700. private static int[] GetWindowList(int[] mag, int extraBits)
  1701. {
  1702. int v = mag[0];
  1703. Debug.Assert(v != 0);
  1704. int leadingBits = BitLen(v);
  1705. int resultSize = (((mag.Length - 1) << 5) + leadingBits) / (1 + extraBits) + 2;
  1706. int[] result = new int[resultSize];
  1707. int resultPos = 0;
  1708. int bitPos = 33 - leadingBits;
  1709. v <<= bitPos;
  1710. int mult = 1, multLimit = 1 << extraBits;
  1711. int zeroes = 0;
  1712. int i = 0;
  1713. for (; ; )
  1714. {
  1715. for (; bitPos < 32; ++bitPos)
  1716. {
  1717. if (mult < multLimit)
  1718. {
  1719. mult = (mult << 1) | (int)((uint)v >> 31);
  1720. }
  1721. else if (v < 0)
  1722. {
  1723. result[resultPos++] = CreateWindowEntry(mult, zeroes);
  1724. mult = 1;
  1725. zeroes = 0;
  1726. }
  1727. else
  1728. {
  1729. ++zeroes;
  1730. }
  1731. v <<= 1;
  1732. }
  1733. if (++i == mag.Length)
  1734. {
  1735. result[resultPos++] = CreateWindowEntry(mult, zeroes);
  1736. break;
  1737. }
  1738. v = mag[i];
  1739. bitPos = 0;
  1740. }
  1741. result[resultPos] = -1;
  1742. return result;
  1743. }
  1744. private static int CreateWindowEntry(int mult, int zeroes)
  1745. {
  1746. while ((mult & 1) == 0)
  1747. {
  1748. mult >>= 1;
  1749. ++zeroes;
  1750. }
  1751. return mult | (zeroes << 8);
  1752. }
  1753. /**
  1754. * return w with w = x * x - w is assumed to have enough space.
  1755. */
  1756. private static int[] Square(
  1757. int[] w,
  1758. int[] x)
  1759. {
  1760. // Note: this method allows w to be only (2 * x.Length - 1) words if result will fit
  1761. // if (w.Length != 2 * x.Length)
  1762. // throw new ArgumentException("no I don't think so...");
  1763. ulong c;
  1764. int wBase = w.Length - 1;
  1765. for (int i = x.Length - 1; i > 0; --i)
  1766. {
  1767. ulong v = (uint)x[i];
  1768. c = v * v + (uint)w[wBase];
  1769. w[wBase] = (int)c;
  1770. c >>= 32;
  1771. for (int j = i - 1; j >= 0; --j)
  1772. {
  1773. ulong prod = v * (uint)x[j];
  1774. c += ((uint)w[--wBase] & UIMASK) + ((uint)prod << 1);
  1775. w[wBase] = (int)c;
  1776. c = (c >> 32) + (prod >> 31);
  1777. }
  1778. c += (uint)w[--wBase];
  1779. w[wBase] = (int)c;
  1780. if (--wBase >= 0)
  1781. {
  1782. w[wBase] = (int)(c >> 32);
  1783. }
  1784. else
  1785. {
  1786. Debug.Assert((c >> 32) == 0);
  1787. }
  1788. wBase += i;
  1789. }
  1790. c = (uint)x[0];
  1791. c = c * c + (uint)w[wBase];
  1792. w[wBase] = (int)c;
  1793. if (--wBase >= 0)
  1794. {
  1795. w[wBase] += (int)(c >> 32);
  1796. }
  1797. else
  1798. {
  1799. Debug.Assert((c >> 32) == 0);
  1800. }
  1801. return w;
  1802. }
  1803. /**
  1804. * return x with x = y * z - x is assumed to have enough space.
  1805. */
  1806. private static int[] Multiply(int[] x, int[] y, int[] z)
  1807. {
  1808. int i = z.Length;
  1809. if (i < 1)
  1810. return x;
  1811. int xBase = x.Length - y.Length;
  1812. do
  1813. {
  1814. long a = z[--i] & IMASK;
  1815. long val = 0;
  1816. if (a != 0)
  1817. {
  1818. for (int j = y.Length - 1; j >= 0; j--)
  1819. {
  1820. val += a * (y[j] & IMASK) + (x[xBase + j] & IMASK);
  1821. x[xBase + j] = (int)val;
  1822. val = (long)((ulong)val >> 32);
  1823. }
  1824. }
  1825. --xBase;
  1826. if (xBase >= 0)
  1827. {
  1828. x[xBase] = (int)val;
  1829. }
  1830. else
  1831. {
  1832. Debug.Assert(val == 0);
  1833. }
  1834. }
  1835. while (i > 0);
  1836. return x;
  1837. }
  1838. /**
  1839. * Calculate mQuote = -m^(-1) mod b with b = 2^32 (32 = word size)
  1840. */
  1841. private int GetMQuote()
  1842. {
  1843. if (mQuote != 0)
  1844. {
  1845. return mQuote; // already calculated
  1846. }
  1847. Debug.Assert(this.sign > 0);
  1848. int d = -magnitude[magnitude.Length - 1];
  1849. Debug.Assert((d & 1) != 0);
  1850. return mQuote = ModInverse32(d);
  1851. }
  1852. private static void MontgomeryReduce(int[] x, int[] m, uint mDash) // mDash = -m^(-1) mod b
  1853. {
  1854. // NOTE: Not a general purpose reduction (which would allow x up to twice the bitlength of m)
  1855. Debug.Assert(x.Length == m.Length);
  1856. int n = m.Length;
  1857. for (int i = n - 1; i >= 0; --i)
  1858. {
  1859. uint x0 = (uint)x[n - 1];
  1860. ulong t = x0 * mDash;
  1861. ulong carry = t * (uint)m[n - 1] + x0;
  1862. Debug.Assert((uint)carry == 0);
  1863. carry >>= 32;
  1864. for (int j = n - 2; j >= 0; --j)
  1865. {
  1866. carry += t * (uint)m[j] + (uint)x[j];
  1867. x[j + 1] = (int)carry;
  1868. carry >>= 32;
  1869. }
  1870. x[0] = (int)carry;
  1871. Debug.Assert(carry >> 32 == 0);
  1872. }
  1873. if (CompareTo(0, x, 0, m) >= 0)
  1874. {
  1875. Subtract(0, x, 0, m);
  1876. }
  1877. }
  1878. /**
  1879. * Montgomery multiplication: a = x * y * R^(-1) mod m
  1880. * <br/>
  1881. * Based algorithm 14.36 of Handbook of Applied Cryptography.
  1882. * <br/>
  1883. * <li> m, x, y should have length n </li>
  1884. * <li> a should have length (n + 1) </li>
  1885. * <li> b = 2^32, R = b^n </li>
  1886. * <br/>
  1887. * The result is put in x
  1888. * <br/>
  1889. * NOTE: the indices of x, y, m, a different in HAC and in Java
  1890. */
  1891. private static void MultiplyMonty(int[] a, int[] x, int[] y, int[] m, uint mDash, bool smallMontyModulus)
  1892. // mDash = -m^(-1) mod b
  1893. {
  1894. int n = m.Length;
  1895. if (n == 1)
  1896. {
  1897. x[0] = (int)MultiplyMontyNIsOne((uint)x[0], (uint)y[0], (uint)m[0], mDash);
  1898. return;
  1899. }
  1900. uint y0 = (uint)y[n - 1];
  1901. int aMax;
  1902. {
  1903. ulong xi = (uint)x[n - 1];
  1904. ulong carry = xi * y0;
  1905. ulong t = (uint)carry * mDash;
  1906. ulong prod2 = t * (uint)m[n - 1];
  1907. carry += (uint)prod2;
  1908. Debug.Assert((uint)carry == 0);
  1909. carry = (carry >> 32) + (prod2 >> 32);
  1910. for (int j = n - 2; j >= 0; --j)
  1911. {
  1912. ulong prod1 = xi * (uint)y[j];
  1913. prod2 = t * (uint)m[j];
  1914. carry += (prod1 & UIMASK) + (uint)prod2;
  1915. a[j + 2] = (int)carry;
  1916. carry = (carry >> 32) + (prod1 >> 32) + (prod2 >> 32);
  1917. }
  1918. a[1] = (int)carry;
  1919. aMax = (int)(carry >> 32);
  1920. }
  1921. for (int i = n - 2; i >= 0; --i)
  1922. {
  1923. uint a0 = (uint)a[n];
  1924. ulong xi = (uint)x[i];
  1925. ulong prod1 = xi * y0;
  1926. ulong carry = (prod1 & UIMASK) + a0;
  1927. ulong t = (uint)carry * mDash;
  1928. ulong prod2 = t * (uint)m[n - 1];
  1929. carry += (uint)prod2;
  1930. Debug.Assert((uint)carry == 0);
  1931. carry = (carry >> 32) + (prod1 >> 32) + (prod2 >> 32);
  1932. for (int j = n - 2; j >= 0; --j)
  1933. {
  1934. prod1 = xi * (uint)y[j];
  1935. prod2 = t * (uint)m[j];
  1936. carry += (prod1 & UIMASK) + (uint)prod2 + (uint)a[j + 1];
  1937. a[j + 2] = (int)carry;
  1938. carry = (carry >> 32) + (prod1 >> 32) + (prod2 >> 32);
  1939. }
  1940. carry += (uint)aMax;
  1941. a[1] = (int)carry;
  1942. aMax = (int)(carry >> 32);
  1943. }
  1944. a[0] = aMax;
  1945. if (!smallMontyModulus && CompareTo(0, a, 0, m) >= 0)
  1946. {
  1947. Subtract(0, a, 0, m);
  1948. }
  1949. Array.Copy(a, 1, x, 0, n);
  1950. }
  1951. private static void SquareMonty(int[] a, int[] x, int[] m, uint mDash, bool smallMontyModulus)
  1952. // mDash = -m^(-1) mod b
  1953. {
  1954. int n = m.Length;
  1955. if (n == 1)
  1956. {
  1957. uint xVal = (uint)x[0];
  1958. x[0] = (int)MultiplyMontyNIsOne(xVal, xVal, (uint)m[0], mDash);
  1959. return;
  1960. }
  1961. ulong x0 = (uint)x[n - 1];
  1962. int aMax;
  1963. {
  1964. ulong carry = x0 * x0;
  1965. ulong t = (uint)carry * mDash;
  1966. ulong prod2 = t * (uint)m[n - 1];
  1967. carry += (uint)prod2;
  1968. Debug.Assert((uint)carry == 0);
  1969. carry = (carry >> 32) + (prod2 >> 32);
  1970. for (int j = n - 2; j >= 0; --j)
  1971. {
  1972. ulong prod1 = x0 * (uint)x[j];
  1973. prod2 = t * (uint)m[j];
  1974. carry += (prod2 & UIMASK) + ((uint)prod1 << 1);
  1975. a[j + 2] = (int)carry;
  1976. carry = (carry >> 32) + (prod1 >> 31) + (prod2 >> 32);
  1977. }
  1978. a[1] = (int)carry;
  1979. aMax = (int)(carry >> 32);
  1980. }
  1981. for (int i = n - 2; i >= 0; --i)
  1982. {
  1983. uint a0 = (uint)a[n];
  1984. ulong t = a0 * mDash;
  1985. ulong carry = t * (uint)m[n - 1] + a0;
  1986. Debug.Assert((uint)carry == 0);
  1987. carry >>= 32;
  1988. for (int j = n - 2; j > i; --j)
  1989. {
  1990. carry += t * (uint)m[j] + (uint)a[j + 1];
  1991. a[j + 2] = (int)carry;
  1992. carry >>= 32;
  1993. }
  1994. ulong xi = (uint)x[i];
  1995. {
  1996. ulong prod1 = xi * xi;
  1997. ulong prod2 = t * (uint)m[i];
  1998. carry += (prod1 & UIMASK) + (uint)prod2 + (uint)a[i + 1];
  1999. a[i + 2] = (int)carry;
  2000. carry = (carry >> 32) + (prod1 >> 32) + (prod2 >> 32);
  2001. }
  2002. for (int j = i - 1; j >= 0; --j)
  2003. {
  2004. ulong prod1 = xi * (uint)x[j];
  2005. ulong prod2 = t * (uint)m[j];
  2006. carry += (prod2 & UIMASK) + ((uint)prod1 << 1) + (uint)a[j + 1];
  2007. a[j + 2] = (int)carry;
  2008. carry = (carry >> 32) + (prod1 >> 31) + (prod2 >> 32);
  2009. }
  2010. carry += (uint)aMax;
  2011. a[1] = (int)carry;
  2012. aMax = (int)(carry >> 32);
  2013. }
  2014. a[0] = aMax;
  2015. if (!smallMontyModulus && CompareTo(0, a, 0, m) >= 0)
  2016. {
  2017. Subtract(0, a, 0, m);
  2018. }
  2019. Array.Copy(a, 1, x, 0, n);
  2020. }
  2021. private static uint MultiplyMontyNIsOne(uint x, uint y, uint m, uint mDash)
  2022. {
  2023. ulong carry = (ulong)x * y;
  2024. uint t = (uint)carry * mDash;
  2025. ulong um = m;
  2026. ulong prod2 = um * t;
  2027. carry += (uint)prod2;
  2028. Debug.Assert((uint)carry == 0);
  2029. carry = (carry >> 32) + (prod2 >> 32);
  2030. if (carry > um)
  2031. {
  2032. carry -= um;
  2033. }
  2034. Debug.Assert(carry < um);
  2035. return (uint)carry;
  2036. }
  2037. public BigInteger Multiply(
  2038. BigInteger val)
  2039. {
  2040. if (val == this)
  2041. return Square();
  2042. if ((sign & val.sign) == 0)
  2043. return Zero;
  2044. if (val.QuickPow2Check()) // val is power of two
  2045. {
  2046. BigInteger result = this.ShiftLeft(val.Abs().BitLength - 1);
  2047. return val.sign > 0 ? result : result.Negate();
  2048. }
  2049. if (this.QuickPow2Check()) // this is power of two
  2050. {
  2051. BigInteger result = val.ShiftLeft(this.Abs().BitLength - 1);
  2052. return this.sign > 0 ? result : result.Negate();
  2053. }
  2054. int resLength = magnitude.Length + val.magnitude.Length;
  2055. int[] res = new int[resLength];
  2056. Multiply(res, this.magnitude, val.magnitude);
  2057. int resSign = sign ^ val.sign ^ 1;
  2058. return new BigInteger(resSign, res, true);
  2059. }
  2060. public BigInteger Square()
  2061. {
  2062. if (sign == 0)
  2063. return Zero;
  2064. if (this.QuickPow2Check())
  2065. return ShiftLeft(Abs().BitLength - 1);
  2066. int resLength = magnitude.Length << 1;
  2067. if ((uint)magnitude[0] >> 16 == 0)
  2068. --resLength;
  2069. int[] res = new int[resLength];
  2070. Square(res, magnitude);
  2071. return new BigInteger(1, res, false);
  2072. }
  2073. public BigInteger Negate()
  2074. {
  2075. if (sign == 0)
  2076. return this;
  2077. return new BigInteger(-sign, magnitude, false);
  2078. }
  2079. public BigInteger NextProbablePrime()
  2080. {
  2081. if (sign < 0)
  2082. throw new ArithmeticException("Cannot be called on value < 0");
  2083. if (CompareTo(Two) < 0)
  2084. return Two;
  2085. BigInteger n = Inc().SetBit(0);
  2086. while (!n.CheckProbablePrime(100, RandomSource, false))
  2087. {
  2088. n = n.Add(Two);
  2089. }
  2090. return n;
  2091. }
  2092. public BigInteger Not()
  2093. {
  2094. return Inc().Negate();
  2095. }
  2096. public BigInteger Pow(int exp)
  2097. {
  2098. if (exp <= 0)
  2099. {
  2100. if (exp < 0)
  2101. throw new ArithmeticException("Negative exponent");
  2102. return One;
  2103. }
  2104. if (sign == 0)
  2105. {
  2106. return this;
  2107. }
  2108. if (QuickPow2Check())
  2109. {
  2110. long powOf2 = (long)exp * (BitLength - 1);
  2111. if (powOf2 > Int32.MaxValue)
  2112. {
  2113. throw new ArithmeticException("Result too large");
  2114. }
  2115. return One.ShiftLeft((int)powOf2);
  2116. }
  2117. BigInteger y = One;
  2118. BigInteger z = this;
  2119. for (;;)
  2120. {
  2121. if ((exp & 0x1) == 1)
  2122. {
  2123. y = y.Multiply(z);
  2124. }
  2125. exp >>= 1;
  2126. if (exp == 0) break;
  2127. z = z.Multiply(z);
  2128. }
  2129. return y;
  2130. }
  2131. public static BigInteger ProbablePrime(
  2132. int bitLength,
  2133. Random random)
  2134. {
  2135. return new BigInteger(bitLength, 100, random);
  2136. }
  2137. private int Remainder(
  2138. int m)
  2139. {
  2140. Debug.Assert(m > 0);
  2141. long acc = 0;
  2142. for (int pos = 0; pos < magnitude.Length; ++pos)
  2143. {
  2144. long posVal = (uint) magnitude[pos];
  2145. acc = (acc << 32 | posVal) % m;
  2146. }
  2147. return (int) acc;
  2148. }
  2149. /**
  2150. * return x = x % y - done in place (y value preserved)
  2151. */
  2152. private static int[] Remainder(
  2153. int[] x,
  2154. int[] y)
  2155. {
  2156. int xStart = 0;
  2157. while (xStart < x.Length && x[xStart] == 0)
  2158. {
  2159. ++xStart;
  2160. }
  2161. int yStart = 0;
  2162. while (yStart < y.Length && y[yStart] == 0)
  2163. {
  2164. ++yStart;
  2165. }
  2166. Debug.Assert(yStart < y.Length);
  2167. int xyCmp = CompareNoLeadingZeroes(xStart, x, yStart, y);
  2168. if (xyCmp > 0)
  2169. {
  2170. int yBitLength = CalcBitLength(1, yStart, y);
  2171. int xBitLength = CalcBitLength(1, xStart, x);
  2172. int shift = xBitLength - yBitLength;
  2173. int[] c;
  2174. int cStart = 0;
  2175. int cBitLength = yBitLength;
  2176. if (shift > 0)
  2177. {
  2178. c = ShiftLeft(y, shift);
  2179. cBitLength += shift;
  2180. Debug.Assert(c[0] != 0);
  2181. }
  2182. else
  2183. {
  2184. int len = y.Length - yStart;
  2185. c = new int[len];
  2186. Array.Copy(y, yStart, c, 0, len);
  2187. }
  2188. for (;;)
  2189. {
  2190. if (cBitLength < xBitLength
  2191. || CompareNoLeadingZeroes(xStart, x, cStart, c) >= 0)
  2192. {
  2193. Subtract(xStart, x, cStart, c);
  2194. while (x[xStart] == 0)
  2195. {
  2196. if (++xStart == x.Length)
  2197. return x;
  2198. }
  2199. //xBitLength = CalcBitLength(xStart, x);
  2200. xBitLength = 32 * (x.Length - xStart - 1) + BitLen(x[xStart]);
  2201. if (xBitLength <= yBitLength)
  2202. {
  2203. if (xBitLength < yBitLength)
  2204. return x;
  2205. xyCmp = CompareNoLeadingZeroes(xStart, x, yStart, y);
  2206. if (xyCmp <= 0)
  2207. break;
  2208. }
  2209. }
  2210. shift = cBitLength - xBitLength;
  2211. // NB: The case where c[cStart] is 1-bit is harmless
  2212. if (shift == 1)
  2213. {
  2214. uint firstC = (uint) c[cStart] >> 1;
  2215. uint firstX = (uint) x[xStart];
  2216. if (firstC > firstX)
  2217. ++shift;
  2218. }
  2219. if (shift < 2)
  2220. {
  2221. ShiftRightOneInPlace(cStart, c);
  2222. --cBitLength;
  2223. }
  2224. else
  2225. {
  2226. ShiftRightInPlace(cStart, c, shift);
  2227. cBitLength -= shift;
  2228. }
  2229. //cStart = c.Length - ((cBitLength + 31) / 32);
  2230. while (c[cStart] == 0)
  2231. {
  2232. ++cStart;
  2233. }
  2234. }
  2235. }
  2236. if (xyCmp == 0)
  2237. {
  2238. Array.Clear(x, xStart, x.Length - xStart);
  2239. }
  2240. return x;
  2241. }
  2242. public BigInteger Remainder(
  2243. BigInteger n)
  2244. {
  2245. if (n.sign == 0)
  2246. throw new ArithmeticException("Division by zero error");
  2247. if (this.sign == 0)
  2248. return Zero;
  2249. // For small values, use fast remainder method
  2250. if (n.magnitude.Length == 1)
  2251. {
  2252. int val = n.magnitude[0];
  2253. if (val > 0)
  2254. {
  2255. if (val == 1)
  2256. return Zero;
  2257. // TODO Make this func work on uint, and handle val == 1?
  2258. int rem = Remainder(val);
  2259. return rem == 0
  2260. ? Zero
  2261. : new BigInteger(sign, new int[]{ rem }, false);
  2262. }
  2263. }
  2264. if (CompareNoLeadingZeroes(0, magnitude, 0, n.magnitude) < 0)
  2265. return this;
  2266. int[] result;
  2267. if (n.QuickPow2Check()) // n is power of two
  2268. {
  2269. // TODO Move before small values branch above?
  2270. result = LastNBits(n.Abs().BitLength - 1);
  2271. }
  2272. else
  2273. {
  2274. result = (int[]) this.magnitude.Clone();
  2275. result = Remainder(result, n.magnitude);
  2276. }
  2277. return new BigInteger(sign, result, true);
  2278. }
  2279. private int[] LastNBits(
  2280. int n)
  2281. {
  2282. if (n < 1)
  2283. return ZeroMagnitude;
  2284. int numWords = (n + BitsPerInt - 1) / BitsPerInt;
  2285. numWords = System.Math.Min(numWords, this.magnitude.Length);
  2286. int[] result = new int[numWords];
  2287. Array.Copy(this.magnitude, this.magnitude.Length - numWords, result, 0, numWords);
  2288. int excessBits = (numWords << 5) - n;
  2289. if (excessBits > 0)
  2290. {
  2291. result[0] &= (int)(UInt32.MaxValue >> excessBits);
  2292. }
  2293. return result;
  2294. }
  2295. private BigInteger DivideWords(int w)
  2296. {
  2297. Debug.Assert(w >= 0);
  2298. int n = magnitude.Length;
  2299. if (w >= n)
  2300. return Zero;
  2301. int[] mag = new int[n - w];
  2302. Array.Copy(magnitude, 0, mag, 0, n - w);
  2303. return new BigInteger(sign, mag, false);
  2304. }
  2305. private BigInteger RemainderWords(int w)
  2306. {
  2307. Debug.Assert(w >= 0);
  2308. int n = magnitude.Length;
  2309. if (w >= n)
  2310. return this;
  2311. int[] mag = new int[w];
  2312. Array.Copy(magnitude, n - w, mag, 0, w);
  2313. return new BigInteger(sign, mag, false);
  2314. }
  2315. /**
  2316. * do a left shift - this returns a new array.
  2317. */
  2318. private static int[] ShiftLeft(
  2319. int[] mag,
  2320. int n)
  2321. {
  2322. int nInts = (int)((uint)n >> 5);
  2323. int nBits = n & 0x1f;
  2324. int magLen = mag.Length;
  2325. int[] newMag;
  2326. if (nBits == 0)
  2327. {
  2328. newMag = new int[magLen + nInts];
  2329. mag.CopyTo(newMag, 0);
  2330. }
  2331. else
  2332. {
  2333. int i = 0;
  2334. int nBits2 = 32 - nBits;
  2335. int highBits = (int)((uint)mag[0] >> nBits2);
  2336. if (highBits != 0)
  2337. {
  2338. newMag = new int[magLen + nInts + 1];
  2339. newMag[i++] = highBits;
  2340. }
  2341. else
  2342. {
  2343. newMag = new int[magLen + nInts];
  2344. }
  2345. int m = mag[0];
  2346. for (int j = 0; j < magLen - 1; j++)
  2347. {
  2348. int next = mag[j + 1];
  2349. newMag[i++] = (m << nBits) | (int)((uint)next >> nBits2);
  2350. m = next;
  2351. }
  2352. newMag[i] = mag[magLen - 1] << nBits;
  2353. }
  2354. return newMag;
  2355. }
  2356. private static int ShiftLeftOneInPlace(int[] x, int carry)
  2357. {
  2358. Debug.Assert(carry == 0 || carry == 1);
  2359. int pos = x.Length;
  2360. while (--pos >= 0)
  2361. {
  2362. uint val = (uint)x[pos];
  2363. x[pos] = (int)(val << 1) | carry;
  2364. carry = (int)(val >> 31);
  2365. }
  2366. return carry;
  2367. }
  2368. public BigInteger ShiftLeft(
  2369. int n)
  2370. {
  2371. if (sign == 0 || magnitude.Length == 0)
  2372. return Zero;
  2373. if (n == 0)
  2374. return this;
  2375. if (n < 0)
  2376. return ShiftRight(-n);
  2377. BigInteger result = new BigInteger(sign, ShiftLeft(magnitude, n), true);
  2378. if (this.nBits != -1)
  2379. {
  2380. result.nBits = sign > 0
  2381. ? this.nBits
  2382. : this.nBits + n;
  2383. }
  2384. if (this.nBitLength != -1)
  2385. {
  2386. result.nBitLength = this.nBitLength + n;
  2387. }
  2388. return result;
  2389. }
  2390. /**
  2391. * do a right shift - this does it in place.
  2392. */
  2393. private static void ShiftRightInPlace(
  2394. int start,
  2395. int[] mag,
  2396. int n)
  2397. {
  2398. int nInts = (int)((uint)n >> 5) + start;
  2399. int nBits = n & 0x1f;
  2400. int magEnd = mag.Length - 1;
  2401. if (nInts != start)
  2402. {
  2403. int delta = (nInts - start);
  2404. for (int i = magEnd; i >= nInts; i--)
  2405. {
  2406. mag[i] = mag[i - delta];
  2407. }
  2408. for (int i = nInts - 1; i >= start; i--)
  2409. {
  2410. mag[i] = 0;
  2411. }
  2412. }
  2413. if (nBits != 0)
  2414. {
  2415. int nBits2 = 32 - nBits;
  2416. int m = mag[magEnd];
  2417. for (int i = magEnd; i > nInts; --i)
  2418. {
  2419. int next = mag[i - 1];
  2420. mag[i] = (int)((uint)m >> nBits) | (next << nBits2);
  2421. m = next;
  2422. }
  2423. mag[nInts] = (int)((uint)mag[nInts] >> nBits);
  2424. }
  2425. }
  2426. /**
  2427. * do a right shift by one - this does it in place.
  2428. */
  2429. private static void ShiftRightOneInPlace(
  2430. int start,
  2431. int[] mag)
  2432. {
  2433. int i = mag.Length;
  2434. int m = mag[i - 1];
  2435. while (--i > start)
  2436. {
  2437. int next = mag[i - 1];
  2438. mag[i] = ((int)((uint)m >> 1)) | (next << 31);
  2439. m = next;
  2440. }
  2441. mag[start] = (int)((uint)mag[start] >> 1);
  2442. }
  2443. public BigInteger ShiftRight(
  2444. int n)
  2445. {
  2446. if (n == 0)
  2447. return this;
  2448. if (n < 0)
  2449. return ShiftLeft(-n);
  2450. if (n >= BitLength)
  2451. return (this.sign < 0 ? One.Negate() : Zero);
  2452. // int[] res = (int[]) this.magnitude.Clone();
  2453. //
  2454. // ShiftRightInPlace(0, res, n);
  2455. //
  2456. // return new BigInteger(this.sign, res, true);
  2457. int resultLength = (BitLength - n + 31) >> 5;
  2458. int[] res = new int[resultLength];
  2459. int numInts = n >> 5;
  2460. int numBits = n & 31;
  2461. if (numBits == 0)
  2462. {
  2463. Array.Copy(this.magnitude, 0, res, 0, res.Length);
  2464. }
  2465. else
  2466. {
  2467. int numBits2 = 32 - numBits;
  2468. int magPos = this.magnitude.Length - 1 - numInts;
  2469. for (int i = resultLength - 1; i >= 0; --i)
  2470. {
  2471. res[i] = (int)((uint) this.magnitude[magPos--] >> numBits);
  2472. if (magPos >= 0)
  2473. {
  2474. res[i] |= this.magnitude[magPos] << numBits2;
  2475. }
  2476. }
  2477. }
  2478. Debug.Assert(res[0] != 0);
  2479. return new BigInteger(this.sign, res, false);
  2480. }
  2481. public int SignValue
  2482. {
  2483. get { return sign; }
  2484. }
  2485. /**
  2486. * returns x = x - y - we assume x is >= y
  2487. */
  2488. private static int[] Subtract(
  2489. int xStart,
  2490. int[] x,
  2491. int yStart,
  2492. int[] y)
  2493. {
  2494. Debug.Assert(yStart < y.Length);
  2495. Debug.Assert(x.Length - xStart >= y.Length - yStart);
  2496. int iT = x.Length;
  2497. int iV = y.Length;
  2498. long m;
  2499. int borrow = 0;
  2500. do
  2501. {
  2502. m = (x[--iT] & IMASK) - (y[--iV] & IMASK) + borrow;
  2503. x[iT] = (int) m;
  2504. // borrow = (m < 0) ? -1 : 0;
  2505. borrow = (int)(m >> 63);
  2506. }
  2507. while (iV > yStart);
  2508. if (borrow != 0)
  2509. {
  2510. while (--x[--iT] == -1)
  2511. {
  2512. }
  2513. }
  2514. return x;
  2515. }
  2516. public BigInteger Subtract(
  2517. BigInteger n)
  2518. {
  2519. if (n.sign == 0)
  2520. return this;
  2521. if (this.sign == 0)
  2522. return n.Negate();
  2523. if (this.sign != n.sign)
  2524. return Add(n.Negate());
  2525. int compare = CompareNoLeadingZeroes(0, magnitude, 0, n.magnitude);
  2526. if (compare == 0)
  2527. return Zero;
  2528. BigInteger bigun, lilun;
  2529. if (compare < 0)
  2530. {
  2531. bigun = n;
  2532. lilun = this;
  2533. }
  2534. else
  2535. {
  2536. bigun = this;
  2537. lilun = n;
  2538. }
  2539. return new BigInteger(this.sign * compare, doSubBigLil(bigun.magnitude, lilun.magnitude), true);
  2540. }
  2541. private static int[] doSubBigLil(
  2542. int[] bigMag,
  2543. int[] lilMag)
  2544. {
  2545. int[] res = (int[]) bigMag.Clone();
  2546. return Subtract(0, res, 0, lilMag);
  2547. }
  2548. public byte[] ToByteArray()
  2549. {
  2550. return ToByteArray(false);
  2551. }
  2552. public byte[] ToByteArrayUnsigned()
  2553. {
  2554. return ToByteArray(true);
  2555. }
  2556. private byte[] ToByteArray(
  2557. bool unsigned)
  2558. {
  2559. if (sign == 0)
  2560. return unsigned ? ZeroEncoding : new byte[1];
  2561. int nBits = (unsigned && sign > 0)
  2562. ? BitLength
  2563. : BitLength + 1;
  2564. int nBytes = GetByteLength(nBits);
  2565. byte[] bytes = new byte[nBytes];
  2566. int magIndex = magnitude.Length;
  2567. int bytesIndex = bytes.Length;
  2568. if (sign > 0)
  2569. {
  2570. while (magIndex > 1)
  2571. {
  2572. uint mag = (uint) magnitude[--magIndex];
  2573. bytes[--bytesIndex] = (byte) mag;
  2574. bytes[--bytesIndex] = (byte)(mag >> 8);
  2575. bytes[--bytesIndex] = (byte)(mag >> 16);
  2576. bytes[--bytesIndex] = (byte)(mag >> 24);
  2577. }
  2578. uint lastMag = (uint) magnitude[0];
  2579. while (lastMag > byte.MaxValue)
  2580. {
  2581. bytes[--bytesIndex] = (byte) lastMag;
  2582. lastMag >>= 8;
  2583. }
  2584. bytes[--bytesIndex] = (byte) lastMag;
  2585. }
  2586. else // sign < 0
  2587. {
  2588. bool carry = true;
  2589. while (magIndex > 1)
  2590. {
  2591. uint mag = ~((uint) magnitude[--magIndex]);
  2592. if (carry)
  2593. {
  2594. carry = (++mag == uint.MinValue);
  2595. }
  2596. bytes[--bytesIndex] = (byte) mag;
  2597. bytes[--bytesIndex] = (byte)(mag >> 8);
  2598. bytes[--bytesIndex] = (byte)(mag >> 16);
  2599. bytes[--bytesIndex] = (byte)(mag >> 24);
  2600. }
  2601. uint lastMag = (uint) magnitude[0];
  2602. if (carry)
  2603. {
  2604. // Never wraps because magnitude[0] != 0
  2605. --lastMag;
  2606. }
  2607. while (lastMag > byte.MaxValue)
  2608. {
  2609. bytes[--bytesIndex] = (byte) ~lastMag;
  2610. lastMag >>= 8;
  2611. }
  2612. bytes[--bytesIndex] = (byte) ~lastMag;
  2613. if (bytesIndex > 0)
  2614. {
  2615. bytes[--bytesIndex] = byte.MaxValue;
  2616. }
  2617. }
  2618. return bytes;
  2619. }
  2620. public override string ToString()
  2621. {
  2622. return ToString(10);
  2623. }
  2624. public string ToString(int radix)
  2625. {
  2626. // TODO Make this method work for other radices (ideally 2 <= radix <= 36 as in Java)
  2627. switch (radix)
  2628. {
  2629. case 2:
  2630. case 8:
  2631. case 10:
  2632. case 16:
  2633. break;
  2634. default:
  2635. throw new FormatException("Only bases 2, 8, 10, 16 are allowed");
  2636. }
  2637. // NB: Can only happen to internally managed instances
  2638. if (magnitude == null)
  2639. return "null";
  2640. if (sign == 0)
  2641. return "0";
  2642. // NOTE: This *should* be unnecessary, since the magnitude *should* never have leading zero digits
  2643. int firstNonZero = 0;
  2644. while (firstNonZero < magnitude.Length)
  2645. {
  2646. if (magnitude[firstNonZero] != 0)
  2647. {
  2648. break;
  2649. }
  2650. ++firstNonZero;
  2651. }
  2652. if (firstNonZero == magnitude.Length)
  2653. {
  2654. return "0";
  2655. }
  2656. StringBuilder sb = new StringBuilder();
  2657. if (sign == -1)
  2658. {
  2659. sb.Append('-');
  2660. }
  2661. switch (radix)
  2662. {
  2663. case 2:
  2664. {
  2665. int pos = firstNonZero;
  2666. sb.Append(Convert.ToString(magnitude[pos], 2));
  2667. while (++pos < magnitude.Length)
  2668. {
  2669. AppendZeroExtendedString(sb, Convert.ToString(magnitude[pos], 2), 32);
  2670. }
  2671. break;
  2672. }
  2673. case 8:
  2674. {
  2675. int mask = (1 << 30) - 1;
  2676. BigInteger u = this.Abs();
  2677. int bits = u.BitLength;
  2678. IList S = BestHTTP.SecureProtocol.Org.BouncyCastle.Utilities.Platform.CreateArrayList();
  2679. while (bits > 30)
  2680. {
  2681. S.Add(Convert.ToString(u.IntValue & mask, 8));
  2682. u = u.ShiftRight(30);
  2683. bits -= 30;
  2684. }
  2685. sb.Append(Convert.ToString(u.IntValue, 8));
  2686. for (int i = S.Count - 1; i >= 0; --i)
  2687. {
  2688. AppendZeroExtendedString(sb, (string)S[i], 10);
  2689. }
  2690. break;
  2691. }
  2692. case 16:
  2693. {
  2694. int pos = firstNonZero;
  2695. sb.Append(Convert.ToString(magnitude[pos], 16));
  2696. while (++pos < magnitude.Length)
  2697. {
  2698. AppendZeroExtendedString(sb, Convert.ToString(magnitude[pos], 16), 8);
  2699. }
  2700. break;
  2701. }
  2702. // TODO This could work for other radices if there is an alternative to Convert.ToString method
  2703. //default:
  2704. case 10:
  2705. {
  2706. BigInteger q = this.Abs();
  2707. if (q.BitLength < 64)
  2708. {
  2709. sb.Append(Convert.ToString(q.LongValue, radix));
  2710. break;
  2711. }
  2712. // TODO Could cache the moduli for each radix (soft reference?)
  2713. IList moduli = BestHTTP.SecureProtocol.Org.BouncyCastle.Utilities.Platform.CreateArrayList();
  2714. BigInteger R = BigInteger.ValueOf(radix);
  2715. while (R.CompareTo(q) <= 0)
  2716. {
  2717. moduli.Add(R);
  2718. R = R.Square();
  2719. }
  2720. int scale = moduli.Count;
  2721. sb.EnsureCapacity(sb.Length + (1 << scale));
  2722. ToString(sb, radix, moduli, scale, q);
  2723. break;
  2724. }
  2725. }
  2726. return sb.ToString();
  2727. }
  2728. private static void ToString(StringBuilder sb, int radix, IList moduli, int scale, BigInteger pos)
  2729. {
  2730. if (pos.BitLength < 64)
  2731. {
  2732. string s = Convert.ToString(pos.LongValue, radix);
  2733. if (sb.Length > 1 || (sb.Length == 1 && sb[0] != '-'))
  2734. {
  2735. AppendZeroExtendedString(sb, s, 1 << scale);
  2736. }
  2737. else if (pos.SignValue != 0)
  2738. {
  2739. sb.Append(s);
  2740. }
  2741. return;
  2742. }
  2743. BigInteger[] qr = pos.DivideAndRemainder((BigInteger)moduli[--scale]);
  2744. ToString(sb, radix, moduli, scale, qr[0]);
  2745. ToString(sb, radix, moduli, scale, qr[1]);
  2746. }
  2747. private static void AppendZeroExtendedString(StringBuilder sb, string s, int minLength)
  2748. {
  2749. for (int len = s.Length; len < minLength; ++len)
  2750. {
  2751. sb.Append('0');
  2752. }
  2753. sb.Append(s);
  2754. }
  2755. private static BigInteger CreateUValueOf(
  2756. ulong value)
  2757. {
  2758. int msw = (int)(value >> 32);
  2759. int lsw = (int)value;
  2760. if (msw != 0)
  2761. return new BigInteger(1, new int[] { msw, lsw }, false);
  2762. if (lsw != 0)
  2763. {
  2764. BigInteger n = new BigInteger(1, new int[] { lsw }, false);
  2765. // Check for a power of two
  2766. if ((lsw & -lsw) == lsw)
  2767. {
  2768. n.nBits = 1;
  2769. }
  2770. return n;
  2771. }
  2772. return Zero;
  2773. }
  2774. private static BigInteger CreateValueOf(
  2775. long value)
  2776. {
  2777. if (value < 0)
  2778. {
  2779. if (value == long.MinValue)
  2780. return CreateValueOf(~value).Not();
  2781. return CreateValueOf(-value).Negate();
  2782. }
  2783. return CreateUValueOf((ulong)value);
  2784. }
  2785. public static BigInteger ValueOf(
  2786. long value)
  2787. {
  2788. if (value >= 0 && value < SMALL_CONSTANTS.Length)
  2789. {
  2790. return SMALL_CONSTANTS[value];
  2791. }
  2792. return CreateValueOf(value);
  2793. }
  2794. public int GetLowestSetBit()
  2795. {
  2796. if (this.sign == 0)
  2797. return -1;
  2798. return GetLowestSetBitMaskFirst(-1);
  2799. }
  2800. private int GetLowestSetBitMaskFirst(int firstWordMask)
  2801. {
  2802. int w = magnitude.Length, offset = 0;
  2803. uint word = (uint)(magnitude[--w] & firstWordMask);
  2804. Debug.Assert(magnitude[0] != 0);
  2805. while (word == 0)
  2806. {
  2807. word = (uint)magnitude[--w];
  2808. offset += 32;
  2809. }
  2810. while ((word & 0xFF) == 0)
  2811. {
  2812. word >>= 8;
  2813. offset += 8;
  2814. }
  2815. while ((word & 1) == 0)
  2816. {
  2817. word >>= 1;
  2818. ++offset;
  2819. }
  2820. return offset;
  2821. }
  2822. public bool TestBit(
  2823. int n)
  2824. {
  2825. if (n < 0)
  2826. throw new ArithmeticException("Bit position must not be negative");
  2827. if (sign < 0)
  2828. return !Not().TestBit(n);
  2829. int wordNum = n / 32;
  2830. if (wordNum >= magnitude.Length)
  2831. return false;
  2832. int word = magnitude[magnitude.Length - 1 - wordNum];
  2833. return ((word >> (n % 32)) & 1) > 0;
  2834. }
  2835. public BigInteger Or(
  2836. BigInteger value)
  2837. {
  2838. if (this.sign == 0)
  2839. return value;
  2840. if (value.sign == 0)
  2841. return this;
  2842. int[] aMag = this.sign > 0
  2843. ? this.magnitude
  2844. : Add(One).magnitude;
  2845. int[] bMag = value.sign > 0
  2846. ? value.magnitude
  2847. : value.Add(One).magnitude;
  2848. bool resultNeg = sign < 0 || value.sign < 0;
  2849. int resultLength = System.Math.Max(aMag.Length, bMag.Length);
  2850. int[] resultMag = new int[resultLength];
  2851. int aStart = resultMag.Length - aMag.Length;
  2852. int bStart = resultMag.Length - bMag.Length;
  2853. for (int i = 0; i < resultMag.Length; ++i)
  2854. {
  2855. int aWord = i >= aStart ? aMag[i - aStart] : 0;
  2856. int bWord = i >= bStart ? bMag[i - bStart] : 0;
  2857. if (this.sign < 0)
  2858. {
  2859. aWord = ~aWord;
  2860. }
  2861. if (value.sign < 0)
  2862. {
  2863. bWord = ~bWord;
  2864. }
  2865. resultMag[i] = aWord | bWord;
  2866. if (resultNeg)
  2867. {
  2868. resultMag[i] = ~resultMag[i];
  2869. }
  2870. }
  2871. BigInteger result = new BigInteger(1, resultMag, true);
  2872. // TODO Optimise this case
  2873. if (resultNeg)
  2874. {
  2875. result = result.Not();
  2876. }
  2877. return result;
  2878. }
  2879. public BigInteger Xor(
  2880. BigInteger value)
  2881. {
  2882. if (this.sign == 0)
  2883. return value;
  2884. if (value.sign == 0)
  2885. return this;
  2886. int[] aMag = this.sign > 0
  2887. ? this.magnitude
  2888. : Add(One).magnitude;
  2889. int[] bMag = value.sign > 0
  2890. ? value.magnitude
  2891. : value.Add(One).magnitude;
  2892. // TODO Can just replace with sign != value.sign?
  2893. bool resultNeg = (sign < 0 && value.sign >= 0) || (sign >= 0 && value.sign < 0);
  2894. int resultLength = System.Math.Max(aMag.Length, bMag.Length);
  2895. int[] resultMag = new int[resultLength];
  2896. int aStart = resultMag.Length - aMag.Length;
  2897. int bStart = resultMag.Length - bMag.Length;
  2898. for (int i = 0; i < resultMag.Length; ++i)
  2899. {
  2900. int aWord = i >= aStart ? aMag[i - aStart] : 0;
  2901. int bWord = i >= bStart ? bMag[i - bStart] : 0;
  2902. if (this.sign < 0)
  2903. {
  2904. aWord = ~aWord;
  2905. }
  2906. if (value.sign < 0)
  2907. {
  2908. bWord = ~bWord;
  2909. }
  2910. resultMag[i] = aWord ^ bWord;
  2911. if (resultNeg)
  2912. {
  2913. resultMag[i] = ~resultMag[i];
  2914. }
  2915. }
  2916. BigInteger result = new BigInteger(1, resultMag, true);
  2917. // TODO Optimise this case
  2918. if (resultNeg)
  2919. {
  2920. result = result.Not();
  2921. }
  2922. return result;
  2923. }
  2924. public BigInteger SetBit(
  2925. int n)
  2926. {
  2927. if (n < 0)
  2928. throw new ArithmeticException("Bit address less than zero");
  2929. if (TestBit(n))
  2930. return this;
  2931. // TODO Handle negative values and zero
  2932. if (sign > 0 && n < (BitLength - 1))
  2933. return FlipExistingBit(n);
  2934. return Or(One.ShiftLeft(n));
  2935. }
  2936. public BigInteger ClearBit(
  2937. int n)
  2938. {
  2939. if (n < 0)
  2940. throw new ArithmeticException("Bit address less than zero");
  2941. if (!TestBit(n))
  2942. return this;
  2943. // TODO Handle negative values
  2944. if (sign > 0 && n < (BitLength - 1))
  2945. return FlipExistingBit(n);
  2946. return AndNot(One.ShiftLeft(n));
  2947. }
  2948. public BigInteger FlipBit(
  2949. int n)
  2950. {
  2951. if (n < 0)
  2952. throw new ArithmeticException("Bit address less than zero");
  2953. // TODO Handle negative values and zero
  2954. if (sign > 0 && n < (BitLength - 1))
  2955. return FlipExistingBit(n);
  2956. return Xor(One.ShiftLeft(n));
  2957. }
  2958. private BigInteger FlipExistingBit(
  2959. int n)
  2960. {
  2961. Debug.Assert(sign > 0);
  2962. Debug.Assert(n >= 0);
  2963. Debug.Assert(n < BitLength - 1);
  2964. int[] mag = (int[]) this.magnitude.Clone();
  2965. mag[mag.Length - 1 - (n >> 5)] ^= (1 << (n & 31)); // Flip bit
  2966. //mag[mag.Length - 1 - (n / 32)] ^= (1 << (n % 32));
  2967. return new BigInteger(this.sign, mag, false);
  2968. }
  2969. }
  2970. }
  2971. #pragma warning restore
  2972. #endif