arith03.c 52 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692
  1. /* arith03.c Copyright (C) 1990-96 Codemist Ltd */
  2. /*
  3. * Arithmetic functions.
  4. * quotient.
  5. *
  6. */
  7. /* Signature: 7f5d27b7 31-May-1997 */
  8. #include <stdarg.h>
  9. #include <string.h>
  10. #include <ctype.h>
  11. #include <math.h>
  12. #include "machine.h"
  13. #include "tags.h"
  14. #include "cslerror.h"
  15. #include "externs.h"
  16. #include "arith.h"
  17. #ifdef TIMEOUT
  18. #include "timeout.h"
  19. #endif
  20. /*
  21. * Division
  22. */
  23. #ifndef IDIVIDE
  24. #ifdef MULDIV64
  25. unsigned32 Idivide(unsigned32 *qp, unsigned32 a, unsigned32 b, unsigned32 c)
  26. /*
  27. * *qp = (a,b) / c, return the remainder
  28. *
  29. * The double-length value (a,b) is given as a 62-bit positive number. Its
  30. * low order 31 bits are in b, and the 0x80000000 bit of b is zero. The
  31. * high order part (in a) is also represented with the 0x80000000 bit zero.
  32. * The divisor c is a positive value that is at most 0x7fffffff, and the
  33. * procedure will only be called when the correct quotient is in the
  34. * (inclusive) range 0 to 0x7fffffff. The above constraints can be thought
  35. * of in two ways: (a) Idivide used a 31-bit word and is working on
  36. * unsigned values. The 31-bit quantities happen to be being passed around
  37. * in 32 bit registers, but the top bit of those registers will never be used
  38. * and will contain zero, or (b) the range of values used is such that
  39. * a 64 by 32-bit division can be performed, and by constraining the range
  40. * of values this 64 by 32-bit division only has to handle positive numbers,
  41. * but depending on the hardware of your computer you are entitled to use
  42. * either signed or unsigned arithmetic.
  43. */
  44. {
  45. unsigned64 p = ((unsigned64)a << 31) | (unsigned64)b;
  46. *qp = (unsigned32)(p / (unsigned64)c);
  47. return (unsigned32)(p % (unsigned64)c);
  48. }
  49. #else
  50. unsigned32 Idivide(unsigned32 *qp, unsigned32 a, unsigned32 b, unsigned32 c)
  51. /*
  52. * *qp = (a,b) / c, return the remainder
  53. *
  54. * The double-length value (a,b) is given as a 62-bit positive number. Its
  55. * low order 31 bits are in b, and the 0x80000000 bit of b is zero. The
  56. * high order part (in a) is also represented with the 0x80000000 bit zero.
  57. * The divisor c is a positive value that is at most 0x7fffffff, and the
  58. * procedure will only be called when the correct quotient is in the
  59. * (inclusive) range 0 to 0x7fffffff. The above constraints can be thought
  60. * of in two ways: (a) Idivide used a 31-bit word and is working on
  61. * unsigned values. The 31-bit quantities happen to be being passed around
  62. * in 32 bit registers, but the top bit of those registers will never be used
  63. * and will contain zero, or (b) the range of values used is such that
  64. * a 64 by 32-bit division can be performed, and by constraining the range
  65. * of values this 64 by 32-bit division only has to handle positive numbers,
  66. * but depending on the hardware of your computer you are entitled to use
  67. * either signed or unsigned arithmetic.
  68. */
  69. {
  70. int i;
  71. unsigned32 q = 0;
  72. /*
  73. * I have a loop that needs to be executed 32 times, with just
  74. * slightly different behaviour the last time around. Since it is
  75. * fairly short, I unwind it three-fold into a loop executed 10 times
  76. * plus a postlude. I also do a test at the start that decides if the
  77. * quotient will be small (less than about 16 bits) an in that case
  78. * loop rather less often - my reasoning is that a plausible distribution
  79. * of quotients is exponential so the short route will be taken about
  80. * half of the time, and will save almost half of the work done here at the
  81. * cost of a not-too-expensive extra test to begin with.
  82. */
  83. if (a < (c >> 15))
  84. { a = (a << 15) | (b >> 16);
  85. b = (b << 15) & 0x7fffffff;
  86. i = 5;
  87. }
  88. else i = 0;
  89. do
  90. { q = q << 3; /* accumulate quotient 3 bits at a time */
  91. if (a >= c) /* bit 1 of 3... */
  92. { a = a - c;
  93. q = q + 4;
  94. }
  95. a = a << 1; /* shift (a,b) doubleword left 1 bit */
  96. b = b << 1;
  97. if (((int32)b) < 0) a = a + 1;
  98. if (a >= c) /* bit 2 of 3... */
  99. { a = a - c;
  100. q = q + 2;
  101. }
  102. a = a << 1;
  103. b = b << 1;
  104. if (((int32)b) < 0) a = a + 1;
  105. if (a >= c) /* bit 3 of 3... */
  106. { a = a - c;
  107. q = q + 1;
  108. }
  109. a = a << 1;
  110. b = b << 1;
  111. if (((int32)b) < 0) a = a + 1;
  112. i++;
  113. } while (i < 10);
  114. q = q << 2; /* 2 more bits to be included */
  115. if (a >= c)
  116. { a = a - c;
  117. q = q + 2;
  118. }
  119. a = a << 1;
  120. b = b << 1;
  121. if (((int32)b) < 0) a = a + 1;
  122. if (a >= c) /* the final bit of the quotient */
  123. { a = a - c; /* leave remainder in a, b should be zero */
  124. q = q + 1;
  125. }
  126. *qp = q;
  127. return a;
  128. }
  129. #endif
  130. #endif /* IDIVIDE */
  131. #ifdef COMMON
  132. static Lisp_Object quotis(Lisp_Object a, Lisp_Object b)
  133. {
  134. Float_union bb;
  135. bb.i = b - TAG_SFLOAT;
  136. if (bb.f == 0.0) return aerror2("bad arg for quotient", a, b);
  137. bb.f = (float) ((float)int_of_fixnum(a) / bb.f);
  138. return (bb.i & ~(int32)0xf) + TAG_SFLOAT;
  139. }
  140. #endif
  141. static Lisp_Object quotib(Lisp_Object a, Lisp_Object b)
  142. {
  143. /*
  144. * a fixnum divided by a bignum always gives 0, regardless of signs etc.,
  145. * save in the one case of (-#x8000000)/(#x8000000) which must yield -1
  146. */
  147. if (int_of_fixnum(a) == fix_mask && bignum_length(b) == 8 &&
  148. bignum_digits(b)[0] == 0x08000000) return fixnum_of_int(-1);
  149. else return fixnum_of_int(0);
  150. }
  151. #ifdef COMMON
  152. static Lisp_Object CLquotib(Lisp_Object a, Lisp_Object b)
  153. {
  154. Lisp_Object g, nil = C_nil;
  155. bool w;
  156. push2(a, b);
  157. w = minusp(b);
  158. errexitn(2);
  159. g = gcd(stack[0], stack[-1]);
  160. errexitn(2);
  161. if (w)
  162. { g = negate(g);
  163. errexitn(2);
  164. }
  165. a = stack[-1];
  166. push(g);
  167. a = quot2(a, g);
  168. errexitn(3);
  169. pop2(g, b);
  170. stack[0] = a;
  171. b = quot2(b, g);
  172. pop(a);
  173. errexit();
  174. return make_ratio(a, b);
  175. }
  176. static Lisp_Object CLquotbi(Lisp_Object a, Lisp_Object b)
  177. {
  178. return CLquotib(a, b);
  179. }
  180. static Lisp_Object CLquotbb(Lisp_Object a, Lisp_Object b)
  181. {
  182. return CLquotib(a, b);
  183. }
  184. static Lisp_Object quotir(Lisp_Object a, Lisp_Object b)
  185. {
  186. Lisp_Object w, nil;
  187. if (a == fixnum_of_int(0)) return a;
  188. push3(b, a, C_nil);
  189. #define g stack[0]
  190. #define a stack[-1]
  191. #define b stack[-2]
  192. g = gcd(a, numerator(b));
  193. /*
  194. * the gcd() function guarantees to hand back a positive result.
  195. */
  196. nil = C_nil;
  197. if (exception_pending()) goto fail;
  198. a = quot2(a, g);
  199. nil = C_nil;
  200. if (exception_pending()) goto fail;
  201. w = minusp(numerator(b));
  202. nil = C_nil;
  203. if (exception_pending()) goto fail;
  204. if (w)
  205. { g = negate(g);
  206. nil = C_nil;
  207. if (exception_pending()) goto fail;
  208. }
  209. g = quot2(numerator(b), g); /* denominator of result will be +ve */
  210. nil = C_nil;
  211. if (exception_pending()) goto fail;
  212. a = times2(a, denominator(b));
  213. nil = C_nil;
  214. if (exception_pending()) goto fail;
  215. w = make_ratio(a, g);
  216. popv(3);
  217. return w;
  218. fail:
  219. popv(3);
  220. return nil;
  221. #undef a
  222. #undef b
  223. #undef g
  224. }
  225. static Lisp_Object quotic(Lisp_Object a, Lisp_Object b)
  226. /*
  227. * Used for all cases of xxx/<p+iq>. This is coded in a fairly naive
  228. * way, multiplying both numerator and denominator of what will end up
  229. * as the result by the complex conjugate of b. If floating point
  230. * arithmetic is used this can lead to grossly premature overflow. For
  231. * /* the moment I will ignore that miserable fact
  232. */
  233. {
  234. Lisp_Object u, v, nil;
  235. push2(a, b);
  236. #define b stack[0]
  237. #define a stack[-1]
  238. /*
  239. * a / (p + iq) is computed as follows:
  240. * (a * (p - iq)) / (p^2 + q^2)
  241. */
  242. u = negate(imag_part(b));
  243. nil = C_nil;
  244. if (exception_pending()) goto fail;
  245. u = make_complex(real_part(b), u);
  246. nil = C_nil;
  247. if (exception_pending()) goto fail;
  248. a = times2(a, u);
  249. u = real_part(b);
  250. u = times2(u, u);
  251. nil = C_nil;
  252. if (exception_pending()) goto fail;
  253. v = imag_part(b);
  254. b = u;
  255. u = times2(v, v);
  256. nil = C_nil;
  257. if (exception_pending()) goto fail;
  258. u = plus2(u, b);
  259. nil = C_nil;
  260. if (exception_pending()) goto fail;
  261. v = a;
  262. popv(2);
  263. return quot2(v, u);
  264. #undef a
  265. #undef b
  266. fail:
  267. popv(2);
  268. return nil;
  269. }
  270. #endif
  271. static Lisp_Object quotif(Lisp_Object a, Lisp_Object b)
  272. {
  273. double d = float_of_number(b);
  274. if (d == 0.0) return aerror2("bad arg for quotient", a, b);
  275. d = (double)int_of_fixnum(a) / d;
  276. return make_boxfloat(d, type_of_header(flthdr(b)));
  277. }
  278. #ifdef COMMON
  279. static Lisp_Object quotsi(Lisp_Object a, Lisp_Object b)
  280. {
  281. Float_union aa;
  282. if (b == fixnum_of_int(0)) return aerror2("bad arg for quotient", a, b);
  283. aa.i = a - TAG_SFLOAT;
  284. aa.f = (float) (aa.f / (float)int_of_fixnum(b));
  285. return (aa.i & ~(int32)0xf) + TAG_SFLOAT;
  286. }
  287. static Lisp_Object quotsb(Lisp_Object a, Lisp_Object b)
  288. {
  289. double d = float_of_number(b);
  290. if (d == 0.0) return aerror2("bad arg for quotient", a, b);
  291. d = float_of_number(a) / d;
  292. return make_sfloat(d);
  293. }
  294. #define quotsr(a, b) quotsb(a, b)
  295. #define quotsc(a, b) quotic(a, b)
  296. #endif
  297. static Lisp_Object quotsf(Lisp_Object a, Lisp_Object b)
  298. {
  299. double d = float_of_number(b);
  300. if (d == 0.0) return aerror2("bad arg for quotient", a, b);
  301. d = float_of_number(a) / d;
  302. return make_boxfloat(d, type_of_header(flthdr(b)));
  303. }
  304. Lisp_Object quotbn(Lisp_Object a, int32 n)
  305. /*
  306. * Divide a bignum by an integer, where the integer is (by now)
  307. * a natural C int32 but limited to 31 not 32 bits active. I.e.
  308. * this can be used when dividing by a fixnum or by dividing by
  309. * a one-word bignum. I will not call this with n=-1, which would
  310. * be the one case where it could cause the quotient to be bigger
  311. * than the dividend. Leaves nwork set to the remainder.
  312. */
  313. {
  314. Lisp_Object nil;
  315. int sign;
  316. int32 lena = (bignum_length(a)>>2)-2, i, lenc, lenx;
  317. unsigned32 carry;
  318. if (lena == 0) /* one-word bignum as numerator */
  319. { int32 p = (int32)bignum_digits(a)[0], w;
  320. nil = C_nil; /* So I can access nwork */
  321. nwork = p % n;
  322. /*
  323. * C does not define what happens on non-exact division involving
  324. * negative quantities, so I adjust things here so that the remainder
  325. * has the sign that I want and the division I do is exact.
  326. */
  327. if (p < 0)
  328. { if (nwork > 0) nwork -= n;
  329. }
  330. else if (nwork < 0) nwork += n;
  331. /* Remainder should be correct (with correct sign) by now, regardless */
  332. p = (p - nwork) / n;
  333. w = p & fix_mask;
  334. if (w == 0 || w == fix_mask) return fixnum_of_int(p);
  335. else return make_one_word_bignum(p);
  336. }
  337. else if (lena == 1)
  338. /*
  339. * I treat division of a 2-word bignum by a fixnum or 1-word bignum as
  340. * a special case since it can lead to a fixnum result - if the divisor is
  341. * just one word long and the dividend is 3 or more words I would
  342. * certainly have a bignum result. Thus by separating off the code here
  343. * I (a) remove the need for test relating to big- to fixnum conversion
  344. * later on and (b) avoid allocating heap in this tolerably common case
  345. * when sometimes I will not need to.
  346. */
  347. { unsigned32 a0 = bignum_digits(a)[0];
  348. int32 a1 = (int32)bignum_digits(a)[1];
  349. if (a1 < 0)
  350. { sign = 3;
  351. if (a0 == 0) a1 = -a1;
  352. else
  353. { a0 = clear_top_bit(-(int32)a0);
  354. a1 = ~a1;
  355. }
  356. }
  357. else sign = 0;
  358. if (n < 0) sign ^= 1, n = -n;
  359. if (a1 < n) /* result can be calculated in one word - good! */
  360. { unsigned32 q, r;
  361. Ddivide(r, q, (unsigned32)a1, a0, n);
  362. if ((sign & 2) != 0) r = -(int32)r;
  363. nil = C_nil;
  364. nwork = r;
  365. a0 = q;
  366. if (a0 == 0) return fixnum_of_int(0);
  367. if ((sign & 1) == 0)
  368. { if ((a0 & fix_mask) == 0) return fixnum_of_int(a0);
  369. else if ((a0 & 0x40000000) == 0)
  370. return make_one_word_bignum(a0);
  371. else return make_two_word_bignum(0, a0);
  372. }
  373. a0 = -(int32)a0;
  374. if ((a0 & fix_mask) == fix_mask) return fixnum_of_int(a0);
  375. else if ((a0 & 0xc0000000U) == 0xc0000000U)
  376. return make_one_word_bignum(a0);
  377. else return make_two_word_bignum(-1, clear_top_bit(a0));
  378. }
  379. /*
  380. * here the quotient will involve two digits.
  381. */
  382. else
  383. { unsigned32 q1, q0, r;
  384. Ddivide(r, q1, 0, (unsigned32)a1, n);
  385. Ddivide(r, q0, r, a0, n);
  386. if ((sign & 2) != 0) r = -(int32)r;
  387. nil = C_nil;
  388. nwork = r;
  389. if ((sign & 1) != 0)
  390. { if (q0 == 0) q1 = -(int32)q1;
  391. else
  392. { q1 = ~q1;
  393. q0 = clear_top_bit(-(int32)q0);
  394. }
  395. }
  396. if (q1 == 0 && (q0 & 0x40000000) == 0)
  397. return make_one_word_bignum(q0);
  398. return make_two_word_bignum(q1, q0);
  399. }
  400. }
  401. /*
  402. * That has dealt with the special cases - now the input is a bignum
  403. * with at least 3 digits and the quotient will certainly be a bignum.
  404. * Start by allocating a workspace copy of the dividend. negateb will
  405. * leave a a bignum, although it may change its length.
  406. */
  407. if ((int32)bignum_digits(a)[lena] < 0) a = negateb(a), sign = 3;
  408. else a = copyb(a), sign = 0;
  409. errexit();
  410. if (n < 0)
  411. { sign ^= 1;
  412. n = -n;
  413. }
  414. /*
  415. * Now both a and n have been made positive, and their original signs
  416. * have been recorded so that I can tidy up at the end. The 1 bit in sign
  417. * tells me what sign the quotient will have, the 2 bit gives the sign
  418. * for the remainder.
  419. */
  420. lena = (bignum_length(a)>>2)-1;
  421. carry = 0;
  422. for (i=lena-1; i>=0; i--)
  423. Ddivide(carry, bignum_digits(a)[i], carry, bignum_digits(a)[i], n);
  424. if ((sign & 2) != 0) carry = -(int32)carry;
  425. nil = C_nil;
  426. nwork = carry; /* leave remainder available to caller */
  427. lena--;
  428. while (bignum_digits(a)[lena] == 0) lena--;
  429. if ((bignum_digits(a)[lena] & 0x40000000) != 0) lena++;
  430. if ((sign & 1) != 0) /* quotient will be negative */
  431. { carry = 0xffffffffU;
  432. for (i=0; i<lena; i++)
  433. { carry = clear_top_bit(~bignum_digits(a)[i]) + top_bit(carry);
  434. bignum_digits(a)[i] = clear_top_bit(carry);
  435. }
  436. carry = ~bignum_digits(a)[i] + top_bit(carry);
  437. if (carry == -1 && (bignum_digits(a)[i-1] & 0x40000000) != 0)
  438. { bignum_digits(a)[lena] = 0;
  439. lena--;
  440. bignum_digits(a)[i-1] |= ~0x7fffffff;
  441. }
  442. else bignum_digits(a)[i] = carry;
  443. }
  444. /* I need to free up some space here, I guess */
  445. lenx = (bignum_length(a)>>2)-1;
  446. lenx |= 1; /* round up to allow for allocation in doublewords */
  447. lenc = (lena+1) | 1;
  448. if (lenc != lenx) /* space to discard? */
  449. bignum_digits(a)[lenc] = make_bighdr(lenx-lenc);
  450. numhdr(a) = make_bighdr(lena+2);
  451. return a;
  452. }
  453. Lisp_Object quotbn1(Lisp_Object a, int32 n)
  454. /*
  455. * Divide a bignum by an integer, where the integer is (by now)
  456. * a natural C int32 but limited to 31 not 32 bits active. I.e.
  457. * this can be used when dividing by a fixnum or by dividing by
  458. * a one-word bignum. I will not call this with n=-1, which would
  459. * be the one case where it could cause the quotient to be bigger
  460. * than the dividend. Leaves nwork set to the remainder.
  461. * This version is JUST the same as quotbn() except that it does not
  462. * hand back a useful quotient - it is intended for use when only
  463. * the remainder is wanted. For consistency I leave that where quotbn() did.
  464. */
  465. {
  466. Lisp_Object nil;
  467. int sign;
  468. int32 lena = (bignum_length(a)>>2)-2, i;
  469. unsigned32 carry;
  470. if (lena == 0) /* one-word bignum as numerator */
  471. { int32 p = (int32)bignum_digits(a)[0];
  472. nil = C_nil; /* So I can access nwork */
  473. nwork = p % n;
  474. /*
  475. * C does not define what happens on non-exact division involving
  476. * negative quantities, so I adjust things here so that the remainder
  477. * has the sign that I want and the division I do is exact.
  478. */
  479. if (p < 0)
  480. { if (nwork > 0) nwork -= n;
  481. }
  482. else if (nwork < 0) nwork += n;
  483. /* Remainder should be correct (with correct sign) by now, regardless */
  484. else return nil;
  485. }
  486. else if (lena == 1)
  487. /*
  488. * I treat division of a 2-word bignum by a fixnum or 1-word bignum as
  489. * a special case since it can lead to a fixnum result - if the divisor is
  490. * just one word long and the dividend is 3 or more words I would
  491. * certainly have a bignum result. Thus by separating off the code here
  492. * I (a) remove the need for test relating to big- to fixnum conversion
  493. * later on and (b) avoid allocating heap in this tolerably common case
  494. * when sometimes I will not need to.
  495. */
  496. { unsigned32 a0 = bignum_digits(a)[0];
  497. int32 a1 = (int32)bignum_digits(a)[1];
  498. if (a1 < 0)
  499. { sign = 3;
  500. if (a0 == 0) a1 = -a1;
  501. else
  502. { a0 = clear_top_bit(-(int32)a0);
  503. a1 = ~a1;
  504. }
  505. }
  506. else sign = 0;
  507. if (n < 0) sign ^= 1, n = -n;
  508. if (a1 < n) /* result can be calculated in one word - good! */
  509. { unsigned32 q, r;
  510. Ddivide(r, q, (unsigned32)a1, a0, n);
  511. if ((sign & 2) != 0) r = -(int32)r;
  512. nil = C_nil;
  513. nwork = r;
  514. return nil;
  515. }
  516. /*
  517. * here the quotient will involve two digits.
  518. */
  519. else
  520. { unsigned32 q1, q0, r;
  521. Ddivide(r, q1, 0, (unsigned32)a1, n);
  522. Ddivide(r, q0, r, a0, n);
  523. if ((sign & 2) != 0) r = -(int32)r;
  524. nil = C_nil;
  525. nwork = r;
  526. return nil;
  527. }
  528. }
  529. /*
  530. * That has dealt with the special cases - now the input is a bignum
  531. * with at least 3 digits and the quotient will certainly be a bignum.
  532. * Start by allocating a workspace copy of the dividend. negateb will
  533. * leave a a bignum, although it may change its length.
  534. */
  535. /*
  536. * Also note that I could fold the negateb() in with the division, and
  537. * thereby save allocating a big hunk of extra memory.
  538. */
  539. if ((int32)bignum_digits(a)[lena] < 0) a = negateb(a), sign = 3;
  540. else a = copyb(a), sign = 0;
  541. errexit();
  542. if (n < 0)
  543. { sign ^= 1;
  544. n = -n;
  545. }
  546. /*
  547. * Now both a and n have been made positive, and their original signs
  548. * have been recorded so that I can tidy up at the end. The 1 bit in sign
  549. * tells me what sign the quotient will have, the 2 bit gives the sign
  550. * for the remainder.
  551. */
  552. lena = (bignum_length(a)>>2)-1;
  553. carry = 0;
  554. for (i=lena-1; i>=0; i--)
  555. Ddivide(carry, bignum_digits(a)[i], carry, bignum_digits(a)[i], n);
  556. if ((sign & 2) != 0) carry = -(int32)carry;
  557. nil = C_nil;
  558. nwork = carry; /* leave remainder available to caller */
  559. return nil;
  560. }
  561. static Lisp_Object quotbi(Lisp_Object a, Lisp_Object b)
  562. {
  563. /*
  564. * dividing by 1 or -1 seems worth optimising.
  565. */
  566. if (b == fixnum_of_int(1)) return a;
  567. else if (b == fixnum_of_int(-1)) return negateb(a);
  568. else if (b == fixnum_of_int(0))
  569. return aerror2("bad arg for quotient", a, b);
  570. else return quotbn(a, int_of_fixnum(b));
  571. }
  572. #ifdef COMMON
  573. #define quotbs(a, b) quotsb(a, b)
  574. #endif
  575. #ifdef __powerc
  576. /* If you have trouble compiling this just comment it out, please */
  577. #pragma options(!global_optimizer)
  578. #endif
  579. Lisp_Object quotbb(Lisp_Object a, Lisp_Object b)
  580. /*
  581. * Divide one bignum by another. This can compute the
  582. * remainder at the same time as the quotient, and leaves same around
  583. * in mv_2. If it is not needed then setting up the remainder is
  584. * a cost - but usually a modest one (in context), so I think that
  585. * the simplification is worth-while.
  586. */
  587. {
  588. Lisp_Object nil, olda;
  589. int sign;
  590. int32 lena, lenb, lenc, lenx, i;
  591. unsigned32 carry, scale, v1;
  592. Lisp_Object q, r;
  593. /*
  594. * If I am dividing by a one-word bignum I can go a bit faster...
  595. */
  596. lenb = (bignum_length(b)>>2)-2;
  597. if (lenb == 0)
  598. { q = quotbn(a, bignum_digits(b)[0]);
  599. errexit();
  600. /*
  601. * Now lots of frivolity packing up the remainder...
  602. */
  603. r = (Lisp_Object)(nwork & fix_mask);
  604. if ((int32)r == 0 || (int32)r == fix_mask)
  605. mv_2 = fixnum_of_int(nwork);
  606. else
  607. { push(q);
  608. a = make_one_word_bignum(nwork);
  609. pop(q);
  610. errexit();
  611. mv_2 = a;
  612. }
  613. return q;
  614. }
  615. /*
  616. * Convert to sign and magnitude representation
  617. */
  618. push2(a, b);
  619. lena = (bignum_length(a)>>2)-2;
  620. if ((int32)bignum_digits(a)[lena] < 0)
  621. { a = negateb(a);
  622. sign = 3;
  623. }
  624. else
  625. { a = copyb(a);
  626. sign = 0;
  627. }
  628. pop(b);
  629. errexit();
  630. lena = (bignum_length(a)>>2)-2;
  631. push(a);
  632. if ((int32)bignum_digits(b)[lenb] < 0)
  633. { b = negateb(b);
  634. /*
  635. * I MUST NOT do lenb = (bignum_length(b)>>2)-2; here because the
  636. * negateb may have failed and therefore handed back junk rather than
  637. * a bignum. E.g. an interrupt provoked by the user or a store-jam might
  638. * lead to a failure report here.
  639. */
  640. sign ^= 1;
  641. }
  642. else b = copyb(b);
  643. pop2(a, olda); /* original a, with original sign */
  644. errexit();
  645. lenb = (bignum_length(b)>>2)-2;
  646. /*
  647. * Now that the numbers are unsigned it could be that I can drop
  648. * a leading zero that had previously been necessary. That could reveal
  649. * that I have a one-word divisor after all...
  650. */
  651. if (bignum_digits(b)[lenb] == 0) lenb--;
  652. if (lenb == 0)
  653. { a = quotbn(a, bignum_digits(b)[0]);
  654. errexit();
  655. if ((sign & 2) != 0) nwork = -nwork;
  656. if ((sign & 1) != 0) a = negate(a);
  657. errexit();
  658. r = (Lisp_Object)(nwork & fix_mask);
  659. if (r == 0 || (int32)r == fix_mask) mv_2 = fixnum_of_int(nwork);
  660. else
  661. { push(a);
  662. if (signed_overflow(nwork))
  663. { if ((int32)nwork < 0)
  664. b = make_two_word_bignum(-1, clear_top_bit(nwork));
  665. else b = make_two_word_bignum(0, nwork);
  666. }
  667. else b = make_one_word_bignum(nwork);
  668. pop(a);
  669. errexit();
  670. mv_2 = b;
  671. }
  672. return a;
  673. }
  674. /*
  675. * Now the divisor is at least two digits long.
  676. */
  677. if (bignum_digits(a)[lena] == 0) lena--;
  678. /*
  679. * I will take a special case when the lengths of the two numbers
  680. * match since in that case the quotient will be only one digit long.
  681. * Also after I have filtered the lena<=lenb cases out, and have treated
  682. * one-word cases of b specially I will know that the numerator a has
  683. * at least two digits.
  684. */
  685. if (lena < lenb)
  686. { mv_2 = olda;
  687. return fixnum_of_int(0);
  688. }
  689. else if (lenb == lena)
  690. { unsigned32 p0 = bignum_digits(a)[lena-1], p1 = bignum_digits(a)[lena],
  691. q0 = bignum_digits(b)[lenb-1], q1 = bignum_digits(b)[lena],
  692. r0, r1, q, w, carry1;
  693. /*
  694. * If the dividend is smaller than the divisor I can return zero right now.
  695. */
  696. if (p1 < q1 || (p1 == q1 && p0 < q0))
  697. { mv_2 = olda;
  698. return fixnum_of_int(0);
  699. }
  700. /*
  701. * I scale the divisor so that it will have its top bit set (top wrt the
  702. * 31 bit field that is in use, that is), and scale the top two digits
  703. * of the dividend to match. The resulting values can then be divided to get
  704. * a pretty good estimate for the true quotient. Note that the division on
  705. * the next line is UNSIGNED arithmetic.
  706. */
  707. scale = 0x80000000U / ((unsigned32)1 + q1);
  708. Dmultiply(carry, w, q0, scale, 0);
  709. Dmultiply(q, w, q1, scale, carry);
  710. q = w;
  711. Dmultiply(carry, w, p0, scale, 0);
  712. Dmultiply(carry, w, p1, scale, carry);
  713. if (carry == q) q = 0x7fffffff;
  714. else
  715. { Ddivide(q, w, carry, w, q);
  716. q = w;
  717. }
  718. /*
  719. * q is now my estimated quotient, based on a quick look at high digits.
  720. */
  721. Dmultiply(carry, w, q0, q, 0);
  722. r0 = w;
  723. Dmultiply(carry, w, q1, q, carry);
  724. r1 = w;
  725. r0 = r0 - p0;
  726. if ((int32)r0 < 0)
  727. { r0 = clear_top_bit(r0);
  728. r1 = r1 - p1 - 1;
  729. }
  730. else r1 = r1 - p1;
  731. /* /*
  732. * the next line is a cop-out for now - if my estimated quotient
  733. * was close enough to the true value than the residual I get here
  734. * ought to be fairly small - if it is not I have bungled. Over a year
  735. * or so of tests I have not seen the disaster message triggered, but the
  736. * code should stay in until I can write the paragraph of comment that
  737. * should go here explaing why all is well...
  738. */
  739. if (carry != 0 && (int32)r1 < 0 &&
  740. carry != 1 && r1 != ~0x7fffffff)
  741. {
  742. err_printf("\n+++!!!+++ carry needed (line %d of \"%s\")\n",
  743. __LINE__, __FILE__);
  744. my_exit(EXIT_FAILURE);
  745. }
  746. /*
  747. * That obtained the remainder from (p1,p0)/(q1,q0) - now adjust q until
  748. * the remainder has the sign I want it to have.
  749. */
  750. /* /* I do not look at carry here - I have a nasty suspicion that I should.. */
  751. while ((int32)r1 > 0 ||
  752. (r1 == 0 && r0 != 0))
  753. { q--;
  754. r0 = r0 - q0;
  755. if ((int32)r0 < 0)
  756. { r0 = clear_top_bit(r0);
  757. r1 = r1 - q1 - 1;
  758. }
  759. else r1 = r1 - q1;
  760. }
  761. /*
  762. * Now q is (p1,p0)/(q1,q0), which is an over-estimate for the true
  763. * quotient, but by at worst 1. Compute a proper full-length remainder
  764. * now, using the original unscaled input numbers.
  765. */
  766. carry = 0;
  767. carry1 = 0xffffffffU;
  768. for (i=0; i<=lena; i++)
  769. { Dmultiply(carry, w, bignum_digits(b)[i], q, carry);
  770. carry1 = bignum_digits(a)[i] + clear_top_bit(~w) + top_bit(carry1);
  771. bignum_digits(a)[i] = clear_top_bit(carry1);
  772. }
  773. if ((int32)carry1 >= 0) /* need to correct it! */
  774. { q--;
  775. for (i=0; i<=lena; i++)
  776. { carry1 = bignum_digits(a)[i] + bignum_digits(b)[i] +
  777. top_bit(carry1);
  778. bignum_digits(a)[i] = clear_top_bit(carry1);
  779. }
  780. }
  781. /*
  782. * Now the quotient is correct - but since I want to hand back a neat
  783. * remainder I have more to do - I must convert the value stored in a
  784. * into a legitimate result, and store it in mv_2.
  785. */
  786. while (lena > 0 && bignum_digits(a)[lena] == 0) lena--;
  787. if ((bignum_digits(a)[lena] & 0x40000000) != 0) lena++;
  788. if (lena == 0) /* Maybe fixnum remainder? */
  789. { int32 rr = bignum_digits(a)[0];
  790. if (rr != 0 && (sign & 2) != 0)
  791. { rr = -rr;
  792. if ((rr & fix_mask) == fix_mask)
  793. { mv_2 = fixnum_of_int(rr);
  794. goto return_q;
  795. }
  796. }
  797. else if ((rr & fix_mask) == 0)
  798. { mv_2 = fixnum_of_int(rr);
  799. goto return_q;
  800. }
  801. }
  802. if ((sign & 2) != 0)
  803. { carry = 0xffffffffU;
  804. for (i=0; i<lena; i++)
  805. { carry = clear_top_bit(~bignum_digits(a)[i]) + top_bit(carry);
  806. bignum_digits(a)[i] = clear_top_bit(carry);
  807. }
  808. carry = ~bignum_digits(a)[i] + top_bit(carry);
  809. bignum_digits(a)[i] = carry;
  810. /*
  811. * if my remainder is -2^(31*n-1) then I can need to shrink lena here, which
  812. * seems like a messy case to have to consider.
  813. */
  814. if (carry == -1 && (bignum_digits(a)[i-1] & 0x40000000) != 0)
  815. { bignum_digits(a)[lena] = 0;
  816. lena--;
  817. bignum_digits(a)[i-1] |= ~0x7fffffff;
  818. }
  819. }
  820. /*
  821. * Now tidy up the space that I am about to discard, so that it will not
  822. * upset the garbage collector.
  823. */
  824. lenx = (bignum_length(a)>>2)-1;
  825. lenx |= 1; /* round up to allow for allocation in doublewords */
  826. lenc = (lena+1) | 1;
  827. if (lenc != lenx) /* space to discard? */
  828. bignum_digits(a)[lenc] = make_bighdr(lenx-lenc);
  829. numhdr(a) = make_bighdr(lena+2);
  830. mv_2 = a;
  831. return_q:
  832. if (q != 0 && (sign & 1) != 0)
  833. { q = -(int32)q;
  834. if ((q & fix_mask) == fix_mask) return fixnum_of_int(q);
  835. if ((q & 0x40000000) == 0)
  836. return make_two_word_bignum(-1, clear_top_bit(q));
  837. }
  838. else if ((q & fix_mask) == 0) return fixnum_of_int(q);
  839. else if ((q & 0x40000000) != 0) return make_two_word_bignum(0, q);
  840. return make_one_word_bignum(q);
  841. }
  842. /*
  843. * Now both a and b have been made positive, and their original signs
  844. * have been recorded so that I can tidy up at the end. The 1 bit in sign
  845. * tells me what sign the quotient will have, the 2 bit gives the sign
  846. * for the remainder. Both a and b have been copied so it will be OK to
  847. * overwrite them in the process of doing the division. The length of
  848. * a is strictly greater than that of b, which itself is at least 2 digits,
  849. * so this is a real case of long division, although the quotient may still
  850. * end up small (but non-zero).
  851. */
  852. /*
  853. * I find a scale factor to multiply a and b by so that the leading digit of
  854. * the divisor is fairly large. Again note that the arithmetic is UNSIGNED.
  855. */
  856. scale = 0x80000000U / ((unsigned32)1 + (unsigned32)bignum_digits(b)[lenb]);
  857. carry = 0;
  858. for (i=0; i<=lenb; i++)
  859. Dmultiply(carry, bignum_digits(b)[i],
  860. bignum_digits(b)[i], scale, carry);
  861. v1 = bignum_digits(b)[lenb];
  862. carry = 0;
  863. for (i=0; i<=lena; i++)
  864. Dmultiply(carry, bignum_digits(a)[i],
  865. bignum_digits(a)[i], scale, carry);
  866. /*
  867. * The scaled value of the numerator a may involve an extra digit
  868. * beyond those that I have space for in the vector. Hold this digit
  869. * in the variable 'carry'. Note that scaling the divisor did NOT cause
  870. * it to change length.
  871. */
  872. while (lena >= lenb)
  873. { unsigned32 w, h0, h1, v2, u2, carry1, carry2, qq, rr;
  874. if (carry == v1)
  875. { qq = 0x7fffffff;
  876. rr = carry + bignum_digits(a)[lena];
  877. }
  878. else
  879. { Ddivide(rr, w, carry, bignum_digits(a)[lena], v1);
  880. qq = w;
  881. }
  882. v2 = bignum_digits(b)[lenb-1];
  883. Dmultiply(h1, h0, v2, qq, 0);
  884. u2 = bignum_digits(a)[lena-1];
  885. if (h1 > rr ||
  886. (h1 == rr && h0 > u2))
  887. { qq--;
  888. h0 -= v2;
  889. if ((int32)h0 < 0)
  890. { h0 = clear_top_bit(h0);
  891. h1--;
  892. }
  893. rr += v1;
  894. if (h1 > rr ||
  895. (h1 == rr && h0 > u2)) qq--;
  896. }
  897. /*
  898. * Now subtract off q*b from a, up as far as lena.
  899. */
  900. carry1 = 0;
  901. carry2 = 0xffffffffU;
  902. for (i=0; i<=lenb; i++)
  903. { Dmultiply(carry1, w, bignum_digits(b)[i], qq, carry1);
  904. carry2 = bignum_digits(a)[lena-lenb+i] +
  905. clear_top_bit(~w) + top_bit(carry2);
  906. bignum_digits(a)[lena-lenb+i] = clear_top_bit(carry2);
  907. }
  908. carry2 = carry + clear_top_bit(~carry1) + top_bit(carry2);
  909. if ((int32)carry2 >= 0) /* overshot.. */
  910. { qq--;
  911. carry2 = 0;
  912. for (i=0; i<=lenb; i++)
  913. { carry2 = bignum_digits(a)[lena-lenb+i] +
  914. bignum_digits(b)[i] + top_bit(carry2);
  915. bignum_digits(a)[lena-lenb+i] = clear_top_bit(carry2);
  916. }
  917. }
  918. carry = bignum_digits(a)[lena];
  919. bignum_digits(a)[lena] = qq; /* good place to store the quotient */
  920. lena--;
  921. }
  922. lena = (bignum_length(a)>>2)-2;
  923. /*
  924. * Now the quotient is in the top digits of a, and the remainder is left
  925. * (scaled up) in the low bit (plus carry)
  926. */
  927. bignum_digits(b)[lenb] = carry/scale;
  928. carry = carry%scale;
  929. for (i=lenb-1; i>=0; i--)
  930. Ddivide(carry, bignum_digits(b)[i],
  931. carry, bignum_digits(a)[i], scale);
  932. if (carry != 0)
  933. { err_printf("+++!!!+++ Incorrect unscaling line %d file \"%s\" (%ld)\n",
  934. __LINE__, __FILE__, (long)carry);
  935. my_exit(EXIT_FAILURE);
  936. }
  937. for (i=0; i<=lena-lenb; i++)
  938. bignum_digits(a)[i] = bignum_digits(a)[i+lenb];
  939. for (;i < lena; i++) bignum_digits(a)[i] = 0;
  940. lena = lena - lenb;
  941. while (bignum_digits(a)[lena] == 0) lena--;
  942. /*
  943. * the above loop terminates properly since the quotient will be non-zero,
  944. * which in turn is because I was dividing p by q where p had strictly
  945. * more digits than q, and was hence strictly greater than q.
  946. */
  947. if ((bignum_digits(a)[lena] & 0x40000000) != 0) lena++;
  948. if ((sign & 1) != 0) /* quotient will be negative */
  949. { carry = 0xffffffffU;
  950. for (i=0; i<lena; i++)
  951. { carry = clear_top_bit(~bignum_digits(a)[i]) + top_bit(carry);
  952. bignum_digits(a)[i] = clear_top_bit(carry);
  953. }
  954. carry = ~bignum_digits(a)[i] + top_bit(carry);
  955. if (carry == -1 && (bignum_digits(a)[i-1] & 0x40000000) != 0)
  956. { bignum_digits(a)[lena] = 0;
  957. lena--;
  958. bignum_digits(a)[i-1] |= ~0x7fffffff;
  959. }
  960. else bignum_digits(a)[i] = carry;
  961. }
  962. lenx = (bignum_length(a)>>2)-1;
  963. lenx |= 1; /* round up to allow for allocation in doublewords */
  964. lenc = (lena+1) | 1;
  965. if (lenc != lenx) /* space to discard? */
  966. bignum_digits(a)[lenc] = make_bighdr(lenx-lenc);
  967. numhdr(a) = make_bighdr(lena+2);
  968. if (lena == 0)
  969. { int32 d = bignum_digits(a)[0];
  970. int32 x = d & fix_mask;
  971. if (x == 0 || x == fix_mask) a = fixnum_of_int(d);
  972. }
  973. /*
  974. * Here I need to tidy up the discarded space that could otherwise
  975. * kill the poor old garbage collector..
  976. */
  977. while (lenb >= 0 && bignum_digits(b)[lenb] == 0) lenb--;
  978. if (lenb < 0) mv_2 = fixnum_of_int(0); /* hah - zero remainder! */
  979. else
  980. {
  981. if ((bignum_digits(b)[lenb] & 0x40000000) != 0) lenb++;
  982. if ((sign & 2) != 0) /* need to negate the remainder */
  983. { carry = 0xffffffffU;
  984. for (i=0; i<lenb; i++)
  985. { carry = clear_top_bit(~bignum_digits(b)[i]) + top_bit(carry);
  986. bignum_digits(b)[i] = clear_top_bit(carry);
  987. }
  988. carry = ~bignum_digits(b)[i] + top_bit(carry);
  989. if (carry == -1 && (bignum_digits(b)[i-1] & 0x40000000) != 0)
  990. { bignum_digits(b)[lenb] = 0;
  991. lenb--;
  992. bignum_digits(b)[i-1] |= ~0x7fffffff;
  993. }
  994. else bignum_digits(b)[i] = carry;
  995. }
  996. lenx = (bignum_length(b)>>2)-1;
  997. lenx |= 1; /* round up to allow for allocation in doublewords */
  998. lenc = (lenb+1) | 1;
  999. if (lenc != lenx) /* space to discard? */
  1000. bignum_digits(b)[lenc] = make_bighdr(lenx-lenc);
  1001. numhdr(b) = make_bighdr(lenb+2);
  1002. mv_2 = b;
  1003. if (lenb == 0)
  1004. { int32 r1, rr;
  1005. rr = bignum_digits(b)[0];
  1006. r1 = rr & fix_mask;
  1007. if (r1 == 0 || r1 == fix_mask) mv_2 = fixnum_of_int(rr);
  1008. }
  1009. }
  1010. return a;
  1011. }
  1012. #ifdef __powerc
  1013. /* If you have trouble compiling this just comment it out, please */
  1014. #pragma options(global_optimizer)
  1015. #endif
  1016. #ifdef COMMON
  1017. #define quotbr(a, b) quotir(a, b)
  1018. #define quotbc(a, b) quotic(a, b)
  1019. #endif
  1020. #define quotbf(a, b) quotsf(a, b)
  1021. #ifdef COMMON
  1022. static Lisp_Object quotri(Lisp_Object a, Lisp_Object b)
  1023. {
  1024. Lisp_Object w, nil;
  1025. if (b == fixnum_of_int(1)) return a;
  1026. else if (b == fixnum_of_int(0))
  1027. return aerror2("bad arg for quotient", a, b);
  1028. push3(a, b, C_nil);
  1029. #define g stack[0]
  1030. #define b stack[-1]
  1031. #define a stack[-2]
  1032. g = gcd(b, numerator(a));
  1033. nil = C_nil;
  1034. if (exception_pending()) goto fail;
  1035. w = minusp(b);
  1036. nil = C_nil;
  1037. if (exception_pending()) goto fail;
  1038. if (w)
  1039. { g = negate(g); /* ensure denominator is +ve */
  1040. nil = C_nil;
  1041. if (exception_pending()) goto fail;
  1042. }
  1043. b = quot2(b, g);
  1044. nil = C_nil;
  1045. if (exception_pending()) goto fail;
  1046. g = quot2(numerator(a), g);
  1047. nil = C_nil;
  1048. if (exception_pending()) goto fail;
  1049. a = times2(b, denominator(a));
  1050. nil = C_nil;
  1051. if (exception_pending()) goto fail;
  1052. w = make_ratio(g, a);
  1053. popv(3);
  1054. return w;
  1055. fail:
  1056. popv(3);
  1057. return nil;
  1058. #undef a
  1059. #undef b
  1060. #undef g
  1061. }
  1062. #define quotrs(a, b) quotsb(a, b)
  1063. #define quotrb(a, b) quotib(a, b)
  1064. static Lisp_Object quotrr(Lisp_Object a, Lisp_Object b)
  1065. {
  1066. Lisp_Object w, nil;
  1067. push5(numerator(a), denominator(a),
  1068. denominator(b), numerator(b), /* NB switched order */
  1069. C_nil);
  1070. #define g stack[0]
  1071. #define db stack[-1]
  1072. #define nb stack[-2]
  1073. #define da stack[-3]
  1074. #define na stack[-4]
  1075. g = gcd(na, db);
  1076. nil = C_nil;
  1077. if (exception_pending()) goto fail;
  1078. w = minusp(db);
  1079. nil = C_nil;
  1080. if (exception_pending()) goto fail;
  1081. if (w)
  1082. { g = negate(g);
  1083. nil = C_nil;
  1084. if (exception_pending()) goto fail;
  1085. }
  1086. na = quot2(na, g);
  1087. nil = C_nil;
  1088. if (exception_pending()) goto fail;
  1089. db = quot2(db, g);
  1090. nil = C_nil;
  1091. if (exception_pending()) goto fail;
  1092. g = gcd(nb, da);
  1093. if (exception_pending()) goto fail;
  1094. nil = C_nil;
  1095. nb = quot2(nb, g);
  1096. if (exception_pending()) goto fail;
  1097. da = quot2(da, g);
  1098. nil = C_nil;
  1099. if (exception_pending()) goto fail;
  1100. na = times2(na, nb);
  1101. nil = C_nil;
  1102. if (exception_pending()) goto fail;
  1103. da = times2(da, db);
  1104. nil = C_nil;
  1105. if (exception_pending()) goto fail;
  1106. w = make_ratio(na, da);
  1107. popv(5);
  1108. return w;
  1109. fail:
  1110. popv(5);
  1111. return nil;
  1112. #undef g
  1113. #undef db
  1114. #undef nb
  1115. #undef da
  1116. #undef na
  1117. }
  1118. #define quotrc(a, b) quotic(a, b)
  1119. #define quotrf(a, b) quotsf(a, b)
  1120. static Lisp_Object quotci(Lisp_Object a, Lisp_Object b)
  1121. {
  1122. Lisp_Object r = real_part(a), i = imag_part(a), nil;
  1123. if (b == fixnum_of_int(0)) return aerror2("bad arg for quotient", a, b);
  1124. push2(b, r);
  1125. i = quot2(i, b);
  1126. pop2(r, b);
  1127. errexit();
  1128. push(i);
  1129. r = quot2(r, b);
  1130. pop(i);
  1131. errexit();
  1132. return make_complex(r, i);
  1133. }
  1134. #define quotcs(a, b) quotci(a, b)
  1135. #define quotcb(a, b) quotci(a, b)
  1136. #define quotcr(a, b) quotci(a, b)
  1137. #define quotcc(a, b) quotic(a, b)
  1138. #define quotcf(a, b) quotci(a, b)
  1139. #endif
  1140. static Lisp_Object quotfi(Lisp_Object a, Lisp_Object b)
  1141. {
  1142. double d;
  1143. if (b == fixnum_of_int(0)) return aerror2("bad arg for quotient", a, b);
  1144. d = float_of_number(a) / (double)int_of_fixnum(b);
  1145. return make_boxfloat(d, type_of_header(flthdr(a)));
  1146. }
  1147. static Lisp_Object quotfs(Lisp_Object a, Lisp_Object b)
  1148. {
  1149. double d = float_of_number(b);
  1150. if (d == 0.0) return aerror2("bad arg for quotient", a, b);
  1151. d = float_of_number(a) / d;
  1152. return make_boxfloat(d, type_of_header(flthdr(a)));
  1153. }
  1154. #define quotfb(a, b) quotfs(a, b)
  1155. #ifdef COMMON
  1156. #define quotfr(a, b) quotfs(a, b)
  1157. #define quotfc(a, b) quotic(a, b)
  1158. #endif
  1159. static Lisp_Object quotff(Lisp_Object a, Lisp_Object b)
  1160. {
  1161. #ifdef COMMON
  1162. int32 ha = type_of_header(flthdr(a)), hb = type_of_header(flthdr(b));
  1163. #endif
  1164. double d = float_of_number(b);
  1165. if (d == 0.0) return aerror2("bad arg for quotient", a, b);
  1166. d = float_of_number(a) / d;
  1167. #ifdef COMMON
  1168. if (ha == TYPE_LONG_FLOAT || hb == TYPE_LONG_FLOAT)
  1169. ha = TYPE_LONG_FLOAT;
  1170. else if (ha == TYPE_DOUBLE_FLOAT || hb == TYPE_DOUBLE_FLOAT)
  1171. ha = TYPE_DOUBLE_FLOAT;
  1172. else ha = TYPE_SINGLE_FLOAT;
  1173. return make_boxfloat(d, ha);
  1174. #else
  1175. return make_boxfloat(d, TYPE_DOUBLE_FLOAT);
  1176. #endif
  1177. }
  1178. Lisp_Object quot2(Lisp_Object a, Lisp_Object b)
  1179. {
  1180. switch ((int)a & TAG_BITS)
  1181. {
  1182. case TAG_FIXNUM:
  1183. switch ((int)b & TAG_BITS)
  1184. {
  1185. case TAG_FIXNUM:
  1186. /*
  1187. * This is where fixnum / fixnum arithmetic happens - the case I most want to
  1188. * make efficient.
  1189. */
  1190. if (b == fixnum_of_int(0))
  1191. return aerror2("bad arg for quotient", a, b);
  1192. else
  1193. { int32 r, aa, bb;
  1194. aa = int_of_fixnum(a);
  1195. bb = int_of_fixnum(b);
  1196. /* calculate remainder and force its sign to be correct */
  1197. r = aa % bb;
  1198. if (aa < 0)
  1199. { if (r > 0) r -= bb;
  1200. }
  1201. else if (r < 0) r += bb;
  1202. /* then the division can be an exact one, as here */
  1203. r = (aa - r)/bb;
  1204. /*
  1205. * the only case where dividing one fixnum by another can lead to
  1206. * a bignum result is (-0x08000000/(-1)) which JUST overflows.
  1207. */
  1208. if (r != 0x08000000) return fixnum_of_int(r);
  1209. else return make_one_word_bignum(r);
  1210. }
  1211. #ifdef COMMON
  1212. case TAG_SFLOAT:
  1213. return quotis(a, b);
  1214. #endif
  1215. case TAG_NUMBERS:
  1216. { int32 hb = type_of_header(numhdr(b));
  1217. switch (hb)
  1218. {
  1219. case TYPE_BIGNUM:
  1220. return quotib(a, b);
  1221. #ifdef COMMON
  1222. case TYPE_RATNUM:
  1223. return quotir(a, b);
  1224. case TYPE_COMPLEX_NUM:
  1225. return quotic(a, b);
  1226. #endif
  1227. default:
  1228. return aerror1("bad arg for quotient", b);
  1229. }
  1230. }
  1231. case TAG_BOXFLOAT:
  1232. return quotif(a, b);
  1233. default:
  1234. return aerror1("bad arg for quotient", b);
  1235. }
  1236. #ifdef COMMON
  1237. case TAG_SFLOAT:
  1238. switch ((int)b & TAG_BITS)
  1239. {
  1240. case TAG_FIXNUM:
  1241. return quotsi(a, b);
  1242. case TAG_SFLOAT:
  1243. { Float_union aa, bb;
  1244. aa.i = a - TAG_SFLOAT;
  1245. bb.i = b - TAG_SFLOAT;
  1246. aa.f = (float) (aa.f / bb.f);
  1247. return (aa.i & ~(int32)0xf) + TAG_SFLOAT;
  1248. }
  1249. case TAG_NUMBERS:
  1250. { int32 hb = type_of_header(numhdr(b));
  1251. switch (hb)
  1252. {
  1253. case TYPE_BIGNUM:
  1254. return quotsb(a, b);
  1255. case TYPE_RATNUM:
  1256. return quotsr(a, b);
  1257. case TYPE_COMPLEX_NUM:
  1258. return quotsc(a, b);
  1259. default:
  1260. return aerror1("bad arg for quotient", b);
  1261. }
  1262. }
  1263. case TAG_BOXFLOAT:
  1264. return quotsf(a, b);
  1265. default:
  1266. return aerror1("bad arg for quotient", b);
  1267. }
  1268. #endif
  1269. case TAG_NUMBERS:
  1270. { int32 ha = type_of_header(numhdr(a));
  1271. switch (ha)
  1272. {
  1273. case TYPE_BIGNUM:
  1274. switch ((int)b & TAG_BITS)
  1275. {
  1276. case TAG_FIXNUM:
  1277. return quotbi(a, b);
  1278. #ifdef COMMON
  1279. case TAG_SFLOAT:
  1280. return quotbs(a, b);
  1281. #endif
  1282. case TAG_NUMBERS:
  1283. { int32 hb = type_of_header(numhdr(b));
  1284. switch (hb)
  1285. {
  1286. case TYPE_BIGNUM:
  1287. return quotbb(a, b);
  1288. #ifdef COMMON
  1289. case TYPE_RATNUM:
  1290. return quotbr(a, b);
  1291. case TYPE_COMPLEX_NUM:
  1292. return quotbc(a, b);
  1293. #endif
  1294. default:
  1295. return aerror1("bad arg for quotient", b);
  1296. }
  1297. }
  1298. case TAG_BOXFLOAT:
  1299. return quotbf(a, b);
  1300. default:
  1301. return aerror1("bad arg for quotient", b);
  1302. }
  1303. #ifdef COMMON
  1304. case TYPE_RATNUM:
  1305. switch ((int)b & TAG_BITS)
  1306. {
  1307. case TAG_FIXNUM:
  1308. return quotri(a, b);
  1309. case TAG_SFLOAT:
  1310. return quotrs(a, b);
  1311. case TAG_NUMBERS:
  1312. { int32 hb = type_of_header(numhdr(b));
  1313. switch (hb)
  1314. {
  1315. case TYPE_BIGNUM:
  1316. return quotrb(a, b);
  1317. case TYPE_RATNUM:
  1318. return quotrr(a, b);
  1319. case TYPE_COMPLEX_NUM:
  1320. return quotrc(a, b);
  1321. default:
  1322. return aerror1("bad arg for quotient", b);
  1323. }
  1324. }
  1325. case TAG_BOXFLOAT:
  1326. return quotrf(a, b);
  1327. default:
  1328. return aerror1("bad arg for quotient", b);
  1329. }
  1330. case TYPE_COMPLEX_NUM:
  1331. switch ((int)b & TAG_BITS)
  1332. {
  1333. case TAG_FIXNUM:
  1334. return quotci(a, b);
  1335. case TAG_SFLOAT:
  1336. return quotcs(a, b);
  1337. case TAG_NUMBERS:
  1338. { int32 hb = type_of_header(numhdr(b));
  1339. switch (hb)
  1340. {
  1341. case TYPE_BIGNUM:
  1342. return quotcb(a, b);
  1343. case TYPE_RATNUM:
  1344. return quotcr(a, b);
  1345. case TYPE_COMPLEX_NUM:
  1346. return quotcc(a, b);
  1347. default:
  1348. return aerror1("bad arg for quotient", b);
  1349. }
  1350. }
  1351. case TAG_BOXFLOAT:
  1352. return quotcf(a, b);
  1353. default:
  1354. return aerror1("bad arg for quotient", b);
  1355. }
  1356. #endif
  1357. default: return aerror1("bad arg for quotient", a);
  1358. }
  1359. }
  1360. case TAG_BOXFLOAT:
  1361. switch ((int)b & TAG_BITS)
  1362. {
  1363. case TAG_FIXNUM:
  1364. return quotfi(a, b);
  1365. #ifdef COMMON
  1366. case TAG_SFLOAT:
  1367. return quotfs(a, b);
  1368. #endif
  1369. case TAG_NUMBERS:
  1370. { int32 hb = type_of_header(numhdr(b));
  1371. switch (hb)
  1372. {
  1373. case TYPE_BIGNUM:
  1374. return quotfb(a, b);
  1375. #ifdef COMMON
  1376. case TYPE_RATNUM:
  1377. return quotfr(a, b);
  1378. case TYPE_COMPLEX_NUM:
  1379. return quotfc(a, b);
  1380. #endif
  1381. default:
  1382. return aerror1("bad arg for quotient", b);
  1383. }
  1384. }
  1385. case TAG_BOXFLOAT:
  1386. return quotff(a, b);
  1387. default:
  1388. return aerror1("bad arg for quotient", b);
  1389. }
  1390. default:
  1391. return aerror1("bad arg for quotient", a);
  1392. }
  1393. }
  1394. #ifdef COMMON
  1395. /*
  1396. * The above version of quot2 is as required for the Standard Lisp QUOTIENT
  1397. * function, and it is called from various internal places in CSL/CCL, eg
  1398. * from within the code for TRUNCATE. Next I have something that will be very
  1399. * similar in behaviour, but which turns quotients of integers into
  1400. * rational numbers when that is necessary.
  1401. */
  1402. /***** Not reconstructed yet!! */
  1403. Lisp_Object CLquot2(Lisp_Object a, Lisp_Object b)
  1404. {
  1405. switch ((int)a & TAG_BITS)
  1406. {
  1407. case TAG_FIXNUM:
  1408. switch ((int)b & TAG_BITS)
  1409. {
  1410. case TAG_FIXNUM:
  1411. /*
  1412. * This is where fixnum / fixnum arithmetic happens - the case I most want to
  1413. * make efficient.
  1414. */
  1415. if (b == fixnum_of_int(0))
  1416. return aerror2("bad arg for /", a, b);
  1417. else
  1418. { int32 r, aa, bb, w;
  1419. aa = int_of_fixnum(a);
  1420. bb = int_of_fixnum(b);
  1421. if (bb < 0) aa = -aa, bb = -bb;
  1422. r = aa % bb;
  1423. if (r == 0) /* Exact division case */
  1424. { r = aa/bb;
  1425. /*
  1426. * the only case where dividing one fixnum by another can lead to
  1427. * a bignum result is (-0x08000000/(-1)) which JUST overflows.
  1428. */
  1429. if (r != 0x08000000) return fixnum_of_int(r);
  1430. else return make_one_word_bignum(r);
  1431. }
  1432. w = bb;
  1433. if (r < 0) r = -r;
  1434. while (r != 0)
  1435. { int32 w1 = w % r;
  1436. w = r;
  1437. r = w1;
  1438. }
  1439. aa = aa / w;
  1440. bb = bb / w;
  1441. return make_ratio(fixnum_of_int(aa), fixnum_of_int(bb));
  1442. }
  1443. #ifdef COMMON
  1444. case TAG_SFLOAT:
  1445. return quotis(a, b);
  1446. #endif
  1447. case TAG_NUMBERS:
  1448. { int32 hb = type_of_header(numhdr(b));
  1449. switch (hb)
  1450. {
  1451. case TYPE_BIGNUM:
  1452. return CLquotib(a, b);
  1453. #ifdef COMMON
  1454. case TYPE_RATNUM:
  1455. return quotir(a, b);
  1456. case TYPE_COMPLEX_NUM:
  1457. return quotic(a, b);
  1458. #endif
  1459. default:
  1460. return aerror1("bad arg for /", b);
  1461. }
  1462. }
  1463. case TAG_BOXFLOAT:
  1464. return quotif(a, b);
  1465. default:
  1466. return aerror1("bad arg for /", b);
  1467. }
  1468. #ifdef COMMON
  1469. case TAG_SFLOAT:
  1470. switch ((int)b & TAG_BITS)
  1471. {
  1472. case TAG_FIXNUM:
  1473. return quotsi(a, b);
  1474. case TAG_SFLOAT:
  1475. { Float_union aa, bb;
  1476. aa.i = a - TAG_SFLOAT;
  1477. bb.i = b - TAG_SFLOAT;
  1478. aa.f = (float) (aa.f / bb.f);
  1479. return (aa.i & ~(int32)0xf) + TAG_SFLOAT;
  1480. }
  1481. case TAG_NUMBERS:
  1482. { int32 hb = type_of_header(numhdr(b));
  1483. switch (hb)
  1484. {
  1485. case TYPE_BIGNUM:
  1486. return quotsb(a, b);
  1487. case TYPE_RATNUM:
  1488. return quotsr(a, b);
  1489. case TYPE_COMPLEX_NUM:
  1490. return quotsc(a, b);
  1491. default:
  1492. return aerror1("bad arg for /", b);
  1493. }
  1494. }
  1495. case TAG_BOXFLOAT:
  1496. return quotsf(a, b);
  1497. default:
  1498. return aerror1("bad arg for /", b);
  1499. }
  1500. #endif
  1501. case TAG_NUMBERS:
  1502. { int32 ha = type_of_header(numhdr(a));
  1503. switch (ha)
  1504. {
  1505. case TYPE_BIGNUM:
  1506. switch ((int)b & TAG_BITS)
  1507. {
  1508. case TAG_FIXNUM:
  1509. return CLquotbi(a, b);
  1510. #ifdef COMMON
  1511. case TAG_SFLOAT:
  1512. return quotbs(a, b);
  1513. #endif
  1514. case TAG_NUMBERS:
  1515. { int32 hb = type_of_header(numhdr(b));
  1516. switch (hb)
  1517. {
  1518. case TYPE_BIGNUM:
  1519. return CLquotbb(a, b);
  1520. #ifdef COMMON
  1521. case TYPE_RATNUM:
  1522. return quotbr(a, b);
  1523. case TYPE_COMPLEX_NUM:
  1524. return quotbc(a, b);
  1525. #endif
  1526. default:
  1527. return aerror1("bad arg for /", b);
  1528. }
  1529. }
  1530. case TAG_BOXFLOAT:
  1531. return quotbf(a, b);
  1532. default:
  1533. return aerror1("bad arg for /", b);
  1534. }
  1535. #ifdef COMMON
  1536. case TYPE_RATNUM:
  1537. switch ((int)b & TAG_BITS)
  1538. {
  1539. case TAG_FIXNUM:
  1540. return quotri(a, b);
  1541. case TAG_SFLOAT:
  1542. return quotrs(a, b);
  1543. case TAG_NUMBERS:
  1544. { int32 hb = type_of_header(numhdr(b));
  1545. switch (hb)
  1546. {
  1547. case TYPE_BIGNUM:
  1548. return quotrb(a, b);
  1549. case TYPE_RATNUM:
  1550. return quotrr(a, b);
  1551. case TYPE_COMPLEX_NUM:
  1552. return quotrc(a, b);
  1553. default:
  1554. return aerror1("bad arg for /", b);
  1555. }
  1556. }
  1557. case TAG_BOXFLOAT:
  1558. return quotrf(a, b);
  1559. default:
  1560. return aerror1("bad arg for /", b);
  1561. }
  1562. case TYPE_COMPLEX_NUM:
  1563. switch ((int)b & TAG_BITS)
  1564. {
  1565. case TAG_FIXNUM:
  1566. return quotci(a, b);
  1567. case TAG_SFLOAT:
  1568. return quotcs(a, b);
  1569. case TAG_NUMBERS:
  1570. { int32 hb = type_of_header(numhdr(b));
  1571. switch (hb)
  1572. {
  1573. case TYPE_BIGNUM:
  1574. return quotcb(a, b);
  1575. case TYPE_RATNUM:
  1576. return quotcr(a, b);
  1577. case TYPE_COMPLEX_NUM:
  1578. return quotcc(a, b);
  1579. default:
  1580. return aerror1("bad arg for /", b);
  1581. }
  1582. }
  1583. case TAG_BOXFLOAT:
  1584. return quotcf(a, b);
  1585. default:
  1586. return aerror1("bad arg for /", b);
  1587. }
  1588. #endif
  1589. default: return aerror1("bad arg for /", a);
  1590. }
  1591. }
  1592. case TAG_BOXFLOAT:
  1593. switch ((int)b & TAG_BITS)
  1594. {
  1595. case TAG_FIXNUM:
  1596. return quotfi(a, b);
  1597. #ifdef COMMON
  1598. case TAG_SFLOAT:
  1599. return quotfs(a, b);
  1600. #endif
  1601. case TAG_NUMBERS:
  1602. { int32 hb = type_of_header(numhdr(b));
  1603. switch (hb)
  1604. {
  1605. case TYPE_BIGNUM:
  1606. return quotfb(a, b);
  1607. #ifdef COMMON
  1608. case TYPE_RATNUM:
  1609. return quotfr(a, b);
  1610. case TYPE_COMPLEX_NUM:
  1611. return quotfc(a, b);
  1612. #endif
  1613. default:
  1614. return aerror1("bad arg for /", b);
  1615. }
  1616. }
  1617. case TAG_BOXFLOAT:
  1618. return quotff(a, b);
  1619. default:
  1620. return aerror1("bad arg for /", b);
  1621. }
  1622. default:
  1623. return aerror1("bad arg for /", a);
  1624. }
  1625. }
  1626. #endif
  1627. /* end of arith03.c */