Keskustelu: Koodit: C: Jakolasku ilman jakolaskua

Olipahan tylsää joten ajattelin heittää tänne tälläisen.

Eli miten lasketaan jakolasku kun ei ole käytettävissä jakolasku operaatiota. No yksi tapa on laskea se iteraatiolla. Yksi hyvä iteraatiofunktio on Newtonin iteraatiofunktio.

x<n+1>=x<n>-(f(x<n>)/(f'(x<n>))

x<n+1> = lukujonon n+1:s luku ja vastaavasti x<n> n:s luku
f(x<n>) = iteroitavan funktion arvo luvulla x<n>
f'(x<n>) = iteroitavan funktion derivaatan arvo luvulla x<n>

Jakolaskua laskettaessa lasketaan jakajan(nimittäjän) käänteisarvo joka kerrotaan jaettavalla(osoittajalla). Joten funktio f on f(x)=a-(1/x). Tästä ja iteraatio menetelmästä seuraa että itse iteraatiolaskua on muotoa:

x<n+1>=x<n>*(2-a*x<n>)

jossa a on jakajan(nimittäjän) arvo.

Kun iteraation virhe abs(x<n+1>-x<n>) on riittävän pieni lopetetaan iteraatio ja kerrotaan sen tulos jaettavalla(osoittajalla).

Alettaessa laskemaan täytyy tehdä alkuarvaus luvulle x0 jotta iteraatiot voidaan aloittaa. Jotta iteraatio lähestyisi arvoa 1/a täytyy x0 valita väliltä 0-(2/a). Jos alkuarvaukseksi laitetaan hyvin pieni luku kestää iteraatio turhan kauan toisaalta jos se menee yli (2/a):n ei iteraatio suppene koskaan. Joten miten valitaan hyvä alkuarvaus?

Helpoiten luvusta a saadaan 'hyvä' alkuarvaus ottamalla luvun a eksponentti jota korotetaan yhdellä ja vaihdetaan sen etumerkki. Esim 3.45e6 olisi 1e-7 joka on alle (2/a)=5.7971e-7.

Tämän voi suorittaa helposti kun huomioidaan että liukuluvut ovat tallennettu muodossa <SGN><EXP><SIGNIFICAND> joista SGN=1 bit EXP=8 bit ja SIGNIFICAND=23 bit kun käytössä on 32bit(float)-liukuluvut. Koska eksponentti on biasoitu eli se on aina positiivinen helpottuu lasku edelleen. Nythän tarvitsee ainoastaan kasvattaa eksponenttia yhdellä ja ottaa sen negaatio.

Yleensä iteraatiot kestävät n. 6 kierrosta, mutta se vähän riippuu siitä kuinka kauas alkuarvaus menee laskettavasta. Tietenkään tämä ei ole nopeuden puolesta verrattavissa normaaliin jakolaskuun, mutta näyttää miten voidaan laskea ko. lasku kun ei oikeaa jakolaskua ole saatavilla.

#define DIV_ERR 1e-9

inline float divide(float a,float b)
{
	int *p;
	float old,ne,t;

	old=b;
	p=(int*)&old;
	//säilytetään etumerkki ja eksponentti
	*p&=0xff800000;
	//korotetaan eksponenttia yhdellä
	*p+=0x00800000;
	//käännetään eksponentti ts. (255-eksponentti)
	*p^=0x7f800000;

	for(;;){
		//lasketaan iteraatio
		ne=old*(2-old*b);
		//iteraation virhe
		t=old-ne;
		//fabs antaa välillä pienille luvuille arvoksi 1.#INF00
		//joten täytyy tehdä tämä itse
		p=(int *)&t;
		*p&=0x7fffffff;
		//virhe 'hyväksyttävä' - lopetetaan iteraatiot
		if(t<DIV_ERR)break;
		old=ne;
		}
	//palauta a/b
	return(a*ne);
}

Mikäs siinä :D

Kiva koodi. :P
Eipä varmaan moni osaa...

...mutta mitä hyötyä kun c++:ssa on kerran jakolasku-ominaisuus? =)

Ei tätä ole tarkoitettu käytettävän jakolaskun sijasta vaan tämä on yksinkertaisesti esimerkki siitä, kuinka voi kiertää jakolaskua. Kiinnostava matemaattisesti.

Kumpi onpi nopeampi?

Aika kiinnostava systeemi. Ei tosiaan varmaan mitään käytännön hyötyä kuitenkaan :)

Itse asiassa tuosta on hyötyä jos tarvitsee tehdä esim. jollekin mikrokontrollerille jokin ohjelma joka tarvitsee kyseisen operaation. Yleensä kuitenkin valmistajat tarjoavat valmista koodia tuon operaation toimittamiseen esim. MicroChip:n PicMicro piireille löytyy IEEE standardin mukaisille liukuluvuille kaikki perusoperaatiot suorittavia lähdekoodi listauksia. Itse asiassa MicroChip:n listauksissa oleva jakolasku käyttää tuota samaa algoritmia paitsi alkuarvaus otetaan taulukosta.

Tuosta nopeudesta sen verran että tuo ko. listaus ei ole hirveän nopea, mutta toinen versio jossa on parempi alkuarvaus ja joka on kokonaan assemblyllä kirjoitettu suoriutuu n. 5 kertaa fdiv:iä hitaammin minun koneella.