Flyttalslib till 16F887 ?

[Marta]

Allt flöt bra hit, men nu är det tvärnit.

Behöver hjälp med antingen saknade siffror som är bortklippta från tabellen jag behöver, eller att få till en alternativ algoritm där det finns tabeller. Den använda har en enda serie. Den alternativa algoritmen har två serier som sedan divideras med varandra och lite till. Det som därutöver görs är black magic för mig.
Det är värt några Birgit att få det löst.

[persika]

Skulle ascii-aritmetik vara något?
Då skulle man kunna få valfritt antal värdesiffror.

[Marta]

Ser inte hur det skulle hjälpa, problemet är på ett annat plan. Det är väl banker o.dyl. som håller på med sådant för at binärt saknar exakt representation av decimaltal. Här är binärt helt ok.

För log(x) vill jag helt enkelt byta till en annan tabell och lägga till en iteration. Den tabellen är definitivt trasig, raderna är brutalt trunkerade.

För e**x utlovas 12 sirror, den viker efter 7. OCR-fel är en möjlighet.

Hade hoppats boken skulle finnas så någon här har tillgång till den så jag kan få dessa tabeller, De definitivt trasiga 2664 och 2665, samt den möjligen trasiga 1045.

[hummel]

Har du kollat med biblioteket (som någon länkade till tidigare) och helt enkelt låna boken?

[Marta]

Jag har lagt en beställning av kopia av två sidor från KB. Förmodligen blir det OCR. Tillsammans med den kopa jag har är nog sannolikheten stor att kunna få fram korrekta tabeller.

Det är av praktiska skäl omöjligt att beställa boklån o.dyl. Är ju invalid som de flesta här redan vet. Jag får noll hjälp med hobbyrelaterade ärenden.

[Icecap]

Jag satt just o hittade formeln för att räkna ut avståndet mellan två geografiska punkter.

Mycket av precisionen hänger på flyttal och då insåg jag att i vissa kompilers kan det finnas float, double samt long double.

Just den siste blev "officiell" med C99 och KAN vara "mer" än double - och då kan den kanske vara intressant för TS.

[pi314]

...avståndet mellan två geografiska punkter.

Jag gissar att du refrerar till Haversines formel?

https://en.wikipedia.org/wiki/Haversine_formula

Kod: Markera allt

Haversine Formula
R = earth’s radius (mean radius = 6,371km)
Δlat = lat2− lat1
Δlong = long2− long1
a = sin²(Δlat/2) + cos(lat1).•cos(lat2).sin²(Δlong/2)
c = 2•atan2(√a, √(1−a))
d = R•c

Jag har använt den i några olika programmerngsspråk.

Det första exempel jag hittade nu var Octave/Matlab.

[VISA]

Kod: Markera allt

function rad = radians(degree) 
% degrees to radians
    rad = degree .* pi / 180;
end; 

function [a,c,dlat,dlon]=haversine(lat1,lon1,lat2,lon2)
% HAVERSINE_FORMULA.AWK - converted from AWK 
    dlat = radians(lat2-lat1);
    dlon = radians(lon2-lon1);
    lat1 = radians(lat1);
    lat2 = radians(lat2);
    a = (sin(dlat./2)).^2 + cos(lat1) .* cos(lat2) .* (sin(dlon./2)).^2;
    c = 2 .* asin(sqrt(a));
    arrayfun(@(x) printf("distance: %.4f km\n",6372.8 * x), c);
end;

[a,c,dlat,dlon] = haversine(36.12,-86.67,33.94,-118.40); % BNA to LAX

Här finns rätt många språk, inklusive C, C# och C++.

https://rosettacode.org/wiki/Haversine_formula


/π

[Icecap]

Nja... mer Vincenty-formeln.

Jag har dock båda liggande som C-kod.

[rvl]

Vincenty om man insett att jorden inte är klotformat och vill ha litet noggranhet, men sen om inte 0,1 mm noggrannhet räcker till, så finns det FORTRAN-bibliotek som är ännu ett strå vassare.

[Marta]

Nu fungerar den här koden hyfsat, funderar ändå på at förbättra ytterligare. Mantissan är nu 47-bit plus ettan. Lagringsformatet är mantissa+exponent+tecken, där det senare lämnar 7 bits tomma.

Med ytterligare 1 byte mantissa blir denna 55-bit plus ettan. Med ett lagringsformat (ja, det är udda) där teckenbit läggs i mantisans toppbyte som bit 0 blir det fortfarande 8 byte för att lagra ett tal. Det blir 2.5 siffra mer, en bit mer än IEEE double.

I en MCU finns näppeligen något intresse av att utbyta tal i IEEE-format med omgivningen och +/- 10**38 räcker långt. En exponent som håller sig i en byte har fördelar. Ofta förekommer det i beräkningar att något är *2, då är det bara att inkrementera exponenten. Mellanlagring av mellanresultat kan undvikas genom detta.

IEEE kräver ycket skyfflande och pusslande. Nämnda format kräver endast två skift för att spara, set carry och tvåskift för att hämta. Väldigt lite overhead.

Vad säger Ni om detta, är det en förbättring, neutralt eller rent av försämring?

[rvl]

Litet smått relaterat:
https://www.righto.com/2025/01/pentium- ... t-ROM.html
Ken har "skrapat" under ytan på Pentium och hittat 304 flyttalskonstanter a 86 bit (17+68).
Även Remez-algoritmen nämns för att ta fram de optimala konstanterna.

[Marta]

Det där med konstanterna i Pentium är något som svävar högt över mitt vetande.

Har tyvärr aldrig fått den utbildning jag velat ha. Därav följande fråga. Antar den är 100% ROTFL för Er som haft styrkan att palla allt meningslöst blaj som krävs för att få en vettig grundutbildning.

Programmet är mitt tafatta försök att implementera formeln i rutan under. Har testat många varianter, men lyckan sviker. Detta är första varianten. Skall kunna beräkna 2**x för 0<x<0.5.
Skulle vara mycket tacksam om någon som har kunskapen kunde rätta programmet så det fungerar. Det är sista funktionen i aktuellt fp-lib för att få det komplett. Övrigt är nu klart och precisionen ökad till 55 bits. Lyckades komplettera de trasiga tabellerna med de snuttar av sidor Google visar från en ebook.

Kod: Markera allt

#include <stdio.h>

void main(){
  int i;
  double x, p, q, r;

  double M_EXP_P[] ={
            1513.906799054338915944,
            20.202065651286927227,
            0.023093347753750234,
          };

  double M_EXP_Q[] ={
            4368.211662727558498354,
            233.184211427481623793,
            1.000000000000000000,
          };

  x=0.32;  //försöker få fram 2**0.32 som testvärde

//        jsr _TABLE_LD          *03FB   2013               M         CALL    _TABLE_LD
//                               *                      00072         FP_ADD                  ; P += Constant
//        jsr _FP_ADD            *03FC   2292               M         CALL    _FP_ADD
//                               *                      00073         TRarg_LOAD      FParg1          ; P *= FParg1
//        lod a,#FParg1          *03FD   3044               M         movlw   FParg1
//        jsr _TRarg_LOAD        *03FE   21E3               M         CALL    _TRarg_LOAD
//                               *                      00074         FP_MUL
//        jsr _FP_MUL            *03FF   22CB               M         CALL    _FP_MUL

  p=0.0;
  for (i=0; i<3; i++){
    p+=M_EXP_P[i];
    p*=x*x;
  };
  p*=x;

  q=0.0;
  for (i=0; i<3; i++){
    q+=M_EXP_Q[i];
    q*=x*x;
  };

  r=(q+p)/(q-p);

  printf("%5.5f,    1.24833 \n", r);

};//main

[pi314]

Du har kastat om ordningen på koefficienterna och beräknat polynomen fel.

Testa detta.

Kod: Markera allt

#include <math.h>
#include <stdio.h>

// For test:
// b**z = e**(z x ln(b))

int main(){
  int i;
  long double x, p, q, r;

  long double M_EXP_P[] ={
            0.023093347753750234,
            20.202065651286927227,
            1513.906799054338915944,
          };

  long double M_EXP_Q[] ={
            1.000000000000000000,
            233.184211427481623793,
            4368.211662727558498354,
          };

  x=0.32;  //försöker få fram 2**0.32 som testvärde

//     jsr _TABLE_LD       *03FB   2013            M      CALL    _TABLE_LD
//                         *                   00072      FP_ADD                 ; P += Constant
//     jsr _FP_ADD         *03FC   2292            M      CALL    _FP_ADD
//                         *                   00073      TRarg_LOAD      FParg1 ; P *= FParg1
//     lod a,#FParg1       *03FD   3044            M      movlw   FParg1
//     jsr _TRarg_LOAD     *03FE   21E3            M      CALL    _TRarg_LOAD
//                         *                   00074      FP_MUL
//     jsr _FP_MUL         *03FF   22CB            M      CALL    _FP_MUL

  p=M_EXP_P[0];
  for (i=1; i<3; i++){
    p*=x*x;
    p+=M_EXP_P[i];
  };
  p=x*p;

  q=M_EXP_Q[0];
  for (i=1; i<3; i++){
    q*=x*x;
    q+=M_EXP_Q[i];
  };

  r=(q+p)/(q-p);

  printf("%10.10Lf,    1.24833055489 \n", r);

  printf("%10.10f \n", exp(x*log(2.0)));


return 0;

};//main

/π

[Marta]

Stort Tack. Ser det nu, skäms.

Orginalprogrammet är skrivet så att det initierar med noll och sedan avslutar med en extra operation.

Nu återstår att få till skalningen. Den förra klarade 0<x<1, den här är 0<x<0.5. Antar att det görs genom att decrementera exponenten ett steg och sedan kompensera slutresultatet för detta.

[rvl]

Jag förväntade mig typ:

Kod: Markera allt

  pp=p=qq=q=0.0;
  xx = x*x;
  for (i=0; i<3; i++){
    p=M_EXP_P[i];
    p*=xx;
    pp+=p;
    q=M_EXP_Q[i];
    q*=xx;
    qq+=q;
    xx*=x*x;
  };
  pp*=x;

  r=(qq+pp)/(qq-pp);

edit: typo fix