2
votes

Comment pourriez-vous diviser un IEEE754 64 bits par 1 000 sur une plate-forme ne prenant en charge que le flottant 32 bits?

J'ai un compteur électrique connecté à un DCS (système de contrôle distribué) par PROFIBUS. Le compteur (Siemens Sentron PAC3200) fournit son compte comme un IEEE 754 double en Wh (wattheures). De plus, le compteur déborde à 1,0e12 Wh ou 1000 GWh. (Scène en coupe: plusieurs années plus tôt, les laboratoires de développement de Siemens. "Voyons, comment transférer une valeur entière non signée de 40 bits? Utilisons le double!")

Mon objectif est de consigner le décompte de manière cohérente avec une précision en kWh. p>

Cependant, le DCS ne prend en charge que les flotteurs à simple précision. Donc, si je prenais la route directe, c'est-à-dire que je pressais les données dans un flotteur, alors à environ sept chiffres décimaux, des erreurs apparaissaient dans la lecture du kWh, c'est-à-dire au plus tard à partir d'environ 100000000 Wh ou 100 MWh. Le décompte actuel est déjà de 600 MWh, donc ce n'est pas possible.

Donc pour l'instant, je mets la mantisse dans un double entier non signé (UDINT, 32 bits sur cette plateforme) et j'effectue la conversion selon IEEE 754, qui donne la valeur correcte en Wh. Cela entraîne cependant un débordement à 2 ^ 32 Wh soit environ 4,3 GWh, qui nous durera à peine dix ans.

Comme je n'ai besoin que de kWh de précision, j'ai eu l'idée de diviser par 1000 au début de la conversion. Cela placerait le débordement variable à 4 300 GWh, et le compteur interne du compteur déborde déjà à 1 000 GWh. Problème résolu, en théorie.

Comme IEEE 754 est un format binaire à virgule flottante, je ne peux que diviser facilement par 1024 (décalage à droite 10 fois), ce qui introduit une erreur substantielle. La multiplication par un facteur de correction de 1,024 par la suite ne se produirait qu'en simple précision sur cette plate-forme, annulant l'effort précédent.

Une autre option serait de générer un UDINT "haut" et "bas" en Wh à partir du conversion, alors je pourrais au moins en théorie calculer à nouveau le kWh, mais cela semble gênant (et compliqué).

J'ai le sentiment subtil que j'ai peut-être oublié quelque chose (une seule personne parler); Je suis ouvert à toute autre idée de la manière dont je pourrais obtenir le 1 / 1000e de la double valeur transférée.

Merci et meilleures salutations

Björn

PS: Pour votre plaisir visuel, c'est la solution basée sur la réponse de @ EricPostpischil - adaptée aux spécificités de la plateforme et des tâches. Le langage utilisé est SCL (langage de contrôle structuré) selon EN 61131-3, qui est une sorte de dialecte Pascal.

  high byte     low byte  decimal value
=======================================
16#41c558c3, 16#2d3f331e,       716_277
16#41EFFFFF, 16#5E000000,     4_294_966
16#41EFFFFF, 16#DB000000,     4_294_967
16#41F00000, 16#2C000000,     4_294_968
16#426D1A94, 16#A1830000,   999_999_999
16#426D1A94, 16#A2000000, 1_000_000_000
16#426D1A94, 16#A27D0000, 1_000_000_001
16#428F3FFF, 16#FFC18000, 4_294_967_294
16#428F3FFF, 16#FFE0C000, 4_294_967_295
16#428F4000, 16#00000000, 4_294_967_296

Les données de test que j'ai utilisées:

FUNCTION_BLOCK PAC3200KON_P

VAR_INPUT
    INH : DWORD;
    INL : DWORD;
END_VAR

VAR_OUTPUT
    OUT : UDINT;
    SGN : BOOL;
END_VAR

VAR
    significand:              UDINT;
    exponent, i, shift:       INT;
    sign:                     BOOL;
    d0, d1, y0, y1, r1, temp: DWORD;
END_VAR
(*
    Convert the energy count delivered by Siemens Sentron PAC3200
    (IEEE 754 binary64 format, a.k.a. double) into an UDINT.

    Peculiarities:
    - This hardware platform only supports binary32 (a.k.a. float).

    - The Sentron's internal counter overflows at 1.0e12 Wh (1000 GWh).

    - kWh resolution suffices.

    - If you converted the double directly to UDINT and divided by 1000
      afterwards, the range would be reduced to (2^32-1)/1000 GWh or about
      4.295 GWh.

    - This is why this function first divides the significand by 1000
      and then proceeds with conversion to UDINT. This expands the
      range to (2^32-1) GWh or about 4295 GWh, which isn't reachable in
      practice since the device's internal counter overflows before.

    Background:

    IEEE 754 binary64 bit assignment:

               High-Byte                         Low-Byte
    66665555555555444444444433333333 3322222222221111111111
    32109876543210987654321098765432 10987654321098765432109876543210
    GEEEEEEEEEEESSSSSSSSSSSSSSSSSSSS SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS

    G: sign (1: negative)
    E: exponent (biased; subtract 1023) (11 bits)
    S: significand (52 bits)
*)

(*
    significand: Bits 19...0 of high byte und complete low byte

    The significand is initially divided by 1000 using integer division. The
    bits are divided into two parts:

    - d1 contains the 31 most significant bits (plus leading 1)
    - d0 contains the next less significant bits

    In total, we use 48 bits of the original significand.
*)

(* d1: insert significand bits from high byte *)
d1 := INH AND     2#0000_0000_0000_1111_1111_1111_1111_1111;
(* result:        2#0000_0000_0000_HHHH_HHHH_HHHH_HHHH_HHHH *)

(* add the 1 before the binary point *)
d1 := d1 OR       2#0000_0000_0001_0000_0000_0000_0000_0000;
(* result:        2#0000_0000_0001_HHHH_HHHH_HHHH_HHHH_HHHH *)

(* "flush left" shift 11 places *)
d1 := d1 * 2048;
(* result:        2#1HHH_HHHH_HHHH_HHHH_HHHH_H000_0000_0000 *)

(* Insert another 11 bits from low byte (msb ones) *)
d1 := d1 OR (INL / 2097152);
(* result:        2#1HHH_HHHH_HHHH_HHHH_HHHH_HLLL_LLLL_LLLL *)

(* Base-65536 division. Integer divide by 1000 and save remainder *)
y1 := d1 / 1000;
r1 := TO_DW(TO_UD(d1) MOD 1000);

(*
   The significand now has leading zeroes. Shift left to make space
   at the other end.
*)
FOR shift := 1 TO 31 BY 1 DO
    y1 := y1 * 2;
    IF (y1 AND 2#1000_0000_0000_0000_0000_0000_0000_0000) <> 0 THEN
        EXIT;
    END_IF;
END_FOR;

(*
   d0: insert next 16 bits from the low byte
   (right shift five times and zero out the leading places)
*)
(* bits:             2#xxxx_xxxx_xxxL_LLLL_LLLL_LLLL_LLLx_xxxx *)
d0 := (INL / 32) AND 2#0000_0000_0000_0000_1111_1111_1111_1111;
(* result:           2#0000_0000_0000_0000_LLLL_LLLL_LLLL_LLLL *)

(* Now divide by 1000, factoring in remainder from before *)
y0 := ((r1 * 65536) OR d0) / 1000;

(*
   y1 and y0 contain results from division by 1000. We'll now build a 32 bit
   significand from these.

   y1 = 2#1HHH_HHHH_HHHH_HHHH_HHHH_HHxx_xxxx_xxxx
   y0 = 2#0000_0000_0000_0000_LLLL_LLLL_LLLL_LLLL

   y1 has an uncertain number of zeroes at its end, resulting from the above
   left shifting (number of steps inside variable "shift"). Fill those with the
   most significant bits from y0.

   y0 has 16 valid bits (0..15). Shift right so that the "highest place zero"
   in y1 corresponds with the MSB from y0. (shift by 16-shift)

   y1 = 2#1HHH_HHHH_HHHH_HHHH_HHHH_HHxx_xxxx_xxxx (ex.: shift=10)
   y0 = 2#0000_0000_0000_0000_0000_00LL_LLLL_LLLL
                              ------>^
*)

FOR i := 1 TO 16 - shift BY 1 DO
    y0 := y0 / 2;
END_FOR;

significand := TO_UD(y1 OR y0);
(* Result: 32-bit significand *)

(*
    Exponent: bits (62-32)...(59-32) or bits 30...20 of high byte, respectively

    Coded with bias of 1023 (needs to be subtracted).

    Special cases as per standard:
    - 16#000: signed zero or underflow (map to zero)
    - 16#7FF: inifinite or NaN (map to overflow)
*)
temp := 2#0111_1111_1111_0000_0000_0000_0000_0000 AND INH;
temp := temp / 1048576 ; (* right shift 20 places (2^20) *)
exponent := TO_IN(TO_DI(temp));
exponent := exponent - 1023; (* remove bias *)

(*
   Above, we already left shifted "shift" times, which needs to be taken into
   account here by shifting less.
*)
exponent := exponent - shift;

(*
    The significand will be output as UDINT, but was initially a binary64 with
    binary point behind the leading 1, after which the coded exponent must be
    "executed".

    temp = 2#1.HHH_HHHH_HHHH_HHHH_HHHH_HLLL_LLLL_LLLL

    As UDINT, this already corresponds to a 31-fold left shift.

    Exponent cases as per IEEE 754:

    - exponent < 0:            result < 1
    - exponent = 0:       1 <= result < 2
    - exponent = x > 0: 2^x <= result < 2^(x+1)

    The UDINT output (32 bit) allows us to represent exponents right up to 31.
    Everything above is mapped to UDINT's maximum value.

    Now determine, after the de facto 31-fold left shift, what shifts remain
    "to do".
*)

IF exponent < 0 THEN
    (* underflow: < 2^0 *)
    significand := 0;
ELSIF exponent > 31 THEN
    (* overflow: > 2^32 - 1 *)
    significand := 4294967295;
ELSE
    (*
        result is significand * 2^exponent or here, as mentioned above,
        significand * 2^(31-exponent).

        The loop index i is the "shift target" after loop execution, which is
        why it starts at 31-1.

        Example: exponent = 27, but de facto we've already got a shift of 31.
        So we'll shift back four times to place the binary point at the right
        position (30, 29, 28, 27):

        before: temp = 2#1HHH_HHHH_HHHH_HHHH_HHHH_HLLL_LLLL_LLLL.

        after:  temp = 2#1HHH_HHHH_HHHH_HHHH_HHHH_HLLL_LLLL.LLLL
                                                           ^<---|
    *)
    FOR i := 30 TO exponent BY -1 DO
        significand := significand / 2;
    END_FOR;
END_IF;

(*
    sign: bit 63 of high byte
*)
sign := (2#1000_0000_0000_0000_0000_0000_0000_0000 AND INH) <> 0;

OUT := significand;
SGN := sign;

END_FUNCTION_BLOCK

BTW, les littéraux entiers de la forme b # 1234 en SCL signifient fondamentalement "le nombre 1234 en base b". Les soulignements sont ignorés (ce sont des séparateurs de chiffres pour une meilleure lisibilité, comme par exemple Python les a).


6 commentaires

Pour commencer, vous pouvez diviser le compteur du compteur par 8, puis le convertir en un entier 32 bits non signé, puis le diviser par 125. Cela augmente votre délai de dépassement à quatre-vingts ans. Ce qui est suffisamment long pour que vous n'ayez pas à vous en occuper personnellement, mais suffisamment court pour que nous ne puissions pas être sûrs que les robots et les IA qui existeront alors résoudront le problème pour nous.


Quelques notes: Il n'y a pas vraiment de problème de virgule flottante ici, car, à part extraire un entier de 40 bits du double , vous n'avez pas d'opérations double et ne pouvez pas utilisez de manière productive les opérations float , donc vous n'avez vraiment que des opérations sur des entiers. • Le terme préféré pour la partie fractionnaire d'un nombre à virgule flottante est «significande». «Mantissa» est un ancien terme pour la partie fractionnaire d'un logarithme. Les signifiants sont linéaires; les mantisses sont logarithmiques. • En ce qui concerne le choix de Siemens, sur un système avec double 64 bits et pas d’entier 64 bits, c’est un moyen raisonnable de gérer les entiers 40 bits.


Je ne pense pas que produire deux UDINT à partir de la conversion soit horrible. Vous pouvez le traiter comme trois chiffres radix-65536 et le diviser par 1000 en utilisant la division longue de l'école primaire: prenez les deux premiers «chiffres» comme un entier de 32 bits, divisez par 1000, multipliez le reste par 1000 et ajoutez le troisième, puis divisez par 1000.


Salut @Eric, merci beaucoup pour vos commentaires détaillés et le code source! La solution des 80 ans pourrait fonctionner - je doute qu'un compteur de type Ferrari vive aussi longtemps, je serais surpris si les électroniques atteignent 40. Mais je vais essayer votre source, cela me semble plus complet Solution. Cela prendra cependant quelques jours à cause d'autres projets.


Hein? FLOAT peut gérer plus de 1e38. Quelle unité l'appareil vous donne-t-il, les microwatts?


Comme mentionné, @RickJames, Watt-heures (unité d'énergie, pas de puissance). Il mesure l'énergie «wattful» dans un comptoir. Ce compteur augmente avec le temps, mais les erreurs d'arrondi doivent rester inférieures à 1 kWh même si le compteur atteint son maximum de 1e12 Wh.


3 Réponses :


3
votes

                                

1 commentaires

Salut @Eric, votre solution "à part entière" LongDivision fonctionne comme un charme. J'ai dû l'ajuster cependant, après avoir divisé l'octet haut, je fais toujours y1 justifié à gauche en décalant et remplis les bits inférieurs en utilisant y0 pour obtenir le significande final. Les décalages déjà exécutés doivent alors être soustraits des décalages à faire restants codés dans l'exposant pour obtenir la valeur finale. Je publierai la source dans les prochains jours. Merci encore! :)



1
votes

J'aborderais le problème d'une manière légèrement différente. Puisque l'OP n'a mentionné aucun langage de programmation utilisé, j'écris ici un pseudocode . Je suppose que le nombre à virgule flottante binary64 est passé au compilateur sous la forme d'une séquence de 8 octets. Je suppose que l'OP s'occupera de l'endianité là où c'est nécessaire.

1. Divisez le binaire64 en trois nombres à virgule flottante binaire32:

A nombre binaire64 à virgule flottante est représenté par un seul bit de signe, 11 bits d'exposant et 52 bits représentant le significatif:

 représentation binaire 64 bits

et est calculé comme suit:

# rescale double-precision d by a single-precision a
function res = Rescale(d,a)
   # first term
   f = (2^23 + Extract(d,1,23)) * 2^(Exponent(d) - 23)
   [p,s] = TwoProduct(f,a)
   # second term
   f = Extract(d,24,46) * 2^(Exponent(d) - 46)
   [h,r] = TwoProduct(f,a)
   [p,q] = TwoSum(p,h)
   s = s + (q + r)       # the error term
   # third term
   f = Extract(d,47,52) * 2^(Exponent(d) - 52)
   [h,r] = TwoProduct(f,a)
   [p,q] = TwoSum(p,h)
   s = s + (q + r)       # the error term
   # the final result
   res = p + s

A le nombre à virgule flottante binary32 est représenté par un bit de signe unique, 8 bits d'exposant et 32 ​​bits représentant le significatif:

 représentation binaire 32 bits

et est calculé comme:

# error-free transformation of the sum of two floating-point numbers
function [x,y] = TwoSum(a,b)
   x = a + b
   z = x - a
   y = ((a - (x - z)) + (b - z))
# Error-free split of a lfoating point number in two parts
function [x,y] Split(a)
   c = (2^12 - 1) * a
   x = (c - (c - a))
   y = a - x
# error-free transformation of the product of two floating-point numbers
function [x,y] = TwoProduct(a,b)
   x = a * b
   [a1,a2] = Split(a); [b1,b2] = Split(b)
   y = (a2*b2 - (((x - a1*b1) - a2*b1) - a1*b2))

L'idée est maintenant de créer t trois nombres à virgule flottante binaire32 f {1,2,3} tels que, lors de l'utilisation de l'arithmétique réelle (pas d'approximations à virgule flottante), le nombre à virgule flottante binaire64 d est donné par:

f1 * 1E-3 + f2 * 1E-3 + f3 * 1E-3

Supposons que la fonction EXTRACT (d, n, m) renvoie un entier extrait des bits n jusqu'à m depuis la représentation binaire 64 bits d:

f1 = (2^23 + Extract(d,1,23)) * 2^(Exponent(d) - 23)
f2 = Extract(d,24,46) * 2^(Exponent(d) - 46)
f3 = Extract(d,47,52) * 2^(Exponent(d) - 52)

et la fonction Exposant ( d) renvoie la valeur e-1023 de la représentation binaire à 64 bits d.

Alors nous savons que

XXX

2. Divisez les valeurs par 1000:

C'est malheureusement plus facile à dire qu'à faire. Il est bien connu que le calcul avec une précision finie implique des erreurs d'arrondi, conduisant à des résultats inexacts pour un calcul. C'est exactement ce que nous essayons d'éviter ici. Si nous calculions simplement

function val Extract(d,n,m)
   val = Sum(b52−i 2n−i;i = m → n )

nous introduirions des erreurs d'arrondi.

Supposons que a et b code > sont 2 nombres à virgule flottante, la fonction fl (x) renvoie le nombre à virgule flottante de la valeur réelle x et a OP b représente le nombre réel complet en arithmétique réelle des opérations de base + , - et * . Avec cela, nous savons que a OP b! = Fl (a OP b) en tant que nombre réel ne peut pas toujours être entièrement représenté par un nombre à virgule flottante. Cependant, on peut montrer que a OP b = fl (a OP b) + y avec y un nombre à virgule flottante. Ce y est l'erreur que nous manquerions dans le calcul ci-dessus en calculant simplement f1 * fl (1E-3) .

Donc, pour calculer d * fl (1E-3) avec précision, nous aurons besoin de garder une trace des termes d'erreur. Pour cela, nous utiliserons quelques transformations sans erreur qui sont examinées dans l'article Somme précise, produit scalaire et évaluation polynomiale en arithmétique flottante complexe :

d = f1 + f2 + f3

3. La fonction complète:

Donc, si nous voulons redimensionner le nombre binaire64 avec une représentation binaire d en utilisant l'arithmétique à virgule flottante binary32, nous devrions utiliser la fonction:

(−1)b31 (1 + Sum(b23−i 2−i;i = 1 → 23 )) × 2e−127

Cela aura gardé une trace de toutes les erreurs numériques dans les mathématiques à virgule flottante et compensé le résultat en conséquence. En conséquence, la valeur res renvoyée par Rescale représentera la valeur simple précision la plus précise de d / 1000 .


7 commentaires

Merci pour votre message, comme mentionné, il me faudra un peu de temps pour l'examiner. En outre, le langage utilisé est appelé "Texte structuré" (ST) et est défini dans la CEI 61131-3 comme l'un des multiples langages standard pour les systèmes d'automatisation industrielle.


@ chr0n0ss Je suppose que vous passez une chaîne de caractères à cette routine (puisque vous ne pouvez pas passer un double). Faites attention à la façon dont vous gérez l'endianness.


Merci! J'utilise plusieurs primitives DWORD (double mot), 32 bits chacune (dans CEI 61131, un mot est défini comme 16 bits). Heureusement, je suis capable de "jouer" avec la configuration, comme je l'ai en effet remarqué lors des premiers essais que l'ordre des octets est inversé entre les deux systèmes.


La déclaration «Cependant, on peut montrer que a OP b = fl (a OP b) + y avec y un nombre à virgule flottante» est fausse. Si, par exemple, a vaut 1 et b vaut 1e − 3, alors a OP b se répète indéfiniment en binaire et ne peut pas être la somme de deux chiffres qui se terminent en binaire. Je n'ai pas cherché à voir si cette réponse pouvait être corrigée.


@EricPostpischil Si a et b sont deux nombres à virgule flottante (pas des réels), alors on peut montrer que a OP b = fl (a OP b) + y_. La valeur 1E-3 n'est pas un nombre à virgule flottante mais un nombre réel qui ne fait pas partie de l'ensemble des nombres à virgule flottante.


Mais le but est de calculer a • 1e − 3 (ou floor (a • 1e − 3)), comme la réponse a tenté avant l'édition il y a quelques minutes, et non un • fl (1e − 3).


@EricPostpischil L'OP déclare Je suis ouvert à toute autre idée sur la façon dont je pourrais obtenir le 1 / 1000e de la valeur double transférée . Ce qui précède fournit un mécanisme sans erreur pour utiliser la valeur double précision complète et la diviser par 1000 (virgule flottante)



0
votes

1e12 Wh / 1kWh = 1e9.

  • Un INT de 4 octets, 32 bits (signé ou non signé) vous donne un peu plus de 9 chiffres significatifs. Mais vous devez vous rappeler que c'est en unités de KWh, pas en Wh. Et chaque fois que vous y ajoutez, vous obtenez potentiellement une autre erreur d'arrondi.

  • FLOAT n'a que ~ 7 chiffres de résolution; vous avez besoin de 9. Le capteur peut envoyer un FLOAT , mais il n'est pas acceptable de accumuler dans FLOAT .

  • DOUBLE est de 8 octets, ~ 16 chiffres significatifs. Cela vous permettrait de stocker Wh.

Un compromis consiste à accumuler en utilisant DOUBLE , mais diviser par 1000 et stocker dans un INT de 4 octets.

Y a-t-il un problème de stockage en DOUBLE ? Outre la prise d'espace supplémentaire, il résout essentiellement tous les problèmes - plus qu'une résolution adéquate et une protection contre les erreurs d'arrondi; capacité à stocker l'unité «naturelle» de Wh; etc. J'utiliserais DOUBLE si possible.

(Un entier 32 bits non signé serait le deuxième choix.)

Je ne considérerais même pas un entier de 40 bits (ou deux entiers de 32 bits) car il est difficile de travailler avec et cause probablement des difficultés de portage.

Sortir des sentiers battus ...

Stockez les sous-totaux pour chaque année. Ensuite, lorsque vous avez besoin du total général, additionnez les sous-totaux. (C'est ce que je pourrais faire de toute façon pour une solution orientée base de données.)


1 commentaires

Comme indiqué initialement, la plate-forme matérielle ne prend en charge que la simple précision ou "float" (format binary32 actuel de IEEE754). Le CPU n'a simplement aucune opération pour binary64 (ou double). Le compteur lui-même n'est pas sous mon contrôle, il est géré par le compteur d'énergie et simplement transféré sur le bus de terrain. Pour reformuler: La valeur binaire64 transférée doit être compressée dans un type 32 bits avec une plage d'au moins 1e9 (car le compteur déborde à 1e12 Wh mais nous n'avons besoin que d'une précision en kWh). Comme mentionné ci-dessous la réponse de @Eric, sa solution fonctionne.