tls_pstm.c 49 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263
  1. /*
  2. * Copyright (C) 2017 Denys Vlasenko
  3. *
  4. * Licensed under GPLv2, see file LICENSE in this source tree.
  5. */
  6. #include "tls.h"
  7. /* The file is taken almost verbatim from matrixssl-3-7-2b-open/crypto/math/.
  8. * Changes are flagged with //bbox
  9. */
  10. /**
  11. * @file pstm.c
  12. * @version 33ef80f (HEAD, tag: MATRIXSSL-3-7-2-OPEN, tag: MATRIXSSL-3-7-2-COMM, origin/master, origin/HEAD, master)
  13. *
  14. * Multiprecision number implementation.
  15. */
  16. /*
  17. * Copyright (c) 2013-2015 INSIDE Secure Corporation
  18. * Copyright (c) PeerSec Networks, 2002-2011
  19. * All Rights Reserved
  20. *
  21. * The latest version of this code is available at http://www.matrixssl.org
  22. *
  23. * This software is open source; you can redistribute it and/or modify
  24. * it under the terms of the GNU General Public License as published by
  25. * the Free Software Foundation; either version 2 of the License, or
  26. * (at your option) any later version.
  27. *
  28. * This General Public License does NOT permit incorporating this software
  29. * into proprietary programs. If you are unable to comply with the GPL, a
  30. * commercial license for this software may be purchased from INSIDE at
  31. * http://www.insidesecure.com/eng/Company/Locations
  32. *
  33. * This program is distributed in WITHOUT ANY WARRANTY; without even the
  34. * implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  35. * See the GNU General Public License for more details.
  36. *
  37. * You should have received a copy of the GNU General Public License
  38. * along with this program; if not, write to the Free Software
  39. * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  40. * http://www.gnu.org/copyleft/gpl.html
  41. */
  42. /******************************************************************************/
  43. //bbox
  44. //#include "../cryptoApi.h"
  45. #ifndef DISABLE_PSTM
  46. static int32 pstm_mul_2d(pstm_int *a, int b, pstm_int *c); //bbox: was int16 b
  47. /******************************************************************************/
  48. /*
  49. init an pstm_int for a given size
  50. */
  51. int32 pstm_init_size(psPool_t *pool, pstm_int * a, uint32 size)
  52. {
  53. //bbox
  54. // uint16 x;
  55. /*
  56. alloc mem
  57. */
  58. a->dp = xzalloc(sizeof (pstm_digit) * size);//bbox
  59. //bbox a->pool = pool;
  60. a->used = 0;
  61. a->alloc = size;
  62. a->sign = PSTM_ZPOS;
  63. /*
  64. zero the digits
  65. */
  66. //bbox
  67. // for (x = 0; x < size; x++) {
  68. // a->dp[x] = 0;
  69. // }
  70. return PSTM_OKAY;
  71. }
  72. /******************************************************************************/
  73. /*
  74. Init a new pstm_int.
  75. */
  76. int32 pstm_init(psPool_t *pool, pstm_int * a)
  77. {
  78. //bbox
  79. // int32 i;
  80. /*
  81. allocate memory required and clear it
  82. */
  83. a->dp = xzalloc(sizeof (pstm_digit) * PSTM_DEFAULT_INIT);//bbox
  84. /*
  85. set the digits to zero
  86. */
  87. //bbox
  88. // for (i = 0; i < PSTM_DEFAULT_INIT; i++) {
  89. // a->dp[i] = 0;
  90. // }
  91. /*
  92. set the used to zero, allocated digits to the default precision and sign
  93. to positive
  94. */
  95. //bbox a->pool = pool;
  96. a->used = 0;
  97. a->alloc = PSTM_DEFAULT_INIT;
  98. a->sign = PSTM_ZPOS;
  99. return PSTM_OKAY;
  100. }
  101. /******************************************************************************/
  102. /*
  103. Grow as required
  104. */
  105. int32 pstm_grow(pstm_int * a, int size)
  106. {
  107. int i; //bbox: was int16
  108. pstm_digit *tmp;
  109. /*
  110. If the alloc size is smaller alloc more ram.
  111. */
  112. if (a->alloc < size) {
  113. /*
  114. Reallocate the array a->dp
  115. We store the return in a temporary variable in case the operation
  116. failed we don't want to overwrite the dp member of a.
  117. */
  118. tmp = xrealloc(a->dp, sizeof (pstm_digit) * size);//bbox
  119. /*
  120. reallocation succeeded so set a->dp
  121. */
  122. a->dp = tmp;
  123. /*
  124. zero excess digits
  125. */
  126. i = a->alloc;
  127. a->alloc = size;
  128. for (; i < a->alloc; i++) {
  129. a->dp[i] = 0;
  130. }
  131. }
  132. return PSTM_OKAY;
  133. }
  134. /******************************************************************************/
  135. /*
  136. copy, b = a (b must be pre-allocated)
  137. */
  138. int32 pstm_copy(pstm_int * a, pstm_int * b)
  139. {
  140. int32 res, n;
  141. /*
  142. If dst == src do nothing
  143. */
  144. if (a == b) {
  145. return PSTM_OKAY;
  146. }
  147. /*
  148. Grow dest
  149. */
  150. if (b->alloc < a->used) {
  151. if ((res = pstm_grow (b, a->used)) != PSTM_OKAY) {
  152. return res;
  153. }
  154. }
  155. /*
  156. Zero b and copy the parameters over
  157. */
  158. {
  159. register pstm_digit *tmpa, *tmpb;
  160. /* pointer aliases */
  161. /* source */
  162. tmpa = a->dp;
  163. /* destination */
  164. tmpb = b->dp;
  165. /* copy all the digits */
  166. for (n = 0; n < a->used; n++) {
  167. *tmpb++ = *tmpa++;
  168. }
  169. /* clear high digits */
  170. for (; n < b->used; n++) {
  171. *tmpb++ = 0;
  172. }
  173. }
  174. /*
  175. copy used count and sign
  176. */
  177. b->used = a->used;
  178. b->sign = a->sign;
  179. return PSTM_OKAY;
  180. }
  181. /******************************************************************************/
  182. /*
  183. Trim unused digits
  184. This is used to ensure that leading zero digits are trimed and the
  185. leading "used" digit will be non-zero. Typically very fast. Also fixes
  186. the sign if there are no more leading digits
  187. */
  188. void pstm_clamp(pstm_int * a)
  189. {
  190. /* decrease used while the most significant digit is zero. */
  191. while (a->used > 0 && a->dp[a->used - 1] == 0) {
  192. --(a->used);
  193. }
  194. /* reset the sign flag if used == 0 */
  195. if (a->used == 0) {
  196. a->sign = PSTM_ZPOS;
  197. }
  198. }
  199. /******************************************************************************/
  200. /*
  201. clear one (frees).
  202. */
  203. void pstm_clear(pstm_int * a)
  204. {
  205. int32 i;
  206. /*
  207. only do anything if a hasn't been freed previously
  208. */
  209. if (a != NULL && a->dp != NULL) {
  210. /*
  211. first zero the digits
  212. */
  213. for (i = 0; i < a->used; i++) {
  214. a->dp[i] = 0;
  215. }
  216. psFree (a->dp, a->pool);
  217. /*
  218. reset members to make debugging easier
  219. */
  220. a->dp = NULL;
  221. a->alloc = a->used = 0;
  222. a->sign = PSTM_ZPOS;
  223. }
  224. }
  225. /******************************************************************************/
  226. /*
  227. clear many (frees).
  228. */
  229. void pstm_clear_multi(pstm_int *mp0, pstm_int *mp1, pstm_int *mp2,
  230. pstm_int *mp3, pstm_int *mp4, pstm_int *mp5,
  231. pstm_int *mp6, pstm_int *mp7)
  232. {
  233. int32 n; /* Number of ok inits */
  234. pstm_int *tempArray[9];
  235. tempArray[0] = mp0;
  236. tempArray[1] = mp1;
  237. tempArray[2] = mp2;
  238. tempArray[3] = mp3;
  239. tempArray[4] = mp4;
  240. tempArray[5] = mp5;
  241. tempArray[6] = mp6;
  242. tempArray[7] = mp7;
  243. tempArray[8] = NULL;
  244. for (n = 0; tempArray[n] != NULL; n++) {
  245. if ((tempArray[n] != NULL) && (tempArray[n]->dp != NULL)) {
  246. pstm_clear(tempArray[n]);
  247. }
  248. }
  249. }
  250. /******************************************************************************/
  251. /*
  252. Set to zero.
  253. */
  254. void pstm_zero(pstm_int * a)
  255. {
  256. int32 n;
  257. pstm_digit *tmp;
  258. a->sign = PSTM_ZPOS;
  259. a->used = 0;
  260. tmp = a->dp;
  261. for (n = 0; n < a->alloc; n++) {
  262. *tmp++ = 0;
  263. }
  264. }
  265. /******************************************************************************/
  266. /*
  267. Compare maginitude of two ints (unsigned).
  268. */
  269. int32 pstm_cmp_mag(pstm_int * a, pstm_int * b)
  270. {
  271. int n; //bbox: was int16
  272. pstm_digit *tmpa, *tmpb;
  273. /*
  274. compare based on # of non-zero digits
  275. */
  276. if (a->used > b->used) {
  277. return PSTM_GT;
  278. }
  279. if (a->used < b->used) {
  280. return PSTM_LT;
  281. }
  282. /* alias for a */
  283. tmpa = a->dp + (a->used - 1);
  284. /* alias for b */
  285. tmpb = b->dp + (a->used - 1);
  286. /*
  287. compare based on digits
  288. */
  289. for (n = 0; n < a->used; ++n, --tmpa, --tmpb) {
  290. if (*tmpa > *tmpb) {
  291. return PSTM_GT;
  292. }
  293. if (*tmpa < *tmpb) {
  294. return PSTM_LT;
  295. }
  296. }
  297. return PSTM_EQ;
  298. }
  299. /******************************************************************************/
  300. /*
  301. Compare two ints (signed)
  302. */
  303. int32 pstm_cmp(pstm_int * a, pstm_int * b)
  304. {
  305. /*
  306. compare based on sign
  307. */
  308. if (a->sign != b->sign) {
  309. if (a->sign == PSTM_NEG) {
  310. return PSTM_LT;
  311. } else {
  312. return PSTM_GT;
  313. }
  314. }
  315. /*
  316. compare digits
  317. */
  318. if (a->sign == PSTM_NEG) {
  319. /* if negative compare opposite direction */
  320. return pstm_cmp_mag(b, a);
  321. } else {
  322. return pstm_cmp_mag(a, b);
  323. }
  324. }
  325. /******************************************************************************/
  326. /*
  327. pstm_ints can be initialized more precisely when they will populated
  328. using pstm_read_unsigned_bin since the length of the byte stream is known
  329. */
  330. int32 pstm_init_for_read_unsigned_bin(psPool_t *pool, pstm_int *a, uint32 len)
  331. {
  332. int32 size;
  333. /*
  334. Need to set this based on how many words max it will take to store the bin.
  335. The magic + 2:
  336. 1 to round up for the remainder of this integer math
  337. 1 for the initial carry of '1' bits that fall between DIGIT_BIT and 8
  338. */
  339. size = (((len / sizeof(pstm_digit)) * (sizeof(pstm_digit) * CHAR_BIT))
  340. / DIGIT_BIT) + 2;
  341. return pstm_init_size(pool, a, size);
  342. }
  343. /******************************************************************************/
  344. /*
  345. Reads a unsigned char array into pstm_int format. User should have
  346. called pstm_init_for_read_unsigned_bin first. There is some grow logic
  347. here if the default pstm_init was used but we don't really want to hit it.
  348. */
  349. int32 pstm_read_unsigned_bin(pstm_int *a, unsigned char *b, int32 c)
  350. {
  351. /* zero the int */
  352. pstm_zero (a);
  353. /*
  354. If we know the endianness of this architecture, and we're using
  355. 32-bit pstm_digits, we can optimize this
  356. */
  357. #if (defined(ENDIAN_LITTLE) || defined(ENDIAN_BIG)) && !defined(PSTM_64BIT)
  358. /* But not for both simultaneously */
  359. #if defined(ENDIAN_LITTLE) && defined(ENDIAN_BIG)
  360. #error Both ENDIAN_LITTLE and ENDIAN_BIG defined.
  361. #endif
  362. {
  363. unsigned char *pd;
  364. if ((unsigned)c > (PSTM_MAX_SIZE * sizeof(pstm_digit))) {
  365. uint32 excess = c - (PSTM_MAX_SIZE * sizeof(pstm_digit));
  366. c -= excess;
  367. b += excess;
  368. }
  369. a->used = ((c + sizeof(pstm_digit) - 1)/sizeof(pstm_digit));
  370. if (a->alloc < a->used) {
  371. if (pstm_grow(a, a->used) != PSTM_OKAY) {
  372. return PSTM_MEM;
  373. }
  374. }
  375. pd = (unsigned char *)a->dp;
  376. /* read the bytes in */
  377. #ifdef ENDIAN_BIG
  378. {
  379. /* Use Duff's device to unroll the loop. */
  380. int32 idx = (c - 1) & ~3;
  381. switch (c % 4) {
  382. case 0: do { pd[idx+0] = *b++;
  383. case 3: pd[idx+1] = *b++;
  384. case 2: pd[idx+2] = *b++;
  385. case 1: pd[idx+3] = *b++;
  386. idx -= 4;
  387. } while ((c -= 4) > 0);
  388. }
  389. }
  390. #else
  391. for (c -= 1; c >= 0; c -= 1) {
  392. pd[c] = *b++;
  393. }
  394. #endif
  395. }
  396. #else
  397. /* Big enough based on the len? */
  398. a->used = (((c / sizeof(pstm_digit)) * (sizeof(pstm_digit) * CHAR_BIT))
  399. / DIGIT_BIT) + 2;
  400. if (a->alloc < a->used) {
  401. if (pstm_grow(a, a->used) != PSTM_OKAY) {
  402. return PSTM_MEM;
  403. }
  404. }
  405. /* read the bytes in */
  406. for (; c > 0; c--) {
  407. if (pstm_mul_2d (a, 8, a) != PSTM_OKAY) {
  408. return PS_MEM_FAIL;
  409. }
  410. a->dp[0] |= *b++;
  411. a->used += 1;
  412. }
  413. #endif
  414. pstm_clamp (a);
  415. return PS_SUCCESS;
  416. }
  417. /******************************************************************************/
  418. /*
  419. */
  420. int pstm_count_bits (pstm_int * a)
  421. {
  422. int r; //bbox: was int16
  423. pstm_digit q;
  424. if (a->used == 0) {
  425. return 0;
  426. }
  427. /* get number of digits and add that */
  428. r = (a->used - 1) * DIGIT_BIT;
  429. /* take the last digit and count the bits in it */
  430. q = a->dp[a->used - 1];
  431. while (q > ((pstm_digit) 0)) {
  432. ++r;
  433. q >>= ((pstm_digit) 1);
  434. }
  435. return r;
  436. }
  437. /******************************************************************************/
  438. int32 pstm_unsigned_bin_size(pstm_int *a)
  439. {
  440. int32 size = pstm_count_bits (a);
  441. return (size / 8 + ((size & 7) != 0 ? 1 : 0));
  442. }
  443. /******************************************************************************/
  444. void pstm_set(pstm_int *a, pstm_digit b)
  445. {
  446. pstm_zero(a);
  447. a->dp[0] = b;
  448. a->used = a->dp[0] ? 1 : 0;
  449. }
  450. /******************************************************************************/
  451. /*
  452. Right shift
  453. */
  454. void pstm_rshd(pstm_int *a, int x)
  455. {
  456. int y; //bbox: was int16
  457. /* too many digits just zero and return */
  458. if (x >= a->used) {
  459. pstm_zero(a);
  460. return;
  461. }
  462. /* shift */
  463. for (y = 0; y < a->used - x; y++) {
  464. a->dp[y] = a->dp[y+x];
  465. }
  466. /* zero rest */
  467. for (; y < a->used; y++) {
  468. a->dp[y] = 0;
  469. }
  470. /* decrement count */
  471. a->used -= x;
  472. pstm_clamp(a);
  473. }
  474. /******************************************************************************/
  475. /*
  476. Shift left a certain amount of digits.
  477. */
  478. int32 pstm_lshd(pstm_int * a, int b)
  479. {
  480. int x; //bbox: was int16
  481. int32 res;
  482. /*
  483. If its less than zero return.
  484. */
  485. if (b <= 0) {
  486. return PSTM_OKAY;
  487. }
  488. /*
  489. Grow to fit the new digits.
  490. */
  491. if (a->alloc < a->used + b) {
  492. if ((res = pstm_grow (a, a->used + b)) != PSTM_OKAY) {
  493. return res;
  494. }
  495. }
  496. {
  497. register pstm_digit *top, *bottom;
  498. /*
  499. Increment the used by the shift amount then copy upwards.
  500. */
  501. a->used += b;
  502. /* top */
  503. top = a->dp + a->used - 1;
  504. /* base */
  505. bottom = a->dp + a->used - 1 - b;
  506. /*
  507. This is implemented using a sliding window except the window goes the
  508. other way around. Copying from the bottom to the top.
  509. */
  510. for (x = a->used - 1; x >= b; x--) {
  511. *top-- = *bottom--;
  512. }
  513. /* zero the lower digits */
  514. top = a->dp;
  515. for (x = 0; x < b; x++) {
  516. *top++ = 0;
  517. }
  518. }
  519. return PSTM_OKAY;
  520. }
  521. /******************************************************************************/
  522. /*
  523. computes a = 2**b
  524. */
  525. int32 pstm_2expt(pstm_int *a, int b)
  526. {
  527. int z; //bbox: was int16
  528. /* zero a as per default */
  529. pstm_zero (a);
  530. if (b < 0) {
  531. return PSTM_OKAY;
  532. }
  533. z = b / DIGIT_BIT;
  534. if (z >= PSTM_MAX_SIZE) {
  535. return PS_LIMIT_FAIL;
  536. }
  537. /* set the used count of where the bit will go */
  538. a->used = z + 1;
  539. if (a->used > a->alloc) {
  540. if (pstm_grow(a, a->used) != PSTM_OKAY) {
  541. return PS_MEM_FAIL;
  542. }
  543. }
  544. /* put the single bit in its place */
  545. a->dp[z] = ((pstm_digit)1) << (b % DIGIT_BIT);
  546. return PSTM_OKAY;
  547. }
  548. /******************************************************************************/
  549. /*
  550. */
  551. int32 pstm_mul_2(pstm_int * a, pstm_int * b)
  552. {
  553. int32 res;
  554. int x, oldused; //bbox: was int16
  555. /*
  556. grow to accomodate result
  557. */
  558. if (b->alloc < a->used + 1) {
  559. if ((res = pstm_grow (b, a->used + 1)) != PSTM_OKAY) {
  560. return res;
  561. }
  562. }
  563. oldused = b->used;
  564. b->used = a->used;
  565. {
  566. register pstm_digit r, rr, *tmpa, *tmpb;
  567. /* alias for source */
  568. tmpa = a->dp;
  569. /* alias for dest */
  570. tmpb = b->dp;
  571. /* carry */
  572. r = 0;
  573. for (x = 0; x < a->used; x++) {
  574. /*
  575. get what will be the *next* carry bit from the
  576. MSB of the current digit
  577. */
  578. rr = *tmpa >> ((pstm_digit)(DIGIT_BIT - 1));
  579. /*
  580. now shift up this digit, add in the carry [from the previous]
  581. */
  582. *tmpb++ = ((*tmpa++ << ((pstm_digit)1)) | r);
  583. /*
  584. copy the carry that would be from the source
  585. digit into the next iteration
  586. */
  587. r = rr;
  588. }
  589. /* new leading digit? */
  590. if (r != 0 && b->used != (PSTM_MAX_SIZE-1)) {
  591. /* add a MSB which is always 1 at this point */
  592. *tmpb = 1;
  593. ++(b->used);
  594. }
  595. /*
  596. now zero any excess digits on the destination that we didn't write to
  597. */
  598. tmpb = b->dp + b->used;
  599. for (x = b->used; x < oldused; x++) {
  600. *tmpb++ = 0;
  601. }
  602. }
  603. b->sign = a->sign;
  604. return PSTM_OKAY;
  605. }
  606. /******************************************************************************/
  607. /*
  608. unsigned subtraction ||a|| >= ||b|| ALWAYS!
  609. */
  610. int32 s_pstm_sub(pstm_int *a, pstm_int *b, pstm_int *c)
  611. {
  612. int oldbused, oldused; //bbox: was int16
  613. int32 x;
  614. pstm_word t;
  615. if (b->used > a->used) {
  616. return PS_LIMIT_FAIL;
  617. }
  618. if (c->alloc < a->used) {
  619. if ((x = pstm_grow (c, a->used)) != PSTM_OKAY) {
  620. return x;
  621. }
  622. }
  623. oldused = c->used;
  624. oldbused = b->used;
  625. c->used = a->used;
  626. t = 0;
  627. for (x = 0; x < oldbused; x++) {
  628. t = ((pstm_word)a->dp[x]) - (((pstm_word)b->dp[x]) + t);
  629. c->dp[x] = (pstm_digit)t;
  630. t = (t >> DIGIT_BIT)&1;
  631. }
  632. for (; x < a->used; x++) {
  633. t = ((pstm_word)a->dp[x]) - t;
  634. c->dp[x] = (pstm_digit)t;
  635. t = (t >> DIGIT_BIT);
  636. }
  637. for (; x < oldused; x++) {
  638. c->dp[x] = 0;
  639. }
  640. pstm_clamp(c);
  641. return PSTM_OKAY;
  642. }
  643. /******************************************************************************/
  644. /*
  645. unsigned addition
  646. */
  647. static int32 s_pstm_add(pstm_int *a, pstm_int *b, pstm_int *c)
  648. {
  649. int x, y, oldused; //bbox: was int16
  650. register pstm_word t, adp, bdp;
  651. y = a->used;
  652. if (b->used > y) {
  653. y = b->used;
  654. }
  655. oldused = c->used;
  656. c->used = y;
  657. if (c->used > c->alloc) {
  658. if (pstm_grow(c, c->used) != PSTM_OKAY) {
  659. return PS_MEM_FAIL;
  660. }
  661. }
  662. t = 0;
  663. for (x = 0; x < y; x++) {
  664. if (a->used < x) {
  665. adp = 0;
  666. } else {
  667. adp = (pstm_word)a->dp[x];
  668. }
  669. if (b->used < x) {
  670. bdp = 0;
  671. } else {
  672. bdp = (pstm_word)b->dp[x];
  673. }
  674. t += (adp) + (bdp);
  675. c->dp[x] = (pstm_digit)t;
  676. t >>= DIGIT_BIT;
  677. }
  678. if (t != 0 && x < PSTM_MAX_SIZE) {
  679. if (c->used == c->alloc) {
  680. if (pstm_grow(c, c->alloc + 1) != PSTM_OKAY) {
  681. return PS_MEM_FAIL;
  682. }
  683. }
  684. c->dp[c->used++] = (pstm_digit)t;
  685. ++x;
  686. }
  687. c->used = x;
  688. for (; x < oldused; x++) {
  689. c->dp[x] = 0;
  690. }
  691. pstm_clamp(c);
  692. return PSTM_OKAY;
  693. }
  694. /******************************************************************************/
  695. /*
  696. */
  697. int32 pstm_sub(pstm_int *a, pstm_int *b, pstm_int *c)
  698. {
  699. int32 res;
  700. int sa, sb; //bbox: was int16
  701. sa = a->sign;
  702. sb = b->sign;
  703. if (sa != sb) {
  704. /*
  705. subtract a negative from a positive, OR a positive from a negative.
  706. For both, ADD their magnitudes, and use the sign of the first number.
  707. */
  708. c->sign = sa;
  709. if ((res = s_pstm_add (a, b, c)) != PSTM_OKAY) {
  710. return res;
  711. }
  712. } else {
  713. /*
  714. subtract a positive from a positive, OR a negative from a negative.
  715. First, take the difference between their magnitudes, then...
  716. */
  717. if (pstm_cmp_mag (a, b) != PSTM_LT) {
  718. /* Copy the sign from the first */
  719. c->sign = sa;
  720. /* The first has a larger or equal magnitude */
  721. if ((res = s_pstm_sub (a, b, c)) != PSTM_OKAY) {
  722. return res;
  723. }
  724. } else {
  725. /* The result has the _opposite_ sign from the first number. */
  726. c->sign = (sa == PSTM_ZPOS) ? PSTM_NEG : PSTM_ZPOS;
  727. /* The second has a larger magnitude */
  728. if ((res = s_pstm_sub (b, a, c)) != PSTM_OKAY) {
  729. return res;
  730. }
  731. }
  732. }
  733. return PS_SUCCESS;
  734. }
  735. /******************************************************************************/
  736. /*
  737. c = a - b
  738. */
  739. int32 pstm_sub_d(psPool_t *pool, pstm_int *a, pstm_digit b, pstm_int *c)
  740. {
  741. pstm_int tmp;
  742. int32 res;
  743. if (pstm_init_size(pool, &tmp, sizeof(pstm_digit)) != PSTM_OKAY) {
  744. return PS_MEM_FAIL;
  745. }
  746. pstm_set(&tmp, b);
  747. res = pstm_sub(a, &tmp, c);
  748. pstm_clear(&tmp);
  749. return res;
  750. }
  751. /******************************************************************************/
  752. /*
  753. setups the montgomery reduction
  754. */
  755. int32 pstm_montgomery_setup(pstm_int *a, pstm_digit *rho)
  756. {
  757. pstm_digit x, b;
  758. /*
  759. fast inversion mod 2**k
  760. Based on the fact that
  761. XA = 1 (mod 2**n) => (X(2-XA)) A = 1 (mod 2**2n)
  762. => 2*X*A - X*X*A*A = 1
  763. => 2*(1) - (1) = 1
  764. */
  765. b = a->dp[0];
  766. if ((b & 1) == 0) {
  767. psTraceCrypto("pstm_montogomery_setup failure\n");
  768. return PS_ARG_FAIL;
  769. }
  770. x = (((b + 2) & 4) << 1) + b; /* here x*a==1 mod 2**4 */
  771. x *= 2 - b * x; /* here x*a==1 mod 2**8 */
  772. x *= 2 - b * x; /* here x*a==1 mod 2**16 */
  773. x *= 2 - b * x; /* here x*a==1 mod 2**32 */
  774. #ifdef PSTM_64BIT
  775. x *= 2 - b * x; /* here x*a==1 mod 2**64 */
  776. #endif
  777. /* rho = -1/m mod b */
  778. *rho = (pstm_digit)(((pstm_word) 1 << ((pstm_word) DIGIT_BIT)) -
  779. ((pstm_word)x));
  780. return PSTM_OKAY;
  781. }
  782. /******************************************************************************/
  783. /*
  784. * computes a = B**n mod b without division or multiplication useful for
  785. * normalizing numbers in a Montgomery system.
  786. */
  787. int32 pstm_montgomery_calc_normalization(pstm_int *a, pstm_int *b)
  788. {
  789. int32 x;
  790. int bits; //bbox: was int16
  791. /* how many bits of last digit does b use */
  792. bits = pstm_count_bits (b) % DIGIT_BIT;
  793. if (!bits) bits = DIGIT_BIT;
  794. /* compute A = B^(n-1) * 2^(bits-1) */
  795. if (b->used > 1) {
  796. if ((x = pstm_2expt (a, (b->used - 1) * DIGIT_BIT + bits - 1)) !=
  797. PSTM_OKAY) {
  798. return x;
  799. }
  800. } else {
  801. pstm_set(a, 1);
  802. bits = 1;
  803. }
  804. /* now compute C = A * B mod b */
  805. for (x = bits - 1; x < (int32)DIGIT_BIT; x++) {
  806. if (pstm_mul_2 (a, a) != PSTM_OKAY) {
  807. return PS_MEM_FAIL;
  808. }
  809. if (pstm_cmp_mag (a, b) != PSTM_LT) {
  810. if (s_pstm_sub (a, b, a) != PSTM_OKAY) {
  811. return PS_MEM_FAIL;
  812. }
  813. }
  814. }
  815. return PSTM_OKAY;
  816. }
  817. /******************************************************************************/
  818. /*
  819. c = a * 2**d
  820. */
  821. static int32 pstm_mul_2d(pstm_int *a, int b, pstm_int *c)
  822. {
  823. pstm_digit carry, carrytmp, shift;
  824. int x; //bbox: was int16
  825. /* copy it */
  826. if (pstm_copy(a, c) != PSTM_OKAY) {
  827. return PS_MEM_FAIL;
  828. }
  829. /* handle whole digits */
  830. if (b >= DIGIT_BIT) {
  831. if (pstm_lshd(c, b/DIGIT_BIT) != PSTM_OKAY) {
  832. return PS_MEM_FAIL;
  833. }
  834. }
  835. b %= DIGIT_BIT;
  836. /* shift the digits */
  837. if (b != 0) {
  838. carry = 0;
  839. shift = DIGIT_BIT - b;
  840. for (x = 0; x < c->used; x++) {
  841. carrytmp = c->dp[x] >> shift;
  842. c->dp[x] = (c->dp[x] << b) + carry;
  843. carry = carrytmp;
  844. }
  845. /* store last carry if room */
  846. if (carry && x < PSTM_MAX_SIZE) {
  847. if (c->used == c->alloc) {
  848. if (pstm_grow(c, c->alloc + 1) != PSTM_OKAY) {
  849. return PS_MEM_FAIL;
  850. }
  851. }
  852. c->dp[c->used++] = carry;
  853. }
  854. }
  855. pstm_clamp(c);
  856. return PSTM_OKAY;
  857. }
  858. /******************************************************************************/
  859. /*
  860. c = a mod 2**d
  861. */
  862. static int32 pstm_mod_2d(pstm_int *a, int b, pstm_int *c) //bbox: was int16 b
  863. {
  864. int x; //bbox: was int16
  865. /* zero if count less than or equal to zero */
  866. if (b <= 0) {
  867. pstm_zero(c);
  868. return PSTM_OKAY;
  869. }
  870. /* get copy of input */
  871. if (pstm_copy(a, c) != PSTM_OKAY) {
  872. return PS_MEM_FAIL;
  873. }
  874. /* if 2**d is larger than we just return */
  875. if (b >= (DIGIT_BIT * a->used)) {
  876. return PSTM_OKAY;
  877. }
  878. /* zero digits above the last digit of the modulus */
  879. for (x = (b / DIGIT_BIT) + ((b % DIGIT_BIT) == 0 ? 0 : 1); x < c->used; x++)
  880. {
  881. c->dp[x] = 0;
  882. }
  883. /* clear the digit that is not completely outside/inside the modulus */
  884. c->dp[b / DIGIT_BIT] &= ~((pstm_digit)0) >> (DIGIT_BIT - b);
  885. pstm_clamp (c);
  886. return PSTM_OKAY;
  887. }
  888. /******************************************************************************/
  889. /*
  890. c = a * b
  891. */
  892. int32 pstm_mul_d(pstm_int *a, pstm_digit b, pstm_int *c)
  893. {
  894. pstm_word w;
  895. int32 res;
  896. int x, oldused; //bbox: was int16
  897. if (c->alloc < a->used + 1) {
  898. if ((res = pstm_grow (c, a->used + 1)) != PSTM_OKAY) {
  899. return res;
  900. }
  901. }
  902. oldused = c->used;
  903. c->used = a->used;
  904. c->sign = a->sign;
  905. w = 0;
  906. for (x = 0; x < a->used; x++) {
  907. w = ((pstm_word)a->dp[x]) * ((pstm_word)b) + w;
  908. c->dp[x] = (pstm_digit)w;
  909. w = w >> DIGIT_BIT;
  910. }
  911. if (w != 0 && (a->used != PSTM_MAX_SIZE)) {
  912. c->dp[c->used++] = (pstm_digit)w;
  913. ++x;
  914. }
  915. for (; x < oldused; x++) {
  916. c->dp[x] = 0;
  917. }
  918. pstm_clamp(c);
  919. return PSTM_OKAY;
  920. }
  921. /******************************************************************************/
  922. /*
  923. c = a / 2**b
  924. */
  925. int32 pstm_div_2d(psPool_t *pool, pstm_int *a, int b, pstm_int *c,
  926. pstm_int *d)
  927. {
  928. pstm_digit D, r, rr;
  929. int32 res;
  930. int x; //bbox: was int16
  931. pstm_int t;
  932. /* if the shift count is <= 0 then we do no work */
  933. if (b <= 0) {
  934. if (pstm_copy (a, c) != PSTM_OKAY) {
  935. return PS_MEM_FAIL;
  936. }
  937. if (d != NULL) {
  938. pstm_zero (d);
  939. }
  940. return PSTM_OKAY;
  941. }
  942. /* get the remainder */
  943. if (d != NULL) {
  944. if (pstm_init(pool, &t) != PSTM_OKAY) {
  945. return PS_MEM_FAIL;
  946. }
  947. if (pstm_mod_2d (a, b, &t) != PSTM_OKAY) {
  948. res = PS_MEM_FAIL;
  949. goto LBL_DONE;
  950. }
  951. }
  952. /* copy */
  953. if (pstm_copy(a, c) != PSTM_OKAY) {
  954. res = PS_MEM_FAIL;
  955. goto LBL_DONE;
  956. }
  957. /* shift by as many digits in the bit count */
  958. if (b >= (int32)DIGIT_BIT) {
  959. pstm_rshd (c, b / DIGIT_BIT);
  960. }
  961. /* shift any bit count < DIGIT_BIT */
  962. D = (pstm_digit) (b % DIGIT_BIT);
  963. if (D != 0) {
  964. register pstm_digit *tmpc, mask, shift;
  965. /* mask */
  966. mask = (((pstm_digit)1) << D) - 1;
  967. /* shift for lsb */
  968. shift = DIGIT_BIT - D;
  969. /* alias */
  970. tmpc = c->dp + (c->used - 1);
  971. /* carry */
  972. r = 0;
  973. for (x = c->used - 1; x >= 0; x--) {
  974. /* get the lower bits of this word in a temp */
  975. rr = *tmpc & mask;
  976. /* shift the current word and mix in the carry bits from previous */
  977. *tmpc = (*tmpc >> D) | (r << shift);
  978. --tmpc;
  979. /* set the carry to the carry bits of the current word above */
  980. r = rr;
  981. }
  982. }
  983. pstm_clamp (c);
  984. res = PSTM_OKAY;
  985. LBL_DONE:
  986. if (d != NULL) {
  987. if (pstm_copy(&t, d) != PSTM_OKAY) {
  988. res = PS_MEM_FAIL;
  989. }
  990. pstm_clear(&t);
  991. }
  992. return res;
  993. }
  994. /******************************************************************************/
  995. /*
  996. b = a/2
  997. */
  998. int32 pstm_div_2(pstm_int * a, pstm_int * b)
  999. {
  1000. int x, oldused; //bbox: was int16
  1001. if (b->alloc < a->used) {
  1002. if (pstm_grow(b, a->used) != PSTM_OKAY) {
  1003. return PS_MEM_FAIL;
  1004. }
  1005. }
  1006. oldused = b->used;
  1007. b->used = a->used;
  1008. {
  1009. register pstm_digit r, rr, *tmpa, *tmpb;
  1010. /* source alias */
  1011. tmpa = a->dp + b->used - 1;
  1012. /* dest alias */
  1013. tmpb = b->dp + b->used - 1;
  1014. /* carry */
  1015. r = 0;
  1016. for (x = b->used - 1; x >= 0; x--) {
  1017. /* get the carry for the next iteration */
  1018. rr = *tmpa & 1;
  1019. /* shift the current digit, add in carry and store */
  1020. *tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1));
  1021. /* forward carry to next iteration */
  1022. r = rr;
  1023. }
  1024. /* zero excess digits */
  1025. tmpb = b->dp + b->used;
  1026. for (x = b->used; x < oldused; x++) {
  1027. *tmpb++ = 0;
  1028. }
  1029. }
  1030. b->sign = a->sign;
  1031. pstm_clamp (b);
  1032. return PSTM_OKAY;
  1033. }
  1034. /******************************************************************************/
  1035. /*
  1036. Creates "a" then copies b into it
  1037. */
  1038. int32 pstm_init_copy(psPool_t *pool, pstm_int * a, pstm_int * b, int toSqr)
  1039. {
  1040. int x; //bbox: was int16
  1041. int32 res;
  1042. if (a == b) {
  1043. return PSTM_OKAY;
  1044. }
  1045. x = b->alloc;
  1046. if (toSqr) {
  1047. /*
  1048. Smart-size: Increasing size of a if b->used is roughly half
  1049. of b->alloc because usage has shown that a lot of these copies
  1050. go on to be squared and need these extra digits
  1051. */
  1052. if ((b->used * 2) + 2 >= x) {
  1053. x = (b->used * 2) + 3;
  1054. }
  1055. }
  1056. if ((res = pstm_init_size(pool, a, x)) != PSTM_OKAY) {
  1057. return res;
  1058. }
  1059. return pstm_copy(b, a);
  1060. }
  1061. /******************************************************************************/
  1062. /*
  1063. With some compilers, we have seen issues linking with the builtin
  1064. 64 bit division routine. The issues with either manifest in a failure
  1065. to find 'udivdi3' at link time, or a runtime invalid instruction fault
  1066. during an RSA operation.
  1067. The routine below divides a 64 bit unsigned int by a 32 bit unsigned int
  1068. explicitly, rather than using the division operation
  1069. The 64 bit result is placed in the 'numerator' parameter
  1070. The 32 bit mod (remainder) of the division is the return parameter
  1071. Based on implementations by:
  1072. Copyright (C) 2003 Bernardo Innocenti <bernie@develer.com>
  1073. Copyright (C) 1999 Hewlett-Packard Co
  1074. Copyright (C) 1999 David Mosberger-Tang <davidm@hpl.hp.com>
  1075. */
  1076. #if defined(USE_MATRIX_DIV64) && defined(PSTM_32BIT)
  1077. static uint32 psDiv64(uint64 *numerator, uint32 denominator)
  1078. {
  1079. uint64 rem = *numerator;
  1080. uint64 b = denominator;
  1081. uint64 res = 0;
  1082. uint64 d = 1;
  1083. uint32 high = rem >> 32;
  1084. if (high >= denominator) {
  1085. high /= denominator;
  1086. res = (uint64) high << 32;
  1087. rem -= (uint64) (high * denominator) << 32;
  1088. }
  1089. while ((int64)b > 0 && b < rem) {
  1090. b = b+b;
  1091. d = d+d;
  1092. }
  1093. do {
  1094. if (rem >= b) {
  1095. rem -= b;
  1096. res += d;
  1097. }
  1098. b >>= 1;
  1099. d >>= 1;
  1100. } while (d);
  1101. *numerator = res;
  1102. return rem;
  1103. }
  1104. #endif /* USE_MATRIX_DIV64 */
  1105. #if defined(USE_MATRIX_DIV128) && defined(PSTM_64BIT)
  1106. typedef unsigned long uint128 __attribute__ ((mode(TI)));
  1107. static uint64 psDiv128(uint128 *numerator, uint64 denominator)
  1108. {
  1109. uint128 rem = *numerator;
  1110. uint128 b = denominator;
  1111. uint128 res = 0;
  1112. uint128 d = 1;
  1113. uint64 high = rem >> 64;
  1114. if (high >= denominator) {
  1115. high /= denominator;
  1116. res = (uint128) high << 64;
  1117. rem -= (uint128) (high * denominator) << 64;
  1118. }
  1119. while ((uint128)b > 0 && b < rem) {
  1120. b = b+b;
  1121. d = d+d;
  1122. }
  1123. do {
  1124. if (rem >= b) {
  1125. rem -= b;
  1126. res += d;
  1127. }
  1128. b >>= 1;
  1129. d >>= 1;
  1130. } while (d);
  1131. *numerator = res;
  1132. return rem;
  1133. }
  1134. #endif /* USE_MATRIX_DIV128 */
  1135. /******************************************************************************/
  1136. /*
  1137. a/b => cb + d == a
  1138. */
  1139. int32 pstm_div(psPool_t *pool, pstm_int *a, pstm_int *b, pstm_int *c,
  1140. pstm_int *d)
  1141. {
  1142. pstm_int q, x, y, t1, t2;
  1143. int32 res;
  1144. int n, t, i, norm, neg; //bbox: was int16
  1145. /* is divisor zero ? */
  1146. if (pstm_iszero (b) == 1) {
  1147. return PS_LIMIT_FAIL;
  1148. }
  1149. /* if a < b then q=0, r = a */
  1150. if (pstm_cmp_mag (a, b) == PSTM_LT) {
  1151. if (d != NULL) {
  1152. if (pstm_copy(a, d) != PSTM_OKAY) {
  1153. return PS_MEM_FAIL;
  1154. }
  1155. }
  1156. if (c != NULL) {
  1157. pstm_zero (c);
  1158. }
  1159. return PSTM_OKAY;
  1160. }
  1161. /*
  1162. Smart-size inits
  1163. */
  1164. if ((res = pstm_init_size(pool, &t1, a->alloc)) != PSTM_OKAY) {
  1165. return res;
  1166. }
  1167. if ((res = pstm_init_size(pool, &t2, 3)) != PSTM_OKAY) {
  1168. goto LBL_T1;
  1169. }
  1170. if ((res = pstm_init_copy(pool, &x, a, 0)) != PSTM_OKAY) {
  1171. goto LBL_T2;
  1172. }
  1173. /*
  1174. Used to be an init_copy on b but pstm_grow was always hit with triple size
  1175. */
  1176. if ((res = pstm_init_size(pool, &y, b->used * 3)) != PSTM_OKAY) {
  1177. goto LBL_X;
  1178. }
  1179. if ((res = pstm_copy(b, &y)) != PSTM_OKAY) {
  1180. goto LBL_Y;
  1181. }
  1182. /* fix the sign */
  1183. neg = (a->sign == b->sign) ? PSTM_ZPOS : PSTM_NEG;
  1184. x.sign = y.sign = PSTM_ZPOS;
  1185. /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
  1186. norm = pstm_count_bits(&y) % DIGIT_BIT;
  1187. if (norm < (int32)(DIGIT_BIT-1)) {
  1188. norm = (DIGIT_BIT-1) - norm;
  1189. if ((res = pstm_mul_2d(&x, norm, &x)) != PSTM_OKAY) {
  1190. goto LBL_Y;
  1191. }
  1192. if ((res = pstm_mul_2d(&y, norm, &y)) != PSTM_OKAY) {
  1193. goto LBL_Y;
  1194. }
  1195. } else {
  1196. norm = 0;
  1197. }
  1198. /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
  1199. n = x.used - 1;
  1200. t = y.used - 1;
  1201. if ((res = pstm_init_size(pool, &q, n - t + 1)) != PSTM_OKAY) {
  1202. goto LBL_Y;
  1203. }
  1204. q.used = n - t + 1;
  1205. /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
  1206. if ((res = pstm_lshd(&y, n - t)) != PSTM_OKAY) { /* y = y*b**{n-t} */
  1207. goto LBL_Q;
  1208. }
  1209. while (pstm_cmp (&x, &y) != PSTM_LT) {
  1210. ++(q.dp[n - t]);
  1211. if ((res = pstm_sub(&x, &y, &x)) != PSTM_OKAY) {
  1212. goto LBL_Q;
  1213. }
  1214. }
  1215. /* reset y by shifting it back down */
  1216. pstm_rshd (&y, n - t);
  1217. /* step 3. for i from n down to (t + 1) */
  1218. for (i = n; i >= (t + 1); i--) {
  1219. if (i > x.used) {
  1220. continue;
  1221. }
  1222. /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
  1223. * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
  1224. if (x.dp[i] == y.dp[t]) {
  1225. q.dp[i - t - 1] = (pstm_digit)((((pstm_word)1) << DIGIT_BIT) - 1);
  1226. } else {
  1227. pstm_word tmp;
  1228. tmp = ((pstm_word) x.dp[i]) << ((pstm_word) DIGIT_BIT);
  1229. tmp |= ((pstm_word) x.dp[i - 1]);
  1230. #if defined(USE_MATRIX_DIV64) && defined(PSTM_32BIT)
  1231. psDiv64(&tmp, y.dp[t]);
  1232. #elif defined(USE_MATRIX_DIV128) && defined(PSTM_64BIT)
  1233. psDiv128(&tmp, y.dp[t]);
  1234. #else
  1235. tmp /= ((pstm_word) y.dp[t]);
  1236. #endif /* USE_MATRIX_DIV64 */
  1237. q.dp[i - t - 1] = (pstm_digit) (tmp);
  1238. }
  1239. /* while (q{i-t-1} * (yt * b + y{t-1})) >
  1240. xi * b**2 + xi-1 * b + xi-2
  1241. do q{i-t-1} -= 1;
  1242. */
  1243. q.dp[i - t - 1] = (q.dp[i - t - 1] + 1);
  1244. do {
  1245. q.dp[i - t - 1] = (q.dp[i - t - 1] - 1);
  1246. /* find left hand */
  1247. pstm_zero (&t1);
  1248. t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
  1249. t1.dp[1] = y.dp[t];
  1250. t1.used = 2;
  1251. if ((res = pstm_mul_d (&t1, q.dp[i - t - 1], &t1)) != PSTM_OKAY) {
  1252. goto LBL_Q;
  1253. }
  1254. /* find right hand */
  1255. t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
  1256. t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
  1257. t2.dp[2] = x.dp[i];
  1258. t2.used = 3;
  1259. } while (pstm_cmp_mag(&t1, &t2) == PSTM_GT);
  1260. /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
  1261. if ((res = pstm_mul_d(&y, q.dp[i - t - 1], &t1)) != PSTM_OKAY) {
  1262. goto LBL_Q;
  1263. }
  1264. if ((res = pstm_lshd(&t1, i - t - 1)) != PSTM_OKAY) {
  1265. goto LBL_Q;
  1266. }
  1267. if ((res = pstm_sub(&x, &t1, &x)) != PSTM_OKAY) {
  1268. goto LBL_Q;
  1269. }
  1270. /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
  1271. if (x.sign == PSTM_NEG) {
  1272. if ((res = pstm_copy(&y, &t1)) != PSTM_OKAY) {
  1273. goto LBL_Q;
  1274. }
  1275. if ((res = pstm_lshd (&t1, i - t - 1)) != PSTM_OKAY) {
  1276. goto LBL_Q;
  1277. }
  1278. if ((res = pstm_add (&x, &t1, &x)) != PSTM_OKAY) {
  1279. goto LBL_Q;
  1280. }
  1281. q.dp[i - t - 1] = q.dp[i - t - 1] - 1;
  1282. }
  1283. }
  1284. /*
  1285. now q is the quotient and x is the remainder (which we have to normalize)
  1286. */
  1287. /* get sign before writing to c */
  1288. x.sign = x.used == 0 ? PSTM_ZPOS : a->sign;
  1289. if (c != NULL) {
  1290. pstm_clamp (&q);
  1291. if (pstm_copy (&q, c) != PSTM_OKAY) {
  1292. res = PS_MEM_FAIL;
  1293. goto LBL_Q;
  1294. }
  1295. c->sign = neg;
  1296. }
  1297. if (d != NULL) {
  1298. if ((res = pstm_div_2d (pool, &x, norm, &x, NULL)) != PSTM_OKAY) {
  1299. goto LBL_Q;
  1300. }
  1301. /*
  1302. the following is a kludge, essentially we were seeing the right
  1303. remainder but with excess digits that should have been zero
  1304. */
  1305. for (i = b->used; i < x.used; i++) {
  1306. x.dp[i] = 0;
  1307. }
  1308. pstm_clamp(&x);
  1309. if (pstm_copy (&x, d) != PSTM_OKAY) {
  1310. res = PS_MEM_FAIL;
  1311. goto LBL_Q;
  1312. }
  1313. }
  1314. res = PSTM_OKAY;
  1315. LBL_Q:pstm_clear (&q);
  1316. LBL_Y:pstm_clear (&y);
  1317. LBL_X:pstm_clear (&x);
  1318. LBL_T2:pstm_clear (&t2);
  1319. LBL_T1:pstm_clear (&t1);
  1320. return res;
  1321. }
  1322. /******************************************************************************/
  1323. /*
  1324. Swap the elements of two integers, for cases where you can't simply swap
  1325. the pstm_int pointers around
  1326. */
  1327. void pstm_exch(pstm_int * a, pstm_int * b)
  1328. {
  1329. pstm_int t;
  1330. t = *a;
  1331. *a = *b;
  1332. *b = t;
  1333. }
  1334. /******************************************************************************/
  1335. /*
  1336. c = a mod b, 0 <= c < b
  1337. */
  1338. int32 pstm_mod(psPool_t *pool, pstm_int *a, pstm_int *b, pstm_int *c)
  1339. {
  1340. pstm_int t;
  1341. int32 err;
  1342. /*
  1343. Smart-size
  1344. */
  1345. if ((err = pstm_init_size(pool, &t, b->alloc)) != PSTM_OKAY) {
  1346. return err;
  1347. }
  1348. if ((err = pstm_div(pool, a, b, NULL, &t)) != PSTM_OKAY) {
  1349. pstm_clear (&t);
  1350. return err;
  1351. }
  1352. if (t.sign != b->sign) {
  1353. err = pstm_add(&t, b, c);
  1354. } else {
  1355. pstm_exch (&t, c);
  1356. }
  1357. pstm_clear (&t);
  1358. return err;
  1359. }
  1360. /******************************************************************************/
  1361. /*
  1362. d = a * b (mod c)
  1363. */
  1364. int32 pstm_mulmod(psPool_t *pool, pstm_int *a, pstm_int *b, pstm_int *c,
  1365. pstm_int *d)
  1366. {
  1367. int32 res;
  1368. int size; //bbox: was int16
  1369. pstm_int tmp;
  1370. /*
  1371. Smart-size pstm_inits. d is an output that is influenced by this local 't'
  1372. so don't shrink 'd' if it wants to becuase this will lead to an pstm_grow
  1373. in RSA operations
  1374. */
  1375. size = a->used + b->used + 1;
  1376. if ((a == d) && (size < a->alloc)) {
  1377. size = a->alloc;
  1378. }
  1379. if ((res = pstm_init_size(pool, &tmp, size)) != PSTM_OKAY) {
  1380. return res;
  1381. }
  1382. if ((res = pstm_mul_comba(pool, a, b, &tmp, NULL, 0)) != PSTM_OKAY) {
  1383. pstm_clear(&tmp);
  1384. return res;
  1385. }
  1386. res = pstm_mod(pool, &tmp, c, d);
  1387. pstm_clear(&tmp);
  1388. return res;
  1389. }
  1390. /******************************************************************************/
  1391. /*
  1392. * y = g**x (mod b)
  1393. * Some restrictions... x must be positive and < b
  1394. */
  1395. int32 pstm_exptmod(psPool_t *pool, pstm_int *G, pstm_int *X, pstm_int *P,
  1396. pstm_int *Y)
  1397. {
  1398. pstm_int M[32], res; /* Keep this winsize based: (1 << max_winsize) */
  1399. pstm_digit buf, mp;
  1400. pstm_digit *paD;
  1401. int32 err, bitbuf;
  1402. int bitcpy, bitcnt, mode, digidx, x, y, winsize; //bbox: was int16
  1403. uint32 paDlen;
  1404. /* set window size from what user set as optimization */
  1405. x = pstm_count_bits(X);
  1406. if (x < 50) {
  1407. winsize = 2;
  1408. } else {
  1409. winsize = PS_EXPTMOD_WINSIZE;
  1410. }
  1411. /* now setup montgomery */
  1412. if ((err = pstm_montgomery_setup (P, &mp)) != PSTM_OKAY) {
  1413. return err;
  1414. }
  1415. /* setup result */
  1416. if ((err = pstm_init_size(pool, &res, (P->used * 2) + 1)) != PSTM_OKAY) {
  1417. return err;
  1418. }
  1419. /*
  1420. create M table
  1421. The M table contains powers of the input base, e.g. M[x] = G^x mod P
  1422. The first half of the table is not computed though except for M[0] and M[1]
  1423. */
  1424. /* now we need R mod m */
  1425. if ((err = pstm_montgomery_calc_normalization (&res, P)) != PSTM_OKAY) {
  1426. goto LBL_RES;
  1427. }
  1428. /*
  1429. init M array
  1430. init first cell
  1431. */
  1432. if ((err = pstm_init_size(pool, &M[1], res.used)) != PSTM_OKAY) {
  1433. goto LBL_RES;
  1434. }
  1435. /* now set M[1] to G * R mod m */
  1436. if (pstm_cmp_mag(P, G) != PSTM_GT) {
  1437. /* G > P so we reduce it first */
  1438. if ((err = pstm_mod(pool, G, P, &M[1])) != PSTM_OKAY) {
  1439. goto LBL_M;
  1440. }
  1441. } else {
  1442. if ((err = pstm_copy(G, &M[1])) != PSTM_OKAY) {
  1443. goto LBL_M;
  1444. }
  1445. }
  1446. if ((err = pstm_mulmod (pool, &M[1], &res, P, &M[1])) != PSTM_OKAY) {
  1447. goto LBL_M;
  1448. }
  1449. /*
  1450. Pre-allocated digit. Used for mul, sqr, AND reduce
  1451. */
  1452. paDlen = ((M[1].used + 3) * 2) * sizeof(pstm_digit);
  1453. paD = xzalloc(paDlen);//bbox
  1454. /*
  1455. compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times
  1456. */
  1457. if (pstm_init_copy(pool, &M[1 << (winsize - 1)], &M[1], 1) != PSTM_OKAY) {
  1458. err = PS_MEM_FAIL;
  1459. goto LBL_PAD;
  1460. }
  1461. for (x = 0; x < (winsize - 1); x++) {
  1462. if ((err = pstm_sqr_comba (pool, &M[1 << (winsize - 1)],
  1463. &M[1 << (winsize - 1)], paD, paDlen)) != PSTM_OKAY) {
  1464. goto LBL_PAD;
  1465. }
  1466. if ((err = pstm_montgomery_reduce(pool, &M[1 << (winsize - 1)], P, mp,
  1467. paD, paDlen)) != PSTM_OKAY) {
  1468. goto LBL_PAD;
  1469. }
  1470. }
  1471. /*
  1472. now init the second half of the array
  1473. */
  1474. for (x = (1<<(winsize-1)) + 1; x < (1 << winsize); x++) {
  1475. if ((err = pstm_init_size(pool, &M[x], M[1<<(winsize-1)].alloc + 1))
  1476. != PSTM_OKAY) {
  1477. for (y = 1<<(winsize-1); y < x; y++) {
  1478. pstm_clear(&M[y]);
  1479. }
  1480. goto LBL_PAD;
  1481. }
  1482. }
  1483. /* create upper table */
  1484. for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
  1485. if ((err = pstm_mul_comba(pool, &M[x - 1], &M[1], &M[x], paD, paDlen))
  1486. != PSTM_OKAY) {
  1487. goto LBL_MARRAY;
  1488. }
  1489. if ((err = pstm_montgomery_reduce(pool, &M[x], P, mp, paD, paDlen)) !=
  1490. PSTM_OKAY) {
  1491. goto LBL_MARRAY;
  1492. }
  1493. }
  1494. /* set initial mode and bit cnt */
  1495. mode = 0;
  1496. bitcnt = 1;
  1497. buf = 0;
  1498. digidx = X->used - 1;
  1499. bitcpy = 0;
  1500. bitbuf = 0;
  1501. for (;;) {
  1502. /* grab next digit as required */
  1503. if (--bitcnt == 0) {
  1504. /* if digidx == -1 we are out of digits so break */
  1505. if (digidx == -1) {
  1506. break;
  1507. }
  1508. /* read next digit and reset bitcnt */
  1509. buf = X->dp[digidx--];
  1510. bitcnt = (int32)DIGIT_BIT;
  1511. }
  1512. /* grab the next msb from the exponent */
  1513. y = (pstm_digit)(buf >> (DIGIT_BIT - 1)) & 1;
  1514. buf <<= (pstm_digit)1;
  1515. /*
  1516. If the bit is zero and mode == 0 then we ignore it.
  1517. These represent the leading zero bits before the first 1 bit
  1518. in the exponent. Technically this opt is not required but it
  1519. does lower the # of trivial squaring/reductions used
  1520. */
  1521. if (mode == 0 && y == 0) {
  1522. continue;
  1523. }
  1524. /* if the bit is zero and mode == 1 then we square */
  1525. if (mode == 1 && y == 0) {
  1526. if ((err = pstm_sqr_comba(pool, &res, &res, paD, paDlen)) !=
  1527. PSTM_OKAY) {
  1528. goto LBL_MARRAY;
  1529. }
  1530. if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen))
  1531. != PSTM_OKAY) {
  1532. goto LBL_MARRAY;
  1533. }
  1534. continue;
  1535. }
  1536. /* else we add it to the window */
  1537. bitbuf |= (y << (winsize - ++bitcpy));
  1538. mode = 2;
  1539. if (bitcpy == winsize) {
  1540. /* ok window is filled so square as required and mul square first */
  1541. for (x = 0; x < winsize; x++) {
  1542. if ((err = pstm_sqr_comba(pool, &res, &res, paD, paDlen)) !=
  1543. PSTM_OKAY) {
  1544. goto LBL_MARRAY;
  1545. }
  1546. if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD,
  1547. paDlen)) != PSTM_OKAY) {
  1548. goto LBL_MARRAY;
  1549. }
  1550. }
  1551. /* then multiply */
  1552. if ((err = pstm_mul_comba(pool, &res, &M[bitbuf], &res, paD,
  1553. paDlen)) != PSTM_OKAY) {
  1554. goto LBL_MARRAY;
  1555. }
  1556. if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen))
  1557. != PSTM_OKAY) {
  1558. goto LBL_MARRAY;
  1559. }
  1560. /* empty window and reset */
  1561. bitcpy = 0;
  1562. bitbuf = 0;
  1563. mode = 1;
  1564. }
  1565. }
  1566. /* if bits remain then square/multiply */
  1567. if (mode == 2 && bitcpy > 0) {
  1568. /* square then multiply if the bit is set */
  1569. for (x = 0; x < bitcpy; x++) {
  1570. if ((err = pstm_sqr_comba(pool, &res, &res, paD, paDlen)) !=
  1571. PSTM_OKAY) {
  1572. goto LBL_MARRAY;
  1573. }
  1574. if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen))
  1575. != PSTM_OKAY) {
  1576. goto LBL_MARRAY;
  1577. }
  1578. /* get next bit of the window */
  1579. bitbuf <<= 1;
  1580. if ((bitbuf & (1 << winsize)) != 0) {
  1581. /* then multiply */
  1582. if ((err = pstm_mul_comba(pool, &res, &M[1], &res, paD, paDlen))
  1583. != PSTM_OKAY) {
  1584. goto LBL_MARRAY;
  1585. }
  1586. if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD,
  1587. paDlen)) != PSTM_OKAY) {
  1588. goto LBL_MARRAY;
  1589. }
  1590. }
  1591. }
  1592. }
  1593. /*
  1594. Fix up result if Montgomery reduction is used recall that any value in a
  1595. Montgomery system is actually multiplied by R mod n. So we have to reduce
  1596. one more time to cancel out the factor of R.
  1597. */
  1598. if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen)) !=
  1599. PSTM_OKAY) {
  1600. goto LBL_MARRAY;
  1601. }
  1602. /* swap res with Y */
  1603. if ((err = pstm_copy (&res, Y)) != PSTM_OKAY) {
  1604. goto LBL_MARRAY;
  1605. }
  1606. err = PSTM_OKAY;
  1607. LBL_MARRAY:
  1608. for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
  1609. pstm_clear(&M[x]);
  1610. }
  1611. LBL_PAD:psFree(paD, pool);
  1612. LBL_M: pstm_clear(&M[1]);
  1613. LBL_RES:pstm_clear(&res);
  1614. return err;
  1615. }
  1616. /******************************************************************************/
  1617. /*
  1618. */
  1619. int32 pstm_add(pstm_int *a, pstm_int *b, pstm_int *c)
  1620. {
  1621. int32 res;
  1622. int sa, sb; //bbox: was int16
  1623. /* get sign of both inputs */
  1624. sa = a->sign;
  1625. sb = b->sign;
  1626. /* handle two cases, not four */
  1627. if (sa == sb) {
  1628. /* both positive or both negative, add their mags, copy the sign */
  1629. c->sign = sa;
  1630. if ((res = s_pstm_add (a, b, c)) != PSTM_OKAY) {
  1631. return res;
  1632. }
  1633. } else {
  1634. /*
  1635. one positive, the other negative
  1636. subtract the one with the greater magnitude from the one of the lesser
  1637. magnitude. The result gets the sign of the one with the greater mag.
  1638. */
  1639. if (pstm_cmp_mag (a, b) == PSTM_LT) {
  1640. c->sign = sb;
  1641. if ((res = s_pstm_sub (b, a, c)) != PSTM_OKAY) {
  1642. return res;
  1643. }
  1644. } else {
  1645. c->sign = sa;
  1646. if ((res = s_pstm_sub (a, b, c)) != PSTM_OKAY) {
  1647. return res;
  1648. }
  1649. }
  1650. }
  1651. return PS_SUCCESS;
  1652. }
  1653. /******************************************************************************/
  1654. /*
  1655. reverse an array, used for radix code
  1656. */
  1657. static void pstm_reverse (unsigned char *s, int len) //bbox: was int16 len
  1658. {
  1659. int32 ix, iy;
  1660. unsigned char t;
  1661. ix = 0;
  1662. iy = len - 1;
  1663. while (ix < iy) {
  1664. t = s[ix];
  1665. s[ix] = s[iy];
  1666. s[iy] = t;
  1667. ++ix;
  1668. --iy;
  1669. }
  1670. }
  1671. /******************************************************************************/
  1672. /*
  1673. No reverse. Useful in some of the EIP-154 PKA stuff where special byte
  1674. order seems to come into play more often
  1675. */
  1676. int32 pstm_to_unsigned_bin_nr(psPool_t *pool, pstm_int *a, unsigned char *b)
  1677. {
  1678. int32 res;
  1679. int x; //bbox: was int16
  1680. pstm_int t = { 0 };
  1681. if ((res = pstm_init_copy(pool, &t, a, 0)) != PSTM_OKAY) {
  1682. return res;
  1683. }
  1684. x = 0;
  1685. while (pstm_iszero (&t) == 0) {
  1686. b[x++] = (unsigned char) (t.dp[0] & 255);
  1687. if ((res = pstm_div_2d (pool, &t, 8, &t, NULL)) != PSTM_OKAY) {
  1688. pstm_clear(&t);
  1689. return res;
  1690. }
  1691. }
  1692. pstm_clear(&t);
  1693. return PS_SUCCESS;
  1694. }
  1695. /******************************************************************************/
  1696. /*
  1697. */
  1698. int32 pstm_to_unsigned_bin(psPool_t *pool, pstm_int *a, unsigned char *b)
  1699. {
  1700. int32 res;
  1701. int x; //bbox: was int16
  1702. pstm_int t = { 0 };
  1703. if ((res = pstm_init_copy(pool, &t, a, 0)) != PSTM_OKAY) {
  1704. return res;
  1705. }
  1706. x = 0;
  1707. while (pstm_iszero (&t) == 0) {
  1708. b[x++] = (unsigned char) (t.dp[0] & 255);
  1709. if ((res = pstm_div_2d (pool, &t, 8, &t, NULL)) != PSTM_OKAY) {
  1710. pstm_clear(&t);
  1711. return res;
  1712. }
  1713. }
  1714. pstm_reverse (b, x);
  1715. pstm_clear(&t);
  1716. return PS_SUCCESS;
  1717. }
  1718. /******************************************************************************/
  1719. /*
  1720. compare against a single digit
  1721. */
  1722. int32 pstm_cmp_d(pstm_int *a, pstm_digit b)
  1723. {
  1724. /* compare based on sign */
  1725. if ((b && a->used == 0) || a->sign == PSTM_NEG) {
  1726. return PSTM_LT;
  1727. }
  1728. /* compare based on magnitude */
  1729. if (a->used > 1) {
  1730. return PSTM_GT;
  1731. }
  1732. /* compare the only digit of a to b */
  1733. if (a->dp[0] > b) {
  1734. return PSTM_GT;
  1735. } else if (a->dp[0] < b) {
  1736. return PSTM_LT;
  1737. } else {
  1738. return PSTM_EQ;
  1739. }
  1740. }
  1741. /*
  1742. Need invmod for ECC and also private key loading for hardware crypto
  1743. in cases where dQ > dP. The values must be switched and a new qP must be
  1744. calculated using this function
  1745. */
  1746. //bbox: pool unused
  1747. #define pstm_invmod_slow(pool, a, b, c) \
  1748. pstm_invmod_slow( a, b, c)
  1749. static int32 pstm_invmod_slow(psPool_t *pool, pstm_int * a, pstm_int * b,
  1750. pstm_int * c)
  1751. {
  1752. pstm_int x, y, u, v, A, B, C, D;
  1753. int32 res;
  1754. /* b cannot be negative */
  1755. if (b->sign == PSTM_NEG || pstm_iszero(b) == 1) {
  1756. return PS_LIMIT_FAIL;
  1757. }
  1758. /* init temps */
  1759. if (pstm_init_size(pool, &x, b->used) != PSTM_OKAY) {
  1760. return PS_MEM_FAIL;
  1761. }
  1762. /* x = a, y = b */
  1763. if ((res = pstm_mod(pool, a, b, &x)) != PSTM_OKAY) {
  1764. goto LBL_X;
  1765. }
  1766. if (pstm_init_copy(pool, &y, b, 0) != PSTM_OKAY) {
  1767. goto LBL_X;
  1768. }
  1769. /* 2. [modified] if x,y are both even then return an error! */
  1770. if (pstm_iseven (&x) == 1 && pstm_iseven (&y) == 1) {
  1771. res = PS_FAILURE;
  1772. goto LBL_Y;
  1773. }
  1774. /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
  1775. if ((res = pstm_init_copy(pool, &u, &x, 0)) != PSTM_OKAY) {
  1776. goto LBL_Y;
  1777. }
  1778. if ((res = pstm_init_copy(pool, &v, &y, 0)) != PSTM_OKAY) {
  1779. goto LBL_U;
  1780. }
  1781. if ((res = pstm_init_size(pool, &A, sizeof(pstm_digit))) != PSTM_OKAY) {
  1782. goto LBL_V;
  1783. }
  1784. if ((res = pstm_init_size(pool, &D, sizeof(pstm_digit))) != PSTM_OKAY) {
  1785. goto LBL_A;
  1786. }
  1787. pstm_set (&A, 1);
  1788. pstm_set (&D, 1);
  1789. if ((res = pstm_init(pool, &B)) != PSTM_OKAY) {
  1790. goto LBL_D;
  1791. }
  1792. if ((res = pstm_init(pool, &C)) != PSTM_OKAY) {
  1793. goto LBL_B;
  1794. }
  1795. top:
  1796. /* 4. while u is even do */
  1797. while (pstm_iseven (&u) == 1) {
  1798. /* 4.1 u = u/2 */
  1799. if ((res = pstm_div_2 (&u, &u)) != PSTM_OKAY) {
  1800. goto LBL_C;
  1801. }
  1802. /* 4.2 if A or B is odd then */
  1803. if (pstm_isodd (&A) == 1 || pstm_isodd (&B) == 1) {
  1804. /* A = (A+y)/2, B = (B-x)/2 */
  1805. if ((res = pstm_add (&A, &y, &A)) != PSTM_OKAY) {
  1806. goto LBL_C;
  1807. }
  1808. if ((res = pstm_sub (&B, &x, &B)) != PSTM_OKAY) {
  1809. goto LBL_C;
  1810. }
  1811. }
  1812. /* A = A/2, B = B/2 */
  1813. if ((res = pstm_div_2 (&A, &A)) != PSTM_OKAY) {
  1814. goto LBL_C;
  1815. }
  1816. if ((res = pstm_div_2 (&B, &B)) != PSTM_OKAY) {
  1817. goto LBL_C;
  1818. }
  1819. }
  1820. /* 5. while v is even do */
  1821. while (pstm_iseven (&v) == 1) {
  1822. /* 5.1 v = v/2 */
  1823. if ((res = pstm_div_2 (&v, &v)) != PSTM_OKAY) {
  1824. goto LBL_C;
  1825. }
  1826. /* 5.2 if C or D is odd then */
  1827. if (pstm_isodd (&C) == 1 || pstm_isodd (&D) == 1) {
  1828. /* C = (C+y)/2, D = (D-x)/2 */
  1829. if ((res = pstm_add (&C, &y, &C)) != PSTM_OKAY) {
  1830. goto LBL_C;
  1831. }
  1832. if ((res = pstm_sub (&D, &x, &D)) != PSTM_OKAY) {
  1833. goto LBL_C;
  1834. }
  1835. }
  1836. /* C = C/2, D = D/2 */
  1837. if ((res = pstm_div_2 (&C, &C)) != PSTM_OKAY) {
  1838. goto LBL_C;
  1839. }
  1840. if ((res = pstm_div_2 (&D, &D)) != PSTM_OKAY) {
  1841. goto LBL_C;
  1842. }
  1843. }
  1844. /* 6. if u >= v then */
  1845. if (pstm_cmp (&u, &v) != PSTM_LT) {
  1846. /* u = u - v, A = A - C, B = B - D */
  1847. if ((res = pstm_sub (&u, &v, &u)) != PSTM_OKAY) {
  1848. goto LBL_C;
  1849. }
  1850. if ((res = pstm_sub (&A, &C, &A)) != PSTM_OKAY) {
  1851. goto LBL_C;
  1852. }
  1853. if ((res = pstm_sub (&B, &D, &B)) != PSTM_OKAY) {
  1854. goto LBL_C;
  1855. }
  1856. } else {
  1857. /* v - v - u, C = C - A, D = D - B */
  1858. if ((res = pstm_sub (&v, &u, &v)) != PSTM_OKAY) {
  1859. goto LBL_C;
  1860. }
  1861. if ((res = pstm_sub (&C, &A, &C)) != PSTM_OKAY) {
  1862. goto LBL_C;
  1863. }
  1864. if ((res = pstm_sub (&D, &B, &D)) != PSTM_OKAY) {
  1865. goto LBL_C;
  1866. }
  1867. }
  1868. /* if not zero goto step 4 */
  1869. if (pstm_iszero (&u) == 0)
  1870. goto top;
  1871. /* now a = C, b = D, gcd == g*v */
  1872. /* if v != 1 then there is no inverse */
  1873. if (pstm_cmp_d (&v, 1) != PSTM_EQ) {
  1874. res = PS_FAILURE;
  1875. goto LBL_C;
  1876. }
  1877. /* if its too low */
  1878. while (pstm_cmp_d(&C, 0) == PSTM_LT) {
  1879. if ((res = pstm_add(&C, b, &C)) != PSTM_OKAY) {
  1880. goto LBL_C;
  1881. }
  1882. }
  1883. /* too big */
  1884. while (pstm_cmp_mag(&C, b) != PSTM_LT) {
  1885. if ((res = pstm_sub(&C, b, &C)) != PSTM_OKAY) {
  1886. goto LBL_C;
  1887. }
  1888. }
  1889. /* C is now the inverse */
  1890. if ((res = pstm_copy(&C, c)) != PSTM_OKAY) {
  1891. goto LBL_C;
  1892. }
  1893. res = PSTM_OKAY;
  1894. LBL_C: pstm_clear(&C);
  1895. LBL_D: pstm_clear(&D);
  1896. LBL_B: pstm_clear(&B);
  1897. LBL_A: pstm_clear(&A);
  1898. LBL_V: pstm_clear(&v);
  1899. LBL_U: pstm_clear(&u);
  1900. LBL_Y: pstm_clear(&y);
  1901. LBL_X: pstm_clear(&x);
  1902. return res;
  1903. }
  1904. /* c = 1/a (mod b) for odd b only */
  1905. int32 pstm_invmod(psPool_t *pool, pstm_int *a, pstm_int *b, pstm_int *c)
  1906. {
  1907. pstm_int x, y, u, v, B, D;
  1908. int32 res;
  1909. int neg, sanity; //bbox: was uint16
  1910. /* 2. [modified] b must be odd */
  1911. if (pstm_iseven (b) == 1) {
  1912. return pstm_invmod_slow(pool, a,b,c);
  1913. }
  1914. /* x == modulus, y == value to invert */
  1915. if ((res = pstm_init_copy(pool, &x, b, 0)) != PSTM_OKAY) {
  1916. return res;
  1917. }
  1918. if ((res = pstm_init_size(pool, &y, a->alloc)) != PSTM_OKAY) {
  1919. goto LBL_X;
  1920. }
  1921. /* we need y = |a| */
  1922. pstm_abs(a, &y);
  1923. /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
  1924. if ((res = pstm_init_copy(pool, &u, &x, 0)) != PSTM_OKAY) {
  1925. goto LBL_Y;
  1926. }
  1927. if ((res = pstm_init_copy(pool, &v, &y, 0)) != PSTM_OKAY) {
  1928. goto LBL_U;
  1929. }
  1930. if ((res = pstm_init(pool, &B)) != PSTM_OKAY) {
  1931. goto LBL_V;
  1932. }
  1933. if ((res = pstm_init(pool, &D)) != PSTM_OKAY) {
  1934. goto LBL_B;
  1935. }
  1936. pstm_set (&D, 1);
  1937. sanity = 0;
  1938. top:
  1939. /* 4. while u is even do */
  1940. while (pstm_iseven (&u) == 1) {
  1941. /* 4.1 u = u/2 */
  1942. if ((res = pstm_div_2 (&u, &u)) != PSTM_OKAY) {
  1943. goto LBL_D;
  1944. }
  1945. /* 4.2 if B is odd then */
  1946. if (pstm_isodd (&B) == 1) {
  1947. if ((res = pstm_sub (&B, &x, &B)) != PSTM_OKAY) {
  1948. goto LBL_D;
  1949. }
  1950. }
  1951. /* B = B/2 */
  1952. if ((res = pstm_div_2 (&B, &B)) != PSTM_OKAY) {
  1953. goto LBL_D;
  1954. }
  1955. }
  1956. /* 5. while v is even do */
  1957. while (pstm_iseven (&v) == 1) {
  1958. /* 5.1 v = v/2 */
  1959. if ((res = pstm_div_2 (&v, &v)) != PSTM_OKAY) {
  1960. goto LBL_D;
  1961. }
  1962. /* 5.2 if D is odd then */
  1963. if (pstm_isodd (&D) == 1) {
  1964. /* D = (D-x)/2 */
  1965. if ((res = pstm_sub (&D, &x, &D)) != PSTM_OKAY) {
  1966. goto LBL_D;
  1967. }
  1968. }
  1969. /* D = D/2 */
  1970. if ((res = pstm_div_2 (&D, &D)) != PSTM_OKAY) {
  1971. goto LBL_D;
  1972. }
  1973. }
  1974. /* 6. if u >= v then */
  1975. if (pstm_cmp (&u, &v) != PSTM_LT) {
  1976. /* u = u - v, B = B - D */
  1977. if ((res = pstm_sub (&u, &v, &u)) != PSTM_OKAY) {
  1978. goto LBL_D;
  1979. }
  1980. if ((res = pstm_sub (&B, &D, &B)) != PSTM_OKAY) {
  1981. goto LBL_D;
  1982. }
  1983. } else {
  1984. /* v - v - u, D = D - B */
  1985. if ((res = pstm_sub (&v, &u, &v)) != PSTM_OKAY) {
  1986. goto LBL_D;
  1987. }
  1988. if ((res = pstm_sub (&D, &B, &D)) != PSTM_OKAY) {
  1989. goto LBL_D;
  1990. }
  1991. }
  1992. /* if not zero goto step 4 */
  1993. if (sanity++ > 1000) {
  1994. res = PS_LIMIT_FAIL;
  1995. goto LBL_D;
  1996. }
  1997. if (pstm_iszero (&u) == 0) {
  1998. goto top;
  1999. }
  2000. /* now a = C, b = D, gcd == g*v */
  2001. /* if v != 1 then there is no inverse */
  2002. if (pstm_cmp_d (&v, 1) != PSTM_EQ) {
  2003. res = PS_FAILURE;
  2004. goto LBL_D;
  2005. }
  2006. /* b is now the inverse */
  2007. neg = a->sign;
  2008. while (D.sign == PSTM_NEG) {
  2009. if ((res = pstm_add (&D, b, &D)) != PSTM_OKAY) {
  2010. goto LBL_D;
  2011. }
  2012. }
  2013. if ((res = pstm_copy (&D, c)) != PSTM_OKAY) {
  2014. goto LBL_D;
  2015. }
  2016. c->sign = neg;
  2017. res = PSTM_OKAY;
  2018. LBL_D: pstm_clear(&D);
  2019. LBL_B: pstm_clear(&B);
  2020. LBL_V: pstm_clear(&v);
  2021. LBL_U: pstm_clear(&u);
  2022. LBL_Y: pstm_clear(&y);
  2023. LBL_X: pstm_clear(&x);
  2024. return res;
  2025. }
  2026. #endif /* !DISABLE_PSTM */
  2027. /******************************************************************************/