File: common\mpmath_c.c
1 /* MPMath_c.c (C) 1989, Mark C. Peterson, CompuServe [70441,3353]
2 All rights reserved.
3
4 Code may be used in any program provided the author is credited
5 either during program execution or in the documentation. Source
6 code may be distributed only in combination with public domain or
7 shareware source code. Source code may be modified provided the
8 copyright notice and this message is left unchanged and all
9 modifications are clearly documented.
10
11 I would appreciate a copy of any work which incorporates this code,
12 however this is optional.
13
14 Mark C. Peterson
15 405-C Queen St. Suite #181
16 Southington, CT 06489
17 (203) 276-9721
18 */
19
20
21 /* see Fractint.c for a description of the "include" hierarchy */
22 #include "port.h"
23 #include "prototyp.h"
24
25 #ifndef XFRACT
26 #if (_MSC_VER >= 700)
27 #pragma code_seg ("mpmath1_text") /* place following in an overlay */
28 #endif
29
30 struct MP *MPsub(struct MP x, struct MP y) {
31 y.Exp ^= 0x8000;
32 return(MPadd(x, y));
33 }
34
35 /* added by TW */
36 struct MP *MPsub086(struct MP x, struct MP y) {
37 y.Exp ^= 0x8000;
38 return(MPadd086(x, y));
39 }
40
41 /* added by TW */
42 struct MP *MPsub386(struct MP x, struct MP y) {
43 y.Exp ^= 0x8000;
44 return(MPadd386(x, y));
45 }
46
47 struct MP *MPabs(struct MP x) {
48 Ans = x;
49 Ans.Exp &= 0x7fff;
50 return(&Ans);
51 }
52
53 struct MPC MPCsqr(struct MPC x) {
54 struct MPC z;
55
56 z.x = *pMPsub(*pMPmul(x.x, x.x), *pMPmul(x.y, x.y));
57 z.y = *pMPmul(x.x, x.y);
58 z.y.Exp++;
59 return(z);
60 }
61
62 struct MP MPCmod(struct MPC x) {
63 return(*pMPadd(*pMPmul(x.x, x.x), *pMPmul(x.y, x.y)));
64 }
65
66 struct MPC MPCmul(struct MPC x, struct MPC y) {
67 struct MPC z;
68
69 z.x = *pMPsub(*pMPmul(x.x, y.x), *pMPmul(x.y, y.y));
70 z.y = *pMPadd(*pMPmul(x.x, y.y), *pMPmul(x.y, y.x));
71 return(z);
72 }
73
74 struct MPC MPCdiv(struct MPC x, struct MPC y) {
75 struct MP mod;
76
77 mod = MPCmod(y);
78 y.y.Exp ^= 0x8000;
79 y.x = *pMPdiv(y.x, mod);
80 y.y = *pMPdiv(y.y, mod);
81 return(MPCmul(x, y));
82 }
83
84 struct MPC MPCadd(struct MPC x, struct MPC y) {
85 struct MPC z;
86
87 z.x = *pMPadd(x.x, y.x);
88 z.y = *pMPadd(x.y, y.y);
89 return(z);
90 }
91
92 struct MPC MPCsub(struct MPC x, struct MPC y) {
93 struct MPC z;
94
95 z.x = *pMPsub(x.x, y.x);
96 z.y = *pMPsub(x.y, y.y);
97 return(z);
98 }
99
100 struct MPC MPCone = { {0x3fff, 0x80000000l},
101 {0, 0l}
102 };
103
104 struct MPC MPCpow(struct MPC x, int exp) {
105 struct MPC z;
106 struct MPC zz;
107
108 if(exp & 1)
109 z = x;
110 else
111 z = MPCone;
112 exp >>= 1;
113 while(exp) {
114 zz.x = *pMPsub(*pMPmul(x.x, x.x), *pMPmul(x.y, x.y));
115 zz.y = *pMPmul(x.x, x.y);
116 zz.y.Exp++;
117 x = zz;
118 if(exp & 1) {
119 zz.x = *pMPsub(*pMPmul(z.x, x.x), *pMPmul(z.y, x.y));
120 zz.y = *pMPadd(*pMPmul(z.x, x.y), *pMPmul(z.y, x.x));
121 z = zz;
122 }
123 exp >>= 1;
124 }
125 return(z);
126 }
127
128 int MPCcmp(struct MPC x, struct MPC y) {
129 struct MPC z;
130
131 if(pMPcmp(x.x, y.x) || pMPcmp(x.y, y.y)) {
132 z.x = MPCmod(x);
133 z.y = MPCmod(y);
134 return(pMPcmp(z.x, z.y));
135 }
136 else
137 return(0);
138 }
139
140 _CMPLX MPC2cmplx(struct MPC x) {
141 _CMPLX z;
142
143 z.x = *pMP2d(x.x);
144 z.y = *pMP2d(x.y);
145 return(z);
146 }
147
148 struct MPC cmplx2MPC(_CMPLX z) {
149 struct MPC x;
150
151 x.x = *pd2MP(z.x);
152 x.y = *pd2MP(z.y);
153 return(x);
154 }
155
156 /* function pointer versions added by Tim Wegner 12/07/89 */
157 /* int (*ppMPcmp)() = MPcmp086; */
158 int (*pMPcmp)(struct MP x, struct MP y) = MPcmp086;
159 struct MP *(*pMPmul)(struct MP x, struct MP y)= MPmul086;
160 struct MP *(*pMPdiv)(struct MP x, struct MP y)= MPdiv086;
161 struct MP *(*pMPadd)(struct MP x, struct MP y)= MPadd086;
162 struct MP *(*pMPsub)(struct MP x, struct MP y)= MPsub086;
163 struct MP *(*pd2MP)(double x) = d2MP086 ;
164 double *(*pMP2d)(struct MP m) = MP2d086 ;
165 /* struct MP *(*pfg2MP)(long x, int fg) = fg2MP086; */
166
167 void setMPfunctions(void) {
168 if(cpu >= 386)
169 {
170 pMPmul = MPmul386;
171 pMPdiv = MPdiv386;
172 pMPadd = MPadd386;
173 pMPsub = MPsub386;
174 pMPcmp = MPcmp386;
175 pd2MP = d2MP386 ;
176 pMP2d = MP2d386 ;
177 /* pfg2MP = fg2MP386; */
178 }
179 else
180 {
181 pMPmul = MPmul086;
182 pMPdiv = MPdiv086;
183 pMPadd = MPadd086;
184 pMPsub = MPsub086;
185 pMPcmp = MPcmp086;
186 pd2MP = d2MP086 ;
187 pMP2d = MP2d086 ;
188 /* pfg2MP = fg2MP086; */
189 }
190 }
191 #if (_MSC_VER >= 700)
192 #pragma code_seg () /* back to normal segment */
193 #endif
194 #endif /* XFRACT */
195
196 #ifndef sqr
198 #endif
199
200 _CMPLX ComplexPower(_CMPLX xx, _CMPLX yy) {
201 _CMPLX z, cLog, t;
202 double e2x, siny, cosy;
203
204 /* fixes power bug - if any complaints, backwards compatibility hook
205 goes here TIW 3/95 */
206 if(ldcheck == 0)
207 if(xx.x == 0 && xx.y == 0) {
208 z.x = z.y = 0.0;
209 return(z);
210 }
211
212 FPUcplxlog(&xx, &cLog);
213 FPUcplxmul(&cLog, &yy, &t);
214
215 if(fpu >= 387)
216 FPUcplxexp387(&t, &z);
217 else {
218 if(t.x < -690)
219 e2x = 0;
220 else
221 e2x = exp(t.x);
222 FPUsincos(&t.y, &siny, &cosy);
223 z.x = e2x * cosy;
224 z.y = e2x * siny;
225 }
226 return(z);
227 }
228
229 /*
230
231 The following Complex function routines added by Tim Wegner November 1994.
232
233 */
234
235 #define Sqrtz(z,rz) (*(rz) = ComplexSqrtFloat((z).x, (z).y))
236
237 /* rz=Arcsin(z)=-i*Log{i*z+sqrt(1-z*z)} */
238 void Arcsinz(_CMPLX z,_CMPLX *rz)
239 {
240 _CMPLX tempz1,tempz2;
241
242 FPUcplxmul( &z, &z, &tempz1);
243 tempz1.x = 1 - tempz1.x; tempz1.y = -tempz1.y; /* tempz1 = 1 - tempz1 */
244 Sqrtz( tempz1, &tempz1);
245
246 tempz2.x = -z.y; tempz2.y = z.x; /* tempz2 = i*z */
247 tempz1.x += tempz2.x; tempz1.y += tempz2.y; /* tempz1 += tempz2 */
248 FPUcplxlog( &tempz1, &tempz1);
249 rz->x = tempz1.y; rz->y = -tempz1.x; /* rz = (-i)*tempz1 */
250 } /* end. Arcsinz */
251
252
253 /* rz=Arccos(z)=-i*Log{z+sqrt(z*z-1)} */
254 void Arccosz(_CMPLX z,_CMPLX *rz)
255 {
256 _CMPLX temp;
257
258 FPUcplxmul( &z, &z, &temp);
259 temp.x -= 1; /* temp = temp - 1 */
260 Sqrtz( temp, &temp);
261
262 temp.x += z.x; temp.y += z.y; /* temp = z + temp */
263
264 FPUcplxlog( &temp, &temp);
265 rz->x = temp.y; rz->y = -temp.x; /* rz = (-i)*tempz1 */
266 } /* end. Arccosz */
267
268 void Arcsinhz(_CMPLX z,_CMPLX *rz)
269 {
270 _CMPLX temp;
271
272 FPUcplxmul( &z, &z, &temp);
273 temp.x += 1; /* temp = temp + 1 */
274 Sqrtz( temp, &temp);
275 temp.x += z.x; temp.y += z.y; /* temp = z + temp */
276 FPUcplxlog( &temp, rz);
277 } /* end. Arcsinhz */
278
279 /* rz=Arccosh(z)=Log(z+sqrt(z*z-1)} */
280 void Arccoshz(_CMPLX z,_CMPLX *rz)
281 {
282 _CMPLX tempz;
283 FPUcplxmul( &z, &z, &tempz);
284 tempz.x -= 1; /* tempz = tempz - 1 */
285 Sqrtz( tempz, &tempz);
286 tempz.x = z.x + tempz.x; tempz.y = z.y + tempz.y; /* tempz = z + tempz */
287 FPUcplxlog( &tempz, rz);
288 } /* end. Arccoshz */
289
290 /* rz=Arctanh(z)=1/2*Log{(1+z)/(1-z)} */
291 void Arctanhz(_CMPLX z,_CMPLX *rz)
292 {
293 _CMPLX temp0,temp1,temp2;
294
295 if( z.x == 0.0){
296 rz->x = 0;
297 rz->y = atan( z.y);
298 return;
299 }
300 else{
301 if( fabs(z.x) == 1.0 && z.y == 0.0){
302 return;
303 }
304 else if( fabs( z.x) < 1.0 && z.y == 0.0){
305 rz->x = log((1+z.x)/(1-z.x))/2;
306 rz->y = 0;
307 return;
308 }
309 else{
310 temp0.x = 1 + z.x; temp0.y = z.y; /* temp0 = 1 + z */
311 temp1.x = 1 - z.x; temp1.y = -z.y; /* temp1 = 1 - z */
312 FPUcplxdiv( &temp0, &temp1, &temp2);
313 FPUcplxlog( &temp2, &temp2);
314 rz->x = .5*temp2.x; rz->y = .5*temp2.y; /* rz = .5*temp2 */
315 return;
316 }
317 }
318 } /* end. Arctanhz */
319
320 /* rz=Arctan(z)=i/2*Log{(1-i*z)/(1+i*z)} */
321 void Arctanz(_CMPLX z,_CMPLX *rz)
322 {
323 _CMPLX temp0,temp1,temp2,temp3;
324 if( z.x == 0.0 && z.y == 0.0)
325 rz->x = rz->y = 0;
326 else if( z.x != 0.0 && z.y == 0.0){
327 rz->x = atan( z.x);
328 rz->y = 0;
329 }
330 else if( z.x == 0.0 && z.y != 0.0){
331 temp0.x = z.y; temp0.y = 0.0;
332 Arctanhz( temp0, &temp0);
333 rz->x = -temp0.y; rz->y = temp0.x; /* i*temp0 */
334 }
335 else if( z.x != 0.0 && z.y != 0.0){
336
337 temp0.x = -z.y; temp0.y = z.x; /* i*z */
338 temp1.x = 1 - temp0.x; temp1.y = -temp0.y; /* temp1 = 1 - temp0 */
339 temp2.x = 1 + temp0.x; temp2.y = temp0.y; /* temp2 = 1 + temp0 */
340
341 FPUcplxdiv( &temp1, &temp2, &temp3);
342 FPUcplxlog( &temp3, &temp3);
343 rz->x = -temp3.y*.5; rz->y = .5*temp3.x; /* .5*i*temp0 */
344 }
345 } /* end. Arctanz */
346
347 #define SinCosFudge 0x10000L
348 #ifdef LONGSQRT
349 long lsqrt(long f)
350 {
351 int N;
352 unsigned long y0, z;
353 static long a=0, b=0, c=0; /* constant factors */
354
355 if (f == 0)
356 return f;
357 if (f < 0)
358 return 0;
359
360 if (a==0) /* one-time compute consts */
361 {
362 a = (long)(fudge * .41731);
363 b = (long)(fudge * .59016);
364 c = (long)(fudge * .7071067811);
365 }
366
367 N = 0;
368 while (f & 0xff000000L) /* shift arg f into the */
369 { /* range: 0.5 <= f < 1 */
370 N++;
371 f /= 2;
372 }
373 while (!(f & 0xff800000L))
374 {
375 N--;
376 f *= 2;
377 }
378
379 y0 = a + multiply(b, f, bitshift); /* Newton's approximation */
380
381 z = y0 + divide (f, y0, bitshift);
382 y0 = (z>>2) + divide(f, z, bitshift);
383
384 if (N % 2)
385 {
386 N++;
387 y0 = multiply(c,y0, bitshift);
388 }
389 N /= 2;
390 if (N >= 0)
391 return y0 << N; /* correct for shift above */
392 else
393 return y0 >> -N;
394 } 395 #endif
396 LCMPLX ComplexSqrtLong(long x, long y)
397 {
398 double mag, theta;
399 long maglong, thetalong;
400 LCMPLX result;
401
402 #ifndef LONGSQRT
403 mag = sqrt(sqrt(((double) multiply(x,x,bitshift))/fudge +
404 ((double) multiply(y,y,bitshift))/ fudge));
405 maglong = (long)(mag * fudge);
408 #endif
409 theta = atan2((double) y/fudge, (double) x/fudge)/2;
410 thetalong = (long)(theta * SinCosFudge);
411 SinCos086(thetalong, &result.y, &result.x);
412 result.x = multiply(result.x << (bitshift - 16), maglong, bitshift);
413 result.y = multiply(result.y << (bitshift - 16), maglong, bitshift);
414 return result;
415 }
416
417 _CMPLX ComplexSqrtFloat(double x, double y)
418 {
419 double mag;
420 double theta;
421 _CMPLX result;
422
423 if(x == 0.0 && y == 0.0)
424 result.x = result.y = 0.0;
425 else
426 {
427 mag = sqrt(sqrt(x*x + y*y));
428 theta = atan2(y, x) / 2;
429 FPUsincos(&theta, &result.y, &result.x);
430 result.x *= mag;
431 result.y *= mag;
432 }
433 return result;
434 }
435
436
437 /***** FRACTINT specific routines and variables *****/
438
439 #ifndef TESTING_MATH
440
441 BYTE far *LogTable = (BYTE far *)0;
442 long MaxLTSize;
443 int Log_Calc = 0;
444 static double mlf;
445 static unsigned long lf;
446
447 /* int LogFlag;
448 LogFlag == 1 -- standard log palettes
449 LogFlag == -1 -- 'old' log palettes
450 LogFlag > 1 -- compress counts < LogFlag into color #1
451 LogFlag < -1 -- use quadratic palettes based on square roots && compress
452 */
453
454 void SetupLogTable(void) {
455 float l, f, c, m;
456 unsigned long prev, limit, sptop;
457 unsigned n;
458
459 if (save_release > 1920 || Log_Fly_Calc == 1) { /* set up on-the-fly variables */
460 if (LogFlag > 0) { /* new log function */
461 lf = (LogFlag > 1) ? LogFlag : 0;
462 if (lf >= (unsigned long)MaxLTSize)
463 lf = MaxLTSize - 1;
464 mlf = (colors - (lf?2:1)) / log(MaxLTSize - lf);
465 } else if (LogFlag == -1) { /* old log function */
466 mlf = (colors - 1) / log(MaxLTSize);
467 } else if (LogFlag <= -2) { /* sqrt function */
468 if ((lf = 0 - LogFlag) >= (unsigned long)MaxLTSize)
469 lf = MaxLTSize - 1;
470 mlf = (colors - 2) / sqrt(MaxLTSize - lf);
471 }
472 }
473
474 if (Log_Calc)
475 return; /* LogTable not defined, bail out now */
476
477 if (save_release > 1920 && !Log_Calc) {
478 Log_Calc = 1; /* turn it on */
479 for (prev = 0; prev <= (unsigned long)MaxLTSize; prev++)
480 LogTable[prev] = (BYTE)logtablecalc((long)prev);
481 Log_Calc = 0; /* turn it off, again */
482 return;
483 }
484
485 if (LogFlag > -2) {
486 lf = (LogFlag > 1) ? LogFlag : 0;
487 if (lf >= (unsigned long)MaxLTSize)
488 lf = MaxLTSize - 1;
489 Fg2Float((long)(MaxLTSize-lf), 0, m);
490 fLog14(m, m);
491 Fg2Float((long)(colors-(lf?2:1)), 0, c);
492 fDiv(m, c, m);
493 for (prev = 1; prev <= lf; prev++)
494 LogTable[prev] = 1;
495 for (n = (lf?2:1); n < (unsigned int)colors; n++) {
496 Fg2Float((long)n, 0, f);
497 fMul16(f, m, f);
498 fExp14(f, l);
499 limit = (unsigned long)Float2Fg(l, 0) + lf;
500 if (limit > (unsigned long)MaxLTSize || n == (unsigned int)(colors-1))
501 limit = MaxLTSize;
502 while (prev <= limit)
503 LogTable[prev++] = (BYTE)n;
504 }
505 } else {
506 if ((lf = 0 - LogFlag) >= (unsigned long)MaxLTSize)
507 lf = MaxLTSize - 1;
508 Fg2Float((long)(MaxLTSize-lf), 0, m);
509 fSqrt14(m, m);
510 Fg2Float((long)(colors-2), 0, c);
511 fDiv(m, c, m);
512 for (prev = 1; prev <= lf; prev++)
513 LogTable[prev] = 1;
514 for (n = 2; n < (unsigned int)colors; n++) {
515 Fg2Float((long)n, 0, f);
516 fMul16(f, m, f);
517 fMul16(f, f, l);
518 limit = (unsigned long)(Float2Fg(l, 0) + lf);
519 if (limit > (unsigned long)MaxLTSize || n == (unsigned int)(colors-1))
520 limit = MaxLTSize;
521 while (prev <= limit)
522 LogTable[prev++] = (BYTE)n;
523 }
524 }
525 LogTable[0] = 0;
526 if (LogFlag != -1)
527 for (sptop = 1; sptop < (unsigned long)MaxLTSize; sptop++) /* spread top to incl unused colors */
528 if (LogTable[sptop] > LogTable[sptop-1])
529 LogTable[sptop] = (BYTE)(LogTable[sptop-1]+1);
530 }
531
532 long logtablecalc(long citer) {
533 long ret = 0;
534
535 if (LogFlag == 0 && !rangeslen) /* Oops, how did we get here? */
536 return(citer);
537 if (LogTable && !Log_Calc)
538 return(LogTable[(long)min(citer, MaxLTSize)]);
539
540 if (LogFlag > 0) { /* new log function */
541 if ((unsigned long)citer <= lf + 1)
542 ret = 1;
543 else if((citer - lf) / log(citer - lf) <= mlf) {
544 if (save_release < 2002)
545 ret = (long)(citer - lf + (lf?1:0));
546 else
547 ret = (long)(citer - lf);
548 }
549 else
550 ret = (long)(mlf * log(citer - lf)) + 1;
551 } else if (LogFlag == -1) { /* old log function */
552 if (citer == 0)
553 ret = 1;
554 else
555 ret = (long)(mlf * log(citer)) + 1;
556 } else if (LogFlag <= -2) { /* sqrt function */
557 if ((unsigned long)citer <= lf)
558 ret = 1;
559 else if((unsigned long)(citer - lf) <= (unsigned long)(mlf * mlf))
560 ret = (long)(citer - lf + 1);
561 else
562 ret = (long)(mlf * sqrt(citer - lf)) + 1;
563 }
564 return (ret);
565 }
566
567 long far ExpFloat14(long xx) {
568 static float fLogTwo = (float)0.6931472;
569 int f;
570 long Ans;
571
572 f = 23 - (int)RegFloat2Fg(RegDivFloat(xx, *(long*)&fLogTwo), 0);
573 Ans = ExpFudged(RegFloat2Fg(xx, 16), f);
574 return(RegFg2Float(Ans, (char)f));
575 }
576
577 double TwoPi;
578 _CMPLX temp, BaseLog;
579 _CMPLX cdegree = { 3.0, 0.0 }, croot = { 1.0, 0.0 };
580
581 int ComplexNewtonSetup(void) {
582 threshold = .001;
583 periodicitycheck = 0;
584 if(param[0] != 0.0 || param[1] != 0.0 || param[2] != 0.0 ||
585 param[3] != 0.0) {
586 croot.x = param[2];
587 croot.y = param[3];
588 cdegree.x = param[0];
589 cdegree.y = param[1];
590 FPUcplxlog(&croot, &BaseLog);
591 TwoPi = asin(1.0) * 4;
592 }
593 return(1);
594 }
595
596 int ComplexNewton(void) {
597 _CMPLX cd1;
598
599 /* new = ((cdegree-1) * old**cdegree) + croot
600 ----------------------------------
601 cdegree * old**(cdegree-1) */
602
603 cd1.x = cdegree.x - 1.0;
604 cd1.y = cdegree.y;
605
606 temp = ComplexPower(old, cd1);
607 FPUcplxmul(&temp, &old, &new);
608
609 tmp.x = new.x - croot.x;
610 tmp.y = new.y - croot.y;
611 if((sqr(tmp.x) + sqr(tmp.y)) < threshold)
612 return(1);
613
614 FPUcplxmul(&new, &cd1, &tmp);
615 tmp.x += croot.x;
616 tmp.y += croot.y;
617
618 FPUcplxmul(&temp, &cdegree, &cd1);
619 FPUcplxdiv(&tmp, &cd1, &old);
620 if(overflow)
621 {
622 return(1);
623 }
624 new = old;
625 return(0);
626 }
627
628 int ComplexBasin(void) {
629 _CMPLX cd1;
630 double mod;
631
632 /* new = ((cdegree-1) * old**cdegree) + croot
633 ----------------------------------
634 cdegree * old**(cdegree-1) */
635
636 cd1.x = cdegree.x - 1.0;
637 cd1.y = cdegree.y;
638
639 temp = ComplexPower(old, cd1);
640 FPUcplxmul(&temp, &old, &new);
641
642 tmp.x = new.x - croot.x;
643 tmp.y = new.y - croot.y;
644 if((sqr(tmp.x) + sqr(tmp.y)) < threshold) {
645 if(fabs(old.y) < .01)
646 old.y = 0.0;
647 FPUcplxlog(&old, &temp);
648 FPUcplxmul(&temp, &cdegree, &tmp);
649 mod = tmp.y/TwoPi;
650 coloriter = (long)mod;
651 if(fabs(mod - coloriter) > 0.5) {
652 if(mod < 0.0)
653 coloriter--;
654 else
655 coloriter++;
656 }
657 coloriter += 2;
658 if(coloriter < 0)
659 coloriter += 128;
660 return(1);
661 }
662
663 FPUcplxmul(&new, &cd1, &tmp);
664 tmp.x += croot.x;
665 tmp.y += croot.y;
666
667 FPUcplxmul(&temp, &cdegree, &cd1);
668 FPUcplxdiv(&tmp, &cd1, &old);
669 if(overflow)
670 {
671 return(1);
672 }
673 new = old;
674 return(0);
675 }
676
677 /*
678 * Generate a gaussian distributed number.
679 * The right half of the distribution is folded onto the lower half.
680 * That is, the curve slopes up to the peak and then drops to 0.
681 * The larger slope is, the smaller the standard deviation.
682 * The values vary from 0+offset to range+offset, with the peak
683 * at range+offset.
684 * To make this more complicated, you only have a
685 * 1 in Distribution*(1-Probability/Range*con)+1 chance of getting a
686 * Gaussian; otherwise you just get offset.
687 */
688 int GausianNumber(int Probability, int Range) {
689 int n, r;
690 long Accum = 0, p;
691
692 p = divide((long)Probability << 16, (long)Range << 16, 16);
693 p = multiply(p, con, 16);
694 p = multiply((long)Distribution << 16, p, 16);
695 if(!(rand15() % (Distribution - (int)(p >> 16) + 1))) {
696 for(n = 0; n < Slope; n++)
697 Accum += rand15();
698 Accum /= Slope;
699 r = (int)(multiply((long)Range << 15, Accum, 15) >> 14);
700 r = r - Range;
701 if(r < 0)
702 r = -r;
703 return(Range - r + Offset);
704 }
705 return(Offset);
706 }
707
708 #endif
709