Mathlib.mod 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576
  1. (*
  2. (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
  3. See the copyright notice in the ACK home directory, in the file "Copyright".
  4. *)
  5. (*$R-*)
  6. IMPLEMENTATION MODULE Mathlib;
  7. (*
  8. Module: Mathematical functions
  9. Author: Ceriel J.H. Jacobs
  10. Version: $Id$
  11. *)
  12. FROM EM IMPORT FIF, FEF;
  13. FROM Traps IMPORT Message;
  14. CONST
  15. OneRadianInDegrees = 57.295779513082320876798155D;
  16. OneDegreeInRadians = 0.017453292519943295769237D;
  17. OneOverSqrt2 = 0.70710678118654752440084436210484904D;
  18. (* basic functions *)
  19. PROCEDURE pow(x: REAL; i: INTEGER): REAL;
  20. BEGIN
  21. RETURN SHORT(longpow(LONG(x), i));
  22. END pow;
  23. PROCEDURE longpow(x: LONGREAL; i: INTEGER): LONGREAL;
  24. VAR val: LONGREAL;
  25. ri: LONGREAL;
  26. BEGIN
  27. ri := FLOATD(i);
  28. IF x < 0.0D THEN
  29. val := longexp(longln(-x) * ri);
  30. IF ODD(i) THEN RETURN -val;
  31. ELSE RETURN val;
  32. END;
  33. ELSIF x = 0.0D THEN
  34. RETURN 0.0D;
  35. ELSE
  36. RETURN longexp(longln(x) * ri);
  37. END;
  38. END longpow;
  39. PROCEDURE sqrt(x: REAL): REAL;
  40. BEGIN
  41. RETURN SHORT(longsqrt(LONG(x)));
  42. END sqrt;
  43. PROCEDURE longsqrt(x: LONGREAL): LONGREAL;
  44. VAR
  45. temp: LONGREAL;
  46. exp, i: INTEGER;
  47. BEGIN
  48. IF x <= 0.0D THEN
  49. IF x < 0.0D THEN
  50. Message("sqrt: negative argument");
  51. HALT
  52. END;
  53. RETURN 0.0D;
  54. END;
  55. temp := FEF(x,exp);
  56. (*
  57. * NOTE
  58. * this wont work on 1's comp
  59. *)
  60. IF ODD(exp) THEN
  61. temp := 2.0D * temp;
  62. DEC(exp);
  63. END;
  64. temp := 0.5D*(1.0D + temp);
  65. WHILE exp > 28 DO
  66. temp := temp * 16384.0D;
  67. exp := exp - 28;
  68. END;
  69. WHILE exp < -28 DO
  70. temp := temp / 16384.0D;
  71. exp := exp + 28;
  72. END;
  73. WHILE exp >= 2 DO
  74. temp := temp * 2.0D;
  75. exp := exp - 2;
  76. END;
  77. WHILE exp <= -2 DO
  78. temp := temp / 2.0D;
  79. exp := exp + 2;
  80. END;
  81. FOR i := 0 TO 5 DO
  82. temp := 0.5D*(temp + x/temp);
  83. END;
  84. RETURN temp;
  85. END longsqrt;
  86. PROCEDURE ldexp(x:LONGREAL; n: INTEGER): LONGREAL;
  87. BEGIN
  88. WHILE n >= 16 DO
  89. x := x * 65536.0D;
  90. n := n - 16;
  91. END;
  92. WHILE n > 0 DO
  93. x := x * 2.0D;
  94. DEC(n);
  95. END;
  96. WHILE n <= -16 DO
  97. x := x / 65536.0D;
  98. n := n + 16;
  99. END;
  100. WHILE n < 0 DO
  101. x := x / 2.0D;
  102. INC(n);
  103. END;
  104. RETURN x;
  105. END ldexp;
  106. PROCEDURE exp(x: REAL): REAL;
  107. BEGIN
  108. RETURN SHORT(longexp(LONG(x)));
  109. END exp;
  110. PROCEDURE longexp(x: LONGREAL): LONGREAL;
  111. (* Algorithm and coefficients from:
  112. "Software manual for the elementary functions"
  113. by W.J. Cody and W. Waite, Prentice-Hall, 1980
  114. *)
  115. CONST
  116. p0 = 0.25000000000000000000D+00;
  117. p1 = 0.75753180159422776666D-02;
  118. p2 = 0.31555192765684646356D-04;
  119. q0 = 0.50000000000000000000D+00;
  120. q1 = 0.56817302698551221787D-01;
  121. q2 = 0.63121894374398503557D-03;
  122. q3 = 0.75104028399870046114D-06;
  123. VAR
  124. neg: BOOLEAN;
  125. n: INTEGER;
  126. xn, g, x1, x2: LONGREAL;
  127. BEGIN
  128. neg := x < 0.0D;
  129. IF neg THEN
  130. x := -x;
  131. END;
  132. n := TRUNC(x/longln2 + 0.5D);
  133. xn := FLOATD(n);
  134. x1 := FLOATD(TRUNCD(x));
  135. x2 := x - x1;
  136. g := ((x1 - xn * 0.693359375D)+x2) - xn * (-2.1219444005469058277D-4);
  137. IF neg THEN
  138. g := -g;
  139. n := -n;
  140. END;
  141. xn := g*g;
  142. x := g*((p2*xn+p1)*xn+p0);
  143. INC(n);
  144. RETURN ldexp(0.5D + x/((((q3*xn+q2)*xn+q1)*xn+q0) - x), n);
  145. END longexp;
  146. PROCEDURE ln(x: REAL): REAL; (* natural log *)
  147. BEGIN
  148. RETURN SHORT(longln(LONG(x)));
  149. END ln;
  150. PROCEDURE longln(x: LONGREAL): LONGREAL; (* natural log *)
  151. (* Algorithm and coefficients from:
  152. "Software manual for the elementary functions"
  153. by W.J. Cody and W. Waite, Prentice-Hall, 1980
  154. *)
  155. CONST
  156. p0 = -0.64124943423745581147D+02;
  157. p1 = 0.16383943563021534222D+02;
  158. p2 = -0.78956112887491257267D+00;
  159. q0 = -0.76949932108494879777D+03;
  160. q1 = 0.31203222091924532844D+03;
  161. q2 = -0.35667977739034646171D+02;
  162. q3 = 1.0D;
  163. VAR
  164. exp: INTEGER;
  165. z, znum, zden, w: LONGREAL;
  166. BEGIN
  167. IF x <= 0.0D THEN
  168. Message("ln: argument <= 0");
  169. HALT
  170. END;
  171. x := FEF(x, exp);
  172. IF x > OneOverSqrt2 THEN
  173. znum := (x - 0.5D) - 0.5D;
  174. zden := x * 0.5D + 0.5D;
  175. ELSE
  176. znum := x - 0.5D;
  177. zden := znum * 0.5D + 0.5D;
  178. DEC(exp);
  179. END;
  180. z := znum / zden;
  181. w := z * z;
  182. x := z + z * w * (((p2*w+p1)*w+p0)/(((q3*w+q2)*w+q1)*w+q0));
  183. z := FLOATD(exp);
  184. x := x + z * (-2.121944400546905827679D-4);
  185. RETURN x + z * 0.693359375D;
  186. END longln;
  187. PROCEDURE log(x: REAL): REAL; (* log with base 10 *)
  188. BEGIN
  189. RETURN SHORT(longlog(LONG(x)));
  190. END log;
  191. PROCEDURE longlog(x: LONGREAL): LONGREAL; (* log with base 10 *)
  192. BEGIN
  193. RETURN longln(x)/longln10;
  194. END longlog;
  195. (* trigonometric functions; arguments in radians *)
  196. PROCEDURE sin(x: REAL): REAL;
  197. BEGIN
  198. RETURN SHORT(longsin(LONG(x)));
  199. END sin;
  200. PROCEDURE sinus(x: LONGREAL; cosflag: BOOLEAN) : LONGREAL;
  201. (* Algorithm and coefficients from:
  202. "Software manual for the elementary functions"
  203. by W.J. Cody and W. Waite, Prentice-Hall, 1980
  204. *)
  205. CONST
  206. r0 = -0.16666666666666665052D+00;
  207. r1 = 0.83333333333331650314D-02;
  208. r2 = -0.19841269841201840457D-03;
  209. r3 = 0.27557319210152756119D-05;
  210. r4 = -0.25052106798274584544D-07;
  211. r5 = 0.16058936490371589114D-09;
  212. r6 = -0.76429178068910467734D-12;
  213. r7 = 0.27204790957888846175D-14;
  214. A1 = 3.1416015625D;
  215. A2 = -8.908910206761537356617D-6;
  216. VAR
  217. x1, x2, y : LONGREAL;
  218. neg : BOOLEAN;
  219. BEGIN
  220. IF x < 0.0D THEN
  221. neg := TRUE;
  222. x := -x
  223. ELSE neg := FALSE
  224. END;
  225. IF cosflag THEN
  226. neg := FALSE;
  227. y := longhalfpi + x
  228. ELSE
  229. y := x
  230. END;
  231. y := y / longpi + 0.5D;
  232. IF FIF(y, 1.0D, y) < 0.0D THEN ; END;
  233. IF FIF(y, 0.5D, x1) # 0.0D THEN neg := NOT neg END;
  234. IF cosflag THEN y := y - 0.5D END;
  235. x2 := FIF(x, 1.0, x1);
  236. x := x1 - y * A1;
  237. x := x + x2;
  238. x := x - y * A2;
  239. IF x < 0.0D THEN
  240. neg := NOT neg;
  241. x := -x
  242. END;
  243. y := x * x;
  244. x := x + x * y * (((((((r7*y+r6)*y+r5)*y+r4)*y+r3)*y+r2)*y+r1)*y+r0);
  245. IF neg THEN RETURN -x END;
  246. RETURN x;
  247. END sinus;
  248. PROCEDURE longsin(x: LONGREAL): LONGREAL;
  249. BEGIN
  250. RETURN sinus(x, FALSE);
  251. END longsin;
  252. PROCEDURE cos(x: REAL): REAL;
  253. BEGIN
  254. RETURN SHORT(longcos(LONG(x)));
  255. END cos;
  256. PROCEDURE longcos(x: LONGREAL): LONGREAL;
  257. BEGIN
  258. IF x < 0.0D THEN x := -x; END;
  259. RETURN sinus(x, TRUE);
  260. END longcos;
  261. PROCEDURE tan(x: REAL): REAL;
  262. BEGIN
  263. RETURN SHORT(longtan(LONG(x)));
  264. END tan;
  265. PROCEDURE longtan(x: LONGREAL): LONGREAL;
  266. (* Algorithm and coefficients from:
  267. "Software manual for the elementary functions"
  268. by W.J. Cody and W. Waite, Prentice-Hall, 1980
  269. *)
  270. CONST
  271. p1 = -0.13338350006421960681D+00;
  272. p2 = 0.34248878235890589960D-02;
  273. p3 = -0.17861707342254426711D-04;
  274. q0 = 1.0D;
  275. q1 = -0.46671683339755294240D+00;
  276. q2 = 0.25663832289440112864D-01;
  277. q3 = -0.31181531907010027307D-03;
  278. q4 = 0.49819433993786512270D-06;
  279. A1 = 1.57080078125D;
  280. A2 = -4.454455103380768678308D-06;
  281. VAR y, x1, x2: LONGREAL;
  282. negative: BOOLEAN;
  283. invert: BOOLEAN;
  284. BEGIN
  285. negative := x < 0.0D;
  286. y := x / longhalfpi + 0.5D;
  287. (* Use extended precision to calculate reduced argument.
  288. Here we used 12 bits of the mantissa for a1.
  289. Also split x in integer part x1 and fraction part x2.
  290. *)
  291. IF FIF(y, 1.0D, y) < 0.0D THEN ; END;
  292. invert := FIF(y, 0.5D, x1) # 0.0D;
  293. x2 := FIF(x, 1.0D, x1);
  294. x := x1 - y * A1;
  295. x := x + x2;
  296. x := x - y * A2;
  297. y := x * x;
  298. x := x + x * y * ((p3*y+p2)*y+p1);
  299. y := (((q4*y+q3)*y+q2)*y+q1)*y+q0;
  300. IF negative THEN x := -x END;
  301. IF invert THEN RETURN -y/x END;
  302. RETURN x/y;
  303. END longtan;
  304. PROCEDURE arcsin(x: REAL): REAL;
  305. BEGIN
  306. RETURN SHORT(longarcsin(LONG(x)));
  307. END arcsin;
  308. PROCEDURE arcsincos(x: LONGREAL; cosfl: BOOLEAN): LONGREAL;
  309. CONST
  310. p0 = -0.27368494524164255994D+02;
  311. p1 = 0.57208227877891731407D+02;
  312. p2 = -0.39688862997540877339D+02;
  313. p3 = 0.10152522233806463645D+02;
  314. p4 = -0.69674573447350646411D+00;
  315. q0 = -0.16421096714498560795D+03;
  316. q1 = 0.41714430248260412556D+03;
  317. q2 = -0.38186303361750149284D+03;
  318. q3 = 0.15095270841030604719D+03;
  319. q4 = -0.23823859153670238830D+02;
  320. q5 = 1.0D;
  321. VAR
  322. negative : BOOLEAN;
  323. big: BOOLEAN;
  324. g: LONGREAL;
  325. BEGIN
  326. negative := x < 0.0D;
  327. IF negative THEN x := -x; END;
  328. IF x > 0.5D THEN
  329. big := TRUE;
  330. IF x > 1.0D THEN
  331. Message("arcsin or arccos: argument > 1");
  332. HALT
  333. END;
  334. g := 0.5D - 0.5D * x;
  335. x := -longsqrt(g);
  336. x := x + x;
  337. ELSE
  338. big := FALSE;
  339. g := x * x;
  340. END;
  341. x := x + x * g *
  342. ((((p4*g+p3)*g+p2)*g+p1)*g+p0)/(((((q5*g+q4)*g+q3)*g+q2)*g+q1)*g+q0);
  343. IF cosfl AND NOT negative THEN x := -x END;
  344. IF cosfl = NOT big THEN
  345. x := (x + longquartpi) + longquartpi;
  346. ELSIF cosfl AND negative AND big THEN
  347. x := (x + longhalfpi) + longhalfpi;
  348. END;
  349. IF negative AND NOT cosfl THEN x := -x END;
  350. RETURN x;
  351. END arcsincos;
  352. PROCEDURE longarcsin(x: LONGREAL): LONGREAL;
  353. BEGIN
  354. RETURN arcsincos(x, FALSE);
  355. END longarcsin;
  356. PROCEDURE arccos(x: REAL): REAL;
  357. BEGIN
  358. RETURN SHORT(longarccos(LONG(x)));
  359. END arccos;
  360. PROCEDURE longarccos(x: LONGREAL): LONGREAL;
  361. BEGIN
  362. RETURN arcsincos(x, TRUE);
  363. END longarccos;
  364. PROCEDURE arctan(x: REAL): REAL;
  365. BEGIN
  366. RETURN SHORT(longarctan(LONG(x)));
  367. END arctan;
  368. VAR A: ARRAY[0..3] OF LONGREAL;
  369. arctaninit: BOOLEAN;
  370. PROCEDURE longarctan(x: LONGREAL): LONGREAL;
  371. (* Algorithm and coefficients from:
  372. "Software manual for the elementary functions"
  373. by W.J. Cody and W. Waite, Prentice-Hall, 1980
  374. *)
  375. CONST
  376. p0 = -0.13688768894191926929D+02;
  377. p1 = -0.20505855195861651981D+02;
  378. p2 = -0.84946240351320683534D+01;
  379. p3 = -0.83758299368150059274D+00;
  380. q0 = 0.41066306682575781263D+02;
  381. q1 = 0.86157349597130242515D+02;
  382. q2 = 0.59578436142597344465D+02;
  383. q3 = 0.15024001160028576121D+02;
  384. q4 = 1.0D;
  385. VAR
  386. g: LONGREAL;
  387. neg: BOOLEAN;
  388. n: INTEGER;
  389. BEGIN
  390. IF NOT arctaninit THEN
  391. arctaninit := TRUE;
  392. A[0] := 0.0D;
  393. A[1] := 0.52359877559829887307710723554658381D; (* p1/6 *)
  394. A[2] := longhalfpi;
  395. A[3] := 1.04719755119659774615421446109316763D; (* pi/3 *)
  396. END;
  397. neg := FALSE;
  398. IF x < 0.0D THEN
  399. neg := TRUE;
  400. x := -x;
  401. END;
  402. IF x > 1.0D THEN
  403. x := 1.0D/x;
  404. n := 2
  405. ELSE
  406. n := 0
  407. END;
  408. IF x > 0.26794919243112270647D (* 2-sqrt(3) *) THEN
  409. INC(n);
  410. x := (((0.73205080756887729353D*x-0.5D)-0.5D)+x)/
  411. (1.73205080756887729353D + x);
  412. END;
  413. g := x*x;
  414. x := x + x * g * (((p3*g+p2)*g+p1)*g+p0) / ((((q4*g+q3)*g+q2)*g+q1)*g+q0);
  415. IF n > 1 THEN x := -x END;
  416. x := x + A[n];
  417. IF neg THEN RETURN -x; END;
  418. RETURN x;
  419. END longarctan;
  420. (* hyperbolic functions *)
  421. (* The C math library has better implementations for some of these, but
  422. they depend on some properties of the floating point implementation,
  423. and, for now, we don't want that in the Modula-2 system.
  424. *)
  425. PROCEDURE sinh(x: REAL): REAL;
  426. BEGIN
  427. RETURN SHORT(longsinh(LONG(x)));
  428. END sinh;
  429. PROCEDURE longsinh(x: LONGREAL): LONGREAL;
  430. VAR expx: LONGREAL;
  431. BEGIN
  432. expx := longexp(x);
  433. RETURN (expx - 1.0D/expx)/2.0D;
  434. END longsinh;
  435. PROCEDURE cosh(x: REAL): REAL;
  436. BEGIN
  437. RETURN SHORT(longcosh(LONG(x)));
  438. END cosh;
  439. PROCEDURE longcosh(x: LONGREAL): LONGREAL;
  440. VAR expx: LONGREAL;
  441. BEGIN
  442. expx := longexp(x);
  443. RETURN (expx + 1.0D/expx)/2.0D;
  444. END longcosh;
  445. PROCEDURE tanh(x: REAL): REAL;
  446. BEGIN
  447. RETURN SHORT(longtanh(LONG(x)));
  448. END tanh;
  449. PROCEDURE longtanh(x: LONGREAL): LONGREAL;
  450. VAR expx: LONGREAL;
  451. BEGIN
  452. expx := longexp(x);
  453. RETURN (expx - 1.0D/expx) / (expx + 1.0D/expx);
  454. END longtanh;
  455. PROCEDURE arcsinh(x: REAL): REAL;
  456. BEGIN
  457. RETURN SHORT(longarcsinh(LONG(x)));
  458. END arcsinh;
  459. PROCEDURE longarcsinh(x: LONGREAL): LONGREAL;
  460. VAR neg: BOOLEAN;
  461. BEGIN
  462. neg := FALSE;
  463. IF x < 0.0D THEN
  464. neg := TRUE;
  465. x := -x;
  466. END;
  467. x := longln(x + longsqrt(x*x+1.0D));
  468. IF neg THEN RETURN -x; END;
  469. RETURN x;
  470. END longarcsinh;
  471. PROCEDURE arccosh(x: REAL): REAL;
  472. BEGIN
  473. RETURN SHORT(longarccosh(LONG(x)));
  474. END arccosh;
  475. PROCEDURE longarccosh(x: LONGREAL): LONGREAL;
  476. BEGIN
  477. IF x < 1.0D THEN
  478. Message("arccosh: argument < 1");
  479. HALT
  480. END;
  481. RETURN longln(x + longsqrt(x*x - 1.0D));
  482. END longarccosh;
  483. PROCEDURE arctanh(x: REAL): REAL;
  484. BEGIN
  485. RETURN SHORT(longarctanh(LONG(x)));
  486. END arctanh;
  487. PROCEDURE longarctanh(x: LONGREAL): LONGREAL;
  488. BEGIN
  489. IF (x <= -1.0D) OR (x >= 1.0D) THEN
  490. Message("arctanh: ABS(argument) >= 1");
  491. HALT
  492. END;
  493. RETURN longln((1.0D + x)/(1.0D - x)) / 2.0D;
  494. END longarctanh;
  495. (* conversions *)
  496. PROCEDURE RadianToDegree(x: REAL): REAL;
  497. BEGIN
  498. RETURN SHORT(longRadianToDegree(LONG(x)));
  499. END RadianToDegree;
  500. PROCEDURE longRadianToDegree(x: LONGREAL): LONGREAL;
  501. BEGIN
  502. RETURN x * OneRadianInDegrees;
  503. END longRadianToDegree;
  504. PROCEDURE DegreeToRadian(x: REAL): REAL;
  505. BEGIN
  506. RETURN SHORT(longDegreeToRadian(LONG(x)));
  507. END DegreeToRadian;
  508. PROCEDURE longDegreeToRadian(x: LONGREAL): LONGREAL;
  509. BEGIN
  510. RETURN x * OneDegreeInRadians;
  511. END longDegreeToRadian;
  512. BEGIN
  513. arctaninit := FALSE;
  514. END Mathlib.